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

基于雙層非穩態導熱過程的井筒溫度場半解析模型

2016-04-06 03:51:48趙金洲彭瑀李勇明田植升符東宇
天然氣工業 2016年1期
關鍵詞:模型

趙金洲彭 瑀李勇明田植升符東宇

1.“油氣藏地質及開發工程”國家重點實驗室?西南石油大學 2.中國石化河南石油工程有限公司

趙金洲等.基于雙層非穩態導熱過程的井筒溫度場半解析模型.天然氣工業,2016,36(1):68-75.

?

基于雙層非穩態導熱過程的井筒溫度場半解析模型

趙金洲1彭 瑀1李勇明1田植升2符東宇1

1.“油氣藏地質及開發工程”國家重點實驗室?西南石油大學 2.中國石化河南石油工程有限公司

趙金洲等.基于雙層非穩態導熱過程的井筒溫度場半解析模型.天然氣工業,2016,36(1):68-75.

摘 要在油氣井酸化、壓裂增產施工中,溫度控制著入井工作液的流變性和酸巖反應速度,是影響裂縫幾何形態、酸蝕導流能力和措施后產能的重要因素。但目前的井筒溫度場模型還缺乏對非穩態過程的嚴密推導和精確描述。為此,首先在物性分析的基礎上,建立了軸向上離散、徑向上解析的雙層非穩態導熱井筒溫度場半解析模型,并通過拉普拉斯變換和Stehfest數值反演方法進行了求解;其次,分析了不同水泥環導熱過程處理方法對井筒溫度分布的影響,認為假設水泥環具有穩態熱阻,模擬輸出的誤差過大且難以修正;最后,對地層和施工參數進行了敏感性分析,認為巖性、施工排量和注入溫度變化對井筒溫度分布具有顯著的影響,而油管腐蝕的影響則可以忽略?,F場實際應用效果表明,該數學模型精度較高,模擬結果與監測結果吻合良好。該成果對于提高單井酸化壓裂成功率、增加措施后的產能和油氣田采收率具有重要意義。

關鍵詞井筒 溫度場 非穩態導熱 酸化 壓裂(巖石) 拉普拉斯變換 數學模型

在油氣井酸化、壓裂增產施工中,溫度控制著入井工作液的流變性[1-3]、濾失能力[4]和酸巖反應速度,是影響裂縫幾何形態、單井增產效果和措施后產能的重要因素[5]。自1964年,Ramey[6]引入Van Everdinger等[7]推導的時間函數建立了第一個非穩態井筒溫度場模型以來,對井筒非穩態傳熱的研究就進入了蓬勃發展期。其中最具有代表性的有:1970年,Eickmeier等[8]采用傳熱過程的分析方法建立了非穩態的井筒溫度場差分模型,但在該非穩態模型中油套環空、油套管壁和水泥環都是單層的、熱阻也是恒定的,因此該模型不能體現近井位置的非穩態導熱特征。1990年,勞倫斯伯克利實驗室的吳玉樹等[9]通過非穩態導熱微分方程和拉普拉斯變換建立了井筒溫度場的解析模型,但在該模型中將從油管到水泥環的傳熱考慮為穩態過程,并且沒有對這一假設導致的誤差進行分析。20世紀90年代,Hasan 和Kabir等[10-13]又在SPE上發表了大量有關于井筒溫度場的研究成果,形成了其獨有的溫度場模擬體系和經驗處理方法,對井筒溫度場的發展做出了巨大貢獻。但他們在地層非穩態導熱模型中引入了松弛距離的概念,其本質就是Ramey模型中的時間函數,因此與Ramey模型一樣,存在不能給出早期精確解的問題。2000年以后,Hagoort、Mcspadden和Yoshida等又分別研究了Ramey模型的理論依據[14]、單井模型在井工廠模式下的應用[15]以及井筒溫度場和酸化油藏溫度場的耦合[16],對原有的經典模型進行了擴展。

