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

采用總荷載不變法的非靜水壓隧道攝動拓展解

2023-06-13 09:20:06張常光李宗輝關港輝
哈爾濱工業(yè)大學學報 2023年6期
關鍵詞:圍巖模型

張常光,李宗輝,關港輝,孫 松

(1.長安大學 建筑工程學院,西安 710061;2.地質災害防治與地質環(huán)境保護國家重點實驗室(成都理工大學),成都 610059)

在水利、交通、油氣等設施建設中存在大量隧道工程,實際又以水平地應力與豎向地應力存在明顯差異的非靜水應力場為主,研究非靜水壓力下非圓塑性區(qū)邊界線、塑性區(qū)位移分布及圍巖特征曲線對隧道支護設計具有重要意義。當圍巖塑性區(qū)完全包圍隧道時,非靜水壓力下圓形隧道的塑性區(qū)應力可聯(lián)立圍巖屈服條件和平衡微分方程獲得,進而圍巖彈性區(qū)應力計算成為確定隧道塑性區(qū)半徑和位移的關鍵。基于Mohr-Coulomb準則和理想彈-塑性模型,已有學者采用Kastner法[1-7]、復變函數(shù)法[8-12]、攝動法[13-15]和總荷載不變法[16-18],探討了非靜水壓力下圓形隧道的彈性區(qū)應力和非圓塑性區(qū)邊界線。

Kastner[1]將非靜水壓圓形隧道純彈性狀態(tài)時的應力基爾希公式,直接代入Mohr-Coulomb準則建立隧道非圓塑性區(qū)的邊界線方程。趙志強等[2-6]基于Kastner法提出隧道的蝶形塑性區(qū)理論,于學馥等[7]考慮應力重分布改進了Kastner法。需指出的是,Kastner法不計圍巖彈塑性發(fā)展,也不顧非靜水壓隧道的塑性區(qū)邊界線不是圓形,而使用內(nèi)邊界要求為圓形孔洞的應力基爾希公式。Detournay等[8-12]借助復變函數(shù)法將非圓形的隧道塑性區(qū)邊界保角映射為圓形,難點在于如何有效確定隧道彈性區(qū)的應力復變函數(shù)形式及眾多系數(shù)。魏悅廣等[13-15]依據(jù)攝動法將隧道非圓塑性區(qū)邊界和彈性區(qū)應力求解轉化為泰勒級數(shù)展開式,獲得了精度較高的圍巖彈性區(qū)應力和類橢圓塑性區(qū)邊界。嚴克強等[16-18]選擇基爾希公式作為圍巖彈性區(qū)應力表達式,由總荷載不變法推導隧道水平軸和豎向軸處的塑性區(qū)半徑方程,并由幾何相似原理給定其他方位角處的,但彈性區(qū)應力基爾希公式的內(nèi)邊界非圓問題依然存在。

總荷載不變法從基本的力平衡原理出發(fā),理論基礎簡單,易于被工程接受及實用推廣,同時,以簡潔公式表達的圍巖彈性區(qū)應力攝動解可逐階逼近真實結果達到高精度。此外,峰后巖石強度常有一定跌落,彈-脆-塑性模型相比理想彈-塑性模型更適合描述巖石真實強度變化。因此,針對非靜水壓力作用下的圓形隧道,基于Mohr-Coulomb準則和彈-脆-塑性模型并引入圍巖彈性區(qū)應力的攝動解,首先采用總荷載不變法建立水平軸與豎向軸處的塑性區(qū)半徑方程,繼而利用幾何相似原理拓展至其他方位角處,對比文獻總荷載不變法(以應力基爾希公式為基礎)、Kastner法和復變函數(shù)法進行合理性與正確性驗證,結合非關聯(lián)流動法則推導塑性區(qū)位移解析解,探討側壓力系數(shù)和脆性軟化對隧道塑性區(qū)邊界線、塑性區(qū)位移分布以及圍巖特征曲線的影響規(guī)律。

1 總荷載不變法

假定圍巖為均勻、連續(xù)的各向同性材料,圖1為非靜水壓力下圓形隧道的平面應變力學模型,洞壁處作用均勻支護力pi,無窮遠處作用豎向地應力q、水平地應力λq,λ為側壓力系數(shù)。另外,a為隧道半徑;Rp為塑性區(qū)半徑;r為某點到隧道中心的距離;θ為方位角,以水平右半軸為起點、逆時針旋轉為正。

