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

        ?

        縫洞型介質(zhì)流動模擬的多尺度分解法*

        2019-04-10 02:45:38張慶福黃朝琴姚軍李陽2嚴(yán)俠
        物理學(xué)報 2019年6期
        關(guān)鍵詞:縫洞溶洞尺度

        張慶福 黃朝琴? 姚軍 李陽2) 嚴(yán)俠

        1) (中國石油大學(xué)(華東)石油工程學(xué)院,青島 266580)

        2) (中石化油田勘探開發(fā)部,北京 100728)

        (2018 年8 月31 日收到; 2019 年1 月21 日收到修改稿)

        縫洞型介質(zhì)通常具有非均質(zhì)性強(qiáng)、結(jié)構(gòu)多尺度的特征. 傳統(tǒng)數(shù)值方法在解決此類多尺度流動問題時,難以兼顧計算精度與計算效率,無法實際應(yīng)用. 對此,本文提出了多孔介質(zhì)流體流動的多尺度分解法,并應(yīng)用于縫洞介質(zhì)流動模擬,能夠大幅減小計算的復(fù)雜度,同時可以通過控制均化程度控制計算精度. 該方法將求解空間分為若干個子空間的正交直和,從而獲得一個近線性的計算復(fù)雜度; 以分層計算的方式實現(xiàn)了快速計算,另外這種方法是一種無網(wǎng)格方法,具有較好的地層適應(yīng)性. 同時,采用離散縫洞模型簡化縫洞結(jié)構(gòu),進(jìn)一步提高了計算效率. 詳細(xì)闡述了基于多尺度分解法的多孔介質(zhì)流體流動數(shù)值計算格式的建立,重點介紹了如何在不同的層次上計算基函數(shù). 數(shù)值結(jié)果表明,本文提出的計算方法不僅能夠準(zhǔn)確捕捉多孔介質(zhì)中的精細(xì)流動特征,而且具有很高的計算效率,是一種有效的流動模擬方法.

        1 引 言

        縫洞型介質(zhì)廣泛存在于自然界中,過去幾年間,針對其開展的數(shù)值模擬技術(shù)引起了廣泛關(guān)注,但是縫洞型介質(zhì)流動模擬給常規(guī)數(shù)值方法帶來很大的挑戰(zhàn). 縫洞型介質(zhì)除了本身具有的非均質(zhì)性,其中所包含的裂縫和溶洞橫跨多個尺度. 如果要描述縫洞介質(zhì)所有尺度上的流動特征,則需要很精細(xì)的網(wǎng)格劃分,使用常規(guī)數(shù)值方法求解需要巨大的計算量,即便使用超級計算機(jī)也難以滿足需求. 因此,常規(guī)數(shù)值方法難以描述流體所有尺度上的流動特征. 尺度升級方法被廣泛視為能減小計算量的方法[1,2],但是由于其忽視了小尺度信息導(dǎo)致不能反映精細(xì)的介質(zhì)非均質(zhì)性. 近年來,多尺度方法被應(yīng)用于流動模擬[3?7],這種方法能夠在粗網(wǎng)格上通過計算多尺度基函數(shù)捕捉小尺度特征,從而在保證計算精度的同時減少計算量. 然而,當(dāng)粗網(wǎng)格和細(xì)網(wǎng)格尺度相差較大時,多尺度方法在差值算子及約束算子的構(gòu)造方面存在著困難. 多重網(wǎng)格方法是公認(rèn)的求解偏微分方程的最快的方法之一[8?11],因此,本文首先將多重網(wǎng)格法應(yīng)用于非均質(zhì)介質(zhì)流動問題,并進(jìn)一步構(gòu)建了一種新的縫洞型介質(zhì)流動模擬的多重網(wǎng)格法(多尺度分解法).

        雖然多重網(wǎng)格方法被拓展至求解多種偏微分方程,但是當(dāng)方程系數(shù)粗糙時,會導(dǎo)致該方法收斂性差[12,13]. 因此,對于強(qiáng)非均質(zhì)性油藏,傳統(tǒng)多重網(wǎng)格法也存在收斂性差的問題. 近年來,雖然改進(jìn)的多重網(wǎng)格方法在一定程度上能夠緩和粗糙系數(shù)對收斂性的影響[14?17],但能夠精確適應(yīng)粗糙系數(shù)的多重網(wǎng)格法一直都是數(shù)學(xué)界的難題[18].

        目前,主要有三種方法用于縫洞型碳酸鹽巖數(shù)值模擬: 1)等效介質(zhì)模型[19]; 2)多重介質(zhì)模型[20];3)離散縫洞網(wǎng)絡(luò)(DFVN)模型[21,22]. DFVN 模型能精細(xì)描述縫洞,基巖和裂縫系統(tǒng)視為滲流區(qū)域,溶洞系統(tǒng)視為自由流區(qū)域,是離散裂縫模型的有效延伸. DFVN 模型屬于宏觀尺度的精細(xì)模型,能克服三重介質(zhì)模型和等效介質(zhì)模型的瓶頸,但受目前計算機(jī)硬件和數(shù)值模擬技術(shù)的制約,只能處理小規(guī)模的縫洞型油藏. 針對此問題,對于小尺度裂縫和溶洞,對其進(jìn)行等效處理,用等效后的滲透率場描述其造成的非均質(zhì)性; 使用離散裂縫模型描述大尺度裂縫,準(zhǔn)確刻畫裂縫流動狀態(tài).

        本文基于多尺度分解思想[23,24],構(gòu)建了縫洞型介質(zhì)多重網(wǎng)格計算格式,能夠準(zhǔn)確描述介質(zhì)的非均質(zhì)性,并能準(zhǔn)確模擬大尺度裂縫存在時的壓力差分布. 文中首先介紹了離散縫洞的數(shù)學(xué)模型和多尺度分解法的基本原理,然后構(gòu)建了縫洞介質(zhì)的多尺度分解計算格式,并給出了幾種不同的算例. 數(shù)值結(jié)果表明,多尺度分解法可以在充分反映裂縫溶洞影響的同時大幅減少計算復(fù)雜度,是一種潛力很大的流動模擬數(shù)值方法.

        2 縫洞混合流動數(shù)學(xué)模型

        2.1 縫洞介質(zhì)數(shù)學(xué)模型

        本文通過計算等效滲透率來描述介質(zhì)中的小尺度裂縫和溶洞導(dǎo)致的非均質(zhì)性; 采用離散裂縫模型表征大尺度裂縫,兩者結(jié)合形成多尺度縫洞混合模型(圖1).

        圖1 縫洞型介質(zhì)示意圖Fig. 1. Schematic of fractured-vuggy porous medium.

        將基巖和裂縫系統(tǒng)視為滲流區(qū)域,其運動方程符合Darcy 定律,即

        式中μ為流體黏度(mPa·s);Ki為滲透率( μ m2);vi為滲流速度(m/s); 對于巖塊l=m,對于裂縫l=f;p為多孔介質(zhì)中的平均壓力(Pa);ρ是流體密度(kg/m3);f為單位質(zhì)量力(m/s2).

        溶洞的空間尺度較大,其流動模式視為自由流,為簡單起見本文僅考慮牛頓流體的低雷諾數(shù)運動. 為了避免引入復(fù)雜的耦合邊界條件,降低自由流-滲流耦合問題的復(fù)雜性,采用Brinkman 方程描述溶洞內(nèi)的自由流:

        式中vv是流體真實速度(m/s);Kv是溶洞內(nèi)滲透率( μ m2);?為溶洞內(nèi)孔隙度(無充填時為1);pv為流體壓力(Pa).

        假設(shè)流體流動過程溫度不變,綜合以上方程可以推導(dǎo)出DFVN 模型的控制方程組:

        式中

        v為速度場. 結(jié)合合適的邊界條件可得到完整的離散縫洞網(wǎng)絡(luò)單相穩(wěn)態(tài)滲流數(shù)學(xué)模型.

        2.2 等效滲透率計算

        基于流量等效計算等效滲透率的方法,往往僅考慮目標(biāo)單元內(nèi)的非均質(zhì)性,這樣得到的等效滲透會存在誤差. 對此,本文采用超樣本方法,在包含目標(biāo)單元的更大區(qū)域內(nèi)建立離散縫洞網(wǎng)絡(luò)模型并進(jìn)行有限元數(shù)值求解,從而獲取超樣本單元內(nèi)的壓力和速度場; 最后,基于體積平均法計算目標(biāo)網(wǎng)格內(nèi)的等效滲透率.

        在超樣本單元內(nèi)結(jié)合封閉定壓邊界條件建立離散縫洞網(wǎng)絡(luò)模型,采用有限單元法進(jìn)行求解[21],根據(jù)求得的壓力及速度場,采用體積平均法計算目標(biāo)網(wǎng)格上的速度和壓力梯度的體積平均值[25]:

        式中j=x,y(表示定壓邊界施加的坐標(biāo)軸方向);Vb表示目標(biāo)原始網(wǎng)格的體積;ul和 (?p)l分別表示目標(biāo)原始網(wǎng)格中第l個單元的速度和壓力梯度;Vl表示第l個單元的體積;N表示目標(biāo)原始網(wǎng)格中單元的總數(shù). 在計算得到后,根據(jù)達(dá)西定律,可以求得等效滲透率,詳細(xì)計算過程可以參考文獻(xiàn)[26].

        3 多尺度分解法

        多尺度分解法的目的是盡可能快的、以近線性的計算復(fù)雜度求解偏微分方程. 令q為網(wǎng)格層數(shù),定義I(q)為每一層的細(xì)網(wǎng)格單元索引i= (i1,· · · ,iq) .令i(k)=(i1,··· ,ik),I(k)={i(k):i ∈I(q)} ,其 中1 ≤k≤q,i=(i1,· · · ,iq)∈I(q);i(k,k′) 為滿足j ∈I(k′)和j(k)=i的 單 元,其 中k′ ∈{k,· · · ,q}.為I(k)×I(k+1)矩陣,滿足當(dāng)j∈/i(k,k+1)時同 時其 中I(k)為I(k)×I(k)的單位矩陣.

        不同網(wǎng)格層上的測試函數(shù)定義為

        如 圖2 所 示,?被 分 為 若 干 2?k×2?k的 子 單 元其 中為的體積,為的指示函數(shù). 對于邊界條件問題,使用非零邊界的測試函數(shù). 通過假設(shè)A和B在進(jìn)行一個博弈游戲來定義基函數(shù),具體理論證明見文獻(xiàn)[24]. 根據(jù)博弈理論[27,28]并結(jié)合流動方程,B的最佳選擇為

        圖2 區(qū)域?的網(wǎng)格剖分示意圖Fig. 2. Schematic of grid partition of solution space.

        其中為基函數(shù),定義為

        在實際計算中,可以通過求解以下方程獲得

        可以證明以下定理成立:

        1) 最優(yōu)問題(9)式可以確定唯一一個由方程(8)定義的ψi;

        以均勻介質(zhì)為例,首次給出了流動問題的多重網(wǎng)格法的基函數(shù),如圖3 所示,每一層的基函數(shù)的共同特點是這些基函數(shù)都是局部化的. 本文沒有對基函數(shù)進(jìn)行截取,當(dāng)基函數(shù)包含所有信息時可視為精確的基函數(shù),單項流動中也存在其他高精度基函數(shù)[29]. 局部化的思想是目前的熱點問題之一,該思想有利于信息的集中處理,目前大數(shù)據(jù)等研究中均有涉及.

        定義空間

        其中k∈{1,··· ,q}. 從(6)式可以推出此空間是嵌套的. 因此,對于k∈{1,·· · ,q ?1},

        圖3 分層基函數(shù)示意圖Fig. 3. Illustration of gamblets .

        令A(yù)(k)是剛度矩陣那么至此,以上嵌套公式足以進(jìn)行多重網(wǎng)格計算,但是本文進(jìn)一步對其進(jìn)行多尺度分解. 定義J(k)為k-元組,其中令W(k)為J(k)×I(k)矩 陣,滿 足 當(dāng)j(k?1)=i(k?1)時對于令

        其中i∈J(k),定義空間

        III(k)為II(k?1)關(guān)于II(k)的正交補集. 令⊕λ表示正交直和,那么可以得到空間的多尺度分解為

        進(jìn)一步可以推導(dǎo)出

        此為壓力p的正交分解,p(k+1)?p(k)為流動方程在空間III(k+1)里的解. 流動方程的及多尺度分解如圖4 和圖5 所示,可以理解為在每一層上的基函數(shù).

        當(dāng)介質(zhì)內(nèi)存在大尺度裂縫時,則采用離散裂縫模型表征大裂縫,即對裂縫進(jìn)行降維處理[8]. 假設(shè)裂縫介質(zhì)流動方程為FEQ,離散裂縫模型的計算格式為

        其中di為裂縫開度,?m為基巖區(qū)域,?f為裂縫區(qū)域. 在每一個子空間內(nèi)求解離散裂縫模型,捕捉裂縫對滲流的影響,然后將各個子空間的解加起來獲取最終縫洞介質(zhì)數(shù)值解,即

        4 數(shù)值算例

        本節(jié)首先通過一個數(shù)值算例與參考算例的對比驗證了程序的正確性,把使用傳統(tǒng)模擬有限差分法求得的數(shù)值解作為參考解; 然后通過復(fù)雜離散裂縫模型數(shù)值算例進(jìn)一步驗證了方法的正確性和程序的魯棒性. 對于單相流動而言,著重對比壓力差分布,為了實現(xiàn)壓力差的對比,設(shè)介質(zhì)中一點的壓力為已知,令壓力分布為正值.

        4.1 小尺度縫洞模型算例

        考慮圖6(a)所示的縫洞型介質(zhì)物理模型,模型大小為100 m × 100 m. 模型包含小裂縫和溶洞,不存在大裂縫,模型中的參數(shù)如表1 所列. 溶洞簡化為規(guī)則球形,溶洞內(nèi)不考慮充填(滲透率無充填時為無窮大),其中的流動視為自由流. 圖6(b)是通過均化理論求得的等效滲透率場圖. 令流體自上而下流動,邊界為不流動邊界.

        圖4 基函數(shù)示意圖Fig. 4. Illustration of basis functions

        圖5 區(qū)域?的多尺度解示意圖Fig. 5. Illustration of solution of the multiple spaces.

        令q=4,即網(wǎng)格系統(tǒng)包含4 層網(wǎng)格,如圖6(c)所示. 采用多尺度分解法求解得到每一層上的解,并根據(jù)(17)式獲得最終的壓力分布. 小尺度裂縫與溶洞擁有較高的滲透率,會造成多孔介質(zhì)的非均質(zhì)性,進(jìn)而影響整體壓力分布. 在數(shù)值算例中采用L2范數(shù)計算壓力誤差

        表1 裂縫性介質(zhì)模型基本參數(shù)Table 1. Parameters of fractured porous medium.

        式中pf為參考壓力解,pmg為多尺度分解法所計算的壓力解.

        圖7 給出了參考解和多重網(wǎng)格解的對比. 可以看出,當(dāng)k=3 時多重網(wǎng)格法所得壓力分布與參考解幾乎一致,因此多重網(wǎng)格解能夠精確地反映非均質(zhì)性,從而捕捉到小尺度裂縫和溶洞的影響. 表2 列出了不同k時的計算誤差,其中當(dāng)k=4時為多層分解法的精確解,當(dāng)k=3 和k=2 時為均化計算,計算誤差均很小,體現(xiàn)了該方法的精確性. 小尺度縫洞模型的計算結(jié)果表明,此方法可以有效地進(jìn)行非均質(zhì)地層以及小尺度縫洞介質(zhì)的流動模擬.

        表2 對于小尺度縫洞模型,不同k時的計算誤差Table 2. Relative error in differentkfor a smallscale-fractured vuggy porous medium.

        圖6 (a)小尺度縫洞型油藏幾何模型; (b)等效滲透率場圖; (c) 多層網(wǎng)格系統(tǒng)Fig. 6. (a) Geometrical model of a small-scale-fractured vuggy porous oil reservoir; (b) equivalent permeability of fractured-vuggy medium; (c) nested grid system.

        4.2 長裂縫模型算例

        算例2 考慮存在長裂縫的情況. 考慮如圖8 所示的裂縫型介質(zhì),在小裂縫和溶洞的背景下包含一條長裂縫. 裂縫開度df=1 × 10–3m,裂縫滲透率為/12 μ m2,基巖孔隙度?=0.2 . 邊界均為不流動邊界,流體從上往下流動. 采用離散裂縫模型表征裂縫.

        本算例考慮長裂縫與小裂縫溶洞共存的情況.壓力分布圖9 展示了不同k時多尺度分解法與參考解的對比. 當(dāng)k=3 時計算誤差為0.0131,當(dāng)k=2 時計算誤差為0.0741,在進(jìn)行均化的前提下依然能夠反映裂縫的存在,保證較高的計算精度.參考解與多重網(wǎng)格解的對比結(jié)果說明,本文構(gòu)建的多尺度分解法能夠有效地處理長裂縫的情況,并能通過均化大幅減少計算量,同時保持較高的計算精度.

        圖7 對于小尺度縫洞模型,參考解和多重網(wǎng)格解對比 (a) 參考解; (b)k=3 時的多重網(wǎng)格解Fig. 7. Comparison of reference solution and gamblets solution for a small-scale-fractured vuggy porous medium: (a) Reference solution; (b) gamblets solution withk=3.

        圖8 長裂縫介質(zhì)模型Fig. 8. Geometrical model of a fractured-vuggy porous medium with a long fracture.

        4.3 大尺度縫洞模型算例

        自然介質(zhì)中,除了小裂縫,往往伴隨著壓裂等增產(chǎn)措施產(chǎn)生大裂縫,該算例檢驗本文提出的方法對大裂縫的模擬能力. 圖10 所示為10 m × 10 m縫洞型介質(zhì),其中包含的小尺度裂縫和溶洞的參數(shù)與算例1 相同. 在包含小尺度溶洞的同時,還包含著尺度較大的裂縫網(wǎng)絡(luò)系統(tǒng). 裂縫網(wǎng)絡(luò)由6 條長裂縫組成,長裂縫開度df=1 × 10–3m,裂縫滲透率為/12 μ m2,基巖孔隙度?=0.2 . 流體從上至下流動.

        多孔介質(zhì)中的長裂縫作為導(dǎo)流通道,對壓力分布有顯著的影響,并且裂縫之間也會相互影響. 該方法在每一個子空間內(nèi)求解離散裂縫模型,圖11為離散裂縫模型在各層上的解,可以很明顯地看出,各個層次上的解都能反映出裂縫的存在. 然后將k個層次上的解加起來,即得到最終的解.

        圖9 對于長裂縫模型,參考解和多重網(wǎng)格解對比 (a)參考解; (b)k=3 時的多重網(wǎng)絡(luò)解; (c)k=2 時的多重網(wǎng)格解Fig. 9. Comparison of reference solution and gamblets solution for a fractured-vuggy porous medium with a long fracture: (a) Reference solution; (b) gamblets solution withk=3; (c) gamblets solution withk=2.

        圖10 大尺度縫洞介質(zhì)幾何模型Fig. 10. Geometrical model of a large-scale-fractured vuggy porous medium.

        圖12 給出了參考解和多重網(wǎng)格解的對比,可以看出,本文構(gòu)建的數(shù)值算法在捕捉小尺度裂縫和溶洞的同時,能夠精確模擬長裂縫對壓力場的影響,并反映裂縫間的相互作用. 結(jié)合表3 可知,當(dāng)均化程度較大(k=1)時存在較大誤差. 圖13 給出了沿x=5 m 的參考解與不同k時多重網(wǎng)格解的對比. 裂縫在流動過程中可視為一個等勢體,對應(yīng)于圖11 中較平緩的部分,與實際符合,隨著均化程度加大誤差變大. 該算例展示了本文構(gòu)建的數(shù)值方法對長裂縫的精確模擬能力.

        實際上,對于任意形狀的一個裂縫性介質(zhì),一旦求出剛度矩陣及對應(yīng)的W和 π ,就可以有效地進(jìn)行任意形狀的裂縫性介質(zhì)數(shù)值模擬. 因此,該方法能夠處理形狀復(fù)雜的裂縫系統(tǒng). 但是,對于裂縫分布極其復(fù)雜的裂縫系統(tǒng),尤其裂縫之間距離非常小的情況,多層網(wǎng)格剖分會存在一定困難,造成分層的基函數(shù)不適應(yīng)高滲透的裂縫,進(jìn)而影響該方法的計算精度.

        圖11 大尺度縫洞介質(zhì)模型在各層上的解Fig. 11. Solutions in different levels for a large-scale-fractured vuggy porous medium.

        圖12 對于大尺度縫洞介質(zhì)模型,參考解和多重網(wǎng)格解對比 (a) 參考解; (b)k=3 時的多重網(wǎng)格解; (c)k=2 時的多重網(wǎng)格解;(d)k=1 時的多重網(wǎng)格解Fig. 12. Comparison of reference solution and gamblets solution for a large-scale-fractured vuggy porous medium: (a) Reference solution; (b) gamblets solution withk=3; (c) gamblets solution withk=2; (d) gamblets solution withk=1.

        表3 對于大尺度縫洞介質(zhì)模型,不同k時的計算誤差Table 3. Relative error with differentkfor a large-scalefractured vuggy porous medium.

        圖13 沿x=5 m 的參考解和不同均化程度下多重網(wǎng)格解對比曲線圖Fig. 13. Comparison of reference solution and gamblets solution with differentkwhenx=5 m.

        5 結(jié) 論

        1)隨著現(xiàn)代地質(zhì)建模技術(shù)的發(fā)展,縫洞型地質(zhì)模型可能包含數(shù)百萬甚至數(shù)億個網(wǎng)格,采用傳統(tǒng)的數(shù)值方法對其進(jìn)行求解,計算量巨大,甚至超出當(dāng)今計算機(jī)的計算能力. 本文首次提出縫洞介質(zhì)的多尺度分解計算模式,構(gòu)建局部化的基函數(shù),該方法能夠快速求解流動方程,非常適合對強(qiáng)非均質(zhì)型油藏進(jìn)行精細(xì)流動模擬;

        2)結(jié)合離散裂縫模型和等效介質(zhì)模型,既可以精確反映大尺度裂縫的流體流動,又可以捕捉小尺度溶洞和裂縫的影響,從而可以對縫洞型介質(zhì)進(jìn)行高效流動模擬;

        3)每一層網(wǎng)格的基函數(shù)是獨立計算的,可以采用并行計算得到,進(jìn)一步減少了計算量. 因此,本文提出的方法對于縫洞型油藏數(shù)值模擬有很高的潛在價值. 滲透率場對計算精度有一定影響,因此提高算法的魯棒性是研究重點. 下一步可采用壓力梯度對比的形式驗證方法對多相流動計算的精確性[30].

        猜你喜歡
        縫洞溶洞尺度
        別有洞天
        碳酸鹽巖縫洞儲集體分尺度量化表征
        財產(chǎn)的五大尺度和五重應(yīng)對
        出發(fā)吧,去溶洞
        哈拉哈塘奧陶系縫洞型成巖圈閉及其成因
        神秘的溶洞
        幼兒100(2017年31期)2017-11-27 02:37:45
        宇宙的尺度
        太空探索(2016年5期)2016-07-12 15:17:55
        縫洞型介質(zhì)結(jié)構(gòu)對非混相氣驅(qū)油采收率的影響
        9
        縫洞型碳酸鹽巖油藏數(shù)值模擬技術(shù)與應(yīng)用
        少妇人妻出水中文字幕乱码| 亚洲国产av导航第一福利网| 欧美激情二区| 亚洲精品无人区一区二区三区| 人妻少妇被粗大爽视频| 蜜臀亚洲av无码精品国产午夜.| 亚洲国产精品久久亚洲精品| 亚洲五月婷婷久久综合| av天堂中文亚洲官网| 亚洲国产精品久久艾草| 亚洲 自拍 另类 欧美 综合| 亚洲精品自拍视频在线观看 | 国产精品成人va| 激情五月天俺也去综合网| 国产成人高清在线观看视频| 国产三区在线成人av| 亚洲天堂手机在线| 在线观看人成网站深夜免费| 国产高清成人在线观看视频| av无码天堂一区二区三区| 日韩欧美国产亚洲中文| 久久久精品国产老熟女| 女人被男人爽到呻吟的视频| 97人妻熟女成人免费视频| 亚洲色图视频在线播放| 日本高清一区二区三区在线观看| 精品国产一二三产品区别在哪 | 免费无码又爽又刺激网站直播| 亚洲男人天堂| 亚洲精品国产熟女久久| 极品av一区二区三区| 成人免费看www网址入口| 人妻被猛烈进入中文字幕| 亚洲精品国产综合久久| 成人免费直播| 国产一区a| 亚洲av高清一区三区三区| 91蜜桃精品一区二区三区毛片| av剧情演绎福利对白| 国产色秀视频在线播放| 亚洲香蕉毛片久久网站老妇人|