井筒傳熱是一個比例極不協調的對流換熱、非穩態導熱耦合問題。在軸向上動輒千米,采用簡單的解析模型不能得到精確結果;在徑向上,層數多且厚薄不一,如果要體現每一層的非穩態導熱特征則會導致計算量過大,甚至結果不收斂,降低計算效率。因此,筆者在物性分析的基礎上,建立了軸向上離散,徑向上解析,基于雙層非穩態導熱的井筒溫度場半解析模型,并進行了參數分析。

1 模型假設條件

根據已有的研究和目前的工藝現狀,需要對井筒溫度分布的物理模型做如下假設(圖1):①酸化、壓裂介質均為液體,不考慮焦湯效應,流體不可壓縮,機械能不向內能轉化;②在軸向上按照巖性和精度要求劃分層段,在徑向上僅有油、套管壁具有穩態熱阻;③酸液在油管中流動時,與油管壁在緩蝕劑作用下微量反應;④地層、管串和環空是均質、各向同性的,并且關于井軸對稱;⑤油管內的強迫對流換熱和油套環空腔體內的自然對流換熱物性參數由Hermite插值獲得,定性溫度為平均溫度;⑥忽略徑向上各層間的界面熱阻。

圖1 井筒溫度場物理模型和物性參數圖

導溫系數是衡量溫度在某種介質中傳播速度的指標,從表1可得低碳鋼的導溫系數最高,而水泥環的導溫系數最低。溫度在低碳鋼中的傳播速度分別是在砂巖、白云巖、石灰巖、頁巖、水泥環中傳播速度的4.5~20倍、8.0~15.0倍、10.7~17.0倍、15.1~17.4倍、26.7~57.9倍。因此,油套管內外壁面的溫差不大,可以考慮這兩層具有穩態熱阻;而水泥環和地層中的導熱則是明顯的非穩態過程。

表1 井筒溫度模擬需要的材料熱物性參數

2 模型建立

井筒溫度場可以分為1區、2區和過渡區,如圖2所示。r1、r2、r3、r4、r5分別為油管內壁、油管外壁、套管內壁、套管外壁和井壁到油管中心的距離。在酸壓過程中1區為強迫對流換熱,2區為非穩態導熱,過渡區從內到外依次為油管壁導熱、環空的腔體自然對流換熱和套管壁導熱。

圖2 徑向上的傳熱過程圖

2.1 管內熱平衡方程

通過軸向上的微元分析,可得管內的熱平衡方程式:

上式中q1,j是徑向上的傳熱,表示為:

式(2)中的Rj/rx是距井軸rx遠處的綜合傳熱系數,其中Rj為:

式(3)中的α1,j和αe,j分別為油管內的強迫對流換熱系數和環空內的腔體自然對流換熱系數。根據對流換熱理論[17-19],前者可用式(4)計算,而后者可以參考圖版的取值(圖3)。即

圖3 腔體自然對流換熱系數圖版

從圖3可得,當H/L大于20時,腔體自然對流換熱系數與尺寸比例無關,努賽爾數的表達式退化為:

2.2 地層非穩態導熱

水泥環和地層中的非穩態導熱微分方程為:

初始狀態為熱平衡條件,各段井筒及地層溫度滿足地溫梯度線性函數,即

油管的上邊界條件等于注入溫度或者上一段的計算值:

非穩態導熱的內邊界條件滿足牛頓冷卻定律,外邊界溫度滿足地溫梯度函數,水泥環和地層的界面滿足第一類和第二類邊界條件:

由式(6)~(9)可以構成地層中第j段關于溫度的定解問題。

3 模型求解

式(6)~(9)給出的定解問題是一個初邊值問題,可在進行拉氏變換后,將原問題化簡為常微分方程求解。為了得到0階的貝塞爾方程,以及簡化推導過程,采用如式(10)所示的方案進行初邊值問題的無因次化。

將無因次化后的方程進行拉氏變換:

