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

基于EDEM的糧食立筒倉卸料動(dòng)態(tài)側(cè)壓力研究

2023-12-14 11:11:22劉強(qiáng)張茜茜劉濤丁永剛許啟鏗孫啟帥
科學(xué)技術(shù)與工程 2023年32期

劉強(qiáng), 張茜茜, 劉濤, 丁永剛, 許啟鏗, 孫啟帥

(1.河南工業(yè)大學(xué)土木工程學(xué)院, 鄭州 450001; 2.河南省糧油倉儲(chǔ)建筑與安全重點(diǎn)實(shí)驗(yàn)室, 鄭州 450001;3.機(jī)械工業(yè)第六設(shè)計(jì)研究院有限公司, 鄭州 450001)

立筒倉以其占地面積小、易于機(jī)械化、自動(dòng)化作業(yè)、流通費(fèi)用低、造價(jià)低等特點(diǎn),被廣泛應(yīng)用于糧食儲(chǔ)備。但貯料[1]與倉體存在著復(fù)雜的相互作用,貯料會(huì)受到倉壁的約束作用,在卸料過程中又會(huì)對(duì)倉壁產(chǎn)生一定的反作用力,發(fā)生超壓現(xiàn)象從而加大倉體破壞的風(fēng)險(xiǎn),影響糧食的供給保障。

隨著筒倉需求量越來越大,相關(guān)的研究也在不斷完善。Wang等[2]運(yùn)用細(xì)砂和玻璃微珠兩種散體材料,對(duì)比研究筒倉卸料時(shí)試驗(yàn)和模擬的流動(dòng)規(guī)律及法向壁壓分布情況,得出試驗(yàn)和模擬的結(jié)果基本一致。Wójcik等[3]通過多次裝卸料試驗(yàn),得出在整體流狀態(tài)時(shí),最大水平壓力和豎向壓力位于倉壁底部范圍內(nèi)。Fullard[4]采用PIV對(duì)筒倉的重力驅(qū)動(dòng)卸料試驗(yàn)進(jìn)行錄像分析,表明了在一個(gè)比例時(shí)間內(nèi),遠(yuǎn)離排放孔,流動(dòng)行為是相同的。原方等[5]對(duì)比分析普通筒倉安裝減壓管時(shí)筒倉單側(cè)壁卸料的動(dòng)態(tài)側(cè)壓力,得出減壓管能夠減小卸料時(shí)儲(chǔ)料對(duì)倉壁的沖擊,從而提高卸料時(shí)筒倉的穩(wěn)定性。王振清等[6]進(jìn)行了50%、80%和100%倉容3 種狀態(tài)下的側(cè)壓力試驗(yàn),發(fā)現(xiàn)管狀流動(dòng)的出現(xiàn)位置與初始的儲(chǔ)糧倉容相關(guān),結(jié)果又與規(guī)范值進(jìn)行對(duì)比,可為糧食筒倉設(shè)計(jì)提供參考。徐志軍等[7]分析顆粒在普通筒倉雙側(cè)壁卸料和帶流槽側(cè)壁卸料過程中的力學(xué)行為,驗(yàn)證了溜槽的減壓機(jī)理。程遠(yuǎn)浩[8]將曲線漏斗筒倉與錐形漏斗筒倉對(duì)比,曲線漏斗改善流態(tài),從而降低了筒倉動(dòng)態(tài)側(cè)壓力。石鑫等[9]建立的改進(jìn)顆粒模型與傳統(tǒng)單元模型比較,發(fā)現(xiàn)其倉壁側(cè)壓力峰值顯著減小,力鏈分叉現(xiàn)象明顯。余漢華等[10]采用3種機(jī)器學(xué)習(xí)方法預(yù)測(cè)筒倉動(dòng)態(tài)側(cè)壓力,同時(shí)進(jìn)行隨機(jī)抽樣得到筒倉動(dòng)態(tài)側(cè)壓力的概率分布,為筒倉結(jié)構(gòu)的可靠度提供了理論基礎(chǔ)。李坤由等[11]通過EDEM(engineering discrete element method)離散元軟件,對(duì)比分析不同摩擦情況下力鏈的細(xì)觀參數(shù),得出顆粒間摩擦系數(shù)影響卸料時(shí)間,增大顆粒間滾動(dòng)摩擦?xí)黾庸靶?yīng)。

