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

兩相環(huán)境下EPDM絕熱層多因素耦合燒蝕預估①

2021-11-24 07:16:36馮喜平陳嘉輝盧景旭
固體火箭技術 2021年5期
關鍵詞:實驗模型

王 樂,馮喜平,陳嘉輝,盧景旭,梁 群,侯 曉

(西北工業(yè)大學 航天學院 燃燒、熱結構與內流場重點實驗室,西安 710072)

符號表

j組分擴散率dc燃燒室直徑

ρ密度D組分擴散系數(shù)

u、v流速A指前因子

h對流換熱系數(shù)E活化能

hD傳質系數(shù)R氣體常數(shù)

cp燃氣定壓比容f總氣化分數(shù)

p壓強y絕熱層厚度

M分子量n質量通量

Mgw燃氣平均分子量α顆粒沖刷角度

T溫度ε孔隙率

Ck當量熱容

hs、hw、hr燒蝕壁焓值、氣動壁焓值、恢復焓值

下標

c 碳化層 d 擴散控制機制

vir 基體層 e 附面層外緣

pyr 熱解層 g 氣體

ch 動力學控制機制 w 壁面

i組分 p 顆粒

0 引言

固體火箭發(fā)動機工作過程中,為確保燃燒室結構的完整性,常采用EPDM或改性EPDM絕熱層來保護燃燒室[1-3]。而在燃燒室的高溫、高壓、兩相流動燃氣環(huán)境中,絕熱層的設計與其燒蝕預估密切相關。固體發(fā)動機燃燒室內的流動多為兩相流動,在此環(huán)境下絕熱層的燒蝕由熱化學反應造成的熱化學燒蝕、燃氣流動造成的氣流剝蝕和離散相顆粒撞擊造成的顆粒侵蝕共同組成。

目前,國內外針對EPDM絕熱層不同因素耦合情況下的燒蝕性能進行了一些理論與實驗研究。KOO等[4]研究了有無顆粒侵蝕情況下的燒蝕性能,實驗結果表明,顆粒侵蝕條件下的絕熱層燒蝕率是無顆粒侵蝕情況下的2~12倍,說明燒蝕研究中不能忽略顆粒侵蝕;YANG等[5]數(shù)值模擬了酚醛樹脂等炭化類絕熱材料在熱化學燒蝕及顆粒侵蝕耦合作用下的燒蝕性能,并與實驗值進行了比較。研究表明,熱化學燒蝕會加劇顆粒侵蝕,反過來顆粒的作用也加劇了熱化學燒蝕;DIMITRIENKO[6]將燒蝕分為體積燒蝕和表面燒蝕,其中體燒蝕主要指熱解,將會使材料密度下降;面燒蝕指燒蝕表面的質量損失,將會使結構厚度變化。CURRY等[7]建立了更為復雜的雙區(qū)體燒蝕模型,分為多孔介質區(qū)和固體區(qū),在多孔介質區(qū)發(fā)生了復雜的物理化學過程,其中包含沉積反應。CAUTY等[8]將超聲技術應用于絕熱材料的數(shù)據(jù)收集和估計火箭發(fā)動機表面的熱載荷,并證明了該方法對于熱流估計的準確性。