可以得到拉氏空間中的定解問題:

其中

將貝塞爾方程通解帶入邊界條件[20-22],求解關于待定系數的線性方程組,可分別得到水泥環和地層中的非穩態導熱解析解:

其中

油管內溫度的解析解可以通過求解常微分方程得到,當j=1時為:

當j>1時為:

將求解的管內溫度代入到式(13)、(14),通過Stehfest數值反演[23-25]方法,即式(17)、(18),就可以得到原空間中水泥環和地層中的溫度分布。

4 實例驗證與參數分析

4.1 實例驗證

根據某氣井的地層和施工參數進行了井底溫度模擬,并與實際的井溫監測數據進行了比較,結果如圖4所示。

圖4 模擬溫度與實測溫度對比圖

從圖4可以看出,井筒溫度場半解析模型模擬出的井底溫度與實測井底溫度曲線具有較高的吻合度,特別是在30 min之后兩條曲線基本重合。該模型的精度較高,可以用來解決一些與酸化壓裂井溫計算相關的模擬與優化問題。

4.2 水泥環的處理方法

國內外對水泥環傳熱的處理方法主要有3種。

1)考慮水泥環的熱物性與地層接近,把水泥環視為地層的一部分,將式(15)、(16)中的β5替換為如下形式即可得到該假設條件下的解。

2)考慮水泥環具有穩態熱阻,將綜合傳熱系數的倒數改為如下形式:

3)本模型中,將水泥環也作為非穩態導熱過程考慮。3種假設條件下模擬輸出的曲線分別如圖5、6、7所示。

圖5 忽略水泥環的模擬結果圖

圖6 水泥環穩態導熱的模擬結果圖

圖7 水泥環非穩態導熱的模擬結果圖

由圖5~7可得到假設情況1)和2)的模擬結果比較接近,對圖5的結果進行平移修正,就基本可以滿足精度要求,因此將水泥環歸為地層的一部分比較合理。假設情況2)得到的模擬結果精度較差,井底溫度的變化范圍很窄、修正困難,因此不能認為水泥環具有穩態熱阻。

4.3 厚層分段

圖8中黑色實線和紅色密短線代表將整個厚層作為單段考慮時的模擬結果,而藍色點線和綠色疏短線代表將整個厚層均分為4段的模擬結果。圖8的曲線顯示,當模擬的層厚超過2 500 m時,采用單段和多段的模擬結果就會產生一定偏差。因此在巖性較為均一的巨厚層需要對其進行網格加密,增加模擬精度。另外,隨著注液時間的增加,厚層分段與否對模擬結果的影響逐漸變小。

圖8 單段與多段模擬對比圖

4.4 巖性分段

在井深方向上的巖性差異,也是進行垂向分段的重要依據之一。圖9模擬的施工井在0~3 796 m砂巖比重高,而井深3 796~6 895 m石灰巖是主要層段。

圖9 巖性分層對井筒溫度分布的影響圖

圖9的曲線模擬了巖性變化對井筒溫度的影響。無論是在注液0.5 h還是4 h后,考慮巖性變化進行分段都會在3 796 m處產生明顯拐點。是否考慮井深方向的巖性變化對模擬出的溫度分布具有重大影響,并且該影響隨著施工時間的增加而增大。在具有明顯巖性差異的區域應該考慮在垂向上分段模擬。

4.5 施工排量

圖10指出排量對井筒溫度的分布具有顯著影響,排量提高三分之二,井底溫度會下降近20 ℃,但若繼續提高排量,其降溫效果會變差。此外,隨著施工時間的推移,排量對井筒溫度的影響有減小的趨勢。圖10中還可以看出;排量很低時,在巖性界面會有一明顯拐點;隨著排量增加,拐點逐漸消失,證明在高排量下巖性變化的影響不如在低排量下顯著。

圖10 施工排量對井筒溫度分布的影響圖

4.6 注入溫度

