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

    超聲波電噴推力器羽流中和特性研究

    2018-03-26 22:06:56于博張巖賀偉國杭觀榮康小錄趙青
    物理學(xué)報(bào) 2018年4期
    關(guān)鍵詞:電噴羽流推力器

    于博 張巖 賀偉國 杭觀榮 康小錄 趙青

    1)(電子科技大學(xué)物理電子學(xué)院,成都 611731)

    2)(上海空間推進(jìn)研究所電推進(jìn)事業(yè)部,上海 201112)

    3)(上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院,上海 200240)

    1 引 言

    近些年,隨著器件小型化、集成化技術(shù)的快速發(fā)展[1?3],小衛(wèi)星平臺(tái)的研制已成為航天領(lǐng)域的重點(diǎn)研究方向之一,各類新型微推進(jìn)技術(shù)也隨之受到關(guān)注.其中,電噴推力器是利用帶電流體在高強(qiáng)電場(chǎng)中的極化作用,形成類似泰勒錐幾何結(jié)構(gòu)[4]的發(fā)射液錐,在液滴的不斷供給和發(fā)射的動(dòng)態(tài)平衡下,完成帶電液滴在真空中的發(fā)射,接著,在中和機(jī)理作用后,形成定向運(yùn)動(dòng)的中性液滴,提供推力.電噴推力器的功率一般在1—20 W,推力可以覆蓋1—100μN(yùn)不等[5,6],是應(yīng)用于立方星平臺(tái)的微推進(jìn)技術(shù)[7].現(xiàn)今報(bào)道的絕大多數(shù)都是毛細(xì)孔電噴推力器[8?11]:包括傳統(tǒng)的以離子液滴為工質(zhì)、金屬毛細(xì)孔電噴推力器[8],以液態(tài)碘為工質(zhì)、硅發(fā)射體的毛細(xì)孔電噴推力器[9],外加磁場(chǎng)的鐵磁體液滴電噴推力器[10],以及高液壓通道的膠體電噴[11].以上膠體電噴的發(fā)射機(jī)理都是利用多孔隙中液體的毛細(xì)作用來提供形成泰勒錐的觸發(fā)條件.由此,膠體電噴的推力密度和液滴形成太小會(huì)受限于多孔隙的加工情況,孔隙數(shù)量會(huì)受到加工空間和加工精度的限制,如果這些多孔數(shù)量較多,則加工時(shí)會(huì)影響到旁邊孔隙,或是如果孔隙加工不均勻,發(fā)射液滴也會(huì)明顯存在空間上的分布不均特性[12].因此,在2010年,有學(xué)者提出了一種推力密度更高的電噴推力器——超聲波電噴推力器[12].

    超聲波電噴推力器是基于超聲源對(duì)發(fā)射表面液體提供震動(dòng),形成液體表面規(guī)則的、陣列式的駐波[13,14].這些細(xì)微駐波的形成為帶電液體形成發(fā)射泰勒錐提供了周期性的條件,即液滴的發(fā)射僅在波峰形成的一段極短時(shí)間內(nèi)完成[12,15].近幾年來,本團(tuán)隊(duì)人員對(duì)超聲波電噴的研究都集中在駐波形成機(jī)理、帶電液滴發(fā)射的驗(yàn)證等研究[12?15].然而,在2017年3月,以場(chǎng)發(fā)射陰極作為中和器進(jìn)行性能試驗(yàn)時(shí),發(fā)現(xiàn)了以下兩個(gè)現(xiàn)象:1)推力與理論發(fā)射方向存在明顯的3°—4°的夾角;2)推力器的陽極效率沒有理論計(jì)算值(70%以上)[12,13]那么高,試驗(yàn)測(cè)得只有47%左右.面對(duì)上述問題,需要對(duì)超聲電噴的液滴羽流中和過程進(jìn)行研究,獲得羽流的中和特性機(jī)理,以解答以上兩個(gè)實(shí)際問題,以及對(duì)相應(yīng)的優(yōu)化策略給出參考.

    針對(duì)超聲波電噴推力器的羽流中和過程,采用試驗(yàn)方法的難度較高:1)探針介入測(cè)量會(huì)破壞羽流場(chǎng)的本身特性,影響試驗(yàn)結(jié)果的準(zhǔn)確性;2)羽流的具體中和過程無法通過試驗(yàn)來進(jìn)行細(xì)節(jié)上的捕捉,難以完成對(duì)中和機(jī)理的挖掘.因此,本文考慮使用數(shù)值計(jì)算方法來進(jìn)行研究.首先,由于羽流中和過程是在電子與液滴間完成,這會(huì)涉及電子中和、液滴的流動(dòng)與傳熱、液滴的破碎與重組等過程,故需要建立一種捕捉該過程的新算法——帶電液滴中和模型(neutralization of electrons and charged droplets,NECD模型),具體將在第2節(jié)進(jìn)行詳細(xì)闡述;其次,在第3節(jié),對(duì)該算法進(jìn)行試驗(yàn)驗(yàn)證,介紹算法的計(jì)算精度和不足;第4節(jié)則對(duì)超聲波電噴推力器羽流中和過程進(jìn)行計(jì)算結(jié)果展示和分析,揭示出推力器的羽流中和特性;在第5節(jié),對(duì)上一段中的兩個(gè)問題進(jìn)行解答,以及對(duì)超聲電噴的性能設(shè)計(jì)給出相應(yīng)的優(yōu)化策略.

    2 NECD模型

    電噴推力器羽流的中和模型應(yīng)具備捕捉正負(fù)帶電粒子在中和過程中粒子的輸運(yùn)過程、碰撞過程以及液滴的破碎與重組等過程的能力,并且可以完成中和過程中各粒子的能量變化的追蹤.因此,基于上述要求,輸運(yùn)過程的計(jì)算采用PIC(單元粒子)算法是比較合適的,以捕捉帶電粒子在電場(chǎng)作用下的運(yùn)動(dòng)狀態(tài)以及帶電粒子的空間電勢(shì)分布.PIC算法在大多數(shù)等離子體數(shù)值計(jì)算的文獻(xiàn)中都有介紹[16,17],這里不再贅述.

    除輸運(yùn)過程外,本文模型還應(yīng)包括粒子間碰撞、液滴破碎、傳熱等過程的子模型.關(guān)于這些子模型,本節(jié)進(jìn)行重點(diǎn)介紹.

    2.1 電子-液滴碰撞過程

    在討論電子與液滴的碰撞時(shí),分為兩種情況:1)如果液滴帶正電,電子與液滴的中和過程;2)如果液滴呈中性或負(fù)電,電子被液滴吸附過程.

    需要說明的是,對(duì)電子-液滴碰撞的判斷是基于對(duì)運(yùn)動(dòng)到網(wǎng)格內(nèi)的電子,在當(dāng)前網(wǎng)格內(nèi)各物理參數(shù)制約下的碰撞概率判斷.因而,在判斷一個(gè)電子在液滴群內(nèi)穿梭時(shí),假設(shè)液滴可以作為背景氣體存在,電子能否與液滴發(fā)生碰撞就是指在一個(gè)時(shí)間步長內(nèi),電子移動(dòng)路徑長度上是否包含一個(gè)電子平均自由程.

    首先,一個(gè)時(shí)間步長內(nèi)電子掠過的長度為

    ve為電子運(yùn)動(dòng)速度(m/s),Δt為時(shí)間步長(s).

    其次,電子在當(dāng)前網(wǎng)格內(nèi)運(yùn)動(dòng)的平均自由程由該網(wǎng)格的相關(guān)參數(shù)確定,有

    其中,Rcell為該網(wǎng)格內(nèi)液滴的平均半徑(m),Ncell為該網(wǎng)格內(nèi)液滴的數(shù)密度(m?3).

    接著,對(duì)電子與液滴的碰撞進(jìn)行判斷:

    1)如果電子掠過的路徑大于等于當(dāng)前網(wǎng)格內(nèi)的平均自由程,則認(rèn)為電子在該網(wǎng)格內(nèi)觸發(fā)碰撞,并將該時(shí)間步長內(nèi)在此網(wǎng)格內(nèi)發(fā)生碰撞的電子電荷對(duì)網(wǎng)格內(nèi)所有液滴進(jìn)行平均分配;

    2)如果電子掠過的路徑小于當(dāng)前網(wǎng)格內(nèi)的平均自由程,則認(rèn)為電子有一定概率在該網(wǎng)格內(nèi)發(fā)生碰撞,概率為

    最后,對(duì)電子的能量沉積進(jìn)行判斷:

    1)如果電子與正電液滴碰撞,則認(rèn)為電子與液滴中和,并將全部電離能Eionization沉積給液滴轉(zhuǎn)化為液滴的內(nèi)能,提升液滴溫度ΔTl(K)(ml為液滴質(zhì)量(kg),cp為液滴比熱容(J/kg·K)),

    2)如果電子與中性或負(fù)電液滴碰撞,則認(rèn)為電子被液滴吸附,不做能量變化處理,待到某時(shí)刻正負(fù)液滴發(fā)生中和時(shí),再進(jìn)行電離能的沉積處理,以此保證能量守恒.

    2.2 液滴之間的碰撞過程

    關(guān)于液滴之間的碰撞過程,采用基于蒙特卡羅的直接模擬(DSMC)算法[18]衍生出的算法.這一點(diǎn)與電子-液滴碰撞不同,是對(duì)液滴在網(wǎng)格物理參數(shù)制約下的判斷.這里需要進(jìn)行假設(shè):忽略三個(gè)或以上的液滴碰撞事件,僅關(guān)注兩液滴的碰撞事件.

    首先,任意網(wǎng)格內(nèi),每個(gè)液滴在當(dāng)前網(wǎng)格內(nèi)發(fā)生液滴間碰撞的概率PNTC為

    其中,ncell為當(dāng)前網(wǎng)格內(nèi)液滴的數(shù)量,σT為該液滴在當(dāng)前網(wǎng)格平均液滴尺寸下的碰撞截面(m2),vl為當(dāng)前液滴的運(yùn)動(dòng)速度(m/s).實(shí)際上,該公式是Bird[18]在計(jì)算稀薄等離子體流體時(shí)提出的重粒子間碰撞模擬的算法,通過粒子在一個(gè)時(shí)間步長內(nèi)掃過體積的大小來判斷粒子間的碰撞概率,這是一種基于對(duì)幾何概型的統(tǒng)計(jì)算法,本文將這種思想應(yīng)用在液滴碰撞的計(jì)算上.

    與DSMC算法不同的是,(5)式是基于對(duì)單個(gè)液滴的碰撞概率為前提,而DSMC算法是對(duì)網(wǎng)格內(nèi)每個(gè)碰撞對(duì)逐一進(jìn)行判斷,但是,(5)式從本質(zhì)上與DSMC算法對(duì)碰撞概率的捕捉是一致的.

    接著,碰撞截面σT為

    Rl為該液滴的半徑(m).

    兩液滴碰撞后,可能會(huì)出現(xiàn)融合或反彈.Jayaretne和Mason認(rèn)為[19],液滴的反彈主要發(fā)生在小尺度液滴對(duì)大尺度液面的碰撞中(例如液滴對(duì)水面或小液滴對(duì)大液滴),并且要滿足一定的條件:R>Eboun,其中R為(0,1)之間的隨機(jī)數(shù),Eboun為反彈效率.Eboun除了要考慮兩個(gè)液滴的物理狀態(tài)外,還要考慮兩液滴的相對(duì)速度以及夾角.在推力器的羽流中,液滴的遷移方向幾乎一致,速度矢量夾角較小,因此,Eboun幾乎為0,可認(rèn)為反彈發(fā)生概率極低,考慮到計(jì)算速度,本文予以忽略.

    兩液滴融合為一體,其運(yùn)動(dòng)狀態(tài)由動(dòng)量守恒來確定:

    這里,對(duì)液滴中和過程的能量沉積進(jìn)行判斷:

    1)如果碰撞發(fā)生在正負(fù)液滴之間,則將電荷絕對(duì)值較小液滴的電荷電離能(|qmin|Eionization)沉積給液滴轉(zhuǎn)化為液滴的內(nèi)能,提升液滴溫度,

    2)如果碰撞發(fā)生在同電荷液滴之間,則認(rèn)為兩液滴的電荷都保留給重組后的液滴,不做能量變化處理,待到某時(shí)刻再有正負(fù)液滴發(fā)生中和時(shí),再進(jìn)行電離能的沉積處理,以此保證能量守恒.

    2.3 液滴的破碎過程

    2.3.1 破碎觸發(fā)的判斷

    液滴在氣體中運(yùn)動(dòng)時(shí),由于外力與內(nèi)部張力的作用會(huì)在液滴表面形成一種對(duì)抗,當(dāng)這種對(duì)抗沒有打破維持液滴的平衡時(shí),液滴可以保持原狀態(tài),當(dāng)這種對(duì)抗已經(jīng)足夠打破這種平衡時(shí),液滴就會(huì)發(fā)生破碎,這是一般情況下的液滴破碎.

    然而,當(dāng)液滴帶有電荷,在靜電場(chǎng)和壓力場(chǎng)中運(yùn)動(dòng)時(shí),則會(huì)受到電場(chǎng)力的作用,但是,在初期液滴破碎的判斷中,液滴未形成尖端狀的極化結(jié)構(gòu),故電場(chǎng)對(duì)液滴的破碎作用暫且不考慮,但是,在液滴破碎成小液滴時(shí),電場(chǎng)力的作用就需要進(jìn)行考慮.由此,運(yùn)動(dòng)液滴與背景氣體發(fā)生相對(duì)運(yùn)動(dòng)時(shí)產(chǎn)生的沖擊E1,液滴的表面張力維持作用E2.此刻,若后者大于前者,液滴可以在此刻保持不破碎:

    其中,E2是維持液滴不破碎的能量(J),可描述為如下[17]:

    σ是水液滴的表面張力(N/m),Vl是滴液的當(dāng)前體積(m3),Rl是滴液的當(dāng)前半徑(m).

    接著,E1是氣流沖擊液滴所產(chǎn)生的能量(J),可描述為[20]

    ρg為背景氣體的密度(kg/m3),Δv為液滴與背景氣體的相對(duì)速度(m/s),這里可認(rèn)為是液滴的運(yùn)動(dòng)速度.

    值得注意的是,水滴的表面張力σ是一個(gè)隨溫度變化的物理量,并且這種變化在本文的工況下是不可忽略的,尤其在多次進(jìn)行中和過程的液滴,其溫度的升高會(huì)影響表面張力的作用.因此,針對(duì)液滴的表面張力σ與液滴溫度T,本文進(jìn)行了公式耦合.圖1顯示表面張力隨液滴溫度的變化曲線以及耦合公式.由此,可根據(jù)液滴溫度T直接求得當(dāng)前液滴的表面張力.

    圖1 表面張力隨溫度的變化規(guī)律[21]Fig.1.Surface tension vs the droplet temperature[21].

    2.3.2 液滴破碎后的處理

    為實(shí)現(xiàn)計(jì)算中快速處理液滴破碎后的運(yùn)動(dòng)狀態(tài),本實(shí)驗(yàn)室特意對(duì)帶電液滴在電場(chǎng)中的破碎過程進(jìn)行了數(shù)值計(jì)算.計(jì)算模型可參考文獻(xiàn)[22],其中,在連續(xù)性的控制方程中,源項(xiàng)改為網(wǎng)格節(jié)點(diǎn)極化電荷所受到的電場(chǎng)加速度極化電荷qplr是考慮分子內(nèi)部電荷庫侖作用下的節(jié)點(diǎn)電荷數(shù)量,這部分的處理方法在泰勒錐的模擬中[23]是比較常見的,這里不再贅述.

    需要說明的是,如果電場(chǎng)均勻,帶電液滴僅僅會(huì)表現(xiàn)為整體的遷移運(yùn)動(dòng),只有當(dāng)電場(chǎng)不均勻時(shí),帶電液滴才能夠有被拉扯的作用,并且隨著電場(chǎng)不均勻系數(shù)f(網(wǎng)格內(nèi)最大電場(chǎng)強(qiáng)度Emax與平均電場(chǎng)強(qiáng)度Eave的比值)的增大,破碎將更加徹底,典型的計(jì)算結(jié)果如圖2所示.在液滴破碎過程中,始終會(huì)受到氣流沖擊作用和不均勻電場(chǎng)的拉扯作用.在破碎最開始階段,電場(chǎng)拉扯作用是可以忽略不計(jì)的,但是,隨著液滴幾何形狀的拉伸,電場(chǎng)的拉扯作用將逐漸增強(qiáng).這意味著一旦液滴開始進(jìn)入破碎過程時(shí),液滴將會(huì)在不均勻電場(chǎng)極化作用下加速破碎,迅速裂解為多個(gè)圓液滴或細(xì)長液滴,圓液滴由于表面張力重新占據(jù)優(yōu)勢(shì)而得到了維持,但細(xì)長液滴則會(huì)在電場(chǎng)極化作用下繼續(xù)發(fā)生破碎,直到裂解成能夠維持住的圓液滴為止.

    圖2 液滴破碎過程模擬(液滴初始體積8.29×10?19m3,f=1.4,溫度320 K)Fig.2.Simulation of droplet breaking up(droplet initial volume 8.29×10?19m3,f=1.4,temperature 320 K).

    圖3給出了液滴破碎更為細(xì)節(jié)的一系列數(shù)據(jù),以便我們快速設(shè)定不同情況下的液滴的破碎后的狀態(tài).這里需要說明的是,本文這種快速判斷方法是針對(duì)羽流多物理過程的復(fù)雜性,為提高計(jì)算效率而采取的粗糙、簡(jiǎn)略方法,并沒有實(shí)現(xiàn)逐個(gè)液滴破碎過程的判斷.

    然后,當(dāng)液滴破碎發(fā)生后,對(duì)液滴破碎的能量耗散進(jìn)行統(tǒng)計(jì):首先,原初液滴的表面張力能會(huì)小于各破碎后小液滴的表面張力能之和,這樣才能夠維持小液滴不破碎;接著需將兩者張力能的差值,即小液滴的動(dòng)能損失量之和,平均分?jǐn)偟礁鱾€(gè)液滴.

    其次,根據(jù)數(shù)值計(jì)算,液滴運(yùn)動(dòng)速度的角度變化不會(huì)超過10°,因此,針對(duì)每個(gè)破碎后小液滴的運(yùn)動(dòng)速度角度變化取?10°—10°之間的隨機(jī)值,液滴總質(zhì)量也采取平均分配.

    圖3 不同電場(chǎng)不均勻系數(shù)、不同液滴溫度下的液滴破碎計(jì)算結(jié)果Fig.3.Calculation results of the droplet breakage at different local f and droplet temperature.

    需要補(bǔ)充的是,在數(shù)值計(jì)算中,液滴的破碎是瞬間完成的,破碎后的所有小液滴暫時(shí)共用同一個(gè)空間坐標(biāo),但由于各個(gè)小液滴的隨機(jī)運(yùn)動(dòng)方向不同,會(huì)在下一個(gè)時(shí)間步長分散開.

    2.4 液滴與背景氣體的傳熱和流阻過程

    如前所述,液滴溫度的升高主要源于中和過程帶來的電離能沉積,大量電離能沉積在液滴,必然導(dǎo)致液滴溫度升高,當(dāng)液滴溫度高于背景氣體溫度時(shí)就會(huì)產(chǎn)生與環(huán)境的傳熱,這個(gè)過程涉及到表面熱輻射和對(duì)流換熱.但值得說明的是,周圍環(huán)境溫度和液滴溫度相差很少,輻射換熱的熱流密度將遠(yuǎn)遠(yuǎn)小于對(duì)流換熱產(chǎn)生的熱流密度,兩者會(huì)相差2—3量級(jí),故這里僅考慮背景氣體對(duì)液滴的對(duì)流換熱作用,以此作為液滴降溫的主要途徑.

    液滴在運(yùn)動(dòng)中,受到冷卻作用的降溫功率(W)為

    其中,Al為液滴表面積(m2),hl為液滴與當(dāng)?shù)乇尘皻怏w的表面換熱系數(shù)(W·m?2·K?1),該參數(shù)與液滴的運(yùn)動(dòng)速度有關(guān),T∞為當(dāng)?shù)乇尘皻怏w溫度(K).

    λ∞為背景氣體的導(dǎo)熱系數(shù)(W·m?1·K?1);d為液滴的特征尺度(m),對(duì)于球體,取半徑的1/3;Nul為努賽爾數(shù).

    液滴運(yùn)動(dòng)中由于氣體的沖擊是變形狀態(tài)的,但這里假設(shè)液滴形狀為剛性球體,采用流體外掠圓球的模型[24],那么努賽爾數(shù)Nul可以表示為

    其中,Re為雷諾數(shù),Pr為普朗特?cái)?shù),兩者均屬于無量綱數(shù);η∞和ηl分別為以背景氣體溫度和以液滴溫度來計(jì)算的動(dòng)力黏度(N·s·m?2).

    由于液滴特征尺寸較小,在102m/s量級(jí)的運(yùn)動(dòng)速度下,雷諾數(shù)Re極小,因此,這里考慮采用低雷諾數(shù)下的剛性球體的阻力計(jì)算公式[25]:

    其中,const為表征物體形狀的常數(shù),球體取0.4—0.5,根據(jù)文獻(xiàn)[25],在低雷諾數(shù)下,const取0.5為宜;S為液滴球體的橫截面積(m2).

    3 驗(yàn)證試驗(yàn)

    為驗(yàn)證NECD模型的正確性,開展超聲波電噴推力器的羽流中和試驗(yàn),采用推力測(cè)量?jī)x器進(jìn)行推力的測(cè)量,采用高速相機(jī)對(duì)液滴大小進(jìn)行測(cè)量,并且以上述兩者和計(jì)算結(jié)果進(jìn)行比對(duì)來驗(yàn)證模型的計(jì)算精度.需要說明的是,推力值的驗(yàn)證目的是驗(yàn)證液滴輸運(yùn)過程的計(jì)算精度,即以大量運(yùn)動(dòng)液滴所形成的推力來衡量液滴運(yùn)動(dòng)速度的捕捉精確程度;液滴粒徑大小的驗(yàn)證目的是驗(yàn)證液滴輸運(yùn)過程中破碎與重組的捕捉精確程度.那么,通過以上兩項(xiàng)驗(yàn)證,基本可以實(shí)現(xiàn)NECD模型的計(jì)算精度評(píng)估.

    試驗(yàn)系統(tǒng)圖如圖4所示:超聲波電噴推力器固定于轉(zhuǎn)動(dòng)機(jī)構(gòu)的右側(cè)托板上,托板則由三根支撐柱進(jìn)行固定,并且上下支撐柱可以自由伸縮,中間的支撐柱與托板由轉(zhuǎn)動(dòng)節(jié)點(diǎn)連接,可以令托板在?15°—15°之間進(jìn)行轉(zhuǎn)動(dòng),由此調(diào)整推力器的羽流噴射方向.

    圖4 驗(yàn)證試驗(yàn)系統(tǒng)示意圖Fig.4.Schematic diagram of the experiment system.

    需要說明的是,本文沒有采用傳統(tǒng)的推力架來進(jìn)行推力的測(cè)量,原因如下:超聲波推力器中的超聲波發(fā)生器對(duì)推力測(cè)量系統(tǒng)有很大的震動(dòng)影響,導(dǎo)致測(cè)量結(jié)果波動(dòng)極大,本文曾經(jīng)測(cè)量過推力器正常工作時(shí)的推力,然而,測(cè)量結(jié)果在0.1×10?6N—80×10?6N之間不穩(wěn)定振蕩,已經(jīng)無法獲取準(zhǔn)確數(shù)值.

    所以,本文采用“打靶法”的測(cè)量思路.“打靶法”推力測(cè)量系統(tǒng)由懸臂薄板、電渦流傳感器和示波器組成.首先,提出軸向有效推力的概念,即羽流產(chǎn)生實(shí)際推力矢量在推力器設(shè)計(jì)中心軸線方向上的分量,為該推力器的軸向有效推力.那么,推力測(cè)量系統(tǒng)測(cè)得的推力實(shí)際上是有效推力,而不是推力器的真實(shí)推力,因?yàn)槿魏瓮屏ζ鞯耐屏Ψ较蚺c中心軸線(設(shè)計(jì)方向)都存在微小偏角.于是,有效推力為F0的羽流對(duì)懸臂薄板會(huì)產(chǎn)生一個(gè)沖擊力P,以P來描述F0,有

    其中,sinβ為羽流對(duì)懸臂薄板產(chǎn)生的沖擊形變角度,由于推力較低,β一般可認(rèn)為是90°;沖擊力P會(huì)對(duì)懸臂薄板產(chǎn)生一個(gè)微小位移d,有

    其中,E為懸臂薄板的彈性模量(MPa),w為薄板的橫向?qū)挾?m),t為薄板的厚度(m),L為薄板的縱向長度(m).那么,只要測(cè)得位移d的數(shù)值,就可以推導(dǎo)出有效推力F0.而電渦流傳感器正是用來測(cè)量微小位移的裝置,傳感器利用感應(yīng)線圈來產(chǎn)生感應(yīng)電流,感應(yīng)電流與位移d的大小有關(guān),由此,傳感器可以將位移的大小轉(zhuǎn)化為電流的高低來讀取位移d的數(shù)值.傳感器的性能參數(shù)如表1所列.

    表1 電渦流傳感器技術(shù)指標(biāo)Table 1.Technical parameters of the eddy current sensor.

    需要說明的是,由于超聲波羽流的波動(dòng)以及放電振蕩等因素的影響,由傳感器所讀取的電流數(shù)據(jù)實(shí)際上是一個(gè)范圍,例如,工況2所測(cè)得的電壓值為4.1 mV—8.9 mV,本文取平均值來表示實(shí)際值,即電壓值為6.5 mV,所測(cè)推力為15.2×10?6N.

    圖5 液滴粒徑的測(cè)量方法Fig.5.Measuring method of the drop diameter.

    此外,在正常布置試驗(yàn)系統(tǒng)時(shí),會(huì)將推力器中心軸線與懸臂梁板垂直放置,但由于推力器偏角的存在,所測(cè)得的推力F1一定是偏小的;那么,在不斷調(diào)整轉(zhuǎn)動(dòng)機(jī)構(gòu)時(shí),一定存在一個(gè)轉(zhuǎn)動(dòng)角度θ,使得此刻測(cè)得的推力值Fmax為最大.那么,可認(rèn)為F1/Fmax的反余弦值為推力偏角.事實(shí)上,推力偏角與轉(zhuǎn)動(dòng)角度θ有差別,本文以羽流達(dá)到懸臂薄板時(shí)的推力狀態(tài)為準(zhǔn),所以推力偏角不采用轉(zhuǎn)動(dòng)角度θ.

    最后,液滴的粒徑由高速相機(jī)來進(jìn)行測(cè)量,由于液滴具有一定的運(yùn)動(dòng)速度以及粒徑在的10?6m的尺度,拍攝到清晰照片的難度較高.因此,本文采取一個(gè)對(duì)相片做后處理的方法來獲取液滴的粒徑,如圖5所示.

    3.1 推力的驗(yàn)證

    試驗(yàn)輸入條件設(shè)置見表2.在表2參數(shù)條件下,測(cè)量出三組推力,與相應(yīng)工況的計(jì)算結(jié)果進(jìn)行比對(duì),見圖6.

    根據(jù)圖6對(duì)于推力計(jì)算的相對(duì)誤差在16%左右,推力偏角的計(jì)算誤差在22%左右,并且隨著工況的變化,變化趨勢(shì)上可以保持一致.圖6仿真結(jié)果表現(xiàn)出高于試驗(yàn)結(jié)果,推測(cè)是由于空間電勢(shì)求解模型和破碎模型的誤差放大了羽流電荷分布的不均勻性,導(dǎo)致推力偏角的升高.

    表2 各工況的輸入?yún)?shù)Table 2.Parameters of the test operation in different cases.

    圖6 推力測(cè)量的試驗(yàn)與計(jì)算結(jié)果對(duì)比Fig.6.Comparison of the thrust between the test and calculation results.

    3.2 液滴粒徑的驗(yàn)證

    為驗(yàn)證液滴粒徑的計(jì)算結(jié)果,首先對(duì)不同位置的羽流進(jìn)行了高速相機(jī)的拍照,獲取3組典型位置下的液滴照片,如圖7所示,各照片中,按照液滴從小到大的順序進(jìn)行了編號(hào),圖中參數(shù)l為與羽流來流面的水平距離(該距離最大為3 m).

    圖7 不同羽流位置的液滴照片 (a)l=0.3—0.31 m;(b)l=1.50—1.51 m;(c)l=2.70—2.71 mFig.7.Photos of droplets in different plume position:(a)l=0.3—0.31 m;(b)l=1.50—1.51 m;(c)l=2.70—2.71 m.

    接著,在每一處羽流位置上,從計(jì)算結(jié)果中抽取10個(gè)液滴,以清晰度較高的液滴為優(yōu)先,這樣的抽取方法確實(shí)受到了相機(jī)性能的限制.依然按照從小到大的順序編號(hào),以液滴粒徑的計(jì)算值與試驗(yàn)值進(jìn)行對(duì)比,對(duì)比的結(jié)果見圖8.雖然計(jì)算結(jié)果與試驗(yàn)結(jié)果有誤差不可避免,但在不同羽流位置上的變化趨勢(shì)可以保持一致.

    圖8 液滴粒徑的試驗(yàn)值與計(jì)算值的對(duì)比 (a)l=0.3—0.31 m;(b)l=1.50—1.51 m;(c)l=2.70—2.71 mFig.8.Comparison of droplet diameter between the test and calculation results:(a)l=0.3—0.31 m;(b)l=1.50—1.51 m;(c)l=2.70— 2.71 m.

    至此,通過對(duì)推力、液滴粒徑大小的驗(yàn)證,認(rèn)為該模型的綜合計(jì)算誤差在20%左右,在趨勢(shì)上可以與試驗(yàn)結(jié)果保持一致,可滿足定性分析的需要.但定量計(jì)算精度并不算很高,原因?yàn)槟P捅旧硭紤]的子模型較多,進(jìn)行了多項(xiàng)假設(shè).這里,對(duì)NECD模型影響精度的一些重要假設(shè)進(jìn)行了歸納:

    1)電子與液滴碰撞時(shí),電荷對(duì)網(wǎng)格內(nèi)液滴平均分配的假設(shè);

    2)忽略三個(gè)或三個(gè)以上液滴碰撞的假設(shè);

    3)對(duì)液滴破碎后的狀態(tài)進(jìn)行快速判斷的假設(shè);

    4)對(duì)運(yùn)動(dòng)中液滴的傳熱與流阻計(jì)算采用剛性球體的假設(shè)(液滴運(yùn)動(dòng)時(shí)存在一定的形變).

    上述假設(shè)都是提高計(jì)算效率所必要采取的措施,但即使這樣,采用16線程、32 GB內(nèi)存的服務(wù)器,計(jì)算一個(gè)算例依然需要100 h以上.實(shí)際上,它們確實(shí)對(duì)計(jì)算精度產(chǎn)生了相當(dāng)大的影響,尤其是第3)項(xiàng)假設(shè).然而,若未來計(jì)算機(jī)的計(jì)算速度大幅度提高時(shí),可以對(duì)上述假設(shè)進(jìn)行更精細(xì)化的處理,以此提高此模型的計(jì)算精度.

    4 計(jì)算結(jié)果與分析

    為解釋前文提出的推力偏角和能損耗較高的問題,本節(jié)首先給出液滴的電荷密度、數(shù)密度以及速度分布的云圖,以此說明液滴產(chǎn)生推力偏角的原因;接著,給出說明液滴破碎與重組情況的液滴體積分布和液滴溫度分布云圖,以此說明液滴產(chǎn)生巨大能耗的原因.最后,總結(jié)出超聲波電噴推力器羽流的帶電粒子中和的規(guī)律.計(jì)算的輸入條件見表3.需要說明的是,計(jì)算區(qū)域覆蓋徑向(X方向)0.36 m,軸向(Z方向)3 m的范圍,而羽流區(qū)域僅在計(jì)算域中X=0.09—0.27 m之間.

    超聲波電噴推力器的羽流中和過程比較特殊:液滴的碰撞截面較大,很容易將大多數(shù)的電子攔截下來,因此,正負(fù)電荷的中和過程會(huì)表現(xiàn)出一定的不均勻性.如圖9所示,在左側(cè)最接近陰極的位置,羽流上部分的帶電量會(huì)迅速下降;接著會(huì)有大片的液滴逐漸由正電荷過渡到中性,再過渡到負(fù)電荷,出現(xiàn)了過量中和的特點(diǎn),如圖9的藍(lán)色區(qū);而距離陰極較遠(yuǎn)的羽流下部分會(huì)表現(xiàn)出中和較慢的現(xiàn)象,依然保持一定的正電性.

    表3 計(jì)算輸入條件Table 3.Input conditions of the calculation.

    圖9 液滴電荷密度分布云圖Fig.9.Distribution contour of the droplet charge density.

    在推力器的中和過程中,只要液滴從加速電極運(yùn)動(dòng)出,就會(huì)從加速階段轉(zhuǎn)變?yōu)闇p速階段,開始進(jìn)行動(dòng)能對(duì)電勢(shì)能的返還,但如果盡早與電子中和,則會(huì)不受電場(chǎng)的影響而保留動(dòng)能繼續(xù)運(yùn)動(dòng),形成最終的推力.因此,在圖9的結(jié)果下,由于過度中和效應(yīng),羽流上部分的液滴將具有更高的動(dòng)能,如圖10所示,羽流上部分的液滴速度將比下部分高很多.

    同時(shí),液滴運(yùn)動(dòng)越快的位置,數(shù)密度將越低,反之?dāng)?shù)密度將越高.因此,圖11從側(cè)面顯示了羽流上下部分的液滴速度和數(shù)量的不一致,進(jìn)一步證明了羽流中和的不均勻性.

    圖10 液滴速度分布云圖Fig.10.Distribution contour of the droplet velocity.

    圖11 液滴數(shù)密度分布云圖Fig.11.Distribution contour of the droplet number density.

    因此,在圖9—圖11的計(jì)算結(jié)果中,可以推導(dǎo)出羽流在中和過程中所形成的推力并不與液滴運(yùn)動(dòng)方向完全一致,而是會(huì)形成一個(gè)與陰極中和器背離的偏角.

    接著,給出液滴在羽流輸運(yùn)過程中的能耗統(tǒng)計(jì),見表4.

    表4 羽流輸運(yùn)中的各類能耗比例Table 4.Each energy loss proportion of the plume.

    這里,熱損耗主要是指高溫液滴與環(huán)境氣體之間的對(duì)流換熱損失;液滴破碎能耗是指一個(gè)大液滴破碎成小液滴后,總表面張力能上升而導(dǎo)致動(dòng)能下降的損失;流動(dòng)阻力損耗是指高速液滴受到大氣阻力的能耗;彈性碰撞損耗是指不同速度方向的液滴在碰撞后,由動(dòng)量守恒所導(dǎo)致的動(dòng)能損失.需要說明的是,計(jì)算得到的總功率損耗(42.63%)與試驗(yàn)結(jié)果(53%)依然有偏差,除了計(jì)算模型的誤差外,推測(cè)可能與陰極針通過遂穿效應(yīng)發(fā)射電子和發(fā)射表面通過泰勒錐發(fā)射正電液滴所需要的電場(chǎng)做功并未考慮在本計(jì)算模型中有關(guān).

    進(jìn)一步地,推力器能效較低的主要原因是液滴破碎發(fā)生頻率較高,故討論其內(nèi)在機(jī)理:在羽流中和與輸運(yùn)過程中,液滴主要會(huì)由三方面的因素來影響破碎:1)較高的運(yùn)動(dòng)速度;2)較高的液滴溫度;3)較明顯的電場(chǎng)不均勻程度.

    而在陰極中和位置處羽流區(qū)域的液滴速度較羽流中后部較高(此為因素1),并且液滴發(fā)生電子中和較為集中,電離能沉積導(dǎo)致溫度高(見圖12)(此為因素2),此外,該區(qū)域空間電荷的變化而導(dǎo)致電場(chǎng)變化極快,電場(chǎng)比較不均勻(此為因素3).因此,在陰極出口位置處的羽流區(qū)會(huì)形成破碎最容易、也是最集中的區(qū)域,這也是羽流產(chǎn)生如此多液滴破碎的主要原因,如圖13的羽流前部分,單個(gè)液滴的平均體積出現(xiàn)急劇下降.此外,圖13的羽流中后部的云圖揭示:正負(fù)帶電液滴的重組過程也伴隨著少量的破碎,但多數(shù)液滴的體積是在逐漸增大的,依然以重組過程為主.

    圖12 液滴溫度分布云圖Fig.12.Distribution contour of the droplet temperature.

    需要說明的是,液滴破碎是基于推力器羽流中和過程的特殊性而存在的固有過程:1)存在中和就涉及中和能的沉積,就會(huì)導(dǎo)致液滴溫度的劇增;2)要保證推力性能就必須保證液滴具有相當(dāng)?shù)乃俣?3)中和集中的區(qū)域始終是電場(chǎng)不均勻性最強(qiáng)的區(qū)域.基于這三點(diǎn)因素,液滴破碎的能損很難降低到較小的水平,對(duì)于工質(zhì)水來說,通常都在20%之上.

    綜上所述,可見超聲波電噴的羽流中和過程是非常特殊、復(fù)雜的物理過程,隱含典型的電子-液滴的中和特性,不但涉及帶電粒子的輸運(yùn)、碰撞,還涉及液滴的破碎與重組,以及傳熱過程.圖14對(duì)超聲波電噴推力器的羽流中和特性進(jìn)行了概述:在羽流靠近陰極出口區(qū)域,發(fā)生正負(fù)電荷中和以及電子吸附過程,一方面,由于中和的不均勻性,表現(xiàn)為上部羽流呈負(fù)電性,下部羽流呈正電性,導(dǎo)致羽流上下部分的液滴運(yùn)動(dòng)速度不一致,會(huì)造成推力偏角;另一方面,由于靠近陰極出口的羽流區(qū)液滴的中和能沉積較大,液滴運(yùn)動(dòng)速度較快,是液滴發(fā)生破碎的集中區(qū),大量破碎過程導(dǎo)致了較高的能量損耗,但在隨后的羽流輸運(yùn)過程中,液滴的速度、溫度逐漸過渡到較低的水平,因而在羽流中后部分,主要發(fā)生液滴的重組過程.

    圖13 單個(gè)液滴體積的平均值分布云圖Fig.13.Distribution contour of the mean volume of one droplet.

    圖14 超聲波電噴推力器羽流中和特性示意圖Fig.14.Schematic diagram of the neutralization characteristics in UAET plume.

    5 結(jié) 論

    本文揭示了超聲波電噴推力器的羽流中和過程中的各種物理機(jī)理內(nèi)涵,在工程意義上,為推力器工作性能的優(yōu)化提供了一定的參考;同時(shí),在學(xué)術(shù)意義上,揭示出電子-液滴的中和特性與電子-離子中和特性截然不同,表現(xiàn)出多種物理過程耦合的復(fù)雜特性,為其他與電子-液滴中和相關(guān)的機(jī)理研究提供了參考理論.主要結(jié)論如下:

    1)超聲波電噴推力器羽流的中和過程會(huì)存在較高的電荷分布、速度分布的不均勻性,這是引起推力偏角的直接原因;

    2)超聲波電噴推力器的功率損耗主要集中在液滴破碎和熱流失兩個(gè)方面,其中,以液滴破碎所占比例最大,對(duì)于水工質(zhì)來說,能耗通常在20%以上.

    進(jìn)而,本文對(duì)超聲波電噴推力器的后續(xù)研發(fā)給出一些啟示:

    A)若采用電子中和器,可在其他方向進(jìn)行對(duì)稱布置,以此改善羽流中電荷的中和不均勻性,或者采用兩臺(tái)發(fā)射液滴帶電性相反的推力器進(jìn)行互和,但上述中和策略對(duì)兩種推力器的性能一致性要求較高;

    B)為提高推力器的能效,從降低液滴破碎頻率的角度來說,可降低發(fā)射液滴的粒徑或者更換表面張力系數(shù)較大的工質(zhì),但這樣的設(shè)計(jì)需要對(duì)液滴發(fā)射裝置進(jìn)行相匹配的改進(jìn).

    [1]Zhao Q,Huang X P,Lin E,Jiao J,Liang G F,Chen T 2017Opto-Electronic Engineer.44 140

    [2]Jiao J,Zhao Q,Li X,Liang G F,Huang X P,Luo X 2014Opt.Express22 26277

    [3]Zhao Y,Huang C,Qing A Y,Luo X 2017IEEE Photon.J.99 1

    [4]Taylor G 1964Proc.Roy.Soc.Lond.A280 383

    [5]Romero S,Bocanegra R,Gamero C 2003J.Appl.Phys.94 3599

    [6]Lozano P,Martinez S 200541st Joint Propulsion Conference&ExhibitTuscon,Arizona.July 10–13,2005 p1

    [7]Ober S,Branam R,Huffman R 201149thAIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace ExpositionOrland,Florida,January 4–7,2011 p1

    [8]Legge R,Lozano P 2011J.Propuls.Power27 485

    [9]Reading C,Anderson J,Kubiak C,Greer F,Rouhi N,Wilson D,White V,Dickie M,Mueller R,Singh V,Mackie W,Wirz R,Castano M 2016AIAA Propulsion and Energy ForumSalt Lake City,UT,July 25–27,2016 p1

    [10]Kurt J,Lyon B 201652nd AIAA/SAE/ASEE Joint Propulsion ConferenceSalt Lake City,UT,July 25–27,2016 p1

    [11]Gutierrez E,Castano M 2017J.Propuls.Power33 984

    [12]Song W D,Shumlak U.2010J.Propuls.Power26 353

    [13]Dong L,Song W D,Kang X M,Zhao W S 2012Acta Astron.77 1

    [14]Zhang Y B,Hang G R,Dong L,Kang X M,Zhao W S,Zhang Y,Kang X L 2016Chin.Space Sci.Technol.36 9(in Chinese)[張姚濱,杭觀榮,董磊,康小明,趙萬生,張巖,康小錄2016中國空間科學(xué)技術(shù)36 9]

    [15]Kang X M,Dong L,Zhao W S 2014Acta Astron.98 1

    [16]Passaro A,Nania F,Vicini A 200637th AIAA Plasmadynamics and Lasers ConferenceSan Francisco,California,June 5–8,2006 p1

    [17]Robert S,Eduardo A 200945th Joint Propulsion Conference&ExhibitDenver,Colorado August 2–5,2009 p1

    [18]Bird G 1963Gas.Phys.Fluids6 1518

    [19]Jayaratne O,Mason B 1974Proc.Roy.Soc.Lond.380 218

    [20]Luo T Q,Wang X Y,Zheng J Q,Wang Z T,Mao H M 2007Drainage and Irrigation Machinery25 57(in Chinese)[羅惕乾,王曉英,鄭捷慶,王貞濤,毛惠敏2007排灌機(jī)械25 57]

    [21]Gao S Q,Liu H P 2010Capillary Mechanics(Beijing:Science Press)p60(in Chinese)[高世橋,劉海鵬 2010毛細(xì)力學(xué)(北京:科學(xué)出版社)第60頁]

    [22]Cai B,Lee L,Wang Z L 2003J.Engineer.Thermophys.24 613(in Chinese)[蔡斌,李磊,王照林2003工程熱物理學(xué)報(bào)24 613]

    [23]Higuera F 2003J.Fluid Mech.484 303

    [24]Yang S M,Tao W Q 2006Heat Transfer(4th Ed.)(Beijing:Higher Education Press)p258(in Chinese)[楊世銘,陶文銓2006傳熱學(xué)(第四版)(北京:高等教育出版社)第258頁]

    [25]Landau L(translated by Lee Z)2013Fluid Dynamics(5th Ed.)(Beijing:Higher Education Press)pp201–202(in Chinese)[朗道L著 (李植 譯)2013流體動(dòng)力學(xué)(第五版)(北京:高等教育出版社)第201—202頁]

    猜你喜歡
    電噴羽流推力器
    單組元推力器倒置安裝多余物控制技術(shù)驗(yàn)證
    水下羽流追蹤方法研究進(jìn)展
    電噴汽車發(fā)動(dòng)機(jī)啟動(dòng)故障的診斷與排除探究
    水下管道向下泄漏的羽/射流特性
    用于小行星探測(cè)的離子推力器技術(shù)研究
    離子推力器和霍爾推力器的異同
    太空探索(2015年6期)2015-07-12 12:48:42
    東風(fēng)雪鐵龍C5各電控系統(tǒng)電路圖解析(七):發(fā)動(dòng)機(jī)電噴系統(tǒng)電路
    船用電噴主機(jī)的原理及日常管理淺析
    科技資訊(2014年12期)2014-11-14 07:54:03
    室內(nèi)多股羽流混合運(yùn)動(dòng)機(jī)理模型研究進(jìn)展分析
    固體微型推力器應(yīng)用設(shè)計(jì)
    航天器工程(2012年6期)2012-12-29 04:13:44
    久久久久精品久久久久真实原创| 欧美性猛交╳xxx乱大交人| 2021少妇久久久久久久久久久| 干丝袜人妻中文字幕| av专区在线播放| 插逼视频在线观看| 亚洲怡红院男人天堂| 久久99热这里只频精品6学生| 五月天丁香电影| 亚洲av在线观看美女高潮| 人妻制服诱惑在线中文字幕| av播播在线观看一区| a级毛片免费高清观看在线播放| 色婷婷久久久亚洲欧美| 可以在线观看毛片的网站| 精品熟女少妇av免费看| 国产一区二区三区av在线| 亚洲成人精品中文字幕电影| videos熟女内射| 国产又色又爽无遮挡免| 亚洲欧美一区二区三区黑人 | 免费观看在线日韩| 赤兔流量卡办理| 在线观看一区二区三区激情| 伦精品一区二区三区| 欧美精品人与动牲交sv欧美| 黄片wwwwww| 偷拍熟女少妇极品色| 性色avwww在线观看| 少妇人妻精品综合一区二区| av在线老鸭窝| 国产精品一区二区三区四区免费观看| 亚洲av日韩在线播放| 汤姆久久久久久久影院中文字幕| 亚洲精品成人久久久久久| 亚洲av免费在线观看| 99九九线精品视频在线观看视频| 亚洲精品国产av成人精品| 日本一二三区视频观看| 日本三级黄在线观看| 免费看光身美女| 久久精品人妻少妇| 日韩三级伦理在线观看| 欧美日韩视频高清一区二区三区二| 18禁在线无遮挡免费观看视频| 下体分泌物呈黄色| 久久久久久久国产电影| 热re99久久精品国产66热6| 国产大屁股一区二区在线视频| 亚洲av中文字字幕乱码综合| 黄色一级大片看看| 午夜福利高清视频| 最新中文字幕久久久久| 男女边吃奶边做爰视频| 男男h啪啪无遮挡| 国产 精品1| 亚洲人成网站高清观看| 国产欧美亚洲国产| 两个人的视频大全免费| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 男女边吃奶边做爰视频| av国产久精品久网站免费入址| 免费黄色在线免费观看| 精品99又大又爽又粗少妇毛片| av线在线观看网站| 精品视频人人做人人爽| 国产成人午夜福利电影在线观看| 国产精品不卡视频一区二区| 波多野结衣巨乳人妻| 狠狠精品人妻久久久久久综合| 久久精品夜色国产| 久久6这里有精品| 亚洲内射少妇av| 久久久久久久久久人人人人人人| 大香蕉久久网| 男女啪啪激烈高潮av片| 久久99热这里只有精品18| 青春草亚洲视频在线观看| 久久久久九九精品影院| 亚洲图色成人| 中国美白少妇内射xxxbb| 午夜福利高清视频| 亚洲精品乱久久久久久| 69人妻影院| 国产淫片久久久久久久久| 久久精品夜色国产| 青春草国产在线视频| 又黄又爽又刺激的免费视频.| 人妻 亚洲 视频| 国产免费又黄又爽又色| 亚洲欧美成人精品一区二区| 久久精品国产亚洲av涩爱| 国产高清不卡午夜福利| 成人毛片60女人毛片免费| 日产精品乱码卡一卡2卡三| av国产免费在线观看| 亚洲国产欧美人成| 日本黄色片子视频| 岛国毛片在线播放| 欧美精品国产亚洲| 国产成人aa在线观看| 国产色婷婷99| 亚洲第一区二区三区不卡| 欧美bdsm另类| 精品少妇久久久久久888优播| 国语对白做爰xxxⅹ性视频网站| 青春草国产在线视频| 国产精品久久久久久久久免| 久热久热在线精品观看| 在线观看一区二区三区激情| 一个人看视频在线观看www免费| 夫妻性生交免费视频一级片| 亚洲精品视频女| 亚洲精品成人久久久久久| 久久久久久伊人网av| 噜噜噜噜噜久久久久久91| 亚洲婷婷狠狠爱综合网| 亚洲最大成人手机在线| 精品一区二区三区视频在线| av线在线观看网站| 国产精品99久久久久久久久| 我的老师免费观看完整版| 久久人人爽人人爽人人片va| 欧美日韩亚洲高清精品| 亚洲国产高清在线一区二区三| 2018国产大陆天天弄谢| 亚洲精品国产av蜜桃| 久久久久久久久久人人人人人人| 日韩av免费高清视频| 韩国av在线不卡| 看免费成人av毛片| 精品久久久久久久久av| 三级国产精品片| 亚洲最大成人中文| 少妇人妻久久综合中文| 51国产日韩欧美| 国产成人freesex在线| 菩萨蛮人人尽说江南好唐韦庄| 色婷婷久久久亚洲欧美| 九草在线视频观看| 少妇人妻久久综合中文| 99热这里只有是精品在线观看| 亚洲最大成人中文| 日韩三级伦理在线观看| 国产v大片淫在线免费观看| 国产男女超爽视频在线观看| 亚洲精品国产色婷婷电影| 干丝袜人妻中文字幕| 国产精品一区二区在线观看99| 亚洲av免费在线观看| 久久精品夜色国产| 欧美区成人在线视频| 欧美国产精品一级二级三级 | 国内精品宾馆在线| 亚洲国产精品999| 成人国产麻豆网| 国产老妇伦熟女老妇高清| 91精品一卡2卡3卡4卡| 26uuu在线亚洲综合色| 免费观看的影片在线观看| 午夜视频国产福利| 亚洲av.av天堂| 亚洲成人久久爱视频| 久久久久久久久久成人| 精品午夜福利在线看| 久久久久性生活片| 国产男女超爽视频在线观看| 久久人人爽av亚洲精品天堂 | 夫妻性生交免费视频一级片| 深爱激情五月婷婷| 99热这里只有是精品50| 亚洲精品成人久久久久久| 日本一二三区视频观看| 黄色一级大片看看| 日韩强制内射视频| av在线老鸭窝| 在线观看一区二区三区| 三级国产精品片| .国产精品久久| 少妇的逼好多水| 韩国av在线不卡| 国产黄片美女视频| 久久精品综合一区二区三区| 免费黄频网站在线观看国产| 26uuu在线亚洲综合色| 色婷婷久久久亚洲欧美| 久久久久国产网址| 丰满乱子伦码专区| 伊人久久国产一区二区| 国产一区二区三区av在线| 男插女下体视频免费在线播放| kizo精华| 国产成人精品婷婷| 别揉我奶头 嗯啊视频| 在线观看av片永久免费下载| 午夜老司机福利剧场| 国产精品国产av在线观看| 26uuu在线亚洲综合色| 久久久久精品性色| 国产成人freesex在线| 2021少妇久久久久久久久久久| 美女cb高潮喷水在线观看| 99热这里只有是精品在线观看| 男女无遮挡免费网站观看| 一个人看的www免费观看视频| 麻豆久久精品国产亚洲av| av又黄又爽大尺度在线免费看| 婷婷色综合www| a级毛片免费高清观看在线播放| 亚洲av福利一区| 免费观看的影片在线观看| 国产日韩欧美亚洲二区| 国产毛片a区久久久久| 久久精品人妻少妇| 国产精品三级大全| 青青草视频在线视频观看| 国产在视频线精品| 亚洲久久久久久中文字幕| 国产精品人妻久久久久久| 国产高清不卡午夜福利| 少妇丰满av| 综合色av麻豆| 春色校园在线视频观看| 国产成人精品一,二区| 伦理电影大哥的女人| 亚洲精品,欧美精品| 日本猛色少妇xxxxx猛交久久| 久久6这里有精品| 免费看av在线观看网站| 中文天堂在线官网| 一本久久精品| 看十八女毛片水多多多| 天美传媒精品一区二区| 人人妻人人看人人澡| 国产精品一及| 蜜桃久久精品国产亚洲av| 午夜日本视频在线| 久久久色成人| 久久久久久久久久久免费av| 日韩成人av中文字幕在线观看| 日韩电影二区| 精品国产乱码久久久久久小说| 天美传媒精品一区二区| 国产久久久一区二区三区| 国产老妇女一区| 一级爰片在线观看| 欧美日本视频| 2022亚洲国产成人精品| 国产精品女同一区二区软件| 国产老妇伦熟女老妇高清| 欧美成人午夜免费资源| 日韩,欧美,国产一区二区三区| 亚洲精品aⅴ在线观看| 视频区图区小说| 综合色丁香网| 韩国av在线不卡| 久久精品熟女亚洲av麻豆精品| 99视频精品全部免费 在线| 少妇人妻久久综合中文| 69人妻影院| 亚洲精品第二区| 尾随美女入室| 成人特级av手机在线观看| 中文字幕人妻熟人妻熟丝袜美| 深夜a级毛片| 国产美女午夜福利| 日韩欧美 国产精品| 亚洲国产最新在线播放| 波野结衣二区三区在线| 亚洲色图综合在线观看| 亚洲在久久综合| 麻豆国产97在线/欧美| 亚洲天堂国产精品一区在线| 视频区图区小说| 国产 精品1| 不卡视频在线观看欧美| 少妇人妻一区二区三区视频| 国产精品成人在线| 欧美激情在线99| 国产高清三级在线| 日韩av不卡免费在线播放| 一区二区三区精品91| 久久久久久久午夜电影| 亚洲欧美清纯卡通| 国产成人午夜福利电影在线观看| 精品人妻视频免费看| 大码成人一级视频| 久久99热6这里只有精品| 欧美区成人在线视频| 国产精品三级大全| 在线 av 中文字幕| 午夜视频国产福利| 亚洲综合色惰| 亚洲精品成人av观看孕妇| 网址你懂的国产日韩在线| 啦啦啦在线观看免费高清www| 乱系列少妇在线播放| av在线亚洲专区| 欧美bdsm另类| 在线 av 中文字幕| 看免费成人av毛片| 日韩三级伦理在线观看| 热re99久久精品国产66热6| .国产精品久久| 久久久久精品久久久久真实原创| 午夜福利视频1000在线观看| 一级爰片在线观看| 在线看a的网站| 一级毛片 在线播放| 国产亚洲精品久久久com| 精品久久久噜噜| 天天躁夜夜躁狠狠久久av| 听说在线观看完整版免费高清| 性色av一级| 国产精品国产av在线观看| 中文欧美无线码| 亚洲精品久久午夜乱码| 精品99又大又爽又粗少妇毛片| 白带黄色成豆腐渣| 亚洲av男天堂| 午夜福利视频1000在线观看| 亚洲在线观看片| 亚洲精品日韩在线中文字幕| 亚洲色图综合在线观看| 国产精品不卡视频一区二区| 亚洲av电影在线观看一区二区三区 | 亚洲国产色片| 特级一级黄色大片| 美女国产视频在线观看| 久久久久国产精品人妻一区二区| 日本一本二区三区精品| 黄片无遮挡物在线观看| 久久久久九九精品影院| 熟女av电影| 国产成人免费观看mmmm| 韩国av在线不卡| 亚洲成色77777| 国产亚洲91精品色在线| 狠狠精品人妻久久久久久综合| 大香蕉久久网| 精品久久久久久久久av| av女优亚洲男人天堂| 亚洲精品一二三| av国产免费在线观看| 大香蕉97超碰在线| 欧美97在线视频| 久久99蜜桃精品久久| 欧美zozozo另类| 亚洲国产精品国产精品| 亚洲精品视频女| 精品久久久精品久久久| 日韩在线高清观看一区二区三区| 在线观看一区二区三区| 亚洲欧美中文字幕日韩二区| 国产午夜精品一二区理论片| 亚洲最大成人中文| 日本wwww免费看| 国产乱来视频区| 亚洲精品aⅴ在线观看| 欧美成人a在线观看| 亚洲,欧美,日韩| 自拍欧美九色日韩亚洲蝌蚪91 | 国产精品秋霞免费鲁丝片| 中文乱码字字幕精品一区二区三区| 岛国毛片在线播放| 国产精品三级大全| 久久精品国产亚洲av天美| 欧美亚洲 丝袜 人妻 在线| 欧美成人一区二区免费高清观看| 特级一级黄色大片| 噜噜噜噜噜久久久久久91| 国产淫语在线视频| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 大香蕉久久网| 亚洲国产欧美人成| 在线观看人妻少妇| 婷婷色综合大香蕉| 免费在线观看成人毛片| 亚洲国产高清在线一区二区三| 久久6这里有精品| 亚洲av中文字字幕乱码综合| 亚洲综合色惰| 国产乱来视频区| 亚洲电影在线观看av| 欧美性感艳星| 精品亚洲乱码少妇综合久久| 纵有疾风起免费观看全集完整版| 嫩草影院精品99| 亚洲天堂国产精品一区在线| 久久热精品热| 日韩强制内射视频| 欧美成人精品欧美一级黄| 99久久九九国产精品国产免费| 97超视频在线观看视频| 久久99热6这里只有精品| 老司机影院毛片| 国产精品三级大全| 国产亚洲av片在线观看秒播厂| 亚洲天堂av无毛| 晚上一个人看的免费电影| 五月伊人婷婷丁香| 直男gayav资源| 亚洲精品成人久久久久久| 最近的中文字幕免费完整| 国产免费一区二区三区四区乱码| 成人综合一区亚洲| 亚洲精品日韩av片在线观看| 午夜爱爱视频在线播放| 免费观看av网站的网址| 欧美一区二区亚洲| 男人添女人高潮全过程视频| 97超视频在线观看视频| 成人二区视频| 精品人妻偷拍中文字幕| 一区二区三区乱码不卡18| 看十八女毛片水多多多| 欧美国产精品一级二级三级 | 欧美三级亚洲精品| 大片免费播放器 马上看| 欧美激情久久久久久爽电影| av一本久久久久| 亚洲图色成人| 国产毛片a区久久久久| 亚洲内射少妇av| 亚洲精品国产av蜜桃| 午夜福利在线在线| 在线 av 中文字幕| 久久久久久久久久人人人人人人| 超碰av人人做人人爽久久| 99久久精品一区二区三区| 成人国产av品久久久| 亚洲怡红院男人天堂| 国产一级毛片在线| 神马国产精品三级电影在线观看| 国产亚洲一区二区精品| 国产精品av视频在线免费观看| 中文字幕人妻熟人妻熟丝袜美| 国产男人的电影天堂91| 久久久欧美国产精品| 嫩草影院新地址| 国产精品一区二区在线观看99| 午夜福利视频1000在线观看| 亚洲一区二区三区欧美精品 | 丰满人妻一区二区三区视频av| 夜夜爽夜夜爽视频| 听说在线观看完整版免费高清| 亚洲国产日韩一区二区| 日本一二三区视频观看| freevideosex欧美| 99热这里只有精品一区| 在线观看三级黄色| 国产伦在线观看视频一区| 69av精品久久久久久| 看十八女毛片水多多多| 91在线精品国自产拍蜜月| 欧美日本视频| 舔av片在线| 亚洲无线观看免费| 精品久久久久久久末码| 日韩欧美 国产精品| 成年免费大片在线观看| 亚洲伊人久久精品综合| 秋霞伦理黄片| 一区二区三区四区激情视频| 全区人妻精品视频| 久久精品国产亚洲网站| 久久久久精品久久久久真实原创| 黄色欧美视频在线观看| 青青草视频在线视频观看| 黄色日韩在线| 国产精品99久久99久久久不卡 | 亚洲人成网站高清观看| 久久影院123| 欧美激情久久久久久爽电影| 日本欧美国产在线视频| 国产免费一级a男人的天堂| 少妇人妻精品综合一区二区| 久久久久久久久大av| 99热这里只有是精品在线观看| 亚洲人成网站在线播| 久久女婷五月综合色啪小说 | 成人特级av手机在线观看| 青春草视频在线免费观看| 99久久中文字幕三级久久日本| 成人二区视频| 亚洲精品,欧美精品| 秋霞在线观看毛片| 精品人妻偷拍中文字幕| 成年人午夜在线观看视频| 天堂中文最新版在线下载 | 亚洲人成网站高清观看| 七月丁香在线播放| 国产综合懂色| 精品久久久久久电影网| 成年av动漫网址| 成人国产麻豆网| 中文字幕免费在线视频6| 欧美日韩精品成人综合77777| 免费看日本二区| 一区二区三区精品91| 韩国高清视频一区二区三区| 久久久午夜欧美精品| 亚洲精品第二区| 精品人妻视频免费看| 观看免费一级毛片| 偷拍熟女少妇极品色| .国产精品久久| 国产成人aa在线观看| 精品国产乱码久久久久久小说| 精品久久久噜噜| 色视频在线一区二区三区| 插逼视频在线观看| 天堂俺去俺来也www色官网| 99热网站在线观看| 久久亚洲国产成人精品v| 国模一区二区三区四区视频| 国产黄色免费在线视频| av黄色大香蕉| 欧美性猛交╳xxx乱大交人| 欧美日韩综合久久久久久| 亚洲欧美日韩另类电影网站 | 国产一区二区亚洲精品在线观看| 女人久久www免费人成看片| av女优亚洲男人天堂| 3wmmmm亚洲av在线观看| 免费人成在线观看视频色| 免费观看性生交大片5| 在线看a的网站| 极品教师在线视频| 精品一区二区免费观看| 中文在线观看免费www的网站| 欧美亚洲 丝袜 人妻 在线| 看非洲黑人一级黄片| 美女高潮的动态| 亚洲av中文av极速乱| 在线 av 中文字幕| 永久免费av网站大全| 亚洲最大成人手机在线| 成年女人看的毛片在线观看| 麻豆成人午夜福利视频| 男女边摸边吃奶| 国产精品一二三区在线看| 91aial.com中文字幕在线观看| 成人鲁丝片一二三区免费| 黄片wwwwww| 亚洲av在线观看美女高潮| 我的女老师完整版在线观看| 国产成人精品一,二区| 观看美女的网站| 亚洲av日韩在线播放| 如何舔出高潮| 中文精品一卡2卡3卡4更新| 尤物成人国产欧美一区二区三区| 亚洲av电影在线观看一区二区三区 | 久久精品国产自在天天线| 免费人成在线观看视频色| 国产成人免费无遮挡视频| 91午夜精品亚洲一区二区三区| 日本三级黄在线观看| 国产高清国产精品国产三级 | 我的女老师完整版在线观看| 亚洲,欧美,日韩| 中文乱码字字幕精品一区二区三区| 日韩伦理黄色片| 亚洲人成网站在线观看播放| 亚洲va在线va天堂va国产| 亚洲自拍偷在线| 精品一区二区三卡| 午夜福利视频精品| 日日啪夜夜撸| 午夜福利在线观看免费完整高清在| 国产高清国产精品国产三级 | 国内揄拍国产精品人妻在线| 尤物成人国产欧美一区二区三区| 亚洲欧美成人综合另类久久久| 综合色丁香网| 白带黄色成豆腐渣| 中文乱码字字幕精品一区二区三区| 2021天堂中文幕一二区在线观| 免费黄网站久久成人精品| 爱豆传媒免费全集在线观看| 日韩一区二区视频免费看| 国产视频首页在线观看| 欧美激情久久久久久爽电影| 男插女下体视频免费在线播放| 91午夜精品亚洲一区二区三区| 嘟嘟电影网在线观看| 国产精品久久久久久精品电影| 午夜免费观看性视频| 久久99热这里只频精品6学生| 夜夜看夜夜爽夜夜摸| 麻豆乱淫一区二区| 中文欧美无线码| 亚洲图色成人| 国产精品三级大全| 99久久九九国产精品国产免费| 99久久精品国产国产毛片| 亚洲精品自拍成人| 中文字幕免费在线视频6| 99久久中文字幕三级久久日本| 亚洲三级黄色毛片| 国国产精品蜜臀av免费| 国产成人freesex在线| 久久久久九九精品影院| 久久久亚洲精品成人影院| 国产精品久久久久久精品古装| 免费播放大片免费观看视频在线观看| 国产大屁股一区二区在线视频| 亚洲欧洲日产国产| 伊人久久精品亚洲午夜| 欧美成人一区二区免费高清观看| 亚洲精品aⅴ在线观看| 成人鲁丝片一二三区免费| 久久99热这里只有精品18| 成人一区二区视频在线观看|