國內對EPDM絕熱層的燒蝕問題同樣也進行了深入研究。何洪慶等[9]對炭化型熱防護材料EPDM進行了研究,采用氣動熱化學燒蝕機理,擴散和化學動力學雙控制機制,其中化學表面反應加入了SiO2的反應模型,并計入氣流與粒子的侵蝕,其中粒子侵蝕公式采用經(jīng)驗公式,建立EPDM的燒蝕模型。李江等[10]研究了顆粒濃度和沖刷速率對絕熱材料燒蝕的影響。在實驗的基礎上,建立了沖刷狀態(tài)與炭層厚度的關系,并采用傳統(tǒng)的炭化燒蝕模型進行燒蝕計算。王娟等[11]在顆粒沖刷條件下進行了EPDM燒蝕實驗,分析結果表明,認為某一種顆粒速度為臨界速度,當臨界速度高于沖刷速度時,炭化燒蝕率受粒子濃度、速度和角度影響較小,當臨界速度低于沖刷速度時,炭化燒蝕率隨速度、角度的增加而急劇增大。王書賢等[12]分析了氣相環(huán)境下不同燃氣速度對EPDM絕熱層炭化燒蝕率、炭化層結構特征以及成分分布等的影響。研究發(fā)現(xiàn),燃氣速度較低時,燒蝕以熱解炭化為主;燃氣速度較高時,炭化層的機械破壞顯著。對于EPDM絕熱層燒蝕預估,綜合考慮三因素的研究相對較少,李強等[13]耦合前述三種因素,采用熱化學反應在碳化層內部進行的假設,以碳化層的孔隙率為紐帶,建立了熱化學燒蝕、顆粒侵蝕和氣流剝蝕的耦合燒蝕模型,初步確定熱解氣體沉積是炭化層中致密結構的形成的主要原因,并根據(jù)實驗結果及沉積理論擬合致密層數(shù)值算法。研究表明,建立的模型在燒蝕率預估方面具有較好的精度。劉洋等[4]研究了熱解氣和燃氣在碳化層多孔結構中的化學氣相沉積和熱化學燒蝕過程。以及氣體和顆粒在碳化層上的侵蝕過程,數(shù)值計算結果表明,實驗結果與數(shù)值計算結果的最大誤差小于12%,數(shù)值結果與實驗結果偏差還有待降低。

本文就EPDM在固體發(fā)動機工作環(huán)境下的流固耦合燒蝕工程預估展開研究。基于絕熱層的消耗反應,即炭化層與推進劑燃氣之間的反應僅發(fā)生在炭化層表面的假設,通過炭化層表面孔隙率實現(xiàn)熱化學燒蝕、顆粒侵蝕和氣流剝蝕等的燒蝕耦合,采用動邊界熱傳導差分法進行數(shù)學模型的求解,實現(xiàn)了固體火箭發(fā)動機兩相環(huán)境下EPDM絕熱層炭化燒蝕率的工程預估。

1 物理模型

圖1為燒蝕后絕熱層材料內部的組成和兩相環(huán)境下絕熱層表面外和絕熱層內部發(fā)生的物理現(xiàn)象。

圖1 EPDM絕熱層多因素燒蝕示意圖

發(fā)動機工作過程為:發(fā)動機工作后(點火并燃燒后),首先,推進劑燃燒產(chǎn)生兩相高溫高壓燃氣;然后,高溫高壓燃氣在貼附有絕熱層的燃燒室內流動;最后,燃氣由噴管噴出產(chǎn)生推力。在這一流動過程中,由于高溫高壓燃氣沿貼附于燃燒內表面的絕熱層外表面流動,兩相燃氣與絕熱層外表面之間發(fā)生對流和輻射換熱,這樣使絕熱層表面溫度升高。同時,由于絕熱層內部溫度梯度的存在,絕熱層內部發(fā)生熱傳導,將表面導入的熱量向絕熱層內部傳遞,使得絕熱層內部溫度升高。由于絕熱層為高分子復合材料,當溫度升高時,其內部發(fā)生化學變化,材料物理性能也會發(fā)生變化。

在這一過程中,對流換熱和熱傳導使得絕熱層表面溫度升高,當絕熱層表面溫度達到絕熱層材料的熱解溫度時,表面材料發(fā)生熱解,隨著時間的進行熱界面加厚形成熱解層;隨著熱交換過程的進一步進行,材料表面溫度達到炭化溫度時,材料表面發(fā)生炭化,這一過程中伴隨材料中極大部分高聚物分解成低分子氣體溢出表面,留下多孔焦炭殘渣,并和原先加入基體的無機添加劑一起形成炭化層。