圖1 非靜水壓隧道力學模型(λ<1)

對于非靜水壓力下圓形隧道的彈性區(qū)和塑性區(qū),由地應力的雙軸對稱性知在水平軸和豎向軸處反對稱的切應力恒為0[16-18],進而可對隧道開挖前/后的圍巖環(huán)向正應力進行積分。總荷載不變法依據(jù)最基本的力平衡原理,對非靜水壓圓形隧道分為開挖前/后圍巖的豎向總荷載不變、開挖前/后圍巖的水平向總荷載不變,由方位角θ=0°處(θ面,存在環(huán)向正應力,切應力為零)圍巖豎向總荷載不變可獲得水平軸處塑性區(qū)半徑RH的求解方程,同理,由方位角θ=90°處(θ面,存在環(huán)向正應力,切應力為零)圍巖水平向總荷載不變可獲得豎向軸處塑性區(qū)半徑Rv的求解方程,下文將以沿水平軸的豎向總荷載不變?yōu)槔治觥?/p>

依據(jù)豎向總荷載不變,得式(1a)中隧道開挖前沿水平軸的圍巖豎向總荷載Q0與式(1b)中隧道開挖后的圍巖豎向總荷載Q1相等。

Q0=qL

(1a)

(1b)

2 應力解答

2.1 塑性區(qū)

采用Mohr-Coulomb準則描述圍巖強度的彈-脆-塑性特征,其中式(2)對應圍巖初始強度,式(3)對應圍巖峰后強度。

σ1=Aiσ3+Bi

(2)

σ1=Arσ3+Br

(3)

式中:σ1為第一主應力,σ3為第三主應力,并以壓應力為正;Ai、Bi、Ar、Br為方程參數(shù);ci、φi為圍巖初始的黏聚力和內(nèi)摩擦角(用于彈-塑性交界處彈性區(qū)一側的圍巖強度);cr、φr為圍巖峰后的黏聚力和內(nèi)摩擦角(用于隧道塑性區(qū)內(nèi)的圍巖強度)。

式(4)表示圍巖主應力與各應力分量之間的關系[19],其中,σr、σθ、τrθ分別為徑向正應力、環(huán)向正應力和切應力。

(4a)

(4b)

將式(4)分別代入式(2)和(3)得非靜水壓力下隧道圍巖初始/峰后的強度關系式為

(5)

(6)

平衡微分方程為

(7a)

(7b)

聯(lián)立式(6)和(7),并結合圖1隧道洞壁處的應力邊界條件,即r=a時,σr=pi、τrθ=0,求得圍巖塑性區(qū)的應力為[20]

(8a)

(8b)

(8c)

2.2 彈性區(qū)

(9a)

(9b)

(9c)

需注意的是,非靜水壓圓形隧道的塑性區(qū)并非圓形,這與基爾希公式要求的圓形內(nèi)邊界不一致。

魏悅廣等[13,15]基于Mohr-Coulomb準則和理想彈-塑性模型,設攝動參數(shù)ε=1-λ,并采用攝動法將非圓塑性區(qū)邊界求解轉化為泰勒級數(shù)展開式,以逐階逼近方式得到精度較高的圍巖彈性區(qū)應力攝動解,又因2階攝動解滿足精度要求且表達較3階簡潔,引入式(10)表示的非靜水壓圓形隧道彈性區(qū)應力2階攝動解,其中,(1-λ)、(1-λ)2項分別代表ε1階攝動增量、ε2階攝動增量。

(10a)

(10b)

(10c)

對于彈-脆-塑性模型,彈-塑性交界處彈性區(qū)一側采用圍巖初始強度參數(shù),塑性區(qū)采用圍巖峰后強度參數(shù),重新推導得非靜水壓力下彈-脆-塑性圓形隧道彈性區(qū)應力的2階攝動解為

(11a)

(11b)

(11c)

3 塑性區(qū)半徑

3.1 水平軸和豎向軸

將式(8b)和(11b)代入式(1b),并令θ=0°得

(12)