綜上可知,前人研究中針對(duì)卸料過程中形成的不同流動(dòng)狀態(tài)沒有具體分析其成因;在靜態(tài)、動(dòng)態(tài)側(cè)壓力等方面的研究成果雖有很多,但是對(duì)于超壓系數(shù)沒有固定的取值,且最大值范圍也并不確定;同時(shí)卸料過程中形成的動(dòng)態(tài)拱具有連續(xù)形成和倒塌的特性,會(huì)對(duì)研究帶來較大的困難,有關(guān)動(dòng)態(tài)拱的形成機(jī)理仍缺少相應(yīng)解釋。為此,應(yīng)用離散元軟件EDEM(engineering discrete element method),對(duì)比分析糧食立筒倉卸料時(shí)靜態(tài)、動(dòng)態(tài)側(cè)壓力和超壓系數(shù)差異規(guī)律,探究筒倉卸料時(shí)成拱機(jī)理與倉壁側(cè)壓力的關(guān)系,為立筒倉動(dòng)態(tài)側(cè)壓力提供一種新的研究途徑。

1 離散元模型建立

1.1 立筒倉模型

以某糧庫儲(chǔ)藏大豆的柱承式鋼筋混凝土立筒倉為研究對(duì)象,在本團(tuán)隊(duì)已開展的試驗(yàn)研究[12]基礎(chǔ)上,通過EDEM(engineering discrete element method)離散元軟件建立糧食立筒倉模型,采用有機(jī)玻璃材料屬性,模型總高度為1.62 m,其中倉壁高度1.28 m,內(nèi)徑0.48 m,外徑0.50 m。結(jié)構(gòu)模型物理參數(shù)泊松比為0.33,密度為1 180 kg/m3,彈性模量為3.05×103MPa,結(jié)構(gòu)模型整體與生成貯料區(qū)域重力加速度為9.8 m/s2,糧食立筒倉模型如圖1所示。

圖1 立筒倉模型Fig.1 Silo model

1.2 貯料顆粒模型

所采用的貯料為陶粒,其為球形顆粒形態(tài),在EDEM(Engineering Discrete Element Method)軟件前處理模塊的Bulk Material中創(chuàng)建貯料顆粒,選取Spheres,設(shè)置半徑為0.001 5 m,將顆粒的粒徑分布設(shè)置為正態(tài)分布,平均值設(shè)為1即表示所有顆粒平均半徑大小在長度方向與當(dāng)前所建顆粒模型大小一致,標(biāo)準(zhǔn)差設(shè)為0.05,并通過Calculate Properties自動(dòng)計(jì)算出顆粒質(zhì)量、體積和轉(zhuǎn)動(dòng)慣量X、Y、Z。

考慮到所用散體貯料為陶粒,顆粒數(shù)量多且在卸料過程中顆粒之間、顆粒與筒倉結(jié)構(gòu)之間接觸時(shí)間較長,可通過加速度定律求解顆粒間的接觸碰撞力;同時(shí)顆粒之間、顆粒與筒倉結(jié)構(gòu)之間無黏結(jié),故選擇Hertz-Mindlin(no-slip)模型以模擬其接觸關(guān)系。依據(jù)材性試驗(yàn)與相關(guān)文獻(xiàn)[13-14],顆粒接觸屬性參數(shù)與流動(dòng)參數(shù)如表1所示。滿倉裝糧線為立筒倉頂部下降0.1 m,即貯料計(jì)算高度hn=1.18 m(hn為裝糧線至倉身下部的高度)。

表1 仿真參數(shù)Table 1 Simulation parameters

1.3 測(cè)點(diǎn)選定

