• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    地下熱流固耦合對(duì)EGS熱開(kāi)采的影響*

    2015-06-01 09:25:15曹文炅黃文博蔣方明
    新能源進(jìn)展 2015年6期
    關(guān)鍵詞:熱應(yīng)力應(yīng)力場(chǎng)工質(zhì)

    曹文炅,黃文博,蔣方明

    (中國(guó)科學(xué)院廣州能源研究所,廣州 510640)

    地下熱流固耦合對(duì)EGS熱開(kāi)采的影響*

    曹文炅,黃文博,蔣方明?

    (中國(guó)科學(xué)院廣州能源研究所,廣州 510640)

    人工熱儲(chǔ)的孔隙率及滲透率在增強(qiáng)型地?zé)嵯到y(tǒng)(EGS)地下熱開(kāi)采過(guò)程中受溫度(T)、水力(H)、應(yīng)力(M)的綜合影響。本文建立了EGS熱開(kāi)采過(guò)程THM耦合的三維計(jì)算模型,并采用局部非熱平衡假設(shè)處理液巖對(duì)流換熱。對(duì)一理想的五口井EGS系統(tǒng)采熱過(guò)程進(jìn)行了THM模擬計(jì)算,分析了巖石溫度、孔隙壓力對(duì)巖石應(yīng)力場(chǎng)的作用機(jī)理,進(jìn)一步研究了應(yīng)力場(chǎng)對(duì)EGS采熱性能的影響。結(jié)果表明,開(kāi)采過(guò)程中巖石應(yīng)力場(chǎng)為熱儲(chǔ)內(nèi)孔隙壓力和溫差綜合作用的結(jié)果,由孔隙壓力造成的巖石應(yīng)力為壓應(yīng)力,僅集中于注入井附近,由巖石溫度變化引起的熱應(yīng)力為拉應(yīng)力,隨著熱開(kāi)采區(qū)域的擴(kuò)展而擴(kuò)展。液?巖溫差是觸發(fā)工質(zhì)與巖石熱交換的動(dòng)因,同時(shí)也是產(chǎn)生熱應(yīng)力的根本。

    增強(qiáng)型地?zé)嵯到y(tǒng);局部非熱平衡;熱流固耦合;變物性

    0 引 言

    賦存于地下3~10 km范圍的干熱巖(Hot dry rock, HDR)地?zé)崮埽蚱淝鍧嵖稍偕院涂臻g分布的廣泛性,已成為位居水力、生物質(zhì)能之后的世界第三大可再生能源[1]。增強(qiáng)型地?zé)嵯到y(tǒng)(Enhanced geothermal systems, EGS)旨在將HDR地?zé)崮芴崛≈恋厣喜⒓右岳?。EGS的原理是通過(guò)水力激發(fā)等手段,在低滲透性的HDR內(nèi)產(chǎn)生具有一定連通性的裂隙網(wǎng)絡(luò),形成人工熱儲(chǔ)。然后由注入井灌注冷流體工質(zhì),流體工質(zhì)流過(guò)地下裂隙時(shí)獲取HDR熱量,熱流體經(jīng)由生產(chǎn)井開(kāi)采出來(lái)后用于地面發(fā)電,發(fā)電后的流體工質(zhì)經(jīng)進(jìn)一步梯級(jí)利用降溫后再回注至地下熱儲(chǔ),從而形成循環(huán)生產(chǎn)。地下采熱過(guò)程是EGS的關(guān)鍵,直接影響EGS的產(chǎn)能和壽命,而這一過(guò)程則包含了復(fù)雜的多物理場(chǎng)、多尺度綜合效應(yīng)。研究熱儲(chǔ)內(nèi)熱(Thermal, T),水力(Hydraulic, H),應(yīng)力場(chǎng)(Mechanical, M)耦合作用對(duì)揭示EGS地下熱開(kāi)采機(jī)理,合理高效地獲取HDR熱能等具有重要意義。

    在EGS熱開(kāi)采過(guò)程中,循環(huán)工質(zhì)溫度和壓力的變化將改變其熱物性,從而影響采熱過(guò)程中工質(zhì)流體的輸運(yùn)和巖石?流體換熱。在水力及熱力作用下高溫巖體則發(fā)生有效應(yīng)力場(chǎng)改變和骨架變形,進(jìn)而導(dǎo)致巖石的儲(chǔ)滲能力發(fā)生變化。這些變化反過(guò)來(lái)又影響循環(huán)工質(zhì)的滲流,影響熱量的輸運(yùn),即會(huì)影響EGS的壽命、出力等生產(chǎn)指標(biāo)。KOH等[2]結(jié)合有限元和隨機(jī)生成裂隙構(gòu)建了二維熱儲(chǔ)的THM耦合模型,研究表明高溫巖體所受熱應(yīng)力能夠在長(zhǎng)期開(kāi)采中影響熱儲(chǔ)滲透率,進(jìn)而影響溫度場(chǎng)分布和生產(chǎn)井采出溫度。GHASSEMI等[3]研究了熱應(yīng)力對(duì)裂隙滑移和擴(kuò)張的效應(yīng),并在后續(xù)研究中逐步增加了孔隙壓力效應(yīng)[4]、化學(xué)反應(yīng)及礦物質(zhì)溶解/沉積效應(yīng)[5]。JEANNE等[6]將商用軟件TOUGH與FLAC3d結(jié)合,進(jìn)行了Geysers地?zé)崽锏腡HM計(jì)算,并將位移計(jì)算值與Geysers項(xiàng)目試驗(yàn)觀測(cè)的地表位移進(jìn)行對(duì)比。MCDERMOTT等[7]研究了結(jié)晶巖體裂隙開(kāi)度的熱響應(yīng)問(wèn)題,認(rèn)為裂隙的開(kāi)度直接與裂隙平面的法向張應(yīng)力相關(guān)。然而,在現(xiàn)有的研究中,對(duì)液?巖換熱的處理通常采用局部熱平衡假設(shè),忽略了液?巖溫差,由工質(zhì)物性的變化引起的效應(yīng)也鮮有報(bào)道。

    本文根據(jù)熱儲(chǔ)層巖石中THM耦合的影響機(jī)理,考慮水的熱物性溫度和壓力而變化的情況,從飽和多孔介質(zhì)單相流體角度出發(fā),建立了EGS熱開(kāi)采過(guò)程THM耦合的三維計(jì)算模型。對(duì)一理想的五口井EGS系統(tǒng)采熱過(guò)程進(jìn)行了THM計(jì)算,分析了巖石溫度、孔隙壓力對(duì)巖石應(yīng)力場(chǎng)的作用機(jī)理,并探討液?巖溫差與巖石熱應(yīng)力的相關(guān)性,進(jìn)一步研究了應(yīng)力場(chǎng)對(duì)EGS采熱性能的影響。

    1 EGS熱開(kāi)采過(guò)程的熱流固耦合模型

    1.1 主控方程

    本文在前期EGS采熱過(guò)程數(shù)值模型工作基礎(chǔ)上引入應(yīng)力場(chǎng)計(jì)算部分[8-12]。模型基于局部非熱平衡思想,采用兩個(gè)能量方程來(lái)分別描述熱儲(chǔ)內(nèi)流體和巖石骨架的溫度場(chǎng),可方便地處理采熱過(guò)程中實(shí)際存在的巖石?流體換熱過(guò)程。對(duì)于應(yīng)力場(chǎng)的計(jì)算本文采用HU等[13]的平均總應(yīng)力模型,該模型能夠通過(guò)標(biāo)量控制方程進(jìn)行描述,可實(shí)現(xiàn)與流場(chǎng)及溫度場(chǎng)計(jì)算的強(qiáng)耦合。模型將EGS的地下部分處理為三個(gè)性質(zhì)不同的子區(qū)域:①開(kāi)放流道性質(zhì)的注入井和生產(chǎn)井; ②多孔介質(zhì)性質(zhì)的熱儲(chǔ);③滲透性可忽略不計(jì)的熱儲(chǔ)周圍巖體。模型假設(shè)單相流體流動(dòng),不考慮循環(huán)流體與巖石的化學(xué)作用。采用的控制方程如下。

    循環(huán)工質(zhì)的連續(xù)性方程:

    循環(huán)工質(zhì)的動(dòng)量守恒方程:

    巖石的能量守恒方程:

    循環(huán)工質(zhì)的能量守恒方程:

    巖石的平均總應(yīng)力方程:

    其中,u、p、Tf、Ts及σm分別表示工質(zhì)速度矢量、工質(zhì)壓力、工質(zhì)溫度、巖石溫度及巖石平均總應(yīng)力,為T(mén)HM模型的主要求解變量;ρf、cpf、kf、μ分別表示工質(zhì)密度、比熱容、有效導(dǎo)熱系數(shù)和動(dòng)力粘度,受工質(zhì)的溫度和壓力決定,將通過(guò)后續(xù)變物性模型給出變量隨溫度和壓力變化的關(guān)系;g為重力加速度;ε和K分別為熱儲(chǔ)的孔隙率和滲透率,下文中將通過(guò)等效應(yīng)力模型進(jìn)行描述;ha表示巖石?流體對(duì)流換熱強(qiáng)度,本文取ha=1.0 W·m?3·K?1;下標(biāo)s和f分別代表巖石和流體,上標(biāo)eff表示有效物性。式(5)中υ、α、β、B分別表示巖石的泊松比、Biot數(shù)、線膨脹系數(shù)以及體積模量;F為熱儲(chǔ)所受外力。

    1.2 水工質(zhì)的變物性模型

    本文采用國(guó)際水及蒸汽物性組織(IAPWS)[14]建立的公式進(jìn)行水工質(zhì)的變物性模型建立。IAPWS模型針對(duì)水不同的相建立了相應(yīng)的區(qū)域和方程,考慮到EGS系統(tǒng)中水的溫度和壓力范圍(壓力低于100 MPa,溫度低于623.15 K),本文選用I區(qū)域(液相區(qū)域)進(jìn)行建模。該區(qū)域采用Gibbs自由能描述水的狀態(tài):

    其中,g(p, T)為Gibbs自由能;Rw為水的比氣體常數(shù),Rw=461.526 J·kg?1·K?1;π與τ分別為無(wú)量綱化的壓力和溫度,π=p/p* , τ=T*/T ,參考?jí)毫?p*=16.53 MPa,參考溫度T*=1 386 K;ni、Ii、Ji均為常數(shù),可參考文獻(xiàn)[14]。利用Gibbs自由能的偏導(dǎo)數(shù),就可得到隨溫度和壓力變化的密度及比熱容。

    導(dǎo)熱系數(shù)和粘度系數(shù)則由下述多項(xiàng)式給出[15-16]:

    其中,參考值λ*=1.0 × 10?3W·m?1·K?1,μ*=1.0 × 10?6Pa·s;及為無(wú)量綱溫度及密度, 可表示為=T/ T'及=ρ/ ρ',參考值T'=647.096 K,ρ'=322.0 kg·m?3;方程中下標(biāo)0、1和2的關(guān)聯(lián)項(xiàng)分別為稀釋氣體項(xiàng),有限密度項(xiàng)及臨界點(diǎn)修正項(xiàng)。導(dǎo)熱系數(shù)及粘度系數(shù)的稀釋氣體項(xiàng)和有限密度項(xiàng)可由下列各式計(jì)算:

    式中,Li、Lij、Hi及 Hij均為常系數(shù),取值可由文獻(xiàn)[15-16]查得。由于熱儲(chǔ)內(nèi)的溫度沒(méi)有達(dá)到水的臨界點(diǎn),因此本文計(jì)算中取

    1.3 THM耦合模型的實(shí)現(xiàn)

    本文采用Fluent進(jìn)行控制方程的求解,其中能量守恒方程、平均總應(yīng)力方程均以標(biāo)量控制方程進(jìn)行計(jì)算,THM耦合機(jī)制如圖1所示。熱流兩場(chǎng)之間的耦合由工質(zhì)的變物性模型實(shí)現(xiàn)。熱儲(chǔ)層的孔隙率及滲透率是滲流及傳熱的關(guān)鍵參數(shù),也是應(yīng)力場(chǎng)與熱、流計(jì)算耦合的關(guān)鍵。在巖土力學(xué)中孔隙率及滲透率可表示為巖石有效應(yīng)力的相關(guān)函數(shù),本文采用的計(jì)算模型如下[17-19]:

    其中,σ'為巖石有效應(yīng)力,0ε及rε分別為在零有效應(yīng)力及高有效應(yīng)力狀態(tài)下的孔隙率;K0為初始孔隙率;ξ與ζ均為與材料相關(guān)的常系數(shù)。

    圖1 THM耦合機(jī)制Fig. 1 Mechanism of the THM coupling

    2 算例分析及討論

    2.1 計(jì)算模型

    圖2為某假設(shè)五口井分布EGS,由于對(duì)稱性取1/4進(jìn)行計(jì)算。人工熱儲(chǔ)的體積為500 m × 500 m × 500 m的立方體。熱儲(chǔ)周圍包覆有足夠體積的基巖和蓋巖,避免了人為設(shè)定的邊界條件帶來(lái)的誤差。注入井和生產(chǎn)井均為0.2 m × 0.2 m的方形通道。巖石的初始溫度按照4 K/100 m的地溫梯度隨深度線性增加,地表溫度設(shè)為300 K,熱儲(chǔ)中裂隙流體初始溫度與當(dāng)?shù)貛r石溫度相同。所有與流體接觸的壁面均為非滑移邊界,注入井與生產(chǎn)井均采用定壓力邊界,注入井與生產(chǎn)井壓差取10 MPa,熱儲(chǔ)內(nèi)環(huán)境壓力取40 MPa。注入溫度為343 K,其余計(jì)算參數(shù)列于表1,計(jì)算模型的幾何尺寸設(shè)置如圖2所示。

    表1 模型計(jì)算參數(shù)Table 1 Parameters of the THM model

    圖2 五口井EGS模型幾何尺寸Fig. 2 Geometrical dimensions of the considered quintuplet EGS

    2.2 熱開(kāi)采過(guò)程溫度場(chǎng)及滲流場(chǎng)對(duì)應(yīng)力場(chǎng)的影響

    圖3顯示了采熱過(guò)程中熱儲(chǔ)內(nèi)巖石溫度隨時(shí)間的變化。注入井附近巖石熱量首先被采集,溫度迅速降低。隨著時(shí)間的推移低溫區(qū)域逐漸向生產(chǎn)井一側(cè)擴(kuò)展。由于五口井分布的對(duì)稱特征,熱儲(chǔ)內(nèi)工質(zhì)的流動(dòng)以注入井為中心向四周擴(kuò)展,溫度采集前沿位置為近似的柱面。

    圖4顯示了熱儲(chǔ)內(nèi)孔隙壓力隨時(shí)間的變化。由圖可知在整個(gè)開(kāi)采階段孔隙壓力分布變化較小,高壓區(qū)域集中于注入井附近。注入井附近壓降極大,表明循環(huán)工質(zhì)的動(dòng)量損失主要集中在注入井附近。

    圖5顯示了巖石平均總應(yīng)力空間分布隨熱開(kāi)采進(jìn)行的發(fā)展情況,其量級(jí)和分布與文獻(xiàn)[5,13]具有一致性:在非常接近注入井的區(qū)域平均總應(yīng)力為正值,根據(jù)巖土力學(xué)中壓應(yīng)力為正的約定可知注入井附近的巖石受壓應(yīng)力作用。而在熱儲(chǔ)層內(nèi)還存在著一個(gè)從注入井隨開(kāi)采的進(jìn)行逐漸向生產(chǎn)井發(fā)展的拉應(yīng)力區(qū)域(平均總應(yīng)力為負(fù)值),該區(qū)域隨著熱開(kāi)采的進(jìn)行逐漸向注入井?dāng)U張。

    巖石平均總應(yīng)力是孔隙壓力和巖石熱應(yīng)力綜合作用的結(jié)果,為了說(shuō)明孔隙壓力和溫度場(chǎng)對(duì)應(yīng)力計(jì)算的影響,圖6中提取并對(duì)比了開(kāi)采至第5年時(shí)由孔隙壓力引起的巖石應(yīng)力及溫度變化引起的巖石應(yīng)力,即孔隙彈性應(yīng)力及熱彈性應(yīng)力。

    由圖6可知,由孔隙壓力造成的巖石應(yīng)力為壓應(yīng)力,僅集中于注入井附近,由巖石溫度變化引起的熱應(yīng)力為拉應(yīng)力,隨著熱開(kāi)采區(qū)域的擴(kuò)展而擴(kuò)展。從應(yīng)力幅值來(lái)看,孔隙彈性應(yīng)力幅值大于熱彈性應(yīng)力幅值,說(shuō)明在注入井附近孔隙彈性應(yīng)力起主導(dǎo)作用。熱彈性應(yīng)力作用相對(duì)較弱,但其空間分布大于孔隙彈性應(yīng)力,并隨著熱開(kāi)采進(jìn)行不斷擴(kuò)展。為了進(jìn)一步研究液?巖溫差對(duì)熱應(yīng)力計(jì)算的影響,我們提取了圖6b所示AB路徑的液?巖溫差及熱應(yīng)力進(jìn)行分析,對(duì)比結(jié)果如圖7所示。

    在局部非熱平衡模型下,巖石與循環(huán)工質(zhì)的溫差分布直接反映了熱儲(chǔ)層內(nèi)發(fā)生熱交換的范圍,即開(kāi)采區(qū)域。由圖7可以看出,在距離注入井約0~60 m范圍內(nèi)液?巖溫差為零,結(jié)合圖3可知該區(qū)域巖石熱能已被充分開(kāi)采,巖石溫度已達(dá)到工質(zhì)的注入溫度;在60~250 m范圍存在液?巖溫差,表明該區(qū)域?yàn)楫?dāng)前時(shí)刻熱開(kāi)采進(jìn)行區(qū)域;在遠(yuǎn)離注入井250 m至生產(chǎn)井(距注入井約283 m)范圍內(nèi),液?巖溫差趨于零,工質(zhì)已得到充分加熱并與巖石達(dá)到溫度平衡。值得注意的是,熱儲(chǔ)內(nèi)巖石所受的熱應(yīng)力變化范圍與液?巖溫差曲線很好的吻合:在0~60 m區(qū)域熱應(yīng)力穩(wěn)定于 ?1.6 MPa左右,該區(qū)域熱應(yīng)力是由此時(shí)刻之前的熱開(kāi)采所形成,拉應(yīng)力的最大值為 ?1.62 MPa,其位置吻合于液?巖溫差曲線中低溫平衡區(qū)域與開(kāi)采區(qū)域的交界位置(即距注入井約60 m位置);在60~250 m區(qū)域內(nèi),熱應(yīng)力幅值隨著距注入井距離的增大而逐漸遞減,而熱應(yīng)力變化曲線的拐點(diǎn)則吻合于液?巖溫差曲線的峰值點(diǎn);在遠(yuǎn)離注入井250 m的區(qū)域,熱應(yīng)力值趨于零。從物理角度而言,液?巖溫差是觸發(fā)工質(zhì)與巖石熱交換的動(dòng)因,巖石在消耗熱能的同時(shí)溫度降低并產(chǎn)生體積收縮,從而產(chǎn)生拉應(yīng)力。該拉應(yīng)力隨著巖石溫度的降低而逐漸累積,液?巖溫差最大位置則為當(dāng)前時(shí)刻拉應(yīng)力變化最顯著的位置,當(dāng)巖石溫度降低至工質(zhì)注入溫度時(shí),拉應(yīng)力達(dá)到最大值。以上分析也說(shuō)明了液?巖溫差是造成熱儲(chǔ)內(nèi)巖石熱應(yīng)力變化的根本動(dòng)因。

    圖3 熱儲(chǔ)內(nèi)巖石溫度變化Fig. 3 Rock temperature distribution in the heat reservoir

    圖4 熱儲(chǔ)內(nèi)孔隙壓力變化Fig. 4 Pore pressure distribution in the heat reservoir

    圖5 巖石平均總應(yīng)力變化Fig. 5 Mean total stress distribution in the heat reservoir

    圖6 開(kāi)采至第5年時(shí)孔隙彈性應(yīng)力及熱彈性應(yīng)力對(duì)比Fig. 6 Distribution of poro- and thermo- elastic stress at 5 years into the EGS operation

    圖7 液?巖溫差與巖石熱應(yīng)力對(duì)比Fig. 7 The comparison between the liquid-rock temperature difference and the thermal stress in the reservoir

    2.3 熱開(kāi)采過(guò)程應(yīng)力場(chǎng)對(duì)采熱性能的影響

    由式(15)~式(17)計(jì)算獲得的開(kāi)采至第5年時(shí)熱儲(chǔ)有效應(yīng)力、孔隙率及滲透率分布如圖8所示。在注入井相鄰區(qū)域,雖然平均總應(yīng)力表現(xiàn)為最大壓應(yīng)力,但由于該位置孔隙壓力也為最大值,因此該區(qū)域有效應(yīng)力達(dá)到最大負(fù)值,對(duì)應(yīng)的孔隙率及滲透率均達(dá)到最大值,表明該區(qū)域內(nèi)裂隙開(kāi)度在較大孔隙壓力作用下產(chǎn)生增長(zhǎng)。在距注入井較遠(yuǎn)范圍內(nèi),孔隙率及滲透率同樣大于初始值,該區(qū)域主要是因?yàn)閹r石溫度降低產(chǎn)生收縮,表明該區(qū)域?yàn)闊釕?yīng)力作用區(qū)域。

    圖8 開(kāi)采至第5年時(shí)有效應(yīng)力、孔隙率、滲透率分布Fig. 8 Distributions of effective stress, porosity and permeability at 5 years

    當(dāng)孔隙率及滲透率發(fā)生改變時(shí),熱儲(chǔ)內(nèi)的滲流性能及換熱性能均受到較大影響,從而使EGS采熱性能發(fā)生顯著變化。圖9對(duì)比了有無(wú)應(yīng)力場(chǎng)作用下出口井質(zhì)量流量隨時(shí)間的變化情況??梢钥闯?,在采熱的最初階段,兩種條件下出口井質(zhì)量流量均表現(xiàn)為劇烈降低過(guò)程,這是因?yàn)樽⑷刖浇鼰醿?chǔ)開(kāi)始有低溫工質(zhì)流入,而工質(zhì)在低溫區(qū)域的粘度系數(shù)顯著增大,使得注入井附近工質(zhì)的動(dòng)量損失相應(yīng)增大。而在應(yīng)力場(chǎng)作用條件下,出口井質(zhì)量流量開(kāi)始回升,由圖8可知注入井相鄰區(qū)域滲透率有所增大,降低了該區(qū)域的流動(dòng)阻力。在后續(xù)運(yùn)行中考慮應(yīng)力場(chǎng)效應(yīng)時(shí)出口井質(zhì)量流率則在較高范圍內(nèi)變化(37 kg/s降低至30 kg/s),而不考慮應(yīng)力場(chǎng)效應(yīng)的質(zhì)量流率低于20 kg/s。后續(xù)運(yùn)行中質(zhì)量流率的逐漸降低是因?yàn)闊醿?chǔ)內(nèi)處于低溫的流體逐漸增多,動(dòng)量損失隨著粘度系數(shù)的增大而增大。

    圖9 熱儲(chǔ)層內(nèi)應(yīng)力場(chǎng)對(duì)出口井質(zhì)量流量的影響Fig. 9 Mechanical effects on the mass flow rate in EGS reservoir

    圖10 熱儲(chǔ)層內(nèi)應(yīng)力場(chǎng)對(duì)開(kāi)采壽命及開(kāi)采率的影響Fig. 10 Mechanical effects on the lifetime and the heat extraction ratio of the EGS

    圖10比較了有無(wú)應(yīng)力場(chǎng)影響下的采出溫度及采出熱量隨時(shí)間的變化情況,其中采出熱量的描述由采熱比(Heat Extraction ratio)給出:

    式中,Tg為地表溫度,Ts,ini為熱儲(chǔ)的初始溫度,Vh為熱儲(chǔ)體積,Vb為基巖體積,Vh + Vb構(gòu)成了圖2所示的計(jì)算域體積。等式右端分母為以地表溫度為參考的熱儲(chǔ)內(nèi)總熱能,分子包含了由熱儲(chǔ)內(nèi)與基巖內(nèi)已開(kāi)采出的能量,該參數(shù)的實(shí)質(zhì)是t時(shí)刻由EGS系統(tǒng)采出的無(wú)量綱化的總能量??梢钥闯?,若以采出溫度降低10 K作為系統(tǒng)廢止的條件,即在本文中出口井采出溫度降至450 K的時(shí)刻作為該系統(tǒng)的開(kāi)采壽命,則在應(yīng)力場(chǎng)影響下的開(kāi)采壽命約為6.5年,低于不考慮應(yīng)力場(chǎng)效應(yīng)的開(kāi)采壽命(12年)。在各自廢止的時(shí)刻,應(yīng)力場(chǎng)影響下的EGS熱開(kāi)采率約為0.25,不考慮應(yīng)力場(chǎng)效應(yīng)的熱開(kāi)采率約為0.26,略高于前者,這是因?yàn)殚_(kāi)采壽命較長(zhǎng)條件下周圍巖石能夠依靠熱傳導(dǎo)對(duì)已開(kāi)采的低溫區(qū)域進(jìn)行熱補(bǔ)償,從而提高了熱儲(chǔ)的整體熱開(kāi)采率。以上結(jié)果表明應(yīng)力場(chǎng)效應(yīng)對(duì)EGS的采熱性能具有顯著影響,特別是注入井附近區(qū)域孔隙率滲透率的改變極大影響了后續(xù)熱開(kāi)采過(guò)程。

    3 結(jié) 論

    本文在基于局部非熱平衡假設(shè)的EGS熱開(kāi)采過(guò)程三維計(jì)算模型基礎(chǔ)上,根據(jù)熱儲(chǔ)層巖石中THM耦合的影響機(jī)理,考慮水的熱物性隨溫度和壓力的改變而變化的情況,從飽和多孔介質(zhì)單相流體角度出發(fā),建立了EGS熱開(kāi)采過(guò)程THM耦合的三維計(jì)算模型。對(duì)一理想的五口井EGS系統(tǒng)采熱過(guò)程進(jìn)行了THM計(jì)算,分析了巖石溫度、孔隙壓力對(duì)巖石應(yīng)力場(chǎng)的作用機(jī)理,進(jìn)一步研究了應(yīng)力場(chǎng)對(duì)EGS采熱性能的影響,主要結(jié)論如下:

    (1)由孔隙壓力造成的巖石應(yīng)力為壓應(yīng)力,在距離注入井20 m范圍內(nèi)壓應(yīng)力高于1.0 MPa,由巖石溫度變化引起的熱應(yīng)力為拉應(yīng)力,隨著熱開(kāi)采區(qū)域的擴(kuò)展而擴(kuò)展。

    (2)液?巖溫差是觸發(fā)工質(zhì)與巖石熱交換的動(dòng)因,同時(shí)也是產(chǎn)生熱應(yīng)力的根本。

    (3)應(yīng)力場(chǎng)效應(yīng)對(duì)EGS的采熱性能具有顯著影響,特別是注入井附近區(qū)域孔隙率滲透率的改變極大影響了后續(xù)熱開(kāi)采過(guò)程。

    [1] TESTER J W, ANDERSON B J, BATCHELOR A S, et al. The future of geothermal energy: impact of enhanced geothermal systems (EGS) on the United States in the 21st century[R]. America: Massachussets of Institute of Technology, 2006.

    [2] KOH J, ROSHAN H, RAHMAN S S. A numerical study on the long term thermo-poroelastic effects of cold water injection into naturally fractured geothermal reservoirs[J]. Computers and Geotechnics, 2011, 38(5): 669-682.

    [3] GHASSEMI A, ZHOU X. A three-dimensional thermo- poroelastic model for fracture response to injection/extraction in enhanced geothermal systems[J]. Geothermics, 2011, 40(1): 39-49.

    [4] GHASSEMI A, TARASOVS S, CHENG A H D. A 3-D study of the effects of thermo-mechanical loads on fracture slip in enhanced geothermal reservoirs[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(8): 1132-1148.

    [5] RAWAL C, GHASSEMI A. A reactive thermo-poroelastic analysis of water injection into an enhanced geothermal reservoir[J]. Geothermics, 2014, 50: 10-23.

    [6] JEANNE P, RUTQVIST J, VASCO D, et al. A 3D hydrogeological and geomechanical model of an enhanced geothermal system at the Geysers, California[J]. Geothermics, 2014, 51: 240-252.

    [7] MCDERMOTT C I, RANDRIAMANJATOSOA A R L, TENZER H, et al. Simulation of heat extraction from crystalline rocks: The influence of coupled processes on differential reservoir cooling[J]. Geothermics, 2006, 35(3): 321-344.

    [8] JIANG F M, LUO L, CHEN J L. A novel three-dimensional transient model for subsurface heat exchange in enhanced geothermal systems[J]. International Communications in Heat and Mass Transfer, 2013, 41: 57-62.

    [9] JIANG F M, CHEN J L, HUANG W B, et al. A three-dimensional transient model for EGS subsurface thermo-hydraulic process[J]. Energy 2014, 72: 300-310.

    [10] CHEN J L, JIANG F M. Designing multi-well layout for enhanced geothermal system to better exploit hot dry rock geothermal energy[J]. Renewable Energy, 2015, 74: 37-48.

    [11] 陳繼良, 羅良, 蔣方明. 熱儲(chǔ)周圍巖石熱補(bǔ)償對(duì)增強(qiáng)型地?zé)嵯到y(tǒng)采熱過(guò)程的影響[J]. 計(jì)算物理, 2013, 30(6): 862-870.

    [12] 陳繼良, 蔣方明, 羅良. 增強(qiáng)型地?zé)嵯到y(tǒng)地下滲流場(chǎng)的模擬和分析[J]. 計(jì)算物理, 2013, 30(6): 871-878.

    [13] HU L T, WINTERFELD P H, FAKCHAROENPHOL P, et al. A novel fully-coupled flow and geomechanics model in enhanced geothermal reservoirs[J]. Journal of Petroleum Science and Engineering, 2013, 107: 1-11.

    [14] The International Association for the Properties of Water and Steam. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam[R]. Switzerland, 2007.

    [15] The International Association for the Properties of Water and Steam. lamda Release on the IAPWS Formulation 2011 for the Thermal Conductivity of Ordinary Water Substance[R]. Czech Republic, 2011.

    [16] The International Association for the Properties of Water and Steam. Viscosity Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance[R]. Germany, 2008.

    [17] ARAIRO W, PRUNIER F, DJERAN M I, et al. On the use of effective stress in three-dimensional hydro- mechanical coupled model[J]. Computers and Geotechnics, 2014, 58: 56-68.

    [18] NATHENSON M. The dependence of permeability on effective tress from flow tests at hot dry rock reservoirs at Rosemanowes (Cornwall) and Fenton Hill (New Mexico)[J]. Geothermics, 1999, 28: 315-340.

    [19] RUTQVIST J, WU Y S, TSANG C F. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, anddeformation in fracturedporous rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2002, 39: 429-442.

    The Thermal-Hydraulic-Mechanical Coupling Effects on Heat Extraction of Enhanced Geothermal Systems

    CAO Wen-jiong, HUANG Wen-bo, JIANG Fang-ming
    (Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences, Guangzhou 510640, China)

    During heat extraction in enhanced geothermal systems (EGS), the reservoir porosity and permeability can be greatly affected by the multi-physical coupling of Thermal (T), Hydraulic (H), and Mechanical (M) actions. In the present work we develop a three-dimensional transient model coupling the subsurface THM behaviors during EGS heat extraction process. The local thermal non-equilibrium is assumed when describing the heat exchange between the rock matrix and heat transmission fluid. Case studies with respect to an imaginary quintuplet EGS reveal the involved mechanisms of inter-couplings in-between T-H-M actions, and the results indicate significant mechanical effects on EGS heat extraction performance. The stress of the rock matrix is largely influenced by the pore pressure and the temperature distributions. The stress triggered by fluid pressure is found to be compressive and confined in the very vicinity region of the injection well; the thermal stress is tensile and to some extent also concentrates around the injection well, but its distribution region expands toward the production well with the proceeding of heat extraction process. The temperature difference between rock matrix and heat transmission fluid is not only the driving force of heat extraction from heat reservoir but also significantly affects the formation of thermal stress in the reservoir.

    enhanced geothermal system; local thermal non-equilibrium; THM coupling; variable properties

    TK529

    A

    10.3969/j.issn.2095-560X.2015.06.006

    2095-560X(2015)06-0444-08

    曹文炅(1983-),男,博士,助理研究員,主要從事增強(qiáng)型地?zé)嵯到y(tǒng)的數(shù)值模擬及實(shí)驗(yàn)研究。

    2015-09-18

    2015-10-12

    中科院百人計(jì)劃(FJ);國(guó)家自然科學(xué)基金(51406213);NSFC-廣東聯(lián)合基金(U1401232);廣東省自然科學(xué)基金重大基礎(chǔ)培育項(xiàng)目(2014A030308001)

    ? 通信作者:蔣方明,E-mail:jiangfm@ms.giec.ac.cn

    黃文博(1990-),男,博士研究生,主要從事增強(qiáng)型地?zé)嵯到y(tǒng)地下物理過(guò)程的數(shù)值模擬研究。

    蔣方明(1973-),男,博士,研究員,博士生導(dǎo)師,中國(guó)科學(xué)院“百人計(jì)劃”引進(jìn)海外杰出人才。目前主要從事電化學(xué)能源/動(dòng)力系統(tǒng)、增強(qiáng)型地?zé)嵯到y(tǒng)、微熱流體系統(tǒng)、以及高效節(jié)能技術(shù)/產(chǎn)品等前沿科學(xué)和應(yīng)用技術(shù)研究。

    猜你喜歡
    熱應(yīng)力應(yīng)力場(chǎng)工質(zhì)
    海洋溫差能發(fā)電熱力循環(huán)系統(tǒng)的工質(zhì)優(yōu)選
    WNS型鍋爐煙管管端熱應(yīng)力裂紋原因分析
    采用R1234ze(E)/R245fa的非共沸混合工質(zhì)有機(jī)朗肯循環(huán)系統(tǒng)實(shí)驗(yàn)研究
    采用二元非共沸工質(zhì)的有機(jī)朗肯循環(huán)熱力學(xué)分析
    若干低GWP 純工質(zhì)在空調(diào)系統(tǒng)上的應(yīng)用分析
    采用單元基光滑點(diǎn)插值法的高溫管道熱應(yīng)力分析
    鋁合金多層多道窄間隙TIG焊接頭應(yīng)力場(chǎng)研究
    焊接(2016年9期)2016-02-27 13:05:22
    考慮斷裂破碎帶的丹江口庫(kù)區(qū)地應(yīng)力場(chǎng)與水壓應(yīng)力場(chǎng)耦合反演及地震預(yù)測(cè)
    基于位移相關(guān)法的重復(fù)壓裂裂縫尖端應(yīng)力場(chǎng)研究
    斷塊油氣田(2014年5期)2014-03-11 15:33:49
    基于流熱固耦合的核電蒸汽發(fā)生器傳熱管熱應(yīng)力數(shù)值模擬
    深爱激情五月婷婷| 18禁在线播放成人免费| 天堂√8在线中文| 国产精品福利在线免费观看| 两个人的视频大全免费| 波多野结衣巨乳人妻| 欧美另类亚洲清纯唯美| 亚洲久久久久久中文字幕| 成人性生交大片免费视频hd| 国产精品嫩草影院av在线观看| 日本欧美国产在线视频| 国产精品熟女久久久久浪| 人妻少妇偷人精品九色| 男插女下体视频免费在线播放| 18禁裸乳无遮挡免费网站照片| 最近手机中文字幕大全| 午夜激情福利司机影院| 欧美丝袜亚洲另类| 成人性生交大片免费视频hd| 嫩草影院精品99| 日日摸夜夜添夜夜添av毛片| 人人妻人人看人人澡| 国产精品野战在线观看| 国产在线男女| 男女那种视频在线观看| 黄色欧美视频在线观看| 免费观看精品视频网站| 欧美日本视频| 欧美激情久久久久久爽电影| 成年版毛片免费区| 永久网站在线| 午夜福利在线在线| 久久久久精品久久久久真实原创| 久久久成人免费电影| 久久这里只有精品中国| 黄片wwwwww| 国产精品麻豆人妻色哟哟久久 | 亚洲欧洲日产国产| 天天躁日日操中文字幕| 好男人在线观看高清免费视频| 午夜激情欧美在线| 三级国产精品片| 精品久久久噜噜| 看十八女毛片水多多多| 日产精品乱码卡一卡2卡三| 久久久久久久久久黄片| 国产精品爽爽va在线观看网站| 在现免费观看毛片| 免费看av在线观看网站| 色尼玛亚洲综合影院| kizo精华| 身体一侧抽搐| 亚洲精品日韩在线中文字幕| 五月玫瑰六月丁香| 国产免费视频播放在线视频 | 狂野欧美白嫩少妇大欣赏| 在线免费观看不下载黄p国产| 国产精品国产三级国产专区5o | 亚洲欧洲国产日韩| 色哟哟·www| 欧美激情在线99| 99热6这里只有精品| 亚洲久久久久久中文字幕| 中文在线观看免费www的网站| 免费一级毛片在线播放高清视频| 中文资源天堂在线| 国产精品乱码一区二三区的特点| 欧美最新免费一区二区三区| 男的添女的下面高潮视频| 久久欧美精品欧美久久欧美| 97人妻精品一区二区三区麻豆| 国产精品人妻久久久影院| 大香蕉97超碰在线| 精品久久久久久久人妻蜜臀av| 在线观看av片永久免费下载| 美女xxoo啪啪120秒动态图| 国产精品一二三区在线看| 欧美日韩国产亚洲二区| 国产精品永久免费网站| 国产精品一区二区三区四区免费观看| 国产成人免费观看mmmm| 亚洲欧美一区二区三区国产| 日韩国内少妇激情av| 成人三级黄色视频| 一区二区三区免费毛片| 欧美潮喷喷水| 国产真实伦视频高清在线观看| 综合色丁香网| 中文字幕人妻熟人妻熟丝袜美| 在线天堂最新版资源| 国产精品麻豆人妻色哟哟久久 | 日韩制服骚丝袜av| 国内精品宾馆在线| 最近最新中文字幕大全电影3| 国产av码专区亚洲av| 亚洲综合色惰| 国产精品一区www在线观看| a级一级毛片免费在线观看| 精品人妻偷拍中文字幕| 天美传媒精品一区二区| 中文字幕人妻熟人妻熟丝袜美| 天堂中文最新版在线下载 | 国产精品综合久久久久久久免费| 国产精品精品国产色婷婷| 国产精品一二三区在线看| 不卡视频在线观看欧美| 国产伦精品一区二区三区视频9| 好男人视频免费观看在线| 老师上课跳d突然被开到最大视频| 亚洲第一区二区三区不卡| 99在线人妻在线中文字幕| 精品人妻熟女av久视频| 日韩国内少妇激情av| 国内揄拍国产精品人妻在线| 日韩高清综合在线| 夜夜爽夜夜爽视频| 亚洲精品亚洲一区二区| 免费搜索国产男女视频| av在线老鸭窝| 精品人妻偷拍中文字幕| 亚洲国产精品成人综合色| 亚洲国产精品国产精品| 在线观看av片永久免费下载| 日韩av在线大香蕉| 日韩视频在线欧美| 亚洲自拍偷在线| 两个人视频免费观看高清| 直男gayav资源| 国产精品久久久久久久电影| 狂野欧美激情性xxxx在线观看| 成人二区视频| 熟女电影av网| 嫩草影院精品99| 欧美成人免费av一区二区三区| 人人妻人人澡人人爽人人夜夜 | 乱人视频在线观看| 色尼玛亚洲综合影院| 日本-黄色视频高清免费观看| 欧美一区二区国产精品久久精品| 嘟嘟电影网在线观看| 国产精品一区二区在线观看99 | 天美传媒精品一区二区| 一级爰片在线观看| 久久99蜜桃精品久久| 毛片一级片免费看久久久久| 久久这里有精品视频免费| 久久99蜜桃精品久久| av卡一久久| 亚洲av成人精品一区久久| 日韩一本色道免费dvd| 国产精品女同一区二区软件| 国产成人精品一,二区| 久久久亚洲精品成人影院| 美女xxoo啪啪120秒动态图| av线在线观看网站| 日韩av在线免费看完整版不卡| 校园人妻丝袜中文字幕| 青春草视频在线免费观看| 欧美又色又爽又黄视频| 亚洲综合精品二区| 欧美人与善性xxx| www.色视频.com| 1024手机看黄色片| 亚洲国产日韩欧美精品在线观看| 日韩制服骚丝袜av| 久久99热6这里只有精品| 网址你懂的国产日韩在线| 亚州av有码| 国产精品国产三级国产专区5o | 成人av在线播放网站| 久久久久免费精品人妻一区二区| 国产精品久久久久久精品电影| 99在线视频只有这里精品首页| 国产亚洲精品av在线| 亚洲美女视频黄频| 久久精品国产亚洲网站| 熟女电影av网| 国产精品精品国产色婷婷| 亚洲国产精品久久男人天堂| 国产极品天堂在线| 国产在视频线精品| 国产探花极品一区二区| 免费电影在线观看免费观看| 非洲黑人性xxxx精品又粗又长| 国内精品一区二区在线观看| 九九热线精品视视频播放| 午夜精品国产一区二区电影 | 成人欧美大片| 精品一区二区免费观看| 干丝袜人妻中文字幕| 在线天堂最新版资源| 听说在线观看完整版免费高清| 精品久久久久久成人av| 日本免费在线观看一区| 岛国在线免费视频观看| 99热这里只有是精品50| 欧美bdsm另类| 午夜激情福利司机影院| 欧美日韩综合久久久久久| 国产一区二区亚洲精品在线观看| 内地一区二区视频在线| 国产爱豆传媒在线观看| 天天躁日日操中文字幕| 国产精品99久久久久久久久| 精品久久国产蜜桃| 国产 一区 欧美 日韩| 亚洲精品日韩在线中文字幕| 一本一本综合久久| 白带黄色成豆腐渣| 亚洲av电影在线观看一区二区三区 | 深夜a级毛片| 亚洲欧美清纯卡通| www.av在线官网国产| 国产成人精品婷婷| 国产在线一区二区三区精 | 禁无遮挡网站| 成年版毛片免费区| 观看美女的网站| 久久人妻av系列| 大又大粗又爽又黄少妇毛片口| 国产一区二区三区av在线| 日本一二三区视频观看| 国产黄片美女视频| 国产高清国产精品国产三级 | a级毛片免费高清观看在线播放| 国产91av在线免费观看| 村上凉子中文字幕在线| 能在线免费观看的黄片| 国产熟女欧美一区二区| 两个人视频免费观看高清| 久久这里有精品视频免费| 亚洲美女视频黄频| 秋霞在线观看毛片| 国内精品美女久久久久久| 国产黄片视频在线免费观看| 免费一级毛片在线播放高清视频| 亚洲自偷自拍三级| 韩国av在线不卡| 夜夜看夜夜爽夜夜摸| 亚洲av电影在线观看一区二区三区 | 毛片一级片免费看久久久久| 精品少妇黑人巨大在线播放 | 熟妇人妻久久中文字幕3abv| 丰满乱子伦码专区| 国产日韩欧美在线精品| 亚洲欧美成人综合另类久久久 | 五月伊人婷婷丁香| 三级毛片av免费| 人人妻人人澡人人爽人人夜夜 | 嫩草影院入口| 内射极品少妇av片p| 亚洲内射少妇av| 极品教师在线视频| 只有这里有精品99| 亚洲人与动物交配视频| 国产极品精品免费视频能看的| 美女xxoo啪啪120秒动态图| 十八禁国产超污无遮挡网站| 男女下面进入的视频免费午夜| 国产人妻一区二区三区在| 亚洲经典国产精华液单| av视频在线观看入口| 国产成人午夜福利电影在线观看| 亚州av有码| 99热这里只有精品一区| 波多野结衣高清无吗| 亚洲综合精品二区| 亚洲成人av在线免费| 欧美一区二区亚洲| 久久久精品大字幕| 亚洲av成人精品一二三区| 久久99热6这里只有精品| 纵有疾风起免费观看全集完整版 | 乱系列少妇在线播放| 欧美日韩综合久久久久久| 日本av手机在线免费观看| 精品久久久久久久人妻蜜臀av| 在线免费观看的www视频| 久热久热在线精品观看| 国产三级在线视频| 简卡轻食公司| 日韩欧美 国产精品| 一区二区三区免费毛片| h日本视频在线播放| 最近视频中文字幕2019在线8| 十八禁国产超污无遮挡网站| 国产精品一及| 亚洲国产精品国产精品| 亚洲怡红院男人天堂| 国产精品伦人一区二区| 亚洲精品乱码久久久久久按摩| 日本一二三区视频观看| 国产熟女欧美一区二区| 一级毛片我不卡| 天天一区二区日本电影三级| 日本av手机在线免费观看| 看黄色毛片网站| 国产精品国产三级国产av玫瑰| 国产精品久久电影中文字幕| www日本黄色视频网| 麻豆精品久久久久久蜜桃| 两个人的视频大全免费| 国产片特级美女逼逼视频| 三级男女做爰猛烈吃奶摸视频| 狠狠狠狠99中文字幕| 免费av观看视频| 国产亚洲91精品色在线| 欧美97在线视频| 国产精品久久久久久精品电影| 午夜福利成人在线免费观看| 91av网一区二区| 热99re8久久精品国产| 国产免费福利视频在线观看| 免费观看的影片在线观看| 免费无遮挡裸体视频| 精品免费久久久久久久清纯| 国产精品不卡视频一区二区| 中文字幕熟女人妻在线| 精品久久久久久久人妻蜜臀av| 国产男人的电影天堂91| 久久久a久久爽久久v久久| 久久久久精品久久久久真实原创| 成年av动漫网址| 婷婷色综合大香蕉| 国产精品久久久久久久久免| 日本色播在线视频| 国产精品伦人一区二区| 日本午夜av视频| 久久精品久久久久久久性| 亚洲中文字幕一区二区三区有码在线看| av卡一久久| 黄色配什么色好看| 岛国在线免费视频观看| 色综合色国产| 亚洲国产精品成人综合色| 不卡视频在线观看欧美| 亚洲欧美日韩卡通动漫| 网址你懂的国产日韩在线| 最近中文字幕高清免费大全6| 欧美日韩国产亚洲二区| 我要看日韩黄色一级片| 91精品一卡2卡3卡4卡| 男女边吃奶边做爰视频| 秋霞伦理黄片| 一级毛片电影观看 | 国产伦精品一区二区三区视频9| 老司机福利观看| 欧美又色又爽又黄视频| 青春草国产在线视频| 亚洲精品日韩在线中文字幕| 国产精品一二三区在线看| 2022亚洲国产成人精品| 日韩一区二区三区影片| 亚洲乱码一区二区免费版| 亚洲人与动物交配视频| 汤姆久久久久久久影院中文字幕 | 亚洲美女搞黄在线观看| 菩萨蛮人人尽说江南好唐韦庄 | av免费观看日本| 国产精品伦人一区二区| 国产真实伦视频高清在线观看| 日本黄色片子视频| 成人三级黄色视频| 一卡2卡三卡四卡精品乱码亚洲| 精品人妻视频免费看| 国产单亲对白刺激| 欧美bdsm另类| 亚洲国产精品国产精品| 97超视频在线观看视频| 老女人水多毛片| 久久这里只有精品中国| 老师上课跳d突然被开到最大视频| 美女国产视频在线观看| 真实男女啪啪啪动态图| АⅤ资源中文在线天堂| 小蜜桃在线观看免费完整版高清| 欧美性猛交╳xxx乱大交人| 国产伦理片在线播放av一区| 国产黄片视频在线免费观看| 白带黄色成豆腐渣| 一边摸一边抽搐一进一小说| 亚洲中文字幕一区二区三区有码在线看| 大香蕉97超碰在线| 国产爱豆传媒在线观看| 麻豆乱淫一区二区| 91aial.com中文字幕在线观看| 桃色一区二区三区在线观看| 婷婷色综合大香蕉| 国产在线一区二区三区精 | 欧美成人精品欧美一级黄| av免费观看日本| 高清视频免费观看一区二区 | 久热久热在线精品观看| 久久这里只有精品中国| av在线老鸭窝| 变态另类丝袜制服| 看黄色毛片网站| av在线天堂中文字幕| 亚洲国产日韩欧美精品在线观看| 免费观看性生交大片5| 99久久九九国产精品国产免费| 天天躁日日操中文字幕| 成年女人看的毛片在线观看| 成人av在线播放网站| 又粗又爽又猛毛片免费看| 久久久久性生活片| 美女大奶头视频| 亚洲av电影不卡..在线观看| 国产精品一及| 国产欧美日韩精品一区二区| 听说在线观看完整版免费高清| 欧美人与善性xxx| 久久6这里有精品| 网址你懂的国产日韩在线| 99国产精品一区二区蜜桃av| 色网站视频免费| 毛片女人毛片| 国内少妇人妻偷人精品xxx网站| 国产中年淑女户外野战色| 村上凉子中文字幕在线| 亚洲久久久久久中文字幕| 亚洲av.av天堂| 精品久久久噜噜| 久久久久久国产a免费观看| 欧美高清成人免费视频www| 亚洲av一区综合| 看免费成人av毛片| 97热精品久久久久久| 日本-黄色视频高清免费观看| 在线免费十八禁| 午夜福利在线在线| 午夜福利成人在线免费观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲aⅴ乱码一区二区在线播放| 波野结衣二区三区在线| 99久国产av精品国产电影| 国产高潮美女av| 国产黄色视频一区二区在线观看 | 国产成人福利小说| 99久久精品一区二区三区| 淫秽高清视频在线观看| 黄色配什么色好看| 春色校园在线视频观看| 精品人妻偷拍中文字幕| 久热久热在线精品观看| 色吧在线观看| 赤兔流量卡办理| 少妇裸体淫交视频免费看高清| 看免费成人av毛片| 午夜福利视频1000在线观看| 午夜日本视频在线| 午夜精品国产一区二区电影 | 日本免费在线观看一区| 午夜激情欧美在线| 国产不卡一卡二| 观看免费一级毛片| 国产视频首页在线观看| 国产精品一区二区三区四区免费观看| 国产亚洲精品久久久com| 中文天堂在线官网| 免费看日本二区| 亚洲怡红院男人天堂| 夜夜爽夜夜爽视频| 一级毛片久久久久久久久女| 成人亚洲欧美一区二区av| 精品一区二区三区视频在线| 亚洲成人精品中文字幕电影| av黄色大香蕉| 久久精品综合一区二区三区| 黄色配什么色好看| 午夜精品在线福利| 亚洲乱码一区二区免费版| 日本三级黄在线观看| 国产亚洲91精品色在线| 尾随美女入室| 亚洲高清免费不卡视频| 日本免费一区二区三区高清不卡| 国产一级毛片在线| 国产精品一及| 午夜视频国产福利| 久久久欧美国产精品| 久久精品久久精品一区二区三区| 久久精品91蜜桃| 亚洲成人久久爱视频| 少妇人妻一区二区三区视频| 国产精品久久久久久久久免| eeuss影院久久| 欧美成人精品欧美一级黄| 欧美日韩综合久久久久久| 国产精品av视频在线免费观看| 久久人人爽人人片av| 亚洲欧美日韩无卡精品| 精品99又大又爽又粗少妇毛片| 高清视频免费观看一区二区 | 尾随美女入室| 精品免费久久久久久久清纯| 久久鲁丝午夜福利片| 亚洲欧美日韩卡通动漫| 国产精品99久久久久久久久| 黑人高潮一二区| 国产久久久一区二区三区| 26uuu在线亚洲综合色| 精品99又大又爽又粗少妇毛片| 婷婷色综合大香蕉| 一级av片app| 2022亚洲国产成人精品| 成人av在线播放网站| 亚洲av中文av极速乱| 桃色一区二区三区在线观看| av免费观看日本| 欧美3d第一页| 97人妻精品一区二区三区麻豆| 国产精品一二三区在线看| 亚洲自拍偷在线| 男人狂女人下面高潮的视频| 国产一区二区在线观看日韩| 麻豆成人午夜福利视频| 中文亚洲av片在线观看爽| 天美传媒精品一区二区| av在线亚洲专区| 成人特级av手机在线观看| 成人美女网站在线观看视频| 少妇的逼好多水| 久久亚洲国产成人精品v| 18禁在线播放成人免费| 特级一级黄色大片| 视频中文字幕在线观看| 亚洲无线观看免费| 免费观看性生交大片5| 一个人看的www免费观看视频| 久久久久九九精品影院| 成年女人看的毛片在线观看| 成人性生交大片免费视频hd| 久久欧美精品欧美久久欧美| 亚洲一区高清亚洲精品| 日韩欧美精品v在线| 欧美丝袜亚洲另类| 亚洲真实伦在线观看| 婷婷色av中文字幕| 在线观看一区二区三区| 特大巨黑吊av在线直播| 亚洲精品日韩av片在线观看| 国产精品蜜桃在线观看| 久久精品国产自在天天线| 欧美高清性xxxxhd video| 色综合亚洲欧美另类图片| 日韩欧美国产在线观看| 五月玫瑰六月丁香| 国产精品野战在线观看| 久久久精品欧美日韩精品| 亚洲在久久综合| 国产成人a∨麻豆精品| 精品国产一区二区三区久久久樱花 | 久久精品夜夜夜夜夜久久蜜豆| 色尼玛亚洲综合影院| 日韩欧美三级三区| 久久精品熟女亚洲av麻豆精品 | 亚洲国产精品专区欧美| 99热这里只有是精品50| 精品久久国产蜜桃| 久久久久免费精品人妻一区二区| 搡老妇女老女人老熟妇| 国产成人freesex在线| 成人午夜高清在线视频| 亚洲成色77777| 男人的好看免费观看在线视频| 亚洲电影在线观看av| 99久久成人亚洲精品观看| 99久久精品热视频| 一级爰片在线观看| 国产精品一及| 亚洲在线自拍视频| 中文字幕熟女人妻在线| 黄色一级大片看看| 国产在视频线精品| 久久久亚洲精品成人影院| 少妇熟女欧美另类| 亚洲欧美精品专区久久| 一卡2卡三卡四卡精品乱码亚洲| 欧美bdsm另类| 国产亚洲av片在线观看秒播厂 | 午夜福利在线在线| 99热全是精品| 国产一区二区三区av在线| 亚洲不卡免费看| 黑人高潮一二区| 国产在线男女| 最近中文字幕高清免费大全6| h日本视频在线播放| 99久国产av精品国产电影| 又粗又爽又猛毛片免费看| 干丝袜人妻中文字幕| 久久99热这里只有精品18| 日韩视频在线欧美| 国产精品精品国产色婷婷| 成年版毛片免费区| 2021少妇久久久久久久久久久| 22中文网久久字幕| 麻豆av噜噜一区二区三区| 搡女人真爽免费视频火全软件| 国产成人精品一,二区| 中文字幕熟女人妻在线| 精品酒店卫生间| 日本欧美国产在线视频| 国产极品天堂在线| 精品久久久久久成人av| 日韩强制内射视频| 日本黄大片高清| av线在线观看网站| 观看美女的网站| 久久精品人妻少妇| 国产色爽女视频免费观看| 一二三四中文在线观看免费高清| 欧美又色又爽又黄视频| 国产色婷婷99|