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

    零重力條件下低溫射流抑制大尺寸液氫儲(chǔ)罐熱分層的數(shù)值研究1)

    2021-05-30 02:41:22趙建福胡文瑞
    力學(xué)學(xué)報(bào) 2021年4期
    關(guān)鍵詞:罐體熱帶射流

    郭 斌 趙建福 李 凱,3) 胡文瑞

    ?(中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

    ?(中國(guó)科學(xué)院力學(xué)研究所微重力重點(diǎn)實(shí)驗(yàn)室,北京 100190)

    ??(中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100190)

    引言

    以液氫為代表的低溫推進(jìn)劑在航天領(lǐng)域有著廣泛的應(yīng)用,但是因?yàn)槠浞悬c(diǎn)低、易蒸發(fā)的特性,在地面停放、發(fā)射過程以及在軌運(yùn)行等各個(gè)階段都極易發(fā)生汽化,引起儲(chǔ)罐內(nèi)部熱分層和自增壓,進(jìn)而威脅飛行器系統(tǒng)結(jié)構(gòu)安全.在軌運(yùn)行期間,推進(jìn)劑儲(chǔ)罐會(huì)受到各種空間輻射的影響,比如太陽輻射、地球反照熱外流、地球紅外輻射熱流以及空間黑背景輻射熱流等等.這些熱源的存在都會(huì)使得推進(jìn)劑儲(chǔ)罐因受熱不均勻在局部發(fā)生汽化現(xiàn)象.傳統(tǒng)的排氣方法能有效降低儲(chǔ)罐內(nèi)部壓強(qiáng),但一定程度上造成了推進(jìn)劑的浪費(fèi),使發(fā)射成本增加,同時(shí)制約飛行器運(yùn)載能力,不利于長(zhǎng)期在軌飛行及深空探測(cè)活動(dòng).20 世紀(jì)末,美國(guó)NASA 的研究者提出了零蒸發(fā)存儲(chǔ)(zero boil-off,ZBO)的概念.ZBO 技術(shù)是通過將被動(dòng)絕熱技術(shù)與主動(dòng)制冷技術(shù)有機(jī)結(jié)合起來,實(shí)現(xiàn)低溫推進(jìn)劑無損儲(chǔ)存的技術(shù),其中被動(dòng)絕熱技術(shù)主要是使用絕熱材料減少外部熱量的滲入.Hastings 等[1]發(fā)現(xiàn)對(duì)于需要長(zhǎng)期在軌運(yùn)行的儲(chǔ)罐而言,單純使用被動(dòng)絕熱技術(shù)會(huì)增加儲(chǔ)罐質(zhì)量,進(jìn)而降低航天器的靈活性,因此在使用絕熱材料的基礎(chǔ)上嘗試主動(dòng)制冷技術(shù)的運(yùn)用對(duì)于長(zhǎng)期航天探索工程具有重要意義.主動(dòng)制冷技術(shù)主要有兩種,一種是通過熱傳導(dǎo)裝置將罐內(nèi)熱量傳導(dǎo)到罐外,另一種則是通過機(jī)械攪拌或者低溫射流的方式加強(qiáng)儲(chǔ)罐內(nèi)部流體對(duì)流.熱傳導(dǎo)裝置作用的范圍相對(duì)較小,機(jī)械攪拌的實(shí)現(xiàn)工藝更加復(fù)雜,因此通過施加儲(chǔ)罐內(nèi)部低溫射流的方式進(jìn)行壓強(qiáng)控制成為國(guó)內(nèi)外相關(guān)領(lǐng)域?qū)W者的研究熱點(diǎn).

    通過低溫射流實(shí)現(xiàn)熱分層和自增壓控制的研究可以追溯到20 世紀(jì)70 年代,Poth 等[2]對(duì)比分析了各種熱分層消除裝置的特點(diǎn),認(rèn)為軸向低溫射流具有很高的流體混合性能并且輕便靈活,增壓控制易于實(shí)現(xiàn);Merte 等[3]研究了軸對(duì)稱的圓柱形儲(chǔ)罐中單組分工質(zhì)氣相和液相之間的傳熱和傳質(zhì)過程,其中氣相采用理想氣體假設(shè),將控制方程轉(zhuǎn)換為有限差分形式,進(jìn)行數(shù)值計(jì)算得到了多種計(jì)算工況下的壓力演化曲線.Audelott[4]在Lewis Research Center 5-10 s零重力裝置中實(shí)驗(yàn)研究了向直徑為10 cm的圓柱形儲(chǔ)罐中軸向噴射乙醇時(shí)產(chǎn)生的液體流動(dòng)模式,發(fā)現(xiàn)儲(chǔ)罐內(nèi)部流場(chǎng)分布與儲(chǔ)罐幾何形狀、低溫射流速度、儲(chǔ)罐填充比以及低溫射流位置有關(guān);Lin 等[5]先后實(shí)驗(yàn)研究了體積為0.144 m3的小型儲(chǔ)罐和體積為4.89 m3大型儲(chǔ)罐使用低溫射流混合裝置實(shí)現(xiàn)增壓控制的效果,實(shí)驗(yàn)表明低溫射流可以有效控制罐內(nèi)增壓,但是小罐模型得出的低溫射流時(shí)間和增壓速率之間的關(guān)系無法很好預(yù)測(cè)大罐模型中兩者之間的相對(duì)關(guān)系;Panzarella 等[6]研究了常重力條件下小型儲(chǔ)罐內(nèi)部氣液兩相的自增壓情況,證實(shí)了增強(qiáng)儲(chǔ)罐內(nèi)部對(duì)流可以抑制熱分層現(xiàn)象進(jìn)而實(shí)現(xiàn)增壓控制;Mukka 等[7]研究了常重力條件下低溫射流方式對(duì)消除熱分層效果的影響,認(rèn)為不同低溫射流條件導(dǎo)致的儲(chǔ)罐內(nèi)部流場(chǎng)的差異對(duì)儲(chǔ)罐內(nèi)溫度場(chǎng)和速度場(chǎng)的分布有著重要影響.在認(rèn)識(shí)到低溫射流可以有效消除熱分層現(xiàn)象后,科研工作者建立了很多罐體系統(tǒng)研究如何優(yōu)化低溫射流消除熱分層的效果.Ho 等[8]建立了一個(gè)泵--熱管模型,通過泵裝置促進(jìn)儲(chǔ)罐內(nèi)部流體循環(huán),不斷把高溫流體噴向熱管,再由熱管將熱量導(dǎo)出.通過對(duì)該模型的三維數(shù)值模擬發(fā)現(xiàn)通過增大低溫射流速度可以有效消除儲(chǔ)罐內(nèi)部的高溫區(qū)域;隨后Ho 等[9]利用數(shù)值模擬方法研究了軸向圓盤多孔低溫射流噴嘴結(jié)構(gòu),發(fā)現(xiàn)在保持入射速度不變的情況下,加大孔口直徑也可以顯著提升熱分層消除效果,此外射流噴嘴位置對(duì)熱分層消除效果有著顯著影響,并對(duì)圓盤射流噴嘴系統(tǒng)進(jìn)行了參量分析,研究了入口管直徑、射流噴嘴位置和射流噴嘴直徑對(duì)系統(tǒng)制冷效果的影響.Ho 等[10]又對(duì)之前建立的泵--熱管模型進(jìn)行了深入研究,分析了噴射間隔和噴管長(zhǎng)度等參量對(duì)消除熱分層效果的影響.Belmedani 等[11]開展了不同熱通量的液氮儲(chǔ)罐熱分層的實(shí)驗(yàn)研究,熱通量越大,自由面蒸發(fā)速率越大,自由面和液相區(qū)溫差越明顯,較好地揭示了蒸發(fā)機(jī)理,通過研究低溫儲(chǔ)罐中液體的速度和溫度特性,認(rèn)為強(qiáng)熱流誘發(fā)熱分層現(xiàn)象.Zilliac 等[12]提出了一種平衡熱力學(xué)模型來預(yù)測(cè)儲(chǔ)罐增壓過程,該模型適用于具有高揮發(fā)性的低溫推進(jìn)劑,通過與實(shí)測(cè)結(jié)果進(jìn)行對(duì)比,證實(shí)了該模型的準(zhǔn)確性.Grayson 等[13]模擬了常重力情況下通過外部加熱實(shí)現(xiàn)增壓和通過熱力學(xué)排氣實(shí)現(xiàn)減壓的儲(chǔ)罐數(shù)值計(jì)算模型,假設(shè)液相不可壓縮,其密度僅是溫度的函數(shù).該模型的壓力和溫度預(yù)測(cè)結(jié)果與測(cè)試中的傳感器測(cè)量結(jié)果相比較有很高的一致性.Barsi 等[14]采取了與原有的集總熱力學(xué)模型不同的氣液兩相模型,在該模型中,液相和汽相的守恒方程都得到了求解.研究中為了簡(jiǎn)化計(jì)算將兩相都視為不可壓縮的,并使用積分方法來求解交界面的質(zhì)量交換.通過將儲(chǔ)罐壓力上升的數(shù)值計(jì)算與以往計(jì)算模型的結(jié)果進(jìn)行比較,證實(shí)了該模型的準(zhǔn)確性.Lopez 等[15]提出了一種新的計(jì)算流體動(dòng)力學(xué)模型,模擬了重力作用下橢球形液氫儲(chǔ)罐在外部加熱情形下的壓力控制過程.壓力控制的實(shí)現(xiàn)由位于儲(chǔ)罐內(nèi)部的軸向射流熱動(dòng)力排氣系統(tǒng)(TVS)提供,該系統(tǒng)向儲(chǔ)罐內(nèi)注入冷流體,使液體發(fā)生混合以降低儲(chǔ)罐內(nèi)部壓力.該軸對(duì)稱模型的計(jì)算采用商業(yè)軟件FLOW-3D 進(jìn)行計(jì)算,其中定量模型驗(yàn)證采用1999 年在馬歇爾航天飛行中心進(jìn)行的工程檢驗(yàn)測(cè)試結(jié)果,計(jì)算結(jié)果表明模型預(yù)測(cè)的自增壓速率和流體溫度與試驗(yàn)數(shù)據(jù)吻合良好.這項(xiàng)研究提升了當(dāng)時(shí)用于實(shí)現(xiàn)低溫壓力控制的CFD 建模能力,并為開發(fā)基于CFD 的空間硬件設(shè)計(jì)提供了參考.Barsi 等[16]提出了一種兩相CFD 模型,該模型描述了在常重力條件下,部分填充的LH2 儲(chǔ)罐的自增壓行為,并利用已有的不同填充水平下的實(shí)驗(yàn)數(shù)據(jù),對(duì)模型的預(yù)測(cè)能力進(jìn)行了評(píng)估.評(píng)估結(jié)果表明該模型的預(yù)測(cè)結(jié)果與實(shí)驗(yàn)測(cè)得的壓力曲線吻合較好.Kumar 等[17]計(jì)算了不同縱橫比的大型液氫儲(chǔ)罐中蒸發(fā)對(duì)儲(chǔ)罐熱分層的影響.他們采用均勻兩相模型,分別求解了氣相和液相各自的守恒方程,氣液界面的蒸發(fā)是通過傳質(zhì)源項(xiàng)來實(shí)現(xiàn)的,隨著儲(chǔ)罐的縱橫比增加,分層的程度會(huì)逐漸增加.Oliveira 等[18]建立了一個(gè)儲(chǔ)罐熱分層模型,其中包括了平臺(tái)的調(diào)節(jié)旋轉(zhuǎn).該模型可用于評(píng)估軸向加速度、自旋速率、熱通量以及儲(chǔ)罐幾何形狀等因素對(duì)推進(jìn)劑儲(chǔ)罐內(nèi)部熱分層現(xiàn)象的影響.Li 等[19]通過實(shí)驗(yàn)和數(shù)值方法研究了多層絕緣低溫儲(chǔ)罐中真空損耗引起的液氮熱分層的瞬態(tài)過程,他們認(rèn)為真空損耗可以迅速導(dǎo)致熱分層.在實(shí)驗(yàn)中,利用熱電偶得到了實(shí)驗(yàn)儲(chǔ)罐中液體溫度的分布和演變.通過建立二維模型進(jìn)行數(shù)值模擬計(jì)算,模擬了儲(chǔ)罐內(nèi)部熱分層的形成和減弱過程以及液體溫度場(chǎng)分布,并將數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行了比較.研究還表明,罐內(nèi)的兩相流動(dòng)對(duì)熱分層起著重要作用.Wang 等[20]通過數(shù)值模擬系統(tǒng)對(duì)比了微重力條件下由熱管和噴管組成的ZBO 系統(tǒng)中,噴管數(shù)量、噴管出流方向和熱管蒸發(fā)器熱傳遞效率對(duì)增壓控制性能的影響.Liu 等[21]通過在Rahman 模型的基礎(chǔ)上增加了向下的導(dǎo)管,使儲(chǔ)罐底部也發(fā)生強(qiáng)制對(duì)流,提升ZBO 性能.Liu 等[22]進(jìn)一步通過正交試驗(yàn)設(shè)計(jì),研究了導(dǎo)管出口到罐底的距離、噴管的半徑、射流噴嘴深度、導(dǎo)管半徑和環(huán)形射流噴嘴直徑等參量對(duì)儲(chǔ)罐設(shè)計(jì)的影響.Roh 等[23]使用商業(yè)軟件FLUENT 對(duì)壓縮液化天然氣儲(chǔ)罐中的瞬態(tài)自然對(duì)流進(jìn)行了數(shù)值計(jì)算.計(jì)算結(jié)果表明,蒸發(fā)氣體的生成在很大程度上取決于儲(chǔ)罐內(nèi)部垂直溫度的分布,而儲(chǔ)罐內(nèi)部垂直方向的溫度分布受儲(chǔ)罐自增壓過程的影響.同時(shí)作者對(duì)儲(chǔ)罐壓力、罐體尺寸和增壓過程對(duì)蒸發(fā)氣體生成的影響進(jìn)行了量化研究.Wang 等[24]采用計(jì)算流體動(dòng)力學(xué)模型,對(duì)液氫儲(chǔ)罐的加壓排放過程進(jìn)行了數(shù)值模擬,將壁面區(qū)域和流體區(qū)域同時(shí)考慮為計(jì)算域,采用低雷諾數(shù)k-ε 模型來處理流體和壁面的熱交換效應(yīng).該模型還考慮了氣液相變,并將數(shù)值計(jì)算結(jié)果與已有實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,認(rèn)為該CFD 模型在增壓計(jì)算過程中具有良好的適應(yīng)性.通過該模型可以獲得增壓所需氣體量、儲(chǔ)罐內(nèi)部壓力變化過程、儲(chǔ)罐內(nèi)部溫度分布等詳細(xì)特征.作者還分析了相變效應(yīng)和儲(chǔ)罐結(jié)構(gòu)對(duì)增壓性能的影響,計(jì)算結(jié)果表明,氣液相變對(duì)增壓行為影響不大.隨后Wang 等[25]又建立了一種計(jì)算流體力學(xué)模型,該模型可同時(shí)考慮罐內(nèi)的熱交換和外部空氣動(dòng)力加熱,并對(duì)低溫儲(chǔ)罐在排放過程中的瞬態(tài)熱狀態(tài)和增壓性能進(jìn)行了研究.該模型中計(jì)算域不僅包括了流體和罐壁區(qū)域,還包括了泡沫絕熱層區(qū)域.作者將該模型的計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較,兩者具有很高的一致性.然后利用該模型對(duì)儲(chǔ)罐的增壓排放過程進(jìn)行了預(yù)測(cè),得到了其熱力學(xué)行為、增壓行為過程.Daigle 等[26]提出了一種描述液氫低溫儲(chǔ)罐中自然對(duì)流的溫度分層的動(dòng)態(tài)模型,用MATLAB 實(shí)現(xiàn)了一個(gè)通用的低溫儲(chǔ)罐的溫度分層現(xiàn)象模擬程序,該模型可以模擬常重力以及重力增加和減少條件下儲(chǔ)罐內(nèi)部溫度分層結(jié)果.Fu 等[27]采用數(shù)值計(jì)算的方法,對(duì)部分填充的圓柱肋式液氫儲(chǔ)罐在不同肋間距比下的自增壓過程進(jìn)行了研究.利用商業(yè)軟件FLUENT 進(jìn)行計(jì)算,選取了流體體積法和相變模型,并通過自定義函數(shù)對(duì)模型進(jìn)行適當(dāng)修改,建立了用于預(yù)測(cè)儲(chǔ)罐內(nèi)部流體流動(dòng)和傳熱的求解方案.隨后Fu 等[28]對(duì)低溫儲(chǔ)罐在微重力條件下的蒸發(fā)及其對(duì)蒸汽壓的影響進(jìn)行了數(shù)值研究,研究了表面張力、氣泡接觸角和重力等因素的影響.李佳超等[29]以液氮為研究工質(zhì)利用透明玻璃搭建的低溫儲(chǔ)罐自增壓實(shí)驗(yàn)系統(tǒng),研究了自增壓過程壓力和溫度的變化規(guī)律及填充比對(duì)壓力和溫度變化過程的影響.2018 年王夕等[30]使用FLUENT 軟件進(jìn)行數(shù)值模擬,對(duì)比研究了4 種相變模型對(duì)微重力環(huán)境中液氫推進(jìn)劑受熱蒸發(fā)過程的影響,并與國(guó)外探空火箭試驗(yàn)數(shù)據(jù)進(jìn)行比較,分析對(duì)比了4 種相變模型各自的優(yōu)劣.李鵬等[31]以液氫、液氧等低溫推進(jìn)劑為研究對(duì)象,開展控制低溫推進(jìn)劑在儲(chǔ)罐內(nèi)部因受熱而蒸發(fā)的現(xiàn)象,通過對(duì)不同條件下的蒸發(fā)量控制效果的對(duì)比分析,為今后航天器推進(jìn)劑儲(chǔ)罐的設(shè)計(jì)提供參考依據(jù).郭志釩等[32]分析了高壓儲(chǔ)氫、低溫液態(tài)儲(chǔ)氫、金屬氫化物儲(chǔ)氫等3 種儲(chǔ)氫方式各自的優(yōu)缺點(diǎn)與發(fā)展現(xiàn)狀,為未來發(fā)展提供新的思路.馬原等[33]采用CFD 方法建立兩相流模型,對(duì)微重力條件下在液氫儲(chǔ)罐內(nèi)噴射過冷流體實(shí)現(xiàn)儲(chǔ)罐降壓的過程開展數(shù)值模擬研究,對(duì)比計(jì)算了不同噴射區(qū)域、噴射流量、噴射速度等對(duì)罐內(nèi)溫度場(chǎng)分布與壓力變化的影響,認(rèn)為氣--液相區(qū)噴射降壓性能優(yōu)于單獨(dú)區(qū)域噴射,液相區(qū)噴射降壓效果最弱.2019 年王舜浩等[34]利用數(shù)值方法研究了液氫縮比儲(chǔ)罐內(nèi)部流體的蒸發(fā)性質(zhì),通過與已有的實(shí)驗(yàn)數(shù)據(jù)對(duì)比構(gòu)建了基于VOF 兩相流模型以及(level-set)界面跟蹤方法的儲(chǔ)罐內(nèi)部流體流動(dòng)和相變傳熱傳質(zhì)模型框架,為模擬液氫儲(chǔ)罐地面停放階段的熱物理過程提供了參考.Zuo 等[35]建立了一個(gè)噴嘴可旋轉(zhuǎn)的三維液氫儲(chǔ)罐模型,利用CFD 方法研究了低溫射流對(duì)零重力條件下熱分層的影響,結(jié)果表明可旋轉(zhuǎn)的噴嘴可以很有效的抑制罐體內(nèi)部的熱分層現(xiàn)象.Guo 等[36]采用二維縮比模型儲(chǔ)罐研究了低溫射流消除熱分層現(xiàn)象中噴頭形狀、位置以及射流速度對(duì)消除效果的影響,認(rèn)為圓形射流噴嘴相較于半球形射流噴嘴消除效果更好.Zhang 等[37]研究了一種效率高,界面清晰,適用于三維模型的計(jì)算氣液兩相界面遷移特性的歐拉運(yùn)動(dòng)界面追蹤方法,該方法將‘米’狀相鄰單元Youngs 方法用于運(yùn)動(dòng)界面重構(gòu),將Youngs-VOF 和水平集通過幾何方法耦合,提高運(yùn)動(dòng)界面精度,克服了VOF 和水平集方法存在的缺陷,避免了利用高階導(dǎo)數(shù)本身的穩(wěn)定性去求解水平集對(duì)流方程和距離函數(shù)方程,為今后將熱分層現(xiàn)象的研究拓展到三維以及兩相流提供了思路.

    由上可知國(guó)內(nèi)外學(xué)者已經(jīng)在該問題上進(jìn)行了很多有益嘗試并且取得了一定進(jìn)展,但是以往的研究大多假定儲(chǔ)罐壁表面均勻漏熱,因此在整個(gè)罐壁上采用均勻熱量邊界條件,并且在計(jì)算時(shí)認(rèn)為邊界漏熱和低溫射流是同時(shí)進(jìn)行的,通過對(duì)比低溫射流一段時(shí)間后儲(chǔ)罐內(nèi)部的最高溫度判斷低溫射流系統(tǒng)增壓控制性能的優(yōu)劣.儲(chǔ)罐罐體一般為金屬材料并且在表面覆蓋多層絕熱材料,盡管如此,實(shí)際使用中還是無法完全消除熱量滲入.罐體表面因?yàn)榻Y(jié)構(gòu)裝置等的差異,通過儲(chǔ)罐壁滲入液體推進(jìn)劑的熱量并不是均勻分布的.相對(duì)于罐體其他位置,通過易漏熱的區(qū)域滲入罐體的熱量對(duì)推進(jìn)劑的熱分層有著更大影響.目前對(duì)于儲(chǔ)罐局部漏熱導(dǎo)致的熱分層現(xiàn)象的低溫射流消除系統(tǒng)的設(shè)計(jì)分析研究還不夠充分,因此本文主要研究了大尺寸儲(chǔ)罐在局部區(qū)域漏熱情形下儲(chǔ)罐內(nèi)部出現(xiàn)的明顯熱分層,并對(duì)比研究了低溫射流噴嘴的形狀及其在儲(chǔ)罐內(nèi)部的相對(duì)位置等因素對(duì)利用低溫射流消除熱分層效果的影響.

    1 計(jì)算模型

    1.1 物理模型

    熱分層控制系統(tǒng)模型由罐體結(jié)構(gòu)和入射結(jié)構(gòu)組成,如圖1 所示.低溫射流由儲(chǔ)罐內(nèi)部的射流噴嘴射出,與儲(chǔ)罐內(nèi)部流體發(fā)生混合和熱交換;另一方面儲(chǔ)罐內(nèi)部流體從罐體出流口導(dǎo)出,經(jīng)過儲(chǔ)罐外制冷系統(tǒng)(本文忽略)處理成低溫流體重新流回儲(chǔ)罐內(nèi)部,形成一個(gè)消除熱分層的閉環(huán)機(jī)制.本文主要研究罐體局部漏熱的情況,并考慮罐體出流口作為儲(chǔ)罐與外部系統(tǒng)的銜接段易發(fā)生漏熱的工程經(jīng)驗(yàn),假定罐體表面的條狀區(qū)域及出口為漏熱帶(見圖2),罐體其他區(qū)域視為絕熱情況.本文數(shù)值模擬采用二維軸對(duì)稱模型,儲(chǔ)罐各部分具體尺寸見圖2;采用液氫為研究工質(zhì),填充率保持為100%.本文通過改變射流噴嘴的形狀以及射流噴嘴在儲(chǔ)罐內(nèi)部的相對(duì)位置,研究不同低溫射流條件對(duì)儲(chǔ)罐內(nèi)部流體流動(dòng)和溫度分布時(shí)空演化過程的影響.

    圖1 罐體結(jié)構(gòu)及射流噴嘴樣式示意圖Fig.1 Schematic diagram of tank structure and jet nozzle pattern

    圖2 儲(chǔ)罐幾何尺寸及監(jiān)測(cè)線示意圖Fig.2 Schematic diagram of the cryogenic storage tank and inspection line

    1.2 數(shù)學(xué)模型

    1.2.1 控制方程

    假設(shè)流體是不可壓縮的并且具有恒定的熱物理特性.在微重力條件下忽略重力和浮力的影響,因而在圓柱坐標(biāo)系中整個(gè)計(jì)算域內(nèi)的質(zhì)量守恒、動(dòng)量守恒和能量守恒控制方程如下:

    其中,μ為動(dòng)力學(xué)黏性系數(shù),由工質(zhì)物理特性決定,μt為湍流黏性系數(shù),是空間坐標(biāo)的函數(shù),取決于流動(dòng)狀態(tài)而非物性參數(shù),μt為湍流黏性系數(shù),是由湍動(dòng)能k及耗散率確定ε,表達(dá)式如下所示

    其中?;?shù)Cμ=0.09.

    (3)能量方程

    其中prt≈0.85.

    1.2.2 湍流模型

    采用雷諾時(shí)均模型對(duì)湍流進(jìn)行模擬.為了使方程封閉,引入了新的未知量湍動(dòng)能k以及耗散率ε.Abid[38]對(duì)受限沖擊低溫射流以及受限對(duì)沖低溫射流進(jìn)行了數(shù)值模擬實(shí)驗(yàn)研究,認(rèn)為采用AB 型低雷諾數(shù)k-ε 湍流模型可以很好的對(duì)低溫射流進(jìn)行數(shù)值模擬,因此本文采用AB 低雷諾數(shù)k-ε 湍流模型,其kε 模式方程以及渦黏性如下所示

    1.2.3 邊界條件以及初始條件

    在入射口截面上

    在對(duì)稱軸上

    在罐壁漏熱帶上

    在罐體其他壁面以及導(dǎo)管壁面上

    初始時(shí)刻,罐內(nèi)流體溫度

    出流口處邊界為自由出流條件.

    工質(zhì)的物理特性根參考了Guo 等[36]的研究,選取溫度為20 K 時(shí)的各項(xiàng)參數(shù)值如表1 所示.

    表1 工質(zhì)的熱物理特性參數(shù)值Table 1 Thermal physical property values of liquid hydrogen

    2 數(shù)值模擬

    上述控制方程的求解通過國(guó)際通用的CFD 模擬軟件Fluent 17.0 進(jìn)行計(jì)算,整個(gè)計(jì)算域通過ICEM生成網(wǎng)格,在壁面附近網(wǎng)格加密,使用SIMPLE(semi-implicit method for pressure-linked equations)算法求解,迭代過程采用絕對(duì)收斂標(biāo)準(zhǔn)進(jìn)行控制,能量項(xiàng)的殘差值設(shè)定為10?9,其他各項(xiàng)設(shè)為10?5.本文采用3 種網(wǎng)格數(shù)目分別為71 757,110 249,174 124,開展了網(wǎng)格無關(guān)性檢驗(yàn).如圖2 所示,在計(jì)算域上溫度梯度較大的位置選取了水平監(jiān)測(cè)線和豎直監(jiān)測(cè)線.圖3 顯示了網(wǎng)格無關(guān)性檢驗(yàn)的結(jié)果,可以發(fā)現(xiàn)沿豎直監(jiān)測(cè)線3 種網(wǎng)格數(shù)下的溫度分布基本一致;沿水平監(jiān)測(cè)線兩種較細(xì)網(wǎng)格的溫度分布沒有明顯差異,但在靠近低溫射流噴嘴位置處與粗網(wǎng)格下的溫度分布有明顯差異.由此,之后的研究采用網(wǎng)格數(shù)量110 426 以同時(shí)兼顧計(jì)算準(zhǔn)確性和計(jì)算量.

    圖3 網(wǎng)格無關(guān)性檢驗(yàn)Fig.3 Grid independence test

    3 結(jié)果討論

    主要研究了圓形射流噴嘴(C)在儲(chǔ)罐內(nèi)部不同位置時(shí)利用低溫射流消除熱分層的效果,并選取了射流噴嘴所在典型位置與半球形射流噴嘴(H)進(jìn)行了對(duì)比.表2 列出了各工況進(jìn)行數(shù)值計(jì)算時(shí)的具體參數(shù)設(shè)置.

    表2 各計(jì)算工況具體參數(shù)Table 2 Details of computational cases

    3.1 圓形射流噴嘴

    采用圓形射流噴嘴進(jìn)行研究,其中低溫射流速度沿噴嘴表面均勻分布,將無射流情況下漏熱持續(xù)30 d(2 592 000 s)后儲(chǔ)罐內(nèi)部的溫度場(chǎng)分布作為射流開始的初始狀態(tài)(t=0 s).如圖4 所示,可以發(fā)現(xiàn)此時(shí)儲(chǔ)罐內(nèi)部流體以漏熱帶為中心出現(xiàn)了明顯的環(huán)形熱分層現(xiàn)象.射流噴嘴在儲(chǔ)罐內(nèi)部的相對(duì)位置對(duì)初始熱分布沒有明顯影響,4 種工況中的熱分布基本一致.

    圖4 各工況初始時(shí)刻溫度云圖分布Fig.4 Isothermals of the initial state(t=0 s)for the cryogenic jet

    對(duì)于在漏熱帶施加溫度邊界條件模擬儲(chǔ)罐壁面局部漏熱的情形,漏熱帶上的熱量傳輸速率可以表征漏熱帶附近的溫度分布.熱量傳輸速率指單位時(shí)間通過漏熱帶進(jìn)入罐體的熱量,其值越大說明漏熱帶附近區(qū)域溫度越低.圖5 給出了射流時(shí)間持續(xù)約1 200 s 后各工況條狀漏熱帶熱量傳輸速率隨時(shí)間演化圖.可以看出在不同時(shí)間段,低溫射流噴嘴與漏熱帶的相對(duì)位置對(duì)于漏熱帶附近高溫區(qū)域的熱量傳輸速率存在顯著影響.按照各工況漏熱帶熱量傳輸速率的差異大體可以將射流過程分為3 個(gè)階段.第一階段(0~100 s),對(duì)任意一種工況,由于低溫射流時(shí)間較短,通過射流進(jìn)入罐體內(nèi)部的冷流體都無法直接影響條狀漏熱帶附近區(qū)域.冷流體在罐內(nèi)作用區(qū)域非常有限,因此各工況條狀漏熱帶熱量傳輸速率曲線基本重合.圖6 給出了低溫射流進(jìn)行64.08 s 后各工況儲(chǔ)罐內(nèi)部的溫度云圖分布.可以看出低溫流體主要分布在射流噴嘴附近,距離條狀漏熱帶附近的高溫區(qū)域較遠(yuǎn).此階段儲(chǔ)罐內(nèi)抑制熱分層發(fā)展主要依靠罐內(nèi)流體由出流口經(jīng)過罐外冷卻系統(tǒng)實(shí)現(xiàn).此時(shí)出流口附近的的溫度存在明顯下降.圖7 給出了該時(shí)刻各工況的流場(chǎng)圖,可以看出施加低溫射流之后,會(huì)在射流噴嘴附近形成渦流,由于射流噴嘴位置的差異,渦流在儲(chǔ)罐內(nèi)部的位置有所不同,因此儲(chǔ)罐內(nèi)部流場(chǎng)存在差異,但總體而言在低溫射流初期,各工況抑制熱分層的效果差異不大.對(duì)比圖4 初始時(shí)刻各工況溫度云圖分布,可以發(fā)現(xiàn)條狀漏熱帶附近高溫區(qū)域的熱分層形態(tài)沒有發(fā)生明顯變化,而出口位置的熱分層明顯被消除.在第二階段(100~700 s),各工況低溫射流持續(xù)到424.08 s 時(shí)的溫度云圖(圖8)表明工況C1 和C2 中條狀漏熱帶附近熱分層形態(tài)與第一階段相比發(fā)生了輕微的變化.由射流噴嘴進(jìn)入到罐體內(nèi)部的冷流體(溫度介于16 ~18 K 的低溫流體)作用范圍進(jìn)一步擴(kuò)大,在罐體沿著軸線形成了冷流體柱.工況C3 和C4 中條狀漏熱帶附近高溫區(qū)域的熱分層形態(tài)相較于第一階段則發(fā)生了顯著的變化.冷流體已經(jīng)在罐底聚集并沿著儲(chǔ)罐壁面內(nèi)卷.漏熱帶附近高溫區(qū)域熱分層的形態(tài)和該時(shí)刻各工況中渦流的位置有很大的關(guān)系.總體而言,熱分層的延展方向與對(duì)應(yīng)位置的流體流動(dòng)方向是一致的.由于從低溫射流開始C3 和C4 兩種工況形成的渦流相對(duì)于C1 和C2 更靠近罐體底部,長(zhǎng)時(shí)間作用之后,靠近對(duì)稱軸的熱分層隨著順時(shí)針方向流動(dòng)的流體而向左延展.圖9 給出了低溫射流進(jìn)行到424.08 s 時(shí)的流場(chǎng)圖.可以看出,各工況渦流都沿對(duì)稱軸運(yùn)動(dòng)到罐體底部.C1 中形成的渦流在罐體內(nèi)部的流動(dòng)區(qū)域范圍相對(duì)于其他工況更廣,整體性更強(qiáng),能夠更大范圍帶動(dòng)罐內(nèi)流體的混合,結(jié)合圖5 可以發(fā)現(xiàn)這一時(shí)段,C1 中漏熱帶熱量傳輸速率高于其他工況,并且各工況整體上滿足隨著射流噴嘴伸入罐體內(nèi)部長(zhǎng)度的增加,漏熱帶熱量傳輸速率逐漸降低.這和圖10 給出的各工況儲(chǔ)罐內(nèi)部流體平均速度隨低溫射流時(shí)間演化曲線圖變化趨勢(shì)是一致的,射流初期各工況儲(chǔ)罐內(nèi)部流體平均速度都快速上升,但工況間沒有明顯差異.低溫射流持續(xù)100 s 之后,罐體內(nèi)部流體平均速度隨著射流噴嘴伸入罐體內(nèi)部長(zhǎng)度的增加而逐漸減小.圖11 給出了儲(chǔ)罐內(nèi)部流體平均溫度隨低溫射流時(shí)間的演化曲線,可以發(fā)現(xiàn)各工況儲(chǔ)罐內(nèi)部流體平均溫度的差異也主要是在低溫射流持續(xù)100 s 之后形成的,整體上各工況平均溫度都經(jīng)歷了先快速降低然后逐漸回升的過程,同樣由于儲(chǔ)罐內(nèi)部流場(chǎng)整體性的差異,在同一低溫射流時(shí)刻,隨著射流噴嘴伸入罐體內(nèi)部長(zhǎng)度的增大,罐體內(nèi)部流體平均溫度逐漸減小.

    圖5 圓形射流噴嘴漏熱帶熱量傳輸速率隨時(shí)間演化圖Fig.5 Evolution diagram of the heat transfer rate with cryogenic jet time in the circular jet nozzle

    圖6 各工況低溫射流64.08 s 時(shí)溫度云圖分布Fig.6 Isothermals of cases at t=64.08 s for the cryogenic jet

    圖7 各工況低溫射流64.08 s 時(shí)流場(chǎng)圖Fig.7 Flow field diagram at t=64.08 s for the cryogenic jet

    圖8 各工況低溫射流424.08 s 時(shí)溫度云圖分布Fig.8 Isothermals of cases at t=424.08 s for the cryogenic jet

    圖9 各工況低溫射流424.08 s 時(shí)流場(chǎng)圖Fig.9 Flow field diagram at t=424.08 s for the cryogenic jet

    在第三階段(700 ~1 200 s),由低溫射流持續(xù)1 024.08 s 時(shí)的溫度云圖(圖12)可以發(fā)現(xiàn),該時(shí)刻各工況熱分層的形態(tài)都發(fā)生了顯著變化,由最初低溫射流初始時(shí)刻的環(huán)形分層演變成不規(guī)則的帶狀分層.對(duì)比第二階段的溫度云圖分布(圖8),可以發(fā)現(xiàn)由射流噴嘴進(jìn)入罐體內(nèi)部的帶狀冷流體分布范圍并沒有繼續(xù)在罐體底部聚集或沿罐壁攀升,反而相對(duì)于上一階段有所減少,一方面是因?yàn)殡S著罐體內(nèi)部流體流速的增加(圖10),低溫流體很快被輸運(yùn)到罐體其他位置發(fā)生熱交換,無法沿對(duì)稱軸附近積聚延伸發(fā)展,另一方面隨著漏熱時(shí)間的延長(zhǎng)以及漏熱帶熱量傳輸速率的快速升高(圖5),罐內(nèi)流體溫度總體逐漸上升(圖11),射流噴嘴附近的流體溫度也有所上升,通過射流噴嘴進(jìn)入罐體的冷流體在射流噴嘴附近就會(huì)發(fā)生熱交換.事實(shí)上由罐體內(nèi)部流體平均溫度的變化曲線(圖11)表明,在低溫射流持續(xù)400 s 左右時(shí),各工況罐體內(nèi)部平均溫度都達(dá)到極小值,隨后溫度開始上升,從能量角度來看,表明在400 s 左右從儲(chǔ)罐內(nèi)部輸出的熱量和從外界環(huán)境輸入罐體內(nèi)部的熱量達(dá)到了平衡,從儲(chǔ)罐內(nèi)部輸出的熱量包括由射流噴嘴進(jìn)入罐體內(nèi)部的冷流體攜帶的熱量(負(fù)值)和通過儲(chǔ)罐出口出流的流體攜帶的熱量,而從外界環(huán)境輸入罐體內(nèi)部的熱量主要包括通過條狀漏熱帶和出口現(xiàn)階段漏熱帶進(jìn)入罐體內(nèi)部的熱量.在低溫射流初期,輸入罐體內(nèi)部的熱量值小于輸出罐體內(nèi)部的熱量,罐體內(nèi)部流體平均溫度逐漸降低,而低溫射流后期則反之.圖13 給出了各工況低溫射流1024.08 s 時(shí)流場(chǎng)圖,可以發(fā)現(xiàn)各工況形成的渦流流動(dòng)范圍都已經(jīng)擴(kuò)大到足夠覆蓋條狀漏熱帶附近的高溫區(qū)域,結(jié)合圖5 可得,在這一階段各工況條狀漏熱帶熱流傳輸速率先后開始迅速提高,這是因?yàn)楣摅w內(nèi)部的渦流演化使得渦流中心區(qū)向條狀漏熱帶附近移動(dòng),而渦流的中心區(qū)類似于一個(gè)滯留區(qū),從產(chǎn)生開始就攜帶著低溫流體,渦流中心區(qū)與鄰近區(qū)域的熱量交換進(jìn)行的很慢,當(dāng)其移動(dòng)靠近條狀漏熱帶邊界后,條狀漏熱帶因?yàn)楦浇鼌^(qū)域存在冷流體團(tuán),熱量傳輸速率顯著加快.射流噴嘴越靠近罐體底部,由低溫射流引起的渦流中心區(qū)越容易到達(dá)條狀漏熱帶附近,因此各工況中條狀漏熱帶熱量傳輸速率曲線快速增長(zhǎng)的先后順序是C4,C3,C2,C1.在渦流中心區(qū)靠近條狀漏熱帶一段時(shí)間后,隨著低溫射流時(shí)間的進(jìn)行,渦流中心區(qū)溫度也會(huì)因?yàn)闊峤粨Q而上升,此后通過條狀漏熱帶進(jìn)入罐體的熱量傳輸速率就會(huì)減慢,而從罐體內(nèi)部輸出的熱量在增加,因此罐體內(nèi)部流體平均溫度就會(huì)開始下降(圖11 中C3 和C4 曲線).

    圖10 圓形低溫射流噴嘴儲(chǔ)罐平均速度隨低溫射流時(shí)間演化圖Fig.10 Evolution diagram of the average velocity of the tank with the cryogenic jet time in the circular jet nozzle

    圖11 圓形低溫射流噴嘴儲(chǔ)罐平均溫度隨低溫射流時(shí)間演化圖Fig.11 Evolution diagram of the average temperature of the tank with the cryogenic jet time in the circular jet nozzle

    圖12 各工況低溫射流1 024.08 s 時(shí)溫度云圖分布Fig.12 Isothermals of cases at t=1 024.08 s for the cryogenic jet

    圖13 各工況低溫射流1 024.08 s 時(shí)流場(chǎng)圖Fig.13 Flow field diagram at t=1 024.08 s for the cryogenic jet

    通過以上分析可知4 種工況對(duì)條狀漏熱帶附近高溫區(qū)域的作用機(jī)理是一致的,但由于射流噴嘴在罐內(nèi)位置的不同,發(fā)揮作用的時(shí)間有所不同.整體而言,低溫射流抑制熱分層的機(jī)理主要有兩種,一種是在射流早期,通過促進(jìn)儲(chǔ)罐內(nèi)部流體運(yùn)動(dòng),避免熱量在儲(chǔ)罐內(nèi)部局部區(qū)域大量積聚形成高溫區(qū)域,包括上文中第一和第二階段;另一種則是低溫射流持續(xù)發(fā)展一段時(shí)間,因?yàn)閮?chǔ)罐內(nèi)部流場(chǎng)演化,射流引起的攜帶著低溫流體的渦流中心移動(dòng)到條狀漏熱帶附近與高溫區(qū)域進(jìn)行熱量交換.結(jié)合各工況各時(shí)間段的流場(chǎng)圖,可以發(fā)現(xiàn)在射流初期,渦流中心區(qū)位于射流噴嘴附近并沿著中軸線向罐底運(yùn)動(dòng),此時(shí)低溫射流對(duì)條狀漏熱帶附近高溫區(qū)域影響較小,然后罐體內(nèi)部流場(chǎng)逐步發(fā)展,渦流中心區(qū)沿著罐體壁面爬升向條狀漏熱帶移動(dòng),因此條狀漏熱帶熱量傳輸速率有了顯著提高,在渦流中心區(qū)吸收從漏熱帶進(jìn)入的熱量逐漸升溫之后,漏熱帶的熱量傳輸速率會(huì)逐漸穩(wěn)定.通過對(duì)圖5 中各條曲線對(duì)時(shí)間積分得到各工況漏熱帶傳輸熱量隨時(shí)間演化圖(圖14),可以看出工況C4 在整個(gè)低溫射流時(shí)間內(nèi)轉(zhuǎn)移熱量的能力更強(qiáng).結(jié)合圖15 給出的溫度區(qū)間面積占比累計(jì)圖(溫度區(qū)間面積占比是指將計(jì)算域溫度范圍每隔0.5 K 設(shè)置為一個(gè)溫度區(qū)間,統(tǒng)計(jì)位于各溫度區(qū)間的計(jì)算域面積占總計(jì)算域面積的百分比并繪制成累計(jì)圖)發(fā)現(xiàn)在低溫射流結(jié)束時(shí),各工況介于21.5 ~25 K 之間的溫度占比基本一致,但是19 ~21.5 K 之間的溫度占比C4 明顯低于其他工況,也就是說明工況C4 消除高溫區(qū)域的效果更徹底,作用更明顯,這和通過漏熱帶傳輸熱量作為判據(jù)得出的結(jié)論是一致的.

    圖14 圓形射流噴嘴各工況漏熱帶傳輸熱量隨時(shí)間演化圖Fig.14 Evolution diagram of heat flux over time for each circular nozzle case

    圖15 圓形射流噴嘴各工況溫度區(qū)間面積占比累計(jì)圖Fig.15 Comparison of the cumulative curve of the temperature-area-ratio for each circular nozzle case

    3.2 兩種射流噴嘴形狀對(duì)比

    基于小尺寸儲(chǔ)罐(直徑70 mm)的研究[38]表明:圓形射流噴嘴因?yàn)槌隽鞣较蚋?,罐體內(nèi)部流場(chǎng)演化更劇烈,因此消除熱分層的效果優(yōu)于半球形射流噴嘴.在結(jié)果討論中發(fā)現(xiàn)對(duì)于大尺寸液氫儲(chǔ)罐,圓形射流噴嘴在罐內(nèi)位置對(duì)熱分層的消除機(jī)理是一致的,只是作用時(shí)間存在差異.相對(duì)而言,靠近儲(chǔ)罐出口的位置能包含其他工況罐體內(nèi)部流體流場(chǎng)演變過程,因此以兩種射流噴嘴在靠近出口的位置時(shí)消除熱分層的效果為例進(jìn)行對(duì)比分析,即工況C1 和H1.半球形射流噴嘴截面半徑同為5.5 cm,為了保持圓形射流噴嘴相同入射質(zhì)量流量,將射流速度設(shè)為0.05 m/s,其他設(shè)定與圓形射流噴嘴情況一致.計(jì)算并繪制兩種射流噴嘴漏熱帶熱量傳輸速率隨低溫射流時(shí)間演化曲線(圖16)通過對(duì)比可以發(fā)現(xiàn),圓形射流噴嘴消除熱分層的效果明顯高于半球形射流噴嘴.圖17 和圖18 分別給出半球形射流噴嘴在低溫射流進(jìn)行424.08 s 時(shí)的溫度云圖以及流場(chǎng)圖,從溫度云圖可以看出低溫流體都聚集在射流噴嘴附近,對(duì)條狀漏熱帶附近高溫區(qū)域影響很小.從流場(chǎng)圖可以發(fā)現(xiàn)半球形射流噴嘴形成的渦流發(fā)展緩慢,無法將低溫流體輸運(yùn)到條狀漏熱帶附近的高溫區(qū)域,罐體內(nèi)部熱分層抑制主要依靠流體流動(dòng)使得熱量無法聚積,因此消除熱分層的效果不如圓形射流噴嘴.

    圖16 圓形與半球形射流噴嘴熱分層消除效果對(duì)比圖Fig.16 Contrast diagram of temperature stratification of circular and hemispherical jet nozzles

    圖17 工況H1 在低溫射流424.08 s 的溫度云圖Fig.17 Isothermals of H1 at t=424.08 s for the cryogenic jet

    圖18 工況H1 在低溫射流424.08 s 的流場(chǎng)圖Fig.18 Flow field diagram of H1 at t=424.08 s for the cryogenic jet

    4 結(jié)論

    本文利用軸對(duì)稱的具有低溫射流裝置的零蒸發(fā)儲(chǔ)罐模型,通過數(shù)值模擬研究了微重力條件下液氫儲(chǔ)存過程中利用低溫射流消除環(huán)境漏熱引起的熱分層現(xiàn)象的效果.通過分析儲(chǔ)罐內(nèi)部流場(chǎng)流動(dòng)和溫度分布的時(shí)間演化過程,研究了射流噴嘴在罐內(nèi)不同位置對(duì)利用低溫射流消除儲(chǔ)罐內(nèi)部熱分層效果的影響.研究結(jié)果表明:對(duì)于大尺寸儲(chǔ)罐,當(dāng)采用圓形射流噴嘴且低溫射流條件相同時(shí),射流噴嘴的位置對(duì)罐體內(nèi)部熱分層消除效果影響不是很明顯,在本文的入射條件下,當(dāng)?shù)蜏厣淞髦脫Q率達(dá)到2%,即低溫射流時(shí)間持續(xù)700 s 時(shí),罐體內(nèi)部熱分層的消除效果最顯著;同時(shí),當(dāng)射流噴嘴位于儲(chǔ)罐內(nèi)部同一相對(duì)位置且入射流量相同時(shí),圓形射流噴嘴因出流方向更集中,罐內(nèi)流場(chǎng)演變更快,熱分層消除效果比半球形射流噴嘴更好.

    猜你喜歡
    罐體熱帶射流
    深海逃逸艙射流注水均壓過程仿真分析
    低壓天然氣泄漏射流擴(kuò)散特性研究
    煤氣與熱力(2022年4期)2022-05-23 12:45:00
    一種醫(yī)用塑料桶注塑成型裝置
    熱帶風(fēng)情
    女報(bào)(2020年7期)2020-08-17 07:16:05
    基于Dynaform有限元模擬的3104鋁質(zhì)罐體再拉伸工藝優(yōu)化
    模具制造(2019年7期)2019-09-25 07:29:58
    熱帶的鳥兒
    圓滾滾的熱帶“龍”
    熱帶小鳥
    射流齒形噴嘴射流流場(chǎng)與氣動(dòng)聲學(xué)分析
    基于ANSYS的LNG儲(chǔ)罐罐體溫度場(chǎng)的數(shù)值計(jì)算
    最新的欧美精品一区二区| 狠狠狠狠99中文字幕| 国产一级毛片在线| 岛国在线观看网站| 久久国产亚洲av麻豆专区| 亚洲av成人一区二区三| 亚洲精品国产色婷婷电影| 国产成人免费观看mmmm| 国产精品偷伦视频观看了| 国产有黄有色有爽视频| 十八禁网站网址无遮挡| 欧美老熟妇乱子伦牲交| 亚洲精品久久成人aⅴ小说| 天天影视国产精品| 人人妻,人人澡人人爽秒播| 乱人伦中国视频| 免费黄频网站在线观看国产| 男人添女人高潮全过程视频| 久久人妻福利社区极品人妻图片| 亚洲欧美成人综合另类久久久| 国产成人啪精品午夜网站| 亚洲av片天天在线观看| 老熟妇仑乱视频hdxx| h视频一区二区三区| 欧美激情极品国产一区二区三区| 如日韩欧美国产精品一区二区三区| 亚洲第一欧美日韩一区二区三区 | av又黄又爽大尺度在线免费看| 国产99久久九九免费精品| 女性被躁到高潮视频| 久久久久久久大尺度免费视频| 国产av一区二区精品久久| 中文字幕另类日韩欧美亚洲嫩草| 如日韩欧美国产精品一区二区三区| 日本猛色少妇xxxxx猛交久久| 欧美日韩国产mv在线观看视频| 性色av一级| 日本欧美视频一区| 90打野战视频偷拍视频| 视频在线观看一区二区三区| 大码成人一级视频| 50天的宝宝边吃奶边哭怎么回事| 亚洲午夜精品一区,二区,三区| 久久精品国产a三级三级三级| 中文字幕人妻丝袜一区二区| 69av精品久久久久久 | 久久久久久久久免费视频了| 国产老妇伦熟女老妇高清| 亚洲精品国产av蜜桃| 欧美老熟妇乱子伦牲交| 大型av网站在线播放| 黄色片一级片一级黄色片| 久久久久久久久免费视频了| 黄色视频在线播放观看不卡| 十八禁高潮呻吟视频| 国产欧美亚洲国产| 亚洲人成电影观看| 亚洲国产中文字幕在线视频| 在线观看舔阴道视频| 国产欧美日韩一区二区三区在线| 另类亚洲欧美激情| 久久精品国产亚洲av高清一级| 日本wwww免费看| 9色porny在线观看| 欧美精品亚洲一区二区| 欧美精品高潮呻吟av久久| www.av在线官网国产| 久久九九热精品免费| 欧美另类亚洲清纯唯美| 天天躁日日躁夜夜躁夜夜| 老熟妇仑乱视频hdxx| 日本精品一区二区三区蜜桃| 亚洲激情五月婷婷啪啪| 欧美一级毛片孕妇| 大片免费播放器 马上看| 午夜精品国产一区二区电影| 久久天堂一区二区三区四区| 亚洲av欧美aⅴ国产| 成年人黄色毛片网站| 成人国产av品久久久| 中文字幕精品免费在线观看视频| 亚洲 国产 在线| 高清视频免费观看一区二区| 午夜福利在线观看吧| 老司机福利观看| 热99国产精品久久久久久7| 97在线人人人人妻| 精品久久久久久久毛片微露脸 | 午夜激情av网站| 亚洲欧美色中文字幕在线| 一个人免费在线观看的高清视频 | 下体分泌物呈黄色| 美女午夜性视频免费| 精品免费久久久久久久清纯 | 黄片大片在线免费观看| 法律面前人人平等表现在哪些方面 | 男女午夜视频在线观看| 午夜福利,免费看| 大香蕉久久网| 亚洲欧美一区二区三区黑人| 免费观看av网站的网址| 久久久国产精品麻豆| videosex国产| 91av网站免费观看| 亚洲精品中文字幕一二三四区 | 嫁个100分男人电影在线观看| 纯流量卡能插随身wifi吗| 麻豆乱淫一区二区| 精品国产超薄肉色丝袜足j| 精品亚洲成a人片在线观看| 午夜影院在线不卡| 老汉色av国产亚洲站长工具| 一区二区三区精品91| 中文字幕最新亚洲高清| 久久久久国产精品人妻一区二区| 久久亚洲精品不卡| 欧美精品人与动牲交sv欧美| 亚洲国产中文字幕在线视频| 亚洲午夜精品一区,二区,三区| 男女无遮挡免费网站观看| 成人免费观看视频高清| 国产欧美日韩综合在线一区二区| 亚洲 欧美一区二区三区| 久久久久久人人人人人| 欧美黑人欧美精品刺激| 国产精品偷伦视频观看了| 久久人人爽av亚洲精品天堂| 国产亚洲一区二区精品| 一本一本久久a久久精品综合妖精| 天天添夜夜摸| 国产精品成人在线| 黄色怎么调成土黄色| 另类精品久久| 免费黄频网站在线观看国产| 日本精品一区二区三区蜜桃| 欧美在线一区亚洲| 最新的欧美精品一区二区| 成人影院久久| 久久久久久免费高清国产稀缺| 久久中文看片网| 中文欧美无线码| 2018国产大陆天天弄谢| 美女扒开内裤让男人捅视频| 久久女婷五月综合色啪小说| 黄频高清免费视频| 两性夫妻黄色片| 国产视频一区二区在线看| 乱人伦中国视频| 亚洲第一青青草原| 欧美xxⅹ黑人| 99精品久久久久人妻精品| 多毛熟女@视频| 黑人操中国人逼视频| 黄网站色视频无遮挡免费观看| 精品亚洲成a人片在线观看| 美女高潮喷水抽搐中文字幕| 在线看a的网站| 久久精品久久久久久噜噜老黄| 欧美国产精品一级二级三级| 国产在视频线精品| 老鸭窝网址在线观看| 十八禁人妻一区二区| 免费观看a级毛片全部| 亚洲精品国产一区二区精华液| 精品国内亚洲2022精品成人 | 999精品在线视频| 亚洲伊人色综图| 黄色视频不卡| 欧美日韩黄片免| 欧美一级毛片孕妇| 黄色毛片三级朝国网站| 久久久久国产精品人妻一区二区| 在线观看一区二区三区激情| 国产在线免费精品| 亚洲精品久久久久久婷婷小说| 老司机影院成人| 丰满人妻熟妇乱又伦精品不卡| 亚洲国产精品一区三区| 不卡av一区二区三区| 欧美日韩亚洲高清精品| 午夜免费观看性视频| 90打野战视频偷拍视频| 1024视频免费在线观看| 中文欧美无线码| 视频区欧美日本亚洲| 亚洲国产欧美网| 久久久久久久精品精品| 蜜桃国产av成人99| 亚洲欧美色中文字幕在线| 黄片小视频在线播放| 久久人妻熟女aⅴ| 免费在线观看黄色视频的| 国产精品免费大片| www.精华液| 两个人看的免费小视频| 国产精品熟女久久久久浪| 午夜两性在线视频| av欧美777| 欧美黄色片欧美黄色片| 最近最新中文字幕大全免费视频| 桃花免费在线播放| 99九九在线精品视频| 一边摸一边做爽爽视频免费| 久久精品aⅴ一区二区三区四区| 成在线人永久免费视频| 最新在线观看一区二区三区| 欧美日韩av久久| 男女床上黄色一级片免费看| 一区在线观看完整版| 91老司机精品| 秋霞在线观看毛片| 国产精品.久久久| 色综合欧美亚洲国产小说| 一级黄色大片毛片| 女性生殖器流出的白浆| 女人精品久久久久毛片| 丝袜脚勾引网站| 高清视频免费观看一区二区| 两人在一起打扑克的视频| 久久 成人 亚洲| 女人高潮潮喷娇喘18禁视频| 精品人妻熟女毛片av久久网站| 欧美人与性动交α欧美精品济南到| 又黄又粗又硬又大视频| 精品一品国产午夜福利视频| 午夜精品久久久久久毛片777| 91精品伊人久久大香线蕉| avwww免费| 一区二区三区精品91| 亚洲成人手机| 无遮挡黄片免费观看| 电影成人av| www.熟女人妻精品国产| 悠悠久久av| 一本一本久久a久久精品综合妖精| 极品人妻少妇av视频| 在线观看免费高清a一片| 国产真人三级小视频在线观看| 俄罗斯特黄特色一大片| 免费av中文字幕在线| 亚洲国产成人一精品久久久| 精品国产一区二区三区久久久樱花| 免费不卡黄色视频| 免费日韩欧美在线观看| 日本精品一区二区三区蜜桃| 十分钟在线观看高清视频www| 亚洲精品第二区| videosex国产| 夫妻午夜视频| 成人av一区二区三区在线看 | 午夜福利在线观看吧| 欧美国产精品一级二级三级| 亚洲自偷自拍图片 自拍| 12—13女人毛片做爰片一| 老汉色av国产亚洲站长工具| 99热网站在线观看| 欧美精品亚洲一区二区| 国产精品二区激情视频| 亚洲综合色网址| 婷婷色av中文字幕| 国产伦人伦偷精品视频| 国产精品一区二区在线不卡| videos熟女内射| 国产免费现黄频在线看| 国产亚洲午夜精品一区二区久久| 在线观看www视频免费| 欧美激情极品国产一区二区三区| 亚洲人成电影免费在线| 国产精品1区2区在线观看. | 日本wwww免费看| 97人妻天天添夜夜摸| 久久精品成人免费网站| 久久久久国内视频| 菩萨蛮人人尽说江南好唐韦庄| 国产亚洲午夜精品一区二区久久| 法律面前人人平等表现在哪些方面 | 久久久久久久大尺度免费视频| 下体分泌物呈黄色| 成年av动漫网址| 黄色视频不卡| 制服诱惑二区| 欧美少妇被猛烈插入视频| 在线 av 中文字幕| 国产淫语在线视频| 国产精品一区二区在线观看99| 欧美成狂野欧美在线观看| 欧美日韩亚洲国产一区二区在线观看 | 久久人妻福利社区极品人妻图片| 国产成人欧美在线观看 | 国产男人的电影天堂91| 女人被躁到高潮嗷嗷叫费观| 亚洲av日韩在线播放| 欧美久久黑人一区二区| 97在线人人人人妻| 制服诱惑二区| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲精品一区蜜桃| 中亚洲国语对白在线视频| 日日摸夜夜添夜夜添小说| 少妇被粗大的猛进出69影院| 一进一出抽搐动态| av欧美777| 国产成人av教育| 一级a爱视频在线免费观看| 久久亚洲精品不卡| 麻豆国产av国片精品| 欧美精品av麻豆av| 亚洲少妇的诱惑av| av超薄肉色丝袜交足视频| 色综合欧美亚洲国产小说| 久久久久久人人人人人| 交换朋友夫妻互换小说| 大片免费播放器 马上看| 深夜精品福利| 免费久久久久久久精品成人欧美视频| 国产成人欧美| 天天操日日干夜夜撸| 美女福利国产在线| 国产伦理片在线播放av一区| 最新在线观看一区二区三区| 久久久久久免费高清国产稀缺| 欧美xxⅹ黑人| 亚洲色图综合在线观看| 国产精品国产av在线观看| 亚洲国产欧美在线一区| 国产精品影院久久| 国产av又大| 国产精品免费视频内射| 男女免费视频国产| 精品少妇内射三级| 中文字幕人妻丝袜一区二区| 国产精品 欧美亚洲| 午夜福利视频在线观看免费| 国产极品粉嫩免费观看在线| 成年人午夜在线观看视频| 美女脱内裤让男人舔精品视频| 久久久久视频综合| 韩国精品一区二区三区| 国产av又大| 免费日韩欧美在线观看| 欧美激情 高清一区二区三区| 熟女少妇亚洲综合色aaa.| 黑人巨大精品欧美一区二区蜜桃| 后天国语完整版免费观看| 亚洲精品一卡2卡三卡4卡5卡 | 黄色片一级片一级黄色片| 亚洲精品一卡2卡三卡4卡5卡 | 亚洲欧美成人综合另类久久久| 啦啦啦啦在线视频资源| 久久久久国内视频| 欧美精品亚洲一区二区| 成年美女黄网站色视频大全免费| 中文字幕高清在线视频| 久久av网站| 欧美国产精品一级二级三级| 黄色 视频免费看| 日韩大码丰满熟妇| 精品福利观看| 纵有疾风起免费观看全集完整版| 亚洲欧美清纯卡通| 岛国毛片在线播放| 黄片播放在线免费| 精品乱码久久久久久99久播| 伦理电影免费视频| 免费少妇av软件| 首页视频小说图片口味搜索| 亚洲专区中文字幕在线| 后天国语完整版免费观看| 一本大道久久a久久精品| 国产在线视频一区二区| 精品一区二区三区av网在线观看 | 女人高潮潮喷娇喘18禁视频| 久久这里只有精品19| 亚洲欧美精品综合一区二区三区| 深夜精品福利| 中文字幕最新亚洲高清| 久久狼人影院| 欧美日韩亚洲国产一区二区在线观看 | 亚洲国产欧美一区二区综合| 免费在线观看视频国产中文字幕亚洲 | 99久久精品国产亚洲精品| 韩国精品一区二区三区| a级毛片在线看网站| www日本在线高清视频| 精品国产一区二区三区久久久樱花| 中文字幕高清在线视频| 美女高潮到喷水免费观看| 青青草视频在线视频观看| 蜜桃在线观看..| 人人妻人人澡人人看| 亚洲全国av大片| 成人手机av| 国产精品国产三级国产专区5o| 色精品久久人妻99蜜桃| 日日爽夜夜爽网站| 中文字幕人妻丝袜一区二区| 国产日韩一区二区三区精品不卡| 久久精品熟女亚洲av麻豆精品| 大码成人一级视频| 啦啦啦免费观看视频1| 日本vs欧美在线观看视频| 亚洲国产欧美日韩在线播放| 久久香蕉激情| 久久久久久人人人人人| 91国产中文字幕| 水蜜桃什么品种好| 久久天躁狠狠躁夜夜2o2o| a 毛片基地| 欧美xxⅹ黑人| 亚洲va日本ⅴa欧美va伊人久久 | 99精国产麻豆久久婷婷| 亚洲国产精品一区三区| 悠悠久久av| 80岁老熟妇乱子伦牲交| 国产又色又爽无遮挡免| 亚洲一卡2卡3卡4卡5卡精品中文| 欧美激情久久久久久爽电影 | 动漫黄色视频在线观看| 91成年电影在线观看| 搡老乐熟女国产| 欧美日韩亚洲高清精品| 一区福利在线观看| 黄色片一级片一级黄色片| 亚洲精品成人av观看孕妇| 大陆偷拍与自拍| 日韩大片免费观看网站| 黄色视频在线播放观看不卡| 两个人免费观看高清视频| 国产黄色免费在线视频| 搡老熟女国产l中国老女人| 99精品欧美一区二区三区四区| 久久久国产一区二区| 精品亚洲成a人片在线观看| 欧美精品一区二区免费开放| www日本在线高清视频| 亚洲一码二码三码区别大吗| 丝袜脚勾引网站| 中国美女看黄片| 国产精品一二三区在线看| 91成人精品电影| 免费在线观看黄色视频的| 999久久久国产精品视频| 免费人妻精品一区二区三区视频| 亚洲人成77777在线视频| 国产免费视频播放在线视频| 国产亚洲欧美在线一区二区| 国产亚洲精品一区二区www | 又紧又爽又黄一区二区| 久久精品人人爽人人爽视色| 国产成人免费观看mmmm| 日韩中文字幕欧美一区二区| 视频区欧美日本亚洲| 黄色a级毛片大全视频| 99国产精品免费福利视频| 日日摸夜夜添夜夜添小说| 亚洲人成77777在线视频| 亚洲精品国产色婷婷电影| 人人妻人人爽人人添夜夜欢视频| 50天的宝宝边吃奶边哭怎么回事| 午夜福利乱码中文字幕| 制服人妻中文乱码| av有码第一页| 无限看片的www在线观看| 天天操日日干夜夜撸| 在线十欧美十亚洲十日本专区| 电影成人av| av免费在线观看网站| 青草久久国产| 男人舔女人的私密视频| 一区二区av电影网| 国产精品麻豆人妻色哟哟久久| 国产成人啪精品午夜网站| 免费在线观看完整版高清| 一级毛片精品| 99精国产麻豆久久婷婷| 777久久人妻少妇嫩草av网站| 国产欧美日韩一区二区三 | 搡老岳熟女国产| 欧美另类一区| 啦啦啦 在线观看视频| av国产精品久久久久影院| 视频区欧美日本亚洲| 免费在线观看视频国产中文字幕亚洲 | 精品久久蜜臀av无| 亚洲精品久久成人aⅴ小说| 亚洲欧洲精品一区二区精品久久久| 欧美激情 高清一区二区三区| 淫妇啪啪啪对白视频 | 免费看十八禁软件| 婷婷成人精品国产| 国产一区二区在线观看av| 免费少妇av软件| 国产人伦9x9x在线观看| 亚洲第一青青草原| 国产高清国产精品国产三级| 久久国产精品男人的天堂亚洲| 俄罗斯特黄特色一大片| 欧美日韩成人在线一区二区| 99国产精品一区二区三区| 免费观看人在逋| 天天躁日日躁夜夜躁夜夜| av有码第一页| 亚洲成人手机| av欧美777| 亚洲avbb在线观看| 91字幕亚洲| 亚洲天堂av无毛| 成人av一区二区三区在线看 | 久久久久视频综合| 国产精品麻豆人妻色哟哟久久| 桃花免费在线播放| 国产在线一区二区三区精| 欧美日韩视频精品一区| 国产亚洲av高清不卡| 色精品久久人妻99蜜桃| 一本久久精品| 香蕉丝袜av| av在线播放精品| 国产又色又爽无遮挡免| 亚洲欧洲精品一区二区精品久久久| 久久久精品免费免费高清| 欧美日韩亚洲国产一区二区在线观看 | 下体分泌物呈黄色| 久久人妻熟女aⅴ| 99国产精品99久久久久| 成年人免费黄色播放视频| 91老司机精品| 亚洲第一青青草原| 91麻豆av在线| 国产成人系列免费观看| 欧美人与性动交α欧美精品济南到| 韩国高清视频一区二区三区| 国产免费现黄频在线看| 一级毛片女人18水好多| 在线永久观看黄色视频| 狂野欧美激情性xxxx| 人人妻,人人澡人人爽秒播| 大陆偷拍与自拍| 18禁裸乳无遮挡动漫免费视频| 午夜精品国产一区二区电影| 最黄视频免费看| 人人妻,人人澡人人爽秒播| 多毛熟女@视频| 亚洲中文av在线| 国产老妇伦熟女老妇高清| 欧美精品一区二区免费开放| 热99国产精品久久久久久7| 又大又爽又粗| 狠狠精品人妻久久久久久综合| 老司机午夜福利在线观看视频 | 亚洲国产看品久久| 丰满人妻熟妇乱又伦精品不卡| 国产亚洲欧美在线一区二区| 美女国产高潮福利片在线看| 桃红色精品国产亚洲av| 日本91视频免费播放| √禁漫天堂资源中文www| 亚洲av美国av| 欧美人与性动交α欧美软件| 成在线人永久免费视频| 午夜两性在线视频| 最新在线观看一区二区三区| 乱人伦中国视频| 少妇的丰满在线观看| 精品视频人人做人人爽| 国产精品 欧美亚洲| 韩国精品一区二区三区| 人妻 亚洲 视频| 亚洲av男天堂| 久久精品国产a三级三级三级| 国产精品 国内视频| 欧美日韩精品网址| 最近中文字幕2019免费版| 夜夜骑夜夜射夜夜干| 黄色视频,在线免费观看| 少妇猛男粗大的猛烈进出视频| 欧美日韩中文字幕国产精品一区二区三区 | 亚洲成国产人片在线观看| 巨乳人妻的诱惑在线观看| 国产男女内射视频| 日韩 欧美 亚洲 中文字幕| 日本vs欧美在线观看视频| 精品亚洲成a人片在线观看| 成人国语在线视频| 日本欧美视频一区| 成年人黄色毛片网站| 午夜视频精品福利| 国产成人av教育| 新久久久久国产一级毛片| 国产av又大| 亚洲全国av大片| www.熟女人妻精品国产| 俄罗斯特黄特色一大片| 亚洲第一欧美日韩一区二区三区 | 十八禁网站网址无遮挡| av天堂久久9| 国产av精品麻豆| 另类精品久久| 91麻豆av在线| 午夜福利,免费看| 国产成+人综合+亚洲专区| tocl精华| 欧美av亚洲av综合av国产av| 50天的宝宝边吃奶边哭怎么回事| 男女下面插进去视频免费观看| 另类精品久久| 熟女少妇亚洲综合色aaa.| 90打野战视频偷拍视频| 精品第一国产精品| 男女免费视频国产| 好男人电影高清在线观看| 高清av免费在线| 国产在线免费精品| 侵犯人妻中文字幕一二三四区| 美女主播在线视频|