發(fā)動機工作過程中,由于炭化層的多孔性質和氣體擴散性,來流燃氣中氧化性組分和熱解氣體中最終裂解產(chǎn)生的氧化性氣體同炭化層中的碳發(fā)生化學反應,吸收熱量的同時使得炭化層表面不斷退移。此外,高速燃氣對絕熱層表面產(chǎn)生沖刷以及燃氣中Al2O3顆粒會撞擊絕熱層表面,引起絕熱層表面附近的炭化層表面的進一步退移。

為進行燒蝕機理研究,設計了EPDM絕熱層三因素耦合燒蝕的實驗器,實驗原理如圖2所示。

圖2 EPDM絕熱層實驗原理圖

由圖2可見,4片EPDM絕熱層燒蝕試件分別安裝于燒蝕實驗發(fā)動機方形試驗段的4個側面,每片絕熱層試件的y方向厚度為0.016 m,x方向長度為0.18 m。燃燒室高溫兩相燃氣通過收斂段加速進入實驗段,并流過EPDM絕熱層試件表面對試件進行燒蝕的同時,燃氣經(jīng)過試驗段后繼續(xù)加速,最后通過噴管流出。針對該耦合燒蝕模型,粒子撞擊角度為0,因此該模型主要受到氣流剝蝕影響。

根據(jù)設計實驗工況,確定計算工況:發(fā)動機工作壓強7 MPa,工作時間7 s,噴管喉徑10.8 mm,EPDM絕熱層厚度0.016 m,長度0.18 m。實驗采用的推進劑為高能復合推進劑,燃溫3536 K,其中鋁粉含量16%。實驗可得絕熱材料燒蝕后平均體孔隙率為0.562。

2 數(shù)學模型

2.1 絕熱層表面熱化學燒蝕模型

基本假設:

(1)絕熱層表面燃氣流動為二維凍結流;

(2)化學反應僅發(fā)生于絕熱層表面;

(3)氣體流動符合理想狀態(tài)方程;

(4)絕熱材料不發(fā)生膨脹或收縮;

(5)表面遵循質量、動量和能量守恒定律。

絕熱層表面的質量、動量和能量傳遞關系如圖3所示,根據(jù)該傳遞關系建立守恒方程。

圖3 絕熱層表面質量、動量和能量守恒關系

(1)質量守恒方程

熱解氣體流入控制體的質量流率+化學反應產(chǎn)生氣體質量流率=絕熱層表面流出氣體的質量流率,即:

(1)

(2)能量守恒方程

熱解氣體流入熱流+熱化學反應熱效應+凈熱輻射熱流+顆粒帶來的熱流+對流換熱流入熱流=導入絕熱層內部熱流+流出氣體帶走的焓+分子運動擴散帶走熱流+燃氣與絕熱層熱傳導熱流,即

(2)

(3)組分守恒方程

絕熱層表面流出第i種氣體的組分質量率=濃度梯度擴散的第i種組分質量率+化學反應產(chǎn)生的第i種組分質量率+熱解氣體產(chǎn)生的第i種組分質量率,即:

(3)

絕熱層表面的三個化學反應為

C(s)+CO2=2CO

C(s)+H2O=CO+H2

2C(s)+H2=C2H2

將不參加壁面化學反應的組分都折合歸納入N2,將組分按次序排列為H2O、CO2、H2、CO、C2H2和N2,并依次對應表示為i=1,2,…,6。化學表面反應方程組缺少熱解氣體組分數(shù)據(jù),取熱解氣體組分質量分數(shù)f1=0.35,f2=0.32,f3=0.18。

絕熱層的熱化學燒蝕,一方面受化學反應動力學控制;另一方面受燃氣組分通過附面層向壁面的擴散速率控制,即上述三個反應為雙控制機制[15]。雙控制下氧化性組分的質量消耗率為

(i=1,2,3)

(4)