依據(jù)沿水平軸的圍巖豎向總荷載不變,將式(1a)和(12)代入Q1=Q0,化簡得非靜水壓力下圓形隧道水平軸處的塑性區(qū)半徑(RH)方程為

(13)

(14)

(15)

同理,根據(jù)沿豎向軸的圍巖水平向總荷載不變,求得非靜水壓力下圓形隧道豎向軸處的塑性區(qū)半徑(RV)方程為

(16)

(17)

(18)

3.2 塑性區(qū)拓展

蔡曉鴻等[21]基于Mohr-Coulomb準則和理想彈-塑性模型,根據(jù)彈-塑性交界處應力連續(xù)獲得非靜水壓力下圓形隧道的非圓塑性區(qū)邊界線,結果具有一般性,但理想彈-塑性模型難以反映圍巖強度的脆性變化。董海龍等[17]基于Mohr-Coulomb準則和彈-脆-塑性模型,選擇基爾希公式作為圍巖彈性區(qū)應力表達式,依據(jù)總荷載不變法拓展了文獻[21]的塑性區(qū)邊界線方程。

(19a)

(19b)

Rp=αRC

(19c)

式中:RC(θ=0°)、RC(θ=90°)為文獻[21]非靜水壓力下圓形隧道水平軸和豎向軸處的塑性區(qū)半徑。

1)RC(θ=0°)和RC(θ=90°)均不小于隧道半徑a時,可按式(19)拓展,且λ<1的塑性區(qū)短軸為RC(θ=90°),λ>1的塑性區(qū)短軸為RC(θ=0°),對應如下判定條件:

當λ≤1時

(20)

當λ>1時

(21)

3.3 拓展比較

采用總荷載不變法拓展建立的非靜水壓力下圓形隧道塑性區(qū)半徑解析解即式(19),引入了精度較高的圍巖彈性區(qū)應力2階攝動解,同時考慮了圍巖強度的脆性軟化,理想彈-塑性模型解答為其特例,具有一定的理論意義和工程應用前景。此外,所得解析解對側壓力系數(shù)λ<1、λ=1和λ>1均適用,下文將針對λ<1進行討論,并對比λ=1時靜水壓力的結果。

為說明本文解析解和文獻[17]對文獻[21]的拓展不同,以及總荷載不變法下2階攝動解即式(11)與基爾希公式即式(9)作為彈性區(qū)應力表達式時塑性區(qū)邊界線的差異,表1和圖3給出了本文與文獻[17,21]的結果比較,其中,文獻[17]以水平軸和豎向軸處的拓展分別記為文獻[17]a和[17]b。隧道算例的參數(shù)為[17]:a=2.43 m,pi=0 MPa,q=21.78 MPa,L=∞;彈-脆-塑性圍巖ci=4.8 MPa、φi=32°,cr=1.8 MPa、φr=20°;理想彈-塑性圍巖cr=ci=4.8 MPa,φr=φi=32°。

表1 塑性區(qū)半徑的比

圖3 塑性區(qū)邊界線的不同拓展

由表1可知,本文水平軸處的塑性區(qū)半徑是文獻[21]的1.68~1.88倍,豎向軸處為1.68~1.99倍,這是由于文獻[21]基于理想彈-塑性模型未考慮圍巖強度的峰后下降。

3.4 對比驗證

Kastner[1]直接將隧道純彈性狀態(tài)時的應力基爾希公式代入圍巖強度準則,得到了非靜水壓力下圓形隧道塑性區(qū)半徑的近似方程;呂愛鐘等[10]推導了非靜水壓圓形隧道塑性區(qū)半徑的復變函數(shù)解答。這2種方法均基于Mohr-Coulomb準則,且假定圍巖為理想彈-塑性材料,即cr=ci、φr=φi,圖4為不同方法的隧道塑性區(qū)邊界線對比。隧道算例的參數(shù)為[10]a=2 m,pi=0 MPa,q=15 MPa,L=∞,ci=cr=1 MPa,φi=φr=30°。

圖4 不同方法的塑性區(qū)邊界線對比

由圖4可以看出,本文式(19c)與文獻[10]的高精度復變函數(shù)解答吻合良好,驗證了本文結果的正確性;Kastner法未考慮開挖后圍巖的彈塑性發(fā)展與分區(qū),簡單地使用了純彈性狀態(tài)時的應力基爾希公式,使得文獻[1]的塑性區(qū)范圍小1.45~1.84 m。

