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

        ?

        基于PDE的凍土區(qū)輸氣管道熱影響數(shù)值研究

        2022-11-30 10:23:00王翔鵬蘇文濤陳星元黃紹暢商麗鳳
        節(jié)能技術(shù) 2022年5期
        關(guān)鍵詞:模型

        王翔鵬,潘 振,蘇文濤,劉 利,陳星元,黃紹暢,商麗鳳

        (1.遼寧石油化工大學(xué) 石油天然氣工程學(xué)院,遼寧 撫順 113001;2.撫順市特種設(shè)備監(jiān)督檢驗(yàn)所,遼寧 撫順 113000)

        多年凍土是指溫度等于或低于0℃至少連續(xù)兩年的地面(土壤或巖石,包括冰)[1]。北半球四分之一的面積和17%的裸露地表被永久凍土覆蓋,并且有著大量化石燃料儲(chǔ)量[2]。中國(guó)政府已經(jīng)發(fā)布了2030年碳排放高峰指標(biāo)(CEUP)[3]。天然氣資源比煤炭更清潔,比可再生能源更便宜,因此有很大的發(fā)展空間[4]。管道作為一種安全、可靠、環(huán)保的運(yùn)輸方式,已廣泛應(yīng)用于天然氣的長(zhǎng)途運(yùn)輸[5]。近幾十年來(lái),化石能源的開(kāi)發(fā)和運(yùn)輸導(dǎo)致了凍土區(qū)管道系統(tǒng)的廣泛應(yīng)用。然而,不同工況管道建設(shè)的干擾以及管道與周圍凍土之間的熱相互作用對(duì)管道系統(tǒng)的設(shè)計(jì)、施工、安全運(yùn)行和維護(hù)構(gòu)成了重大挑戰(zhàn)[6-7]。根據(jù)管道運(yùn)行的溫度范圍,管道在運(yùn)行后可能會(huì)出現(xiàn)大量?jī)雒浐?或融化沉降。這些問(wèn)題嚴(yán)重影響凍土地區(qū)管道系統(tǒng)的安全運(yùn)行,不僅增加了施工和維護(hù)成本,還增加了管道沿線油氣泄漏的可能性[8]。

        本文設(shè)計(jì)了基于偏微分方程(PDE)的考慮冰水相變、水分遷移的非飽和凍土區(qū)耦合模型,以西部某凍土區(qū)正溫輸氣管道為例進(jìn)行模型驗(yàn)證。對(duì)新建管道或停輸再啟動(dòng)190 d內(nèi),管-土區(qū)溫度、相對(duì)飽和度、體積含水率的分布及變化進(jìn)行分析。討論了輸氣管道埋深、輸送溫度變化對(duì)凍土的熱影響。為了改善這種地面熱狀況,引用了類似道路和鐵路路堤的一種緩解熱影響管道路堤(下文簡(jiǎn)稱:管堤)的設(shè)計(jì)。管堤主要由碎石層組成。碎石層在寒冷季節(jié)通過(guò)自然空氣對(duì)流將管道中的熱量傳遞到寒冷環(huán)境中。該層已廣泛應(yīng)用于道路和鐵路路堤的施工中,以控制氣候變暖條件下多年凍土路基的溫度狀況。在冬季,高孔隙碎石層內(nèi)的自然空氣對(duì)流可能會(huì)增加永久凍土路基的熱損失,或?qū)⒏嗟睦淠軓沫h(huán)境空氣引入凍土路基[9]。為了評(píng)估設(shè)計(jì)的熱性能,對(duì)管道路堤進(jìn)行模擬,以分析其傳熱過(guò)程,為凍土地區(qū)能源管道系統(tǒng)的設(shè)計(jì)和施工提供信息參考。

        1 模型介紹

        1.1 物理模型

        圖1顯示了埋地輸氣管道和擬建的管道路堤的物理模型。管道管頂埋深設(shè)置為天然地面以下2.5 m,管道未設(shè)置保溫層,管徑和壁厚分別為660 mm,10 mm。與研究現(xiàn)場(chǎng)實(shí)際管線參數(shù)相同[10]。土壤分為三層,即草炭亞粘土(Ⅰ)、亞粘土(Ⅱ)和粉質(zhì)粘土(Ⅲ)。

        圖1 埋地輸氣管道和管堤的物理模型

        對(duì)于管堤的物理模型,路堤由0.8 m厚的碎石層(Ⅳ),1.6 m厚的粗粒土層組成。為了方便管道維修期間車輛進(jìn)入,路堤的表面寬度設(shè)計(jì)為3.5 m。輸氣管道放置在對(duì)凍脹和融沉不敏感的粗顆粒層(Ⅴ)內(nèi)。路基土壤的計(jì)算域從路堤坡腳向外延伸8 m,路基土層的地層組成、管徑、壁厚與埋地輸氣管道模型相同。

        1.2 控制方程

        為了建立非飽和凍土中水運(yùn)動(dòng)、熱傳輸和變形的控制方程,進(jìn)行了以下幾個(gè)基本假設(shè):

        (1)土壤顆粒、水和冰是均勻的、各向同性的、不可壓縮的。冰的生長(zhǎng)/膨脹是各向同性的。水冰相變是基于平衡熱力學(xué)的。

        (2)在有輸氣管道的淺層凍土中,假設(shè)始終有未凍水,水分只以液體形式遷移,遵循達(dá)西定律,不考慮冰晶的運(yùn)動(dòng)。

        1.2.1 水分場(chǎng)

        凍土中瞬態(tài)非飽和液體流動(dòng)控制方程[11]

        (1)

        式中θu——土壤中未凍水的體積含量;

        D——水力擴(kuò)散系數(shù);

        k——土體的滲透系數(shù);

        下標(biāo)w,i——表示水,冰。

        水力擴(kuò)散系數(shù)[11]

        (2)

        式中c(θu)——比水容量/m-1,由滯水模型確定,本文選用Van Genuehten(VG)滯水模型;

        I——阻抗因子[12],表達(dá)式如下

        I=10-10θi

        (3)

        引入一個(gè)無(wú)量綱含水率變量S(本文稱為“相對(duì)飽和度”),它可由飽和含水率θs和殘余含水率θr定義成標(biāo)準(zhǔn)化含水率的形式[13]

        (4)

        若殘余含水率θr等于0,相對(duì)飽和度S等于土的飽和度。本文θs和θr分別取0.4和0.08。

        土體的滲透系數(shù)選用Gardner滲透系數(shù)方程[11]

        k(θu)=ks·Sl[1-(1-S1/m)m]2

        (5)

        式中ks——飽和滲透系數(shù),9.62×10-5m/s。

        比水容量[11]

        (6)

        式中a0、m、l——隨土質(zhì)變化而變化的本構(gòu)參數(shù),本文分別取2,1/m;0.5;0.5。

        1.2.2 溫度場(chǎng)

        考慮冰水相變的凍土傳熱方程如下[14]

        (7)

        式中ρ——土壤密度/kg·m-3;

        ρi——冰的密度/kg·m-3;

        C——土壤的體積熱容/J·kg-1·℃-1;

        λ——土壤的導(dǎo)熱系數(shù)/W·m-1·℃-1;

        θ——體積含水量;

        L——相變潛熱/kJ·kg-1,334.56 kJ/kg;

        T——溫度/℃;

        t——時(shí)間/s。

        管壁中的熱傳導(dǎo)方程[14]

        (8)

        式中ρp——鋼管密度/kg·m-3;

        Cp——鋼管的體積熱容/J·kg-1·℃-1;

        λp——鋼管的導(dǎo)熱系數(shù)/W·m-1·℃-1。

        1.2.3 補(bǔ)充方程

        因?yàn)楸捏w積含量θi、未凍水體積含量θu同時(shí)存在于控制方程(1),(7)中,且含有未知變量溫度T。故還需引入一個(gè)聯(lián)系方程實(shí)現(xiàn)耦合模型求解,本文通過(guò)引入凍土孔隙中冰與未凍水之關(guān)系[13]聯(lián)立求解

        (9)

        式中Tf——土體凍結(jié)溫度/-1 ℃;

        B——常數(shù),與土壤種類和含鹽量有關(guān),本文取0.56。

        1.2.4 應(yīng)力場(chǎng)

        平衡方程[15]

        ?{σ}+[F]=0

        (10)

        幾何關(guān)系[15]

        {ε}=?{δ}

        (11)

        本構(gòu)關(guān)系[15]

        {σ}=[C]({ε}-{εv})

        (12)

        式中F——體積力;

        δ——位移量;

        C——彈性模量;

        εv——有水分相變產(chǎn)生的應(yīng)變,表達(dá)式[13]如下

        (13)

        上述數(shù)值程序由有限元程序COMSOL Multiphysics實(shí)現(xiàn),以初始條件作為初始猜測(cè)通過(guò)使用的Newton-Raphson算法來(lái)尋求解決方案。在每次Newton-Raphson迭代中,使用完全耦合方法中的并行稀疏直接解算器MUMPS求解線性化方程組,得到包含控制方程耦合產(chǎn)生的所有非對(duì)角項(xiàng)的系統(tǒng)雅可比矩陣。偏微分方程導(dǎo)入COMSOL的過(guò)程可參考文獻(xiàn)[16]。

        2 邊界條件與初始參數(shù)

        計(jì)算模型上表面的溫度邊界基于在研究現(xiàn)場(chǎng)收集的空氣和地溫測(cè)量數(shù)據(jù)、附面層理論[9,16-17]。使用表1列出的方程式確定計(jì)算域上表面的溫度。輸氣溫度設(shè)為8 ℃,初始地溫設(shè)為年平均地溫-0.9 ℃,下邊界恒溫層溫度設(shè)為-1 ℃,左右邊界設(shè)為絕熱邊界。不考慮水分補(bǔ)給所以水分場(chǎng)邊界設(shè)為零通量。左右邊界設(shè)為輥支撐,地表設(shè)為自由邊界,停輸再啟動(dòng)或新建管線短期內(nèi)不會(huì)有明顯位移,所以管道和底部邊界設(shè)為固定約束。水熱力耦合模型中涉及的參數(shù)[9,17-19]見(jiàn)表2。路基模型中,碎石層設(shè)為多孔介質(zhì),當(dāng)外部溫度波動(dòng)時(shí),空氣在其中自然對(duì)流??紤]到冬季的厚積雪,路基的兩個(gè)邊坡假設(shè)是不透氣。

        表1 計(jì)算域上表面的溫度邊界

        表2 參數(shù)表

        3 模型驗(yàn)證

        為驗(yàn)證本文模型與數(shù)值計(jì)算的準(zhǔn)確性與可靠性,按照西部某凍土區(qū)正溫輸氣管道周圍測(cè)溫?cái)?shù)據(jù)[10],以幾何模型中距管道豎直中心線水平距離4.5 m的鉛垂線上溫度分布為研究對(duì)象。將12月下旬、3月下旬不同深度位置的土壤自然溫度值與計(jì)算值比較。圖2即不同深度位置土壤自然溫度值與數(shù)值模擬值對(duì)比圖。從圖2中可見(jiàn)此地區(qū)土壤溫度場(chǎng)分布的計(jì)算值與實(shí)際測(cè)溫元件的實(shí)測(cè)值吻合良好,測(cè)點(diǎn)位置最大偏差不超過(guò)1.2℃,說(shuō)明本模型構(gòu)建以及參數(shù)選取合理。

        圖2 實(shí)測(cè)地溫與模擬值隨深度變化圖

        4 結(jié)果與分析

        4.1 埋深、輸運(yùn)溫度對(duì)管周凍土溫度場(chǎng)的影響

        將天然氣輸送溫度依次設(shè)為8 ℃,10 ℃;管道頂端埋深依次設(shè)為2.5 m,3.0 m進(jìn)行190 d秋冬季停輸再啟動(dòng)/新建管線的耦合場(chǎng)模擬分析。運(yùn)行44 d、89 d以及運(yùn)行187 d溫度云圖如圖3所示,管道下方土體受管道影響較大,融化范圍(融化溫度0.1℃等值線圖中用粗實(shí)線表示)逐漸擴(kuò)大,融化下限在6 m深處附近,水平融化距離始終在3 m以內(nèi)。受管道影響管道上方凍土先融化未凍水向兩側(cè)凍結(jié)鋒面遷移,兩側(cè)凍土中未凍水結(jié)冰釋放潛熱導(dǎo)致管道上方溫度場(chǎng)不是均勻分布如圖3(c)。距管道中心線4 m以上受輸氣管道影響較小,取距管道中心1 m、1.5 m、2.5 m處鉛垂線上溫度為研究對(duì)象。當(dāng)?shù)乇頊囟茸畹蜁r(shí),2.5 m埋深,不同輸運(yùn)溫度情況下該鉛垂線上土壤溫度如圖4所示。

        圖3 (a)運(yùn)行44天、(b)運(yùn)行89天、(c)運(yùn)行187天時(shí)溫度場(chǎng)

        由圖4可知當(dāng)管溫提升時(shí),管道附近1.5 m內(nèi)受影響較大有較明顯升溫,2.5 m處升溫不明顯。1 m地溫最高處高于天然場(chǎng)地5 ℃,凍土原有穩(wěn)定性受到較大影響。不同埋深、相同輸運(yùn)溫度時(shí),各處溫度分布如圖5。圖5中,當(dāng)管道埋深增加,熱影響區(qū)明顯下移,1 m地溫最高處從3 m下移至3.6 m。距管道中心線2.5 m處鉛垂深度超過(guò)3.6 m時(shí),受地表溫度影響小,埋深變化導(dǎo)致溫度上升,表明埋深增加橫向影響距離變遠(yuǎn)。

        圖4 不同輸運(yùn)溫度時(shí)土壤各處溫度分布

        圖5 不同埋深時(shí)土壤各處溫度分布

        4.2 相對(duì)飽和度S與含水率的變化

        非飽和凍土相對(duì)飽和度、含水率會(huì)隨水分遷移、凍結(jié)而變化。輸氣管道的散熱是否影響水分運(yùn)動(dòng)過(guò)程,目前實(shí)驗(yàn)研究較少,本文以模擬數(shù)據(jù)分析其變化趨勢(shì)。圖6描繪了地表溫度、距管道中心線1.5 m,0.5 m深處溫度及相對(duì)飽和度S隨時(shí)間的變化。顯然圖6中(1.5,-0.5)點(diǎn)溫度隨氣溫(地表溫度T)變化而變化,且存在20 d左右的滯后。相對(duì)飽和度與該點(diǎn)溫度變化趨勢(shì)基本一致,隨地表溫度變化而波動(dòng)。

        圖6 地表溫度、(1.5,-0.5)點(diǎn)溫度相對(duì)飽和度隨時(shí)間變化圖

        圖7為管道上方距地表0.5 m處溫度、相對(duì)飽和度隨時(shí)間圖。降溫初期相對(duì)飽和度與地溫變化趨勢(shì)有明顯滯后,這與土壤中的水分遷移有關(guān)。由于管道周圍土體融化區(qū)有較多未凍水持續(xù)向上方凍結(jié)區(qū)遷移減緩了該點(diǎn)相對(duì)飽和度的下降。同時(shí),通過(guò)對(duì)比圖6、圖7,管道上方(0,-0.5)點(diǎn)相對(duì)飽和度始終大于(1.5,-0.5)。

        圖7 (0,-0.5)點(diǎn)溫度相對(duì)飽和度隨時(shí)間變化圖

        輸氣管道上方凍土是易受管道熱影響的區(qū)域,為了進(jìn)一步分析管道上方水分遷移規(guī)律,取管道上方距地表1.25 m處含水率、該點(diǎn)溫度、地表溫度以及該點(diǎn)下方溫度為因變量,時(shí)間為自變量繪制圖8。

        圖8 (0,-1.25)點(diǎn)含水率溫度隨時(shí)間變化圖

        從圖8可以看出,輸運(yùn)初期淺層土壤中的未凍水逐漸結(jié)冰,該點(diǎn)受管道影響凍土逐漸融化,土壤中的水分由未凍結(jié)區(qū)向凍結(jié)鋒面遷移并發(fā)生聚集[20],導(dǎo)致含水量下降。管道周圍未凍水向上遷移補(bǔ)充該點(diǎn)水分,所以含水率下降緩慢。然后,該點(diǎn)水向上方更低溫處遷移受環(huán)境溫度影響逐漸凍結(jié),下方未凍水向上補(bǔ)給,所以含水率下降緩慢。下方土體中水逐漸凍結(jié),該點(diǎn)水補(bǔ)給降低,同時(shí)氣溫逐漸升高地溫最低處逐漸下降,該點(diǎn)未凍水向更低溫處遷移含水率快速下降。該點(diǎn)溫度成為地溫最低處后,含水率逐漸上升。當(dāng)?shù)乇砩琳郎貢r(shí)冰融化,向下補(bǔ)給,含水率迅速上升至趨于飽和。

        4.3 埋深、輸運(yùn)溫度對(duì)管周凍土位移場(chǎng)的影響

        對(duì)于凍土區(qū)管道工程而言,管道附近土體狀態(tài)對(duì)管道的穩(wěn)定性有較大影響。圖9描繪了地表溫度最低時(shí),不同埋深、不同管溫影響下管道上方地表與距管道中心線4.5 m處地表的位移圖。圖9中,由于工況的不同,輸送溫度低、埋深大,管道對(duì)地表自然凍脹的影響??;管道溫度為8℃埋深2.5 m時(shí),位移差可達(dá)14.5 mm。輸送溫度升高,管道上方位移減少,位移差增大,在本文模擬工況中位移差增加了1.5%。埋深增加位移差明顯減小,減小幅度達(dá)18%。研究還發(fā)現(xiàn),當(dāng)管道埋深增加,溫度場(chǎng)橫向影響距離變遠(yuǎn)(見(jiàn)圖5),遠(yuǎn)處未凍水增加,地表凍結(jié)水分向上遷移,引起遠(yuǎn)處凍脹位移增加。

        圖9 不同條件下位移對(duì)比圖

        4.4 緩解路堤的熱效應(yīng)

        由于路堤目前在實(shí)際管道工程中應(yīng)用較少,而在路堤中修建管道時(shí)地面熱狀況由環(huán)境溫度、路基、輸氣管道和凍土路基之間復(fù)雜的相互作用決定。以埋地管道模型初始條件運(yùn)行1年后熱狀態(tài)為初始狀態(tài),研究第二年的熱的相互作用。使用了第二年初冬、冬季、氣溫回升時(shí),路基土壤的模擬溫度場(chǎng)。第二年初冬氣溫約為-11℃,而管溫高達(dá)8℃。由于這兩種溫度差異較大,管堤中的管堤填土層會(huì)出現(xiàn)水平方向的溫度傳遞,管道下方的碎石層中的空氣會(huì)因受熱不均勻產(chǎn)生對(duì)流。如圖10(a)所示,在初冬時(shí),粗顆粒層與碎石層中的水平方向溫度梯度較為明顯,這說(shuō)明輸氣管道向外擴(kuò)散的熱量正被管道路堤傳遞到外部環(huán)境中。隨著冬季低氣溫的持續(xù)影響,空氣在碎石層孔隙中不斷運(yùn)動(dòng)。隨著熱量傳遞過(guò)程的進(jìn)行,水平方向溫度梯度逐漸消失,粗顆粒層下邊界的溫度與天然地表溫度基本一樣且低于淺層天然凍土路基。管道路堤下部碎石層中有較為明顯的豎直方向溫度梯度如圖10(b)。這表明,管道路堤下方較高溫度優(yōu)先向上傳遞,從而較小對(duì)凍土路基的影響。當(dāng)氣溫回升時(shí),碎石層內(nèi)觀察到明顯的溫度梯度,但此時(shí)溫度梯度與冬季時(shí)圖10(b)相反。如圖10(c)所示碎石層的上部溫度比下部高,碎石層內(nèi)上層空氣的密度比下層小,分層是穩(wěn)定的,所以碎石層中的傳熱以熱傳導(dǎo)為主。由于碎石層中各級(jí)粒徑顆粒的分配不良,熱傳導(dǎo)過(guò)程將很緩慢[21]。因此,管道路堤可以減少管道中的熱量向路基土壤中傳遞,并限制管道周圍和路基上方區(qū)域熱損失。

        圖10 (a)初冬時(shí)管基土壤溫度場(chǎng);(b)冬季管基土壤溫度場(chǎng);(c)氣溫回升時(shí)管基土壤溫度場(chǎng)

        5 結(jié)論

        本文針對(duì)埋地輸氣管道對(duì)周圍凍土的影響作用,基于PDE設(shè)計(jì)了埋地管道周圍土壤的熱水力耦合模型,并對(duì)模型進(jìn)行了驗(yàn)證。使用Comsol Multiphysics模擬了新建管線或停輸再啟動(dòng)時(shí),周圍土壤的溫度場(chǎng)、水分場(chǎng)、位移場(chǎng)的變化。為改善輸氣管道周圍土壤熱狀況,將路堤應(yīng)用在管道鋪設(shè)中,研究了緩解熱影響管堤的熱效應(yīng)。根據(jù)數(shù)值模擬結(jié)果得出以下結(jié)論:

        (1)新建管線或停輸再啟動(dòng)短期內(nèi),融化范圍逐漸擴(kuò)大,融化下限在6 m深處附近,水平融化距離始終在3 m以內(nèi),當(dāng)管溫提升時(shí),管道附近1.5 m內(nèi)受影響較大有較明顯升溫,凍土原有穩(wěn)定性受到較大影響;當(dāng)管道埋深增加,熱影響區(qū)明顯下移,橫向影響距離變遠(yuǎn)。氣溫回升短時(shí)間內(nèi),管道上方溫度場(chǎng)不均勻分布。

        (2)隨著地表溫度變化,管道附近地表下淺層土壤溫度隨之變化,約有20 d的滯后。相對(duì)飽和度與該點(diǎn)溫度變化趨勢(shì)基本一致,同一深度靠近管道處的相對(duì)飽和度大于遠(yuǎn)離管道處。管道上方含水率的變化分為四個(gè)階段:緩慢下降、下降、上升、快速上升。

        (3)輸氣管道的鋪設(shè)對(duì)變飽和凍土的自然凍脹有較大影響,管道溫度為8℃埋深2.5 m時(shí),遠(yuǎn)離管道處地表與管道上方位移差可達(dá)14.5 mm。輸送溫度升高時(shí),位移差增加幅度較少;埋深增加時(shí),位移差較少幅度較大。

        (4)氣溫下降時(shí),管堤可以有效地將輸氣管道向外擴(kuò)散的熱量傳遞到外部環(huán)境中,并在氣溫回升時(shí),減少凍土路基的熱量攝入。從而達(dá)到保護(hù)凍土原有穩(wěn)定性的目的。在管道運(yùn)行的短年限內(nèi)對(duì)凍土沒(méi)有較大影響。緩解熱影響管堤對(duì)凍土路基的熱狀態(tài)控制效果較好。

        猜你喜歡
        模型
        一半模型
        一種去中心化的域名服務(wù)本地化模型
        適用于BDS-3 PPP的隨機(jī)模型
        提煉模型 突破難點(diǎn)
        函數(shù)模型及應(yīng)用
        p150Glued在帕金森病模型中的表達(dá)及分布
        函數(shù)模型及應(yīng)用
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        3D打印中的模型分割與打包
        精品亚洲av一区二区| 亚洲av无码一区二区三区网站| 欧美国产亚洲日韩在线二区| 国产亚洲精品综合99久久 | 女人被狂c躁到高潮视频| 久久精品国产夜色| 情色视频在线观看一区二区三区| 国产亚洲成人精品久久久| 国产偷国产偷精品高清尤物| 亚洲成色在线综合网站| 亚洲av永久无码精品水牛影视| 国产免费网站在线观看不卡| 又湿又紧又大又爽a视频国产| 国产精品福利视频一区| 国产未成女年一区二区| 美女被黑人巨大入侵的的视频| 国产精品爽爽ⅴa在线观看| 国产精品麻豆aⅴ人妻| 中国产无码一区二区三区| 亚洲中文字幕精品久久吃奶| 亚洲精品成人片在线观看精品字幕| 亚洲av成人一区二区三区av| 性色av成人精品久久| 色熟妇人妻久久中文字幕| 精品日产卡一卡二卡国色天香| 亚洲 国产 哟| 亚洲国产女同在线观看| 无码aⅴ精品一区二区三区浪潮 | 亚洲国产精品午夜一区| 国产一区二区三区中文在线| 久久精品娱乐亚洲领先| 日本视频一区二区三区免费观看| 九一精品少妇一区二区三区| 国产无夜激无码av毛片| 亚洲国产欧美在线成人| 蜜桃av一区在线观看| 国产内射爽爽大片| 亚洲永久无码7777kkk| 国产码欧美日韩高清综合一区| 亚洲国产天堂久久综合网| 亚洲av无码一区二区乱孑伦as|