在需要泵送隔離劑、堵劑等試劑的特殊類型施工中,入井液的溫度可能會高于正常值。因此分析了注入溫度對井筒溫度分布的影響。由圖11可知:隨著注入溫度的升高,井筒的加熱效應逐漸減弱;甚至在注入高溫液體時會在上半段產生明顯的冷卻、降溫作用。在正常施工條件下,井口注入的溫度越高,井筒的溫度分布越均勻。

圖11 注入溫度對井筒溫度場的影響圖

4.7 油管腐蝕

酸壓施工中,油管內壁與酸液會在緩蝕劑的作用下微量反應,因此在推導過程中引入了內熱源項以考察其影響。圖12中的曲線是模擬腐蝕速度分別為15、30、45 g/(m2·h)時井筒溫度的分布。從圖12中可以看出,油管腐蝕對井筒溫度分布的影響不大,可以忽略。

圖12 腐蝕速度對井筒溫度分布的影響圖

5 結論

1)進行了物性分析,認為在井筒傳熱過程中水泥環具有不可忽略的非穩態效應,因此建立了軸向上離散、徑向上解析的雙層非穩態導熱井筒溫度場半解析模型。

2)比較了現場施工監測數據和模擬數據,認為半解析模型的精度較高,可以用來解決與酸化壓裂井溫計算相關的模擬與優化問題。

3)對比分析了不同水泥環導熱過程處理方法的井溫分布計算結果,認為假設水泥環具有穩態熱阻時,計算誤差大且難于修正。

4)考察了軸向分段對井溫分布的影響,認為在劃分軸向網格時,應該以巖性界面為準并按照精度要求在某些厚層加密網格。

5)巖性變化、施工排量和注入溫度對井筒溫度分布具有顯著影響,但巖性變化和注入溫度的影響隨著施工時間的增加而增強,施工排量的影響則隨著施工時間的增加而減弱;油管腐蝕對井溫分布的影響可以忽略。

符 號 說 明

r1、r2、r3、r4、r5分別是油管內壁、油管外壁、套管內壁、套管外壁和井壁到油管中心的距離,m;ρ1是油管內流體密度,kg/m3;c1是油管內流體的比熱,J/(kg·℃);v1是油管內流體的流速,m/s;t是時間,s;s是拉氏變量,無量綱;z是軸向坐標,m;r是徑向坐標,m;T1,j、T2,j、T3,j、Tair、Tinj分別是油管內的溫度、水泥環內的溫度、地層中的溫度、常溫層溫度和注入溫度,℃;q1,j、q2,j、q3,j、qi分別是油管內壁的熱流密度、水泥環中的熱流密度、地層中的熱流密度和單位面積單位時間油管腐蝕產生的熱量,W/m2;a是地溫梯度,℃/ m;f(z)是地溫函數,℃;a2,j、a3,j分別是水泥環和地層中的導溫系數,m2/s;λt、λ2、λ3分別是鋼的導熱系數、水泥環的導熱系數和地層的導熱系數,W/(m·℃);α1、αe分別是油管內的強迫對流換熱系數和環空中的腔體自然對流換熱系數,W/(m2·℃);μf、μw分別是管內溫度和管壁溫度下的流體黏度,Pa·s; Nu是努賽爾數,無量綱; Re是雷諾數,無量綱;Gr是格拉曉夫數,無量綱; Pr是普朗特數,無量綱; Ra是瑞利數,無量綱;D是油管直徑,m;Z是井深,m;H是腔體高度,m;L是腔體寬度,m;N是精度系數(必須是偶數),無量綱;T是一個確定的時間,無量綱;Φi,j是無因次溫度,無量綱;rD是無因次徑向距離,無量綱;zD是無因次軸向距離,無量綱;tD是無因次時間,無量綱;K0是0階的第一類修正貝塞爾函數,無量綱;I0是0階的第二類修正貝塞爾函數,無量綱。

參 考 文 獻