根據(jù)立筒倉模型結(jié)構(gòu)的對(duì)稱性和試驗(yàn)條件,仿真模型數(shù)據(jù)測(cè)點(diǎn)的選取與試驗(yàn)[12]所布測(cè)點(diǎn)位置保持一致,其中P1~P5用來采集貯料顆粒對(duì)倉壁產(chǎn)生的側(cè)壓力,分別對(duì)應(yīng)測(cè)點(diǎn)深度0.2、0.5、0.8、1.1、1.256 m,且由于卸料方式為中心卸料,同一高度處的側(cè)壓力數(shù)據(jù)是相同的,故選取一列數(shù)據(jù)進(jìn)行分析。因筒倉為三維對(duì)稱結(jié)構(gòu),故截取筒倉的中部垂直薄片進(jìn)行數(shù)據(jù)分析,該垂直截面薄片設(shè)置為20 mm厚。在與側(cè)壓力測(cè)點(diǎn)P1~P5對(duì)應(yīng)高度處,設(shè)置顆粒網(wǎng)格組(高度30 mm),其中V1~V5測(cè)點(diǎn)以獲取與側(cè)壓力測(cè)點(diǎn)高度一致的不同高度處顆粒速度,如圖2所示。

1.4 卸料工況

在進(jìn)行卸料側(cè)壓力分析時(shí),由于不同流速狀態(tài)下所產(chǎn)生的動(dòng)態(tài)側(cè)壓力[15]與超壓系數(shù)會(huì)存在明顯不同,故將其分為兩種工況,流速慢與流速快工況分別命名為卸料a、卸料b,進(jìn)而分析兩種流速卸料工況下糧食立筒倉靜態(tài)側(cè)壓力、動(dòng)態(tài)側(cè)壓力與超壓系數(shù)。

2 模擬結(jié)果與分析

2.1 貯料流動(dòng)狀態(tài)與速度分布

根據(jù)EDEM數(shù)值仿真和試驗(yàn)結(jié)果,卸料流速不同時(shí),兩種工況下的流動(dòng)狀態(tài)性質(zhì)相同,故進(jìn)行流動(dòng)狀態(tài)與速度分布分析時(shí),選取卸料a(流速慢)工況為例,分別采用直接觀察和整體流量指數(shù)法進(jìn)行卸料流動(dòng)狀態(tài)分析。

2.1.1 直觀觀察卸料流態(tài)

通過截取不同時(shí)刻所選取的筒倉垂直截面薄片,直接觀察筒倉卸料流態(tài),如圖3(a)~圖3(g)所示;同時(shí)根據(jù)筒倉卸料的實(shí)際工作情況,選取貯料頂部的5個(gè)顆粒進(jìn)行標(biāo)記,隱藏其它顆粒,觀察貯料顆粒的運(yùn)動(dòng)軌跡,如圖3(h)所示。

分圖題為卸料時(shí)間

由圖3(a)可知,當(dāng)糧食達(dá)到穩(wěn)定狀態(tài)后,筒倉開始進(jìn)行卸料;由圖3(b)、圖3(c)可知,貯料為整體流動(dòng)狀態(tài),貯料上部表現(xiàn)為整體下降,下部呈漏斗流,卸料流動(dòng)下降較緩慢;由圖3(d)、圖3(e)可知,當(dāng)貯料下降至0.8 m測(cè)點(diǎn)深度處時(shí),開始出現(xiàn)管狀流動(dòng)趨勢(shì),中間部位貯料流速明顯更快于倉壁兩側(cè)的貯料流速,倉壁附近的貯料顆粒向中心移動(dòng),并逐漸變?yōu)楣軤盍鲃?dòng)階段,貯料下降明顯,流動(dòng)均勻;由圖3(f)、圖3(g)可知,在卸料的末期,貯料表現(xiàn)出漏斗狀,流動(dòng)加快,以卸料口正上方范圍內(nèi)的顆粒速度最快先行流出,而靠近漏斗壁的顆粒由于阻礙作用形成死料區(qū)導(dǎo)致流出較晚,直至卸料完成。

隨著貯料顆粒的流動(dòng),在筒倉下部的顆粒之間、顆粒與倉壁之間的摩擦和擠壓更顯著,導(dǎo)致流動(dòng)阻力增大,并在筒倉下部范圍內(nèi)達(dá)到最大,但是此時(shí)筒倉上部的貯料顆粒不受此影響,阻力較小且流動(dòng)較快,使得筒倉下部貯料更加密實(shí),流動(dòng)更加困難,相對(duì)應(yīng)使得倉壁側(cè)壓力增加。

2.1.2 速度分布與整體流量指數(shù)