式中Ai為指前因子,A1和A2均取0.011 357 kg/(m2·s·Pa),A3=0.000 711 22 kg/(m2·s·Pa);Ei為反應活化能,E1和E2均取236 500 J/mol,E3=284 600 J/mol;pw,i為i組元在壁面處的分壓;R為通用氣體常數(shù);Tw為氣動壁溫度;pw,i為氣動壁面不同組分分壓;ρe為燃氣密度;ue為燃氣速度;CH為傳熱傳質系數(shù);Kie為燃氣中不同組分濃度[16]。

表面化學反應中組分的質量生成率為

(5)

其中,M為相對分子質量。將式(5)代入到絕熱層表面控制方程式(1)和式(3),得絕熱層表面單位面積的碳消耗率為

(6)

則絕熱層線燒蝕率為

(7)

因為實驗后絕熱層燒蝕試件的炭化層與基體層之間存在明顯間隙,同時炭化層會出現(xiàn)膨脹等現(xiàn)象,如采用線燒蝕率來表征絕熱層的燒蝕性能將存在較大誤差。因此,本文采用炭化燒蝕率rc表征絕熱層的燒蝕性能,其定義為:炭化燒蝕率=(試件初始厚度-實驗后試件基體層厚度)/實驗發(fā)動機工作時間。

(8)

根據(jù)熱化學燒蝕實驗結果[18],炭化層內部孔隙率逐漸增大,與炭化層厚度δ成線性關系,t+Δt時刻炭化層內部各點孔隙率采用經(jīng)驗關系式:

(9)

2.2 顆粒侵蝕計算模型

根據(jù)多孔材料強度的概念,孔隙率大小與炭化層強度存在一一對應的關系。由于顆粒侵蝕強度與燃氣流速、顆粒撞擊角度、炭化層強度等均有關系,作用機制復雜。因此,顆粒侵蝕計算模型也常采用基于實驗的經(jīng)驗關系式。

根據(jù)現(xiàn)有實驗結果,顆粒侵蝕炭化層燒蝕表面的臨界孔隙率εc1可表示為[19]

(1)臨界速度以下(v<32 m/s)

(10)

(2)臨界速度以上(v≥32 m/s)

(11)

式中np為顆粒沖刷質量通量;vp為顆粒撞擊速度;α為顆粒撞擊角度。

2.3 氣流剝蝕模型

影響絕熱層氣流剝蝕的主要因素有氣流速度、環(huán)境壓強、炭化層高溫結構強度、熱解氣體逸出產(chǎn)生的粘性力和內壓等,這些影響氣流剝蝕因素的詳細作用機制目前還不清楚。因此,氣流剝蝕計算模型只能采用經(jīng)驗關系式,即以實驗測得的表面孔隙率為氣流剝蝕后炭化層表面的臨界孔隙率。

在認為氣流主要通過剪切力實現(xiàn)剝蝕的前提下,根據(jù)實驗結果擬合得到炭化層臨界孔隙率和炭化層表面剪切力的關系式為[10]

εc2=0.2+0.07×0.896τ

(12)

其中,εc2為氣流剝蝕對炭化層表面的臨界孔隙率。且剪切力與燃氣參數(shù)的實驗關系式如下:

(13)

2.4 絕熱材料傳熱分析數(shù)學模型

由物理模型可知,EPDM絕熱材料燒蝕后可分為三層,針對不同層區(qū)有不同熱傳導方程:

(1)炭化層

(14)

(2)熱解層

(15)

式中 ΔHp為材料的熱解吸熱率。

熱解層的密度、熱導率和比熱可近似認為是隨溫度線性變化的[20]。即:

(3)基體層:

(16)

絕熱層表面溫度不同時,傳入內部熱流量不同,定義絕熱材料開始熱解溫度為Tf=600 K,絕熱材料開始炭化溫度為Tc=1200 K。

(4)初始條件

t=0時,T0=Ts。

(5)邊界條件

1)T

(17)

2)Tf

(18)

3)T>Tc時,絕熱材料表面開始炭化,絕熱材料開始形成三層結構,此時:

(19)

在絕熱材料背面假設為絕熱,則:

(20)

2.4.1 壁面溫度求解