[1] 米強波, 伊向藝, 羅攀登, 任嵐. 塔中超深致密砂巖油藏水平井分段壓裂技術[J]. 西南石油大學學報: 自然科學版, 2015, 37(2): 114-118. Mi Qiangbo, Yi Xiangyi, Luo Pandeng, Ren Lan. Horizontal well staged fracturing technology of tight sandstone reservoirs with super depth in Tazhong Area[J]. Journal of Southwest Petroleum University: Science & Technology Edition, 2015, 37(2): 114-118.

[2] 韓烈祥, 朱麗華, 孫海芳, 譙抗逆. LPG無水壓裂技術[J]. 天然氣工業, 2014, 34(6): 48-54. Han Liexiang, Zhu Lihua, Sun Haifang, Qiao Kangni. LPG waterless fracturing technology[J]. Natural Gas Industry, 2014, 34(6): 48-54.

[3] 宋振云, 蘇偉東, 楊延增, 李勇, 李志航, 汪小宇, 等. CO2干法加砂壓裂技術研究與實踐[J]. 天然氣工業, 2014, 34(6): 55-59. Song Zhenyun, Su Weidong, Yang Yanzeng, Li Yong, Li Zhihang, Wang Xiaoyu, et al. Experimental studies of CO2/ sand dry-frac process[J]. Natural Gas Industry, 2014, 34(6): 55-59.

[4] 張燁, 楊勝來, 焦克波. 塔河油田超大型復合酸壓降濾失技術研究[J]. 西南石油大學學報: 自然科學版, 2014, 36(3): 121-126. Zhang Ye, Yang Shenglai, Jiao Kebo. Technology research about super large composite acid fracturing filtration in Tahe Oilfield[J]. Journal of Southwest Petroleum University: Science & Technology Edition, 2014, 36(3): 121-126.

[5] 曾凡輝, 王樹義, 郭建春, 江啟峰, 張柟喬. 裂縫面非均勻流入的氣藏壓裂水平井產量計算[J]. 天然氣工業, 2014, 34(5): 100-105. Zeng Fanhui, Wang Shuyi, Guo Jianchun, Jiang Qifeng, Zhang Nanqiao. Yield calculation of a fractured horizontal well with a non-uniform gas flow on fracture surface[J]. Natural Gas Industry, 2014, 34(5): 100-105.

[6] Ramey Jr HJ. Wellbore heat transmission[J]. Journal of Petro leum Technology, 1962, 14(4): 427-435.

[7] Van Everdingen AF, Hurst W. The application of the Laplace transformation to flow problems in reservoirs[J]. Journal of Petroleum Technology, 1949, 1(12): 305-324.

[8] Eickmeier JR, Ersoy D, Ramey HJ. Wellbore temperatures and heat losses during production or injection operations[J]. Journal of Canadian Petroleum Technology, 1970, 9(2): 115-121.

[9] Wu YS, Pruess K. An analytical solution for wellbore heat transmission in layered formations[J]. SPE Reservoir Engineering, 1990, 5(4): 531-538.

[10] Hasan AR, Kabir CS. Heat transfer during two-phase flow in wellbores: Part I—formation temperature[C]//SPE Annual Technical Conference and Exhibition, 6-9 October 1991, Dallas, Texas, USA. DOI: http: //dx. doi. org/10. 2118/22866-MS.

[11] Hasan AR, Kabir CS. Aspects of wellbore heat transfer during two-phase flow[J]. SPE Production & Facilities, 1994, 9(3): 211-216.

[12] Hasan AR, Kabir CS. Determination of static reservoir temperature from transient data following mud circulation[J]. SPE Drilling & Completion, 1994, 9(1): 17-24.

[13] Hasan AR, Kabir CS, Lin D. Analytic wellbore temperature model for transient gas-well testing[J]. SPE Reservoir Evaluation & Engineering, 2005, 8(3): 240-247.

[14] Hagoort J. Ramey's wellbore heat transmission revisited[J]. SPE Journal, 2004, 9(4): 465-474.