通過建立的V1~V5顆粒網(wǎng)格組,對(duì)整體卸料的貯料顆粒速度分布情況進(jìn)行統(tǒng)計(jì)分析,根據(jù)模型內(nèi)徑寬度,每個(gè)網(wǎng)格組在水平方向可均分為19個(gè)小網(wǎng)格,用于速度數(shù)據(jù)的采集,所提取速度代表為某一時(shí)刻顆粒經(jīng)過該網(wǎng)格組時(shí)的豎向平均速度,如圖2(b)所示,在不同高度的豎向平均速度如表2所示。

表2 不同高度處豎向平均速度分布Table 2 Vertical average velocity distribution at different heights

通過整體流量指數(shù)[16](global flow index, MFI)對(duì)卸料流態(tài)情況進(jìn)行量化分析,分析每個(gè)顆粒網(wǎng)格組速度分布,計(jì)算公式為

(1)

式(1)中:Vwall為筒倉內(nèi)壁邊界處的貯料顆粒運(yùn)動(dòng)速度;Vcenter為筒倉內(nèi)中心處的貯料顆粒運(yùn)動(dòng)速度。

當(dāng)整體流量指數(shù)值大于0.3時(shí)該貯料流態(tài)為整體流,小于0.3時(shí)該貯料流態(tài)為漏斗流。因此根據(jù)表2計(jì)算不同高度處MFI值,其中V1和V2高度處的貯料表現(xiàn)為整體流,而從V3開始,貯料改變流態(tài),逐漸呈漏斗流,結(jié)合直觀觀察卸料流態(tài),并與試驗(yàn)流態(tài)進(jìn)行對(duì)比,兩者卸料時(shí)的流動(dòng)狀態(tài)大致相同。

由圖4可知,在卸料整體過程中,中間貯料流動(dòng)速度一直大于倉壁貯料流動(dòng)速度,且筒倉倉壁[圖3(h)]豎向范圍內(nèi)由上至下的不同深度位置的速度值逐漸減小,中間部位[圖3(h)]豎向范圍內(nèi)由上至下的不同深度位置的速度值逐漸增大,即MFI值在逐漸減小,說明越靠近漏斗部分,貯料卸料流動(dòng)狀態(tài)越接近漏斗流;在筒倉漏斗部分(V5測(cè)點(diǎn)處),因其接近于卸料口,故一直表現(xiàn)為漏斗流,筒倉在進(jìn)行卸料作業(yè)時(shí),其貯料所表現(xiàn)出的流態(tài)與所處于筒倉的位置有關(guān)。

括號(hào)中數(shù)值為MFI

2.2 立筒倉倉體側(cè)壓力分析

2.2.1 靜態(tài)側(cè)壓力

圖5為立筒倉滿倉靜態(tài)側(cè)壓力模擬值、試驗(yàn)值和規(guī)范值對(duì)比結(jié)果。由圖5可知,靜態(tài)側(cè)壓力的模擬值同試驗(yàn)值、規(guī)范值[17]一致,均隨著測(cè)點(diǎn)深度的增加均呈增大趨勢(shì),最大值出現(xiàn)在倉壁與漏斗連接部位范圍內(nèi),靜態(tài)側(cè)壓力的模擬值較大于規(guī)范值,規(guī)范公式計(jì)算值具有一定的局限性且偏于保守,本次卸料仿真的模擬值與試驗(yàn)值接近,為后續(xù)的動(dòng)態(tài)側(cè)壓力與超壓系數(shù)分析提供可靠的依據(jù)。

圖5 立筒倉滿倉靜態(tài)側(cè)壓力模擬值、試驗(yàn)值和規(guī)范值Fig.5 Simulation, test and specification values of static lateral pressure of silo full

2.2.2 動(dòng)態(tài)側(cè)壓力與超壓系數(shù)