壁面溫度分為氣動壁溫度和燒蝕壁溫度,這是由于化學反應不是在面內進行而是在控制體內進行,控制體上下表面溫度不一致。可用氣膜分析法找到氣動壁溫度、燒蝕壁溫度和氣相附面層外燃氣溫度的關系如下:

(21)

式中a為氣膜厚度,實驗結束后通過電鏡掃描測量獲得。

實際上,當表面沒有炭化時,氣動壁溫度等于燒蝕壁溫度,但燃燒室溫度逐漸升高,開始燒蝕時環(huán)境溫度沒有達到推進劑絕熱燃燒溫度,若仍采取絕熱燃燒溫度計算會造成輻射熱流量和對流換熱量異常大,不符合實際過程。因此,采取上述方法對壁面溫度進行初始計算,以符合實際過程。

2.4.2 輻射熱流計算

可用黑體輻射定律計算輻射熱流量,即

(22)

2.4.3 附面層熱傳導求解

燃氣在絕熱層表面流動可為層流流動,也可為湍流流動。認為是層流流動時可用非定溫繞流平壁的溫度型方程來求解燃氣與絕熱層之間的熱傳導熱流。

本次計算中采用簡化算法,使用對流換熱關系式進行求解,表達式如下:

(23)

2.4.4 焓值計算

焓值計算包括氣動壁處焓值、燒蝕壁處焓值和恢復焓值的計算。

氣動壁處焓值:

Hw=ceTw

(24)

燒蝕壁處焓值:

Hs=ceTs

(25)

恢復焓值:

(26)

式中k為燃氣比熱比;Me為燃氣馬赫數(shù)。

2.4.5 化學反應熱效應計算

化學反應熱效應的計算思路:反應熱×每種反應碳的消耗量。則計算式為

(27)

該反應下,C與H2O、CO2、H2的化學發(fā)應熱分別為-4.688、-7.621、-7.428 kJ/mol。

2.4.6 顆粒沖刷熱增量模型

顆粒沖刷壁面的熱增量計算關系式為

(28)

式中Np為顆粒濃度。

2.5 多因素耦合計算模型

多因素耦合計算模型是以炭化層燒蝕表面的表面孔隙率為耦合參數(shù),對熱化學燒蝕、顆粒侵蝕和氣流剝蝕等各燒蝕模型進行耦合,建立了多因素耦合計算模型。絕熱層多因素耦合計算模型如圖4所示。

圖4 絕熱層多因素耦合計算模型示意圖

絕熱層非穩(wěn)態(tài)多因素耦合燒蝕性能計算步驟:

(3)如此循環(huán)計算,即可得到燒蝕表面各點的炭化線位置和炭化燒蝕率。

3 結果分析與討論

3.1 兩相燒蝕氣流內流場計算結果與分析

為了便于將計算結果與實驗結果進行比較,以燒蝕實驗裝置的實驗段為物理模型,采用二維可壓縮流動假設,設置流動入口邊界條件為壓力入口,出口邊界條件為壓力出口,實驗器殼體為絕熱壁面,兩相流采用DPM模型,粒子含量取燃氣流量的15%。計算得到燒蝕實驗發(fā)動機燒蝕氣流內流場結果如圖5所示。

(a)Distribution of velocity

由圖5可知,燒蝕過程中試驗段(參見圖2,X=0.18~0.36 m)燃氣速度為32 m/s,燃氣壓力為7 MPa,燃氣溫度為3536 K,且氣流的速度、壓力、溫度沿X方向無明顯變化,粒子沖刷質量通量為87 kg/(m2·s),顆粒撞擊絕熱層表面速度為24 m/s。多因素耦合燒蝕性能計算時,選取絕熱層平均炭化燒蝕率作為性能表征,并將其與1-B處的炭化燒蝕率理論計算值rc進行比較,以驗證多因素耦合燒蝕模型。

3.2 絕熱層燒蝕性能計算結果與分析