[15] Mcspadden AR, Coker OD. Multiwell thermal interaction: Predicting wellbore and formation temperatures for closely spaced wells[J]. SPE Drilling & Completion, 2010, 25(4): 448-457.

[16] Yoshida N, Zhu D, Hill AD. Temperature prediction model for a horizontal well with multiple fractures in a shale reservoir[C]//SPE Annual Technical Conference and Exhibition, 30 September - 2 October 2013, New Orleans, Louisiana, USA. DOI: http: //dx. doi. org/10. 2118/166241-MS.

[17] Lienhard JH. A heat transfer textbook[M]. New York: Courier Dover Publications, 2005.

[18] Holman JP. Heat transfer[M]. New York: McGraw-Hill, 1997.

[19] Rohsenow WM. Handbook of heat transfer[M]. New York: M-cGraw-Hill, 1998.

[20] 顧樵. 數學物理方法[M]. 北京: 科學出版社, 2012. Gu Qiao. Mathematical methods for physics[M]. Beijing: Science Press, 2012.

[21] 李順初, 黃炳光. Laplace變換與Bessel函數及試井分析理論基礎[M]. 北京: 石油工業出版社, 2000. Li Shunchu, Huang Bingguang. Basic theory of Laplace transform, Bessel function and well test analysis[M]. Beijing: Petroleum Industry Press, 2000.

[22] 奚定平. 貝塞爾函數[M]. 北京: 高等教育出版社, 1998. Xi Dingping. Bessel function[M]. Beijing: Higher Education Press, 1998.

[23] 同登科, 陳欽雷. 關于Laplace數值反演Stehfest方法的一點注記[J]. 石油學報, 2001, 22(6): 91-92. Tong Dengke, Chen Qinlei. Remark on Stehfest numerical inversion method of Laplace transforms[J]. Acta Petrolei Sinica, 2001, 22(6): 91-92.

[24] Stehfest H. Algorithm368: Numerical inversion of Laplace transforms[J]. Communications of the ACM, 1970, 13(1): 47-49.

[25] Stehfest H. Remark on algorithm 368: Numerical inversion of Laplace transforms[J]. Communications of the ACM, 1970, 13(10): 624.

A semi-analytic model of wellbore temperature field based on double-layer unsteady heat conducting process