圖6為立筒倉兩種卸料工況下動(dòng)態(tài)側(cè)壓力時(shí)程曲線。由圖6(a)可知,13 s時(shí)在0.8 m測(cè)點(diǎn)深度處壓力突增明顯,動(dòng)態(tài)側(cè)壓力達(dá)8.79 kPa,是靜態(tài)側(cè)壓力的2.14倍;由圖6(b)可知,14 s時(shí)在0.8 m測(cè)點(diǎn)深度處壓力突增明顯,動(dòng)態(tài)側(cè)壓力達(dá)8.77 kPa,是靜態(tài)側(cè)壓力的2.27倍。在卸料a與卸料b模擬時(shí)間0~150 s和0~48 s,貯料向下滑移,由于貯料顆粒之間、顆粒與倉壁之間的相互擠壓摩擦,一直對(duì)倉壁產(chǎn)生沖擊,致使倉壁上的P1~P4測(cè)點(diǎn)均出現(xiàn)了壓力突增的現(xiàn)象,并且發(fā)生震蕩現(xiàn)象;在漏斗處即測(cè)點(diǎn)P5處,由于離卸料口處最近,且在此處流態(tài)型式一直為漏斗流,會(huì)相應(yīng)產(chǎn)生一部分死料區(qū),故從卸料0 s開始,動(dòng)態(tài)側(cè)壓力明顯迅速下降至某一段范圍內(nèi)發(fā)生震蕩現(xiàn)象直至試驗(yàn)結(jié)束,且無超壓現(xiàn)象發(fā)生。

圖6 立筒倉兩種卸料工況下動(dòng)態(tài)側(cè)壓力時(shí)程曲線Fig.6 Dynamic lateral pressure time history curve of silo under two discharge conditions

在卸料過程中呈現(xiàn)出的震蕩現(xiàn)象,是由于卸料時(shí)貯料流動(dòng)的過程中在1/3hn范圍內(nèi)形成了超壓動(dòng)力拱,使貯料內(nèi)部反復(fù)的結(jié)拱和塌拱,拱腳對(duì)倉壁不斷進(jìn)行沖擊造成。

由表3、圖7可知,EDEM仿真模擬所得的靜態(tài)、動(dòng)態(tài)側(cè)壓力數(shù)據(jù)值與試驗(yàn)做對(duì)比,兩者相差較小;筒倉倉壁且除漏斗處的各測(cè)點(diǎn)動(dòng)態(tài)側(cè)壓力值相較于靜態(tài)側(cè)壓力值均有不同程度的增加,存在超壓現(xiàn)象,產(chǎn)生這種現(xiàn)象的原因是由于上部貯料重力和下部貯料的持續(xù)流出,導(dǎo)致發(fā)生成拱和破拱的現(xiàn)象,在0.8 m測(cè)點(diǎn)深度處成拱和破拱速度更快。在漏斗處即1.265 m測(cè)點(diǎn)深度發(fā)生動(dòng)態(tài)側(cè)壓力減小的現(xiàn)象。隨著物料的流出,倉內(nèi)貯料不斷地重復(fù)“結(jié)拱起始—結(jié)拱完成—拱塌落”的動(dòng)態(tài)拱過程[18],因而倉壁卸料壓力呈現(xiàn)“較小—局部較大—較小”的震蕩分布。

表3 立筒倉側(cè)壓力模擬值與試驗(yàn)值Table 3 Simulation and test values of silo lateral pressure

圖7 立筒倉靜態(tài)側(cè)壓力和動(dòng)態(tài)側(cè)壓力Fig.7 Silo static lateral pressure and dynamic lateral pressure

由圖8可知,超壓系數(shù)整體呈現(xiàn)先升后降的趨勢(shì),其中在0.8 m測(cè)點(diǎn)深度處最大,模擬卸料a超壓系數(shù)(2.14)和模擬卸料b超壓系數(shù)(2.27)均超過規(guī)范中水平側(cè)壓力修正系數(shù)Ch,其他測(cè)點(diǎn)超壓系數(shù)均在規(guī)范值內(nèi);在0.8 m測(cè)點(diǎn)深度處模擬卸料a、b相較于試驗(yàn)卸料a、b的超壓系數(shù)誤差在3.4%、7.6%,差別較小且在合理誤差范圍內(nèi);且對(duì)比在不同測(cè)點(diǎn)深度處兩個(gè)工況超壓系數(shù)的大小,可發(fā)現(xiàn)不同卸料流速下,流速快工況的超壓系數(shù)均大于流速慢工況的超壓系數(shù)。

3 結(jié)論

通過對(duì)糧食立筒倉卸料進(jìn)行EDEM(Engineering Discrete Element Method)仿真模擬分析,通過直接觀察和整體流量指數(shù)分析貯料流態(tài),對(duì)比研究了不同流速下模擬值與試驗(yàn)結(jié)果的貯料靜態(tài)壓力數(shù)據(jù)、動(dòng)態(tài)側(cè)壓力分布規(guī)律與超壓系數(shù),得出如下結(jié)論。