在前一節(jié)中給出的兩相燒蝕氣流作用下,絕熱層溫度分布如圖6所示。絕熱試件1的1-B處(圖2)典型的燒蝕性能計算結果如圖7~圖10所示。

圖6 EPDM絕熱層溫度分布

圖7 燒蝕線、炭化線和熱解線位置隨燒蝕時間的變化圖

圖8 燒蝕表面溫度隨燒蝕時間的變化

圖9 燒蝕兩相流各熱流率隨燒蝕時間的變化

圖10 炭化燒蝕率隨時間的變化

由圖7和圖8可知,隨著燒蝕時間的增加,絕熱層表面溫度不斷升高、炭化層與熱解層厚度逐漸增加、炭化層內部孔隙率不斷增大,材料燒蝕6.4 s后,炭化層內部各點孔隙率的平均值≥EPDM絕熱層燒蝕實驗測得的炭化層平均孔隙率時而發(fā)生剝離,導致絕熱層表面發(fā)生燒蝕剝離,從而使炭化層厚度驟然減小,燒蝕表面溫度驟然降低,并隨著燒蝕的進行,溫度重新升高。熱解層和原始層厚度所受影響較小。燒蝕繼續(xù)進行到7 s時,燒蝕表面退移了0.69 mm。

由圖9可知,絕熱層表面處不同熱流率隨時間的變化,其中輻射換熱率的量級是106W/m2,對流換熱率與顆粒帶來的熱流率是105W/m2,其余熱流率較小,同時由計算結果可見,當燒蝕進行到6.4 s時,絕熱層剝離會導致表面溫度降低,進而會造成輻射換熱率以及對流換熱率的增加。

圖10為炭化燒蝕率隨時間的變化曲線,隨著燒蝕的進行,炭化燒蝕率逐漸達到穩(wěn)態(tài)為0.21 mm/s;當燒蝕進行到6.68 s時,絕熱層表面發(fā)生剝離,由于氣流剝蝕、顆粒侵蝕帶來的熱流率較輻射換熱率小一個量級,因此其對炭化燒蝕率影響較小,炭化燒蝕率保持穩(wěn)定。燒蝕繼續(xù)進行到7 s時,炭化燒蝕率為0.213 mm/s。

3.3 燒蝕性能數(shù)值解與實驗值對比分析

絕熱層燒蝕實驗的試件表面測點位置如圖11所示。

圖11 絕熱層表面測點位置圖(單位:m)

實驗后,測得絕熱層表面不同位置所對應的燒蝕厚度。之后,對相同X位置的值取平均燒蝕厚度,可得沿X軸不同位置的平均炭化燒蝕率,見圖12。

圖12 實驗后測得炭化燒蝕率

由圖12可得絕熱層平均炭化燒蝕率0.203 mm/s。與理論計算炭化燒蝕率0.22 mm/s比較,理論計算誤差小于10%,誤差主要來源于炭化層比熱容以及熱導率在高溫下變化的影響。由此驗證了本模型在熱化學、氣流剝蝕以及顆粒侵蝕下的燒蝕率預測能力。

本文所提的方法可用于不同發(fā)動機或同一發(fā)動機不同部位的兩相環(huán)境下EPDM絕熱層的燒蝕預估。

3.4 不同燒蝕影響因素下炭化燒蝕率分析

根據(jù)上述內容,對不同燃氣溫度、燃氣壓強、顆粒含量下的耦合燒蝕模型進行了分析計算,計算結果如表1所示。可見,燃氣溫度對絕熱層炭化燒蝕率的影響較大,因為絕熱層表面化學反應速率主要受燃氣溫度影響,燃溫越高,化學反應速率越快,加速燒蝕線、炭化線、熱解線退移,導致炭化燒蝕率變化幅度較大;燃氣壓強以及顆粒含量對炭化燒蝕率的影響相對較小,炭化燒蝕率隨壓強的升高以及顆粒含量的增大呈上升趨勢。

表1 不同燒蝕影響因素下的炭化燒蝕率

4 結論