3.5 實測驗證

文獻[15]對TBM開挖的圓形煤礦巷道擾動范圍開展現(xiàn)場實測,相關參數(shù)為a=2.25 m,pi=0 MPa,q=14.3 MPa,λ=1.33,ci=9 MPa,cr=3 MPa,φi=45°,φr=42°,方位角θ=5°、90°、162°處圍巖的實測擾動區(qū)深度為2.98、3.31和2.89 m。

本文由理想彈-塑性模型所得塑性區(qū)范圍為4.33~4.49 m,對應方位角θ=5°、90°、162°處圍巖的計算塑性區(qū)深度為2.08、2.24和2.09 m,小于實測的擾動區(qū)深度;然而,由彈-脆-塑性模型所得塑性區(qū)范圍為5.28~5.52 m,對應方位角θ=5°、90°、162°處圍巖的計算塑性區(qū)深度為3.03、3.27和3.05 m,與實測的擾動區(qū)深度非常接近,表明彈-脆-塑性模型更適合求解此巷道的擾動范圍。

文獻[22]對淮南礦區(qū)某千米深井巷道松動圈范圍進行了現(xiàn)場實測,相關參數(shù)為a=2.95 m(等效半徑),pi=0.75 MPa,q=21.86 MPa,λ=0.8,ci=5.58 MPa,cr=0.72 MPa,φi=φr=27.83°,圍巖幫部、頂部和底板中部的實測松動圈深度為2.80、2.39和3.84 m。

本文在彈-脆-塑性模型下所得塑性區(qū)范圍為5.56~6.69 m,對應圍巖幫部、頂部和底板中部的計算塑性區(qū)深度為4.15、2.61和4.06 m,可見計算塑性區(qū)深度稍大于實測松動圈深度而在其外側,這符合深井巷道的開挖擾動規(guī)律。相反,由理想彈塑性模型計算所得圍巖幫部、頂部和底板中部的計算塑性區(qū)深度為1.35、0.46和2.16 m,明顯小于實測松動圈深度而不合理。

4 位移解答

4.1 彈性區(qū)

彈性區(qū)圍巖的幾何方程為

(22a)

(22b)

式中:εr、εθ、γrθ為彈性區(qū)的徑向線應變、環(huán)向線應變和切應變,ur、uθ為對應的徑向線位移和環(huán)向線位移。

將圖1中的水平/豎向初始地應力場轉換為徑向/環(huán)向初始地應力場[19],得初始徑向正應力σr0、初始環(huán)向正應力σθ0和初始切應力τrθ0為

(23)

在平面應變狀態(tài)下扣除初始地應力即式(23),由廣義Hooke定律得彈性區(qū)圍巖的物理方程為

(24a)

(24b)

(24c)

式中E、ν為圍巖的彈性模量和泊松比。

聯(lián)立式(11)、(22)~(24)得圍巖彈性區(qū)位移為

(25a)

(25b)

4.2 塑性區(qū)

(26)

在圍巖塑性區(qū)內(nèi),將應變分解為式(27)中的彈性應變分量(對應于上標e)和塑性應變分量(對應于上標p),即

(27)

采用非關聯(lián)流動法則即式(28)描述圍巖的剪脹特性[23]:

(28)

式中β為剪脹參數(shù)。

(29)

將r=a代入式(29),得隧道洞壁位移u0為

(30)

通過式(30)可構建圍巖特征曲線即洞壁位移u0隨支護力pi的變化,進而基于收斂約束法由圍巖特征曲線(GRC)和支護特征曲線(SRC)的交點縱/橫坐標確定支護壓力和圍巖穩(wěn)定變形。

4.3 結果與討論

隧道算例同3.3節(jié),并取λ=0.8。圖5、6給出了本文、文獻[17](豎向軸拓展)和[21]的塑性區(qū)位移分布與收斂約束分析,其中,討論塑性區(qū)位移分布時令支護力pi=0 MPa(下同),以虛線標出塑性區(qū)半徑與隧道半徑的比值。在收斂約束分析時,假定支護為均質等厚的圓環(huán)結構,并在u0/a=2%處施作。

圖5 塑性區(qū)位移分布