Zhao Jinzhou1, Peng Yu1, Li Yongming1, Tian Zhisheng2, Fu Dongyu1
(1. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation//Southwest Petroleum University, Chengdu, Sichuan 610500, China; 2. Sinopec Henan Petroleum Engineering Co., Ltd., Nanyang, Henan 473132, China)
NATUR.GAS IND.VOLUME 36,ISSUE 1,pp.68-75, 1/25/2016.(ISSN 1000-0976;In Chinese)

Abstract:In the process of well acidizing and fracturing, downhole working fluid rheology and acid-rock reaction rate are controlled by temperature, which is an important factor affecting fracture geometry, acid etching conductivity and stimulated productivity. At present, however, the well bore temperature field model lacks strict deduction and fine description of the unsteady process. In this paper, a semianalytic model was established for temperature field of double-layer (axial discretion and radial analysis) unsteady heat conducting wellbores and solved by means of Laplace transformation and Steh fest numerical inversion. Then, analysis was conducted of the effect on wellbore temperature distribution by cement-sheath heat conducting treatment methods. It is believed that the errors of analog outputs are too large to be corrected when cement sheath is assumed to have steady state thermal resistance. Finally, sensitivity analysis was performed on formation and operation parameters. It is regarded that wellbore temperature distribution is remarkably influenced by lithological change, operational discharge capacity and injection temperature, but the effect of tubing corrosion can be neglected. Field application shows that this mathematical model is higher in precision and the simulation results fit monitoring results. Therefore, this research result is of great importance in increasing the success ratio of single-well acid fracturing and improving productivity and recovery factor after stimulation.

Keywords:Wellbore; Temperature field; Unsteady-state heat conducting; Acidizing; fracturing(Rock) ; Laplace transformation; Mathematical model

(收稿日期2015-07-15 編輯 韓曉渝)

通信作者:彭瑀,1988年生,博士研究生;主要從事油氣田增產改造理論與技術方面的研究工作。地址:(610500)四川省成都市新都區新都大道8號。ORCID:0000-0001-5233-4860。E-mail:pengyu_frac@foxmail.com

作者簡介:趙金洲,1962年生,教授,博士生導師,本刊第七屆編委會委員、《Natural Gas Industry B》編委會委員,現任西南石油大學校長;主要從事油氣田開采和增產新技術新理論的研究工作。地址:(610500)四川省成都市新都區新都大道8號。ORCID:0000-0003-1686-1828。E-mail:zhaojz@swpu.edu.cn

基金項目:國家自然科學基金項目“頁巖地層動態隨機裂縫控制機理與無水壓裂理論”(編號:51490653)、 國家重點基礎研究發展計劃(973計劃)項目“頁巖氣儲層壓裂改造機理研究”(編號:2013CB228004)。

DOI:10.3787/j.issn.1000-0976.2016.01.008

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 久久99这里精品8国产| 黄色免费在线网址| 亚洲精品在线91| 小说 亚洲 无码 精品| 色有码无码视频| 国产成人精品优优av| 综合色区亚洲熟妇在线| 亚洲精品在线91| 亚洲精品第一页不卡| 国产精品美人久久久久久AV| 久久久受www免费人成| 在线欧美a| 亚洲激情99| 欧美亚洲国产一区| 午夜一区二区三区| 无码高潮喷水在线观看| 中文字幕一区二区视频| 国产拍揄自揄精品视频网站| 亚洲AV永久无码精品古装片| 国产情侣一区二区三区| 素人激情视频福利| 亚洲女同一区二区| 中文字幕欧美日韩| 99在线视频免费观看| 色噜噜中文网| 国产午夜福利亚洲第一| 欧美日韩中文国产| 一本综合久久| a级毛片免费看| 亚洲日韩AV无码一区二区三区人| 99精品久久精品| 日韩中文欧美| 华人在线亚洲欧美精品| 五月激情综合网| 亚洲中文字幕在线观看| 色偷偷男人的天堂亚洲av| 亚洲欧美另类专区| 欧美性猛交一区二区三区| 亚洲国产精品日韩av专区| 亚洲欧美精品日韩欧美| 一本色道久久88亚洲综合| 精品欧美一区二区三区久久久| 天堂在线亚洲| 狠狠亚洲五月天| 97视频在线观看免费视频| 青青热久免费精品视频6| 国产永久在线视频| 久久毛片网| 欧美亚洲激情| 91精品国产91久无码网站| 亚洲无码高清一区| 国产欧美高清| 久久精品日日躁夜夜躁欧美| 91视频首页| 国产成人在线无码免费视频| 22sihu国产精品视频影视资讯| 亚洲精品卡2卡3卡4卡5卡区| 亚洲天堂久久| 日本不卡免费高清视频| 亚洲九九视频| 福利一区在线| av大片在线无码免费| 亚洲精品国产成人7777| 免费毛片在线| yy6080理论大片一级久久| 欧美综合区自拍亚洲综合天堂| 国产9191精品免费观看| 国产成人精品高清不卡在线| 狠狠色噜噜狠狠狠狠色综合久| 国产亚洲第一页| 中文字幕亚洲乱码熟女1区2区| 精品国产福利在线| 乱码国产乱码精品精在线播放| 午夜啪啪福利| 98超碰在线观看| 亚洲bt欧美bt精品| 亚洲不卡无码av中文字幕| 国产乱人伦精品一区二区| 最新国产成人剧情在线播放| 青青青伊人色综合久久| 亚洲综合天堂网| 最新国产成人剧情在线播放|