(1)隨著燒蝕的進行,絕熱層表面溫度不斷升高、炭化層與熱解層厚度逐漸增加、炭化層內部孔隙率不斷增大;當炭化層內部各點孔隙率的平均值 EPDM絕熱層燒蝕實驗測得的炭化層平均孔隙率時發(fā)生剝離,導致絕熱層表面發(fā)生燒蝕剝離,從而使炭化層厚度驟然減小,燒蝕表面溫度驟然降低并隨著燒蝕的進行,溫度重新升高。熱解層和原始層厚度所受影響較小。

(2)絕熱層理論計算的炭化燒蝕率為0.22 mm/s,與實驗平均炭化燒蝕率0.203 mm/s的誤差為8.37%,且誤差主要來源于炭化層比熱容及熱導率在高溫下變化的影響。由此驗證了本模型在熱化學、氣流剝蝕及顆粒侵蝕下的燒蝕率預測能力。

(3)燃氣溫度相對于燃氣壓強以及顆粒含量對絕熱層炭化燒蝕率的影響較大。燃溫越高,化學反應速率越快,加速燒蝕線、炭化線、熱解線退移,炭化燒蝕率增大;炭化燒蝕率隨壓強的升高以及顆粒含量的增大呈上升趨勢。

猜你喜歡
實驗模型
一半模型
記一次有趣的實驗
微型實驗里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 午夜精品久久久久久久99热下载 | 国产午夜人做人免费视频中文| 91区国产福利在线观看午夜| 九九热免费在线视频| 最新亚洲av女人的天堂| 中文国产成人精品久久一| 国产亚洲视频播放9000| 久久婷婷六月| 欧美精品成人| 久久99国产视频| 1级黄色毛片| 精品国产自| 国产成人精品亚洲77美色| 国产精品成人一区二区不卡| 国产福利免费视频| 午夜a视频| 亚洲国产精品一区二区第一页免 | 久久久亚洲色| 午夜激情福利视频| 欧美午夜理伦三级在线观看| 亚洲日本www| 天堂成人在线| 欧美精品成人| 中文字幕一区二区视频| 国产一区二区三区精品久久呦| 福利视频久久| 国产毛片高清一级国语| 99久久精品视香蕉蕉| 国产精欧美一区二区三区| 亚洲人成日本在线观看| 无码专区国产精品一区| 动漫精品中文字幕无码| AV不卡在线永久免费观看| 乱色熟女综合一区二区| 国产成人亚洲精品无码电影| 国产精品私拍在线爆乳| 福利片91| 亚洲精品无码在线播放网站| 91精品专区国产盗摄| 欧美中日韩在线| 亚洲一区二区视频在线观看| 久久午夜夜伦鲁鲁片无码免费| 国产第一页免费浮力影院| 91久久精品国产| 毛片在线看网站| 特级精品毛片免费观看| 亚洲精品欧美日本中文字幕| 亚洲欧洲日韩综合色天使| 久久www视频| 亚洲色图欧美视频| 国产成人精品男人的天堂下载| 欧美精品成人一区二区视频一| AV天堂资源福利在线观看| 青青草91视频| 激情无码字幕综合| 青青草91视频| 精品久久久无码专区中文字幕| 亚洲美女视频一区| 三上悠亚一区二区| 国产永久无码观看在线| 91蜜芽尤物福利在线观看| 亚洲天堂2014| 久久国产高清视频| 伊人五月丁香综合AⅤ| 中文字幕2区| 91啦中文字幕| 激情無極限的亚洲一区免费| 亚洲无线国产观看| 91色国产在线| 国产爽妇精品| 99偷拍视频精品一区二区| 丝袜久久剧情精品国产| 91年精品国产福利线观看久久| 99久久国产综合精品2020| 亚洲中文字幕97久久精品少妇| 成人综合久久综合| 欧美日韩国产在线播放| 2020国产免费久久精品99| 国产三级精品三级在线观看| 国产小视频免费| 91精品综合| 国产麻豆va精品视频|