由圖5、6可以看出,本文的塑性區(qū)半徑、塑性區(qū)位移、支護壓力和圍巖穩(wěn)定變形均最大,最小的是文獻[21](塑性區(qū)半徑、洞壁位移比本文的分別小1.89~2.26、0.059~0.128 m),文獻[17]的支護壓力和圍巖穩(wěn)定變形稍大于文獻[21],但明顯小于本文的。在圖6(a)中本文對應的支護已到達承載極限,需提高支護承載能力或推后施作位置,而根據(jù)文獻[17]或[21]所做的隧道設計將存在安全隱患。

圖6 收斂約束分析

5 影響因素分析

結合本文采用總荷載不變法建立的隧道塑性攝動拓展解答即式(19c)、(29)和(30),探討側壓力系數(shù)和脆性軟化對塑性區(qū)邊界線、塑性區(qū)位移分布和收斂約束分析的影響規(guī)律。

5.1 側壓力系數(shù)

隧道算例同3.3節(jié),并取E=2 GPa,ν=0.2,β=2,圖7為兩種材料模型的塑性區(qū)邊界線比較,其中,側壓力系數(shù)λ=0.6、0.8和1。可以看出,彈-脆-塑性模型的塑性區(qū)范圍明顯大于理想彈-塑性模型的,水平軸處的塑性區(qū)半徑相差1.41~2.44 m,宜采用彈-脆-塑性模型分析隧道塑性發(fā)展。

圖7 塑性區(qū)邊界線比較

圖8、9給出了λ=0.6、0.8和1時彈-脆-塑性圍巖的塑性區(qū)位移分布與收斂約束分析。

圖8 不同λ時塑性區(qū)位移分布

由圖8(a)、(b)可以看出,λ=0.6、0.8時,在第Ⅰ象限內(nèi)塑性區(qū)位移最大在θ=0°處,最小在θ=90°處,洞壁位移相差0.099 m;在λ=0.6的圖9(a)中,不同方位角處的圍巖特征曲線及收斂約束存在顯著差異,在θ=90°處施作支護時圍巖位移已全部釋放,支護將不再受荷與變形,而在θ=0°和45°處仍存在較大的支護壓力,故非靜水壓力下需針對具體方位角進行收斂約束分析與支護設計;在λ=0.8的圖9(b)中,3個方位角處的支護壓力和圍巖穩(wěn)定變形都明顯不同,在θ=0°處支護已達到承載極限,這與靜水壓力λ=1的圖9(c)不同,依據(jù)理想靜水壓力條件開展的均勻支護設計將影響實際隧道后續(xù)施工和正常運維。

圖9 不同λ時收斂約束分析

5.2 脆性軟化

隧道算例同3.3節(jié),并取λ=0.8,E=2 GPa,ν=0.2,β=2,圖10~12給出了圍巖不同峰后強度下的塑性區(qū)邊界線和收斂約束分析。

圖10 脆性軟化對塑性區(qū)邊界線的影響

由圖10可以看出,脆性軟化幅度越大即圍巖峰后強度越低,塑性區(qū)范圍越大。峰后黏聚力cr每減小0.1 MPa,水平軸處的塑性區(qū)半徑平均增加0.18 m、豎向軸處的增加0.15 m;峰后內(nèi)摩擦角φr每減小1°,水平軸處的塑性區(qū)半徑平均增加0.16 m、豎向軸處的增加0.14 m。

由圖11、12可以看出,支護壓力和圍巖穩(wěn)定變形隨峰后強度的增加明顯減小,相應的支護可由達到承載極限過渡到安全儲備充足,通過注漿、預應力錨桿等主動加固提高圍巖的峰后強度,以改變支護類型或減小支護截面尺寸。

圖11 不同cr時收斂約束分析

圖12 不同φr時收斂約束分析

6 結 論

1)采用總荷載不變法并引入圍巖彈性區(qū)應力的2階攝動解,所推導并由幾何相似原理拓展建立的非靜水壓圓形隧道塑性區(qū)半徑與塑性區(qū)位移的解析解,考慮了圍巖強度的脆性軟化,理想彈-塑性模型解答為其特例,結合收斂約束法可有效指導隧道支護初步設計,具有一定的理論意義和工程應用前景。