(1)筒倉下部1/3hn范圍內(nèi)貯料的流動(dòng)狀態(tài)發(fā)生顯著變化;貯料所表現(xiàn)出的流態(tài)與所處于筒倉的位置有關(guān),越靠近漏斗部分,貯料卸料流動(dòng)狀態(tài)越接近漏斗流;卸料的發(fā)展進(jìn)程為整體流動(dòng)—管狀流動(dòng)—漏斗流。

(2)靜態(tài)側(cè)壓力數(shù)值模擬值與試驗(yàn)值較為接近、動(dòng)態(tài)側(cè)壓力分布規(guī)律一致,EDEM離散元卸料仿真模型動(dòng)態(tài)側(cè)壓力研究具有可行性和合理性。

(3)卸料過程中,倉壁出現(xiàn)側(cè)壓力震蕩增大現(xiàn)象是由于在倉體下部1/3hn范圍內(nèi)出現(xiàn) “成拱—塌拱—再成拱—再塌拱”的動(dòng)態(tài)拱過程。

(4)超壓系數(shù)最大值出現(xiàn)于0.8 m測(cè)點(diǎn)深度處;相同位置范圍內(nèi)不同流速下的超壓系數(shù)均高于規(guī)范中水平側(cè)壓力修正系數(shù)Ch;不同測(cè)點(diǎn)深度處,流速快工況的超壓系數(shù)均大于流速慢工況的超壓系數(shù)。

主站蜘蛛池模板: 一级毛片免费播放视频| 中文毛片无遮挡播放免费| 亚洲中文在线看视频一区| 亚洲精品另类| 国产精品亚洲精品爽爽| 欧美va亚洲va香蕉在线| 免费在线成人网| 亚洲第一色视频| 99re热精品视频国产免费| 亚洲高清在线播放| 日本精品视频| 日本欧美午夜| 久久亚洲国产最新网站| 欧美v在线| 99久久国产精品无码| 欧美日韩久久综合| 久久精品国产免费观看频道| 国产色网站| 亚洲天堂区| 波多野结衣中文字幕一区二区| 久久久久无码国产精品不卡| 狼友av永久网站免费观看| 成人毛片免费观看| 亚洲欧美综合另类图片小说区| 国产亚洲一区二区三区在线| 狠狠色综合网| 欧美视频免费一区二区三区| 亚洲人成色77777在线观看| Jizz国产色系免费| 韩国福利一区| 色视频国产| 白浆免费视频国产精品视频| 国产一区二区三区夜色| 日韩一区精品视频一区二区| 欧美激情视频二区三区| 欧美日在线观看| 国产成人亚洲无吗淙合青草| 国产系列在线| 国产香蕉一区二区在线网站| 婷婷色中文| 亚洲人成日本在线观看| 国产精品福利尤物youwu | 国产女人18水真多毛片18精品| 国产精品美女网站| 久久一色本道亚洲| 在线免费a视频| 日本爱爱精品一区二区| 视频一区亚洲| 热99re99首页精品亚洲五月天| 亚洲日韩精品伊甸| 欧美一区福利| 亚洲精品波多野结衣| 国产传媒一区二区三区四区五区| 亚洲有无码中文网| 久久99国产综合精品1| 欧美第一页在线| 日本人妻一区二区三区不卡影院| 99色亚洲国产精品11p| 亚洲伊人电影| 国产午夜人做人免费视频中文| 日韩一级毛一欧美一国产| 国产精品偷伦视频免费观看国产| 在线中文字幕日韩| 国产精品林美惠子在线观看| 欧美一级黄片一区2区| 国产免费高清无需播放器| 22sihu国产精品视频影视资讯| 国产精品流白浆在线观看| 欧美日韩中文国产va另类| 亚洲一区二区三区香蕉| 一本久道热中字伊人| 久久99热这里只有精品免费看| 国产成人精品男人的天堂下载 | 亚洲大学生视频在线播放| 噜噜噜久久| 国内精品视频在线| 国产成人8x视频一区二区| 亚洲大学生视频在线播放| 在线国产毛片| 国内精品九九久久久精品| 日本高清有码人妻| 亚洲精品国产精品乱码不卞|