2)彈性區(qū)應力表達式對非靜水壓圓形隧道的塑性區(qū)范圍影響顯著,通過對比文獻總荷載不變法、Kastner法、復變函數(shù)法和實測數(shù)據(jù),驗證了本文按插值系數(shù)拓展解答的正確性與彈性區(qū)應力2階攝動解的合理性;Kastner法的塑性區(qū)范圍過小,采用應力基爾希公式所得支護壓力和圍巖穩(wěn)定變形都偏小。

3)非靜水壓力下不同方位角處的塑性區(qū)半徑、塑性區(qū)位移分布與圍巖特征曲線存在較大差異,需根據(jù)收斂約束法確定沿洞周最大的支護壓力和圍巖穩(wěn)定變形,這有別于靜水壓力下圓形塑性區(qū)時的均勻支護設計;圍巖峰后強度全面而明顯地影響隧道塑性發(fā)展,可使支護處于承載極限、安全儲備等狀態(tài)。

猜你喜歡
圍巖模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
隧道開挖圍巖穩(wěn)定性分析
中華建設(2019年12期)2019-12-31 06:47:58
軟弱破碎圍巖隧道初期支護大變形治理技術
江西建材(2018年4期)2018-04-10 12:37:22
3D打印中的模型分割與打包
復雜巖層大斷面硐室群圍巖破壞機理及控制
煤炭學報(2015年10期)2015-12-21 01:55:09
滑動構造帶大斷面弱膠結圍巖控制技術
山西煤炭(2015年4期)2015-12-20 11:36:18
采空側巷道圍巖加固與巷道底臌的防治
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产好痛疼轻点好爽的视频| 亚洲天堂自拍| 无码日韩精品91超碰| 日韩精品专区免费无码aⅴ| 久精品色妇丰满人妻| 亚洲日韩第九十九页| 久久五月视频| 精品一区二区三区自慰喷水| 成人国产免费| 免费看美女自慰的网站| 色九九视频| 99久久免费精品特色大片| 在线观看亚洲精品福利片| 亚洲成A人V欧美综合天堂| 成人福利在线观看| 不卡网亚洲无码| 天堂va亚洲va欧美va国产| 国产激情无码一区二区免费| 国产女人爽到高潮的免费视频 | 欧美 亚洲 日韩 国产| 午夜视频免费试看| 波多野结衣在线se| 国产成人免费手机在线观看视频| 在线人成精品免费视频| 夜夜操天天摸| 青草视频网站在线观看| 亚洲男人的天堂久久香蕉| 99这里只有精品免费视频| 亚洲AⅤ综合在线欧美一区| 国产成人综合久久精品尤物| 亚洲一级毛片在线观| 蜜桃视频一区二区三区| 欧美亚洲另类在线观看| 91精品最新国内在线播放| 欧美视频免费一区二区三区| 亚洲一区网站| 久久精品只有这里有| 亚洲一区二区约美女探花| 亚洲成人网在线播放| 日韩a级毛片| 成人国产免费| 国产情侣一区二区三区| 制服丝袜 91视频| 青草午夜精品视频在线观看| 国产精品尤物在线| 亚国产欧美在线人成| 九九精品在线观看| 国产中文一区二区苍井空| 青青青视频蜜桃一区二区| 欧美日本在线一区二区三区| 狠狠躁天天躁夜夜躁婷婷| 欧美日韩国产在线人| 国产91视频免费观看| 亚洲国产系列| 中国精品久久| 国产九九精品视频| 亚洲成a人片| jizz国产在线| 亚洲综合天堂网| 54pao国产成人免费视频| 欧美性猛交一区二区三区| 精品国产aⅴ一区二区三区| 国产va在线| 成人福利在线视频免费观看| 玩两个丰满老熟女久久网| 欧美激情综合| 国产午夜一级淫片| 亚洲国产成人久久77| 久久精品国产电影| 欧美成人午夜在线全部免费| 欧美在线伊人| 国产中文在线亚洲精品官网| 亚洲第一视频网站| 精品在线免费播放| 精品人妻系列无码专区久久| 国产成人高清精品免费软件| 国产成人资源| 亚洲成a人片77777在线播放| 国产黄在线免费观看| 久久黄色免费电影| 蜜芽一区二区国产精品| 少妇精品在线|