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

固體氧化物燃料電池溫升模擬中入口異常高溫度梯度研究*

2022-08-28 09:27:34申雙林張小坤萬興文鄭克晴凌意瀚王紹榮
物理學報 2022年16期
關鍵詞:模型

申雙林 張小坤 萬興文 鄭克晴 凌意瀚 王紹榮

1)(中國礦業大學,低碳能源與動力工程學院,徐州 221116)

2)(中國礦業大學,材料與物理學院,徐州 221116)

3)(中國礦業大學,化工學院,徐州 221116)

針對固體氧化物燃料電池熱循環失效問題,建立了固體氧化物燃料電池熱氣體預熱動態模型,研究了電池內最大溫度梯度分布規律和入口異常高溫度梯度形成的原因,結果表明: 在熱氣體參數和預熱方式變化時,電池內最大溫度梯度始終處于電池入口邊緣處的電極表面;電池入口處存在異常高的溫度梯度,且在入口一小段區域內,溫度梯度沿流動方向迅速下降;其原因是模型中入口采用均一的平均速度和溫度,“入口效應”強化氣體與電池換熱;采用入口段延長的方式可使入口速度充分發展,降低電池內最大溫度梯度,但由于均一溫度入口并未優化,入口處仍然存在很大的溫度梯度和溫度梯度變化;因此采用數值模擬研究電池預熱升溫安全性時,僅采用最大溫度梯度作為安全性判據會高估電池內熱應力.

1 引言

固體氧化物燃料電池(solid oxide fuel cell,SOFC)是一種清潔高效的發電技術,在分布式供能、備用電源、移動供電、軍事等領域具有廣闊的應用前景[1,2].然而,當前SOFC 的商業化過程中仍然存在成本、壽命、熱循環穩定性、高溫密封等關鍵問題亟需解決[3-5],其中SOFC 在熱循環和啟停過程中熱沖擊引起性能衰減和失效是制約其壽命和運行安全性的重要問題[6-11].

SOFC 工作溫度通常在700 ℃以上,電池啟動時需要首先預熱到啟動溫度(例如: 600 ℃),采用燃燒器和熱氣體預熱SOFC 電堆是一種簡便可行的方法,但熱氣體加熱方法耗時長,氣體進、出口間溫度梯度大,影響電池安全快速啟動[12,13],針對這一問題許多學者做了大量研究工作.Selimovic等[14]通過數值方法研究了預熱過程中金屬和陶瓷連接板對電池內溫度梯度的影響,并指出金屬連接板電池內溫度梯度要小于陶瓷雙極板;Damm和Fedorov[15]建立了SOFC 動態升溫的降階模型,并采用模型研究了電池預熱過程中的傳熱機理和入口參數的影響規律;Colpan等[16]采用數值模型研究了直接內重整SOFC 預熱啟動過程中逆流與順流加熱對電池內熱應力的影響;Chen和Jiang[17]采用數值模型研究了熱氣體不同通入方式(單流道、雙流道順流和雙流道逆流)對電池預熱時間和溫度梯度分布的影響規律;Yuan和Liu[18]采用數值模型研究了熱氣體入口速度分配對電池內溫度梯度分布的影響規律,并指出陽極和陰極氣體入口速度分配對電池內溫度梯度分布和溫升時間有一定影響;Peksen 等[19-21]建立了SOFC 電堆動態模型,研究了不同預熱方式對電堆預熱時間和溫度分布的影響規律;Zheng等[22]采用數值模型,研究了雙極板幾何結構對電池預熱時間和溫度梯度的影響規律,研究發現雙極板吸熱量占電池總預熱熱量90%以上,因此雙極板結構對電池預熱時間和溫度梯度具有重要影響.

基于文獻綜述可知,當前針對SOFC 啟動預熱和熱循環方法的研究大多是采用數值模擬方法,分析電池在預熱過程中的溫度梯度分布,并采用最大溫度梯度作為電池失效判據,這是因為電池內熱應力與溫度梯度大小成正比[23-25].SOFC 的最大安全溫度梯度并沒有明確的結論,Dikwal等[26]采用實驗研究管式SOFC 熱循環過程中溫度梯度影響發現,當在溫度梯度約為3000 K/m 時電池性能下降明顯;Aguiar等[27]基于PEN(陽極、電解質、陰極)強度分析指出電池內最大溫度梯度應小于1000 K/m 時較為安全.然而數值模擬得到的最大溫度梯度大多大于這些安全溫度梯度,例如:Selimovic等[14]給出采用氫氣預熱時最大溫度梯度超過了5000 K/m;Chen和Jiang[17]研究結果表明無論順流還是逆流電池內最大溫度梯度均高于2000 K/m;Zheng等[22]計算得到的垂直于電解質方向最大溫度梯度可達2.5×105K/m.顯然,這些數值方法得到的最大溫度梯度與實際電池失效間存在一定差距,同時本課題組分析數值結果發現,電池預熱過程中在入口處會出現異常高的溫度梯度,因此,本文通過數值方法,研究SOFC 在預熱過程中最大溫度梯度產生位置和形成規律,并針對電池入口處產生的異常高溫度梯度形成機理進行一定分析,研究結果指出采用數值方法優化電池預熱過程時,僅僅采用電池內最大溫度梯度作為安全行判據會得到不合理的結果,而這一結論對指導電池預熱過程優化具有重要意義.

2 SOFC 預熱過程數值模型

2.1 模型假設和幾何結構

SOFC 的預熱過程是將電池從室溫預熱到啟動溫度(例如: 600 ℃),在預熱過程中由于電池溫度遠低于SOFC 工作溫度,電池不會放電,基本不會發生電化學反應,因此可假設電池內電化學反應生成電流為0,電池內只存在流動和傳熱過程.為了簡化模型,本文做出如下假設:

1)模型中SOFC 各組成部分材料屬性均為各向同性,如: 比熱和導熱系數;

2)電極(多孔陽極和陰極)和集流板脊之間的接觸熱阻忽略不計;

3)SOFC 加熱升溫過程中的熱輻射影響忽略不計;

4)由于不考慮電化學反應,加熱升溫過程中電池內部的加熱氣體為穩定狀態,不涉及傳質過程;

5)兩極流道內加熱氣體的流動狀態為層流(以空氣為例,雷諾數約為350).

對于平板式SOFC,電池結構具有重復性,為簡化模型計算量,本文選取電池內單個可重復流道單元作為研究對象[28,29],如圖1 所示.電池各部件的幾何尺寸見表1.

2.2 模型控制方程

根據前文假設可知,SOFC 預熱過程中電化學反應速率為0,電池內不存在電流傳輸和組分傳輸,因此SOFC 的預熱過程只存在氣體流動和傳熱過程,其控制方程分別為質量守恒方程、動量守恒方程和能量守恒方程,具體表達式如下[30]:

1)質量守恒方程(連續性方程)

式中,ρf為加熱氣體密度.

2)動量守恒方程

式中,p為氣體靜壓力;ε為多孔介質孔隙率,對于流道其值取1;為應力張量,可由下式計算:

μ為動力黏度;I為單位張量;SM為動量在多孔介質內的源項,在流道內其值為0,在多孔介質內可由下式計算:

式中,α為多孔介質滲透系數;C2為慣性阻力系數.

3)能量守恒方程

電池預熱過程中輻射傳熱忽略不計,只存在熱傳導和熱對流兩種熱量傳遞方式,因此能量守恒方程可表示為

式中,ρs為多孔介質中固體骨架的密度;E為內能;hi為組分焓值;keff為多孔介質的有效導熱系數,可由下式計算:

式中,kf和ks分別表示加熱氣體和固體域的導熱系數.

2.3 模型求解與驗證

SOFC 預熱過程為動態過程,因此,模型求解需要邊界條件和初始條件.當SOFC 采用熱氣體雙流道逆流預熱時,陽極通入氮氣、陰極通入空氣,陽極與陰極流道入口分別給定氣體入口速度和溫度.基礎情況下,氣體入口速度為5 m/s,入口氣體初始溫度為298 K(25 ℃),氣體溫度以1 K/s 速率增加直到電池啟動溫度(898 K,即625 ℃)后維持不變,通常實驗電池的溫升速率較低(例如:3 K/min),本文選用如此大的溫升速率,是在計算時間和結果適用性間折衷考慮的結果.文獻[15,22]和本文后續研究結果指出: 大的溫升梯度會增大最大溫度梯度的絕對值,但不會改變電池內傳熱規律和溫升規律,而本文的目的主要是研究入口處異常高溫度梯度的產生機理,選用較大的溫升速率不會影響本文研究結論的適用性,且可極大地減少數值模型的計算時間,因此,在研究中主要選用1 K/s溫升速率的結果.氣體流道出口為自然出流邊界條件,即背壓為0 Pa(相對壓力).假設電池從室溫啟動,即電池初始溫度為298 K,流道對稱面為對稱邊界條件,其余表面為絕熱邊界條件.

本文模型采用流體力學計算軟件Fluent?求解.數值模型中的物性參數見表2 所列.圖1 所示求解單元采用結構化網格,選擇不同疏密的網格進行600 s 計算并取電解質內最低溫度結果如圖2(a)所示,當網格單元數為1547792 時,計算結果基本不變.在動態模型中,時間步長對求解穩定性和結果的可靠性至關重要.本文氣體入口速度較大,若依據動量在網格中的傳播來選擇時間步長的話,時間步長會非常小,整個升溫過程的計算時間會非常長.考慮到在整個升溫過程中,氣體入口速度保持不變,電池內流場變化必然比較小,預熱過程中主要是溫度的變化與傳播,因此依據電池內溫度計算的穩定性來選擇時間步長.如圖2 為不同時間步長時,預熱600 s 后電解質內最低溫度的變化曲線,由圖可知當時間步長小于0.1 s 時,電解質最低溫度計算結果基本穩定,因此,本文選擇的時間步長為0.1 s.

圖2 (a)SOFC 數值模型網格無關性和(b)時間步長的驗證Fig.2.(a)Verification of mesh independence and (b)time step size for SOFC numerical model.

表2 SOFC 各部分結構物性參數[15]Table 2.The physical parameters for each component of SOFC.

本文數值模型的驗證通過與文獻[15]中的數值結果進行對比,比較結果如圖3 所示[22].由圖3可知,本文數值模型與文獻中模型結果吻合很好,因此,本文所建立的數值模型可用于SOFC 預熱過程的研究.

圖3 本文預熱溫升數值模型結果與文獻[15]中結果對比圖Fig.3.Comparison of the current preheating model results with model results given in Ref.[15].

3 數值模擬結果與分析

3.1 預熱過程中,最大溫度梯度位置分析

在SOFC 預熱過程研究中,通常只關心最大溫度梯度值的變化,而缺乏對其分布位置的研究.本文為了研究最大溫度梯度形成原因,首先分析不同情況下SOFC 預熱過程中最大溫度梯度的分布位置.且由前期研究結果可知,電池中最大溫度梯度主要集中在垂直于電解質方向(Y方向),因此,本文后面均取Y方向最大溫度梯度進行研究[22].

圖4 為采用雙流道逆流預熱時,熱氣體不同升溫速率(0.2,0.5,1,2,3,4和5 K/s)和不同入口速度(0.1,0.5,1,2,3,4,5,6,7和8 m/s)時,最大溫度梯度在電池內的位置.圖4 中網格為沿流動方向(z方向)陰極入口處第一層網格切面,圖4(a)為1 K/s 升溫速率、5 m/s 入口速度下,預熱600 s 時陰極入口橫截面(垂直于流動方向)處的溫度分布.在研究中,當氣體入口升溫速率變化時,氣體入口速度保持5 m/s 不變;而當氣體入口速度變化時,氣體升溫速率保持1 K/s 不變.由圖4(a)可知,當熱氣體升溫速率變化時,電池預熱過程中最大溫度梯度的位置始終在A 點不變,即在陰極入口處陰極表面第一層網格靠近流道、脊和陰極接觸的角落中.對于雙流道逆流預熱,陰極入口處熱量由陰極流道向陽極流道傳輸,由表2 可知,陰極導熱系數最低,陰極Y方向溫度梯度最大,而陰極表面熱流密度最大,因此,最大溫度梯度在陰極表面第一層網格中,同時對于熱氣體預熱方式,氣體入口處溫差最大,最大溫度梯度在入口第一層網格內.

圖4 加熱氣體不同升溫速率(a)和不同入口速度(b)時SOFC 內最大溫度梯度位置示意圖((a)1 K/s 升溫600 s 后陰極入口橫截面處的溫度分布,圖中溫度單位為K)Fig.4.The location of maximum temperature gradient in SOFC under different temperature rise rate (a)and inlet velocity (b)(Figure (a)is the temperature distribution of the cathode inlet section at 600 s with a 1 K/s increasing rate,temperature unit is K).

由圖4(b)可知,當氣體入口速度大于等于4 m/s 時,最大溫度梯度始終在A 處不變;而當入口速度小于4 m/s 時,最大溫度梯度會隨升溫過程,由A 點到C 點、D 點再到E 點變化,當入口速度小于等于1 m/s 時,最大溫度梯度主要在E 點附近.這主要是由于氣流速度會改變氣體與陰極和脊間的換熱系數,從而改變熱量由熱氣體到脊與陰極的傳輸量.

此外,熱氣體的通入方式對溫度梯度的形成有重要影響[17],本文也研究了分別采用雙流道逆流、雙流道順流、陰極單流道和陽極單流道幾種氣體通入方式對最大溫度梯度位置的影響.研究結果發現當采用雙流道逆流、雙流道順流和陰極單流道氣體通入方式時,最大溫度梯度始終在圖4(a)所示的A 點位置,也即陰極入口處陰極第一層網格內;而當采用陽極單流道預熱時,最大溫度梯度在陽極入口處陽極表面第一層網格內,即陽極內與A 點對應位置.

綜合上述研究結果可知,最大溫度梯度均出現在氣流入口處,這是由于入口處氣體與電池間溫差最大.同時最大溫度梯度均在PEN(陽極、電極質、陰極)表面第一層網格內,對于采用陽極單流道預熱,最大溫度梯度在陽極表面;對于雙流道逆流、雙流道順流和陰極單流道預熱,最大溫度梯度在陰極表面,這是因為PEN 表面熱流密度最大,同時陰極導熱系數最小.

3.2 SOFC 內溫度梯度分布規律

前述結果可知最大溫度梯度在入口處,本文進一步研究了電池內部溫度梯度的分布規律.由前期研究結果可知當氣體入口溫度達到啟動溫度時(即氣體溫度達到898 K 時),電池內溫度梯度最大,因此,后續結果為熱氣體溫度達到898 K時Y方向上的溫度梯度分布結果[22].

圖5 為熱氣體不同升溫速率和不同入口速度時,SOFC內Y方向上溫度梯度沿流道方向分布曲線.這些曲線均是在X=0 mm和Y=1.975 mm這一線上繪制,線的起點即為前文所述的最大溫度梯度點A 點.由圖5 中結果可知,電池內溫度梯度由氣體入口處沿流動方向先急劇下降、后緩慢降低,在超過流道中心處以后再緩慢升高,而到靠近出口處時再迅速上升;在氣體入口段,電池溫度梯度非常高,且隨著溫升速率和進氣速率的降低,入口處溫度梯度曲線越來越陡,也即高溫度梯度段越來越短.文獻研究結果可知,電池內溫度梯度隨入口氣體溫升速率降低和入口速率增大而降低,圖5(a)中溫度梯度的變化表現出相同的規律,不過即使溫升速率低至0.05 K/s(3 K/min)時,入口處仍然有非常高的溫度梯度,顯然如此高的溫度梯度與電池實際工況不符,如果數值模擬僅僅采用最大溫度梯度作為失效判據,會得到不合理的結果;圖5(b)的結果似乎與文獻結果不太一致,圖中結果似乎是隨著入口速率增大溫度梯度變大了,然而仔細比較最大溫度梯度結果可知本文結果與文獻結果一致,也即隨著入口速度增大最大溫度梯度會逐漸降低,而隨著入口速度降低,溫度梯度在入口處越來越陡,高溫度梯度區域主要集中在入口處,這主要是由于氣體質量流量降低,氣體攜帶進入電池內的總熱值較低,熱量會在入口處迅速被電池吸收,顯然在氣體入口速率的優化中僅僅考慮最大溫度梯度是不合適的.圖6 為采用不同預熱氣體通入方式時,SOFC 內部溫度梯度沿流道方向分布曲線.由圖6 可以發現,當采用陽極和陰極單流道預熱時,異常高溫度梯度分別出現在陽極和陰極入口段;雙流道順流預熱方式,異常高溫度梯度同樣出現在入口區域;而雙流道逆流方式由于兩端均有進氣,因此異常高溫度梯度在流道兩端均存在.

圖5 不同氣體入口升溫速率(a)和入口速度(b)時,SOFC 內溫度梯度沿流道方向變化曲線Fig.5.The plotting of temperature gradient in the SOFC on Z direction under different gases temperature rise rate (a)and different inlet velocity (b).

圖6 不同預熱方式時,SOFC 內溫度梯度沿流道方向變化曲線Fig.6.The plotting of temperature gradient in the SOFC on z direction under different preheating method.

由以上數值模擬結果可知: 在流道氣體入口區域存在非常大的溫度梯度,且這一溫度梯度沿流動方向急劇下降;這一溫度梯度明顯與實際電池預熱過程不符,若僅僅采用最大溫度梯度作為數值模擬研究和優化判據,必然會得到不合理的結果,而采用平均溫度梯度作為判據也會一定程度上高估電池內的熱應力;同時無論如何改變氣體預熱參數和預熱方式,這一異常高溫度梯度區域始終存在,因此這一結果并非模型計算錯誤,其形成原因必然與模型方程和邊界條件設置相關,需要進一步深入研究.

3.3 “入口效應”對最大溫度梯度的影響

由流體力學和傳熱學可知,在管道入口處由于邊界層的發展會存在“入口效應”強化流體與管道壁面的換熱,Gamrat 等[31]針對氣體在方形微管道中的流動與換熱研究中也表明: 由于入口效應的存在,管道入口氣體表面換熱系數會大于充分發展段.SOFC 預熱模型中,熱氣體在方形流道內為層流流動,與電池本體進行對流換熱,模擬結果中入口異常高溫度梯度有可能是由于入口效應引起,因此,采用將流道入口延長21 mm 的方式,使流道入口氣流為充分發展,模擬結果如圖7—圖9 所示.

圖7 600 s 時,SOFC 陰極入口處流道對稱面上的速度分布圖Fig.7.The distribution of velocity in symmetrical surface of SOFC cathode channel near inlet at 600 s.

圖7 為原模型和入口延長后數值計算得到的陰極入口處速度分布圖,由圖可以看出: 入口延長后,在入口延長段,氣體流動達到充分發展,在電池流道入口處,氣體流場略有變化,但邊界層和流場分布基本穩定;同時入口延長后,電池流道入口靠近陰極處的速度遠小于入口平均速度,因此在陰極入口邊緣處的對流換熱會明顯降低,有助于緩解“入口效應”導致的入口處較強的對流換熱.

圖8 為600 s 時,入口延長后溫度梯度沿流動方向分布結果與原模型結果的對比圖.入口延長后電池入口處的溫度梯度會大幅降低,也即電池內最大溫度梯度會大幅降低,這主要是由于入口處陰極邊緣的速度降低,從而降低入口處的局部對流換熱,進而降低電池內的溫度梯度.不過由圖8 結果可知,入口延長后,電池的最大溫度梯度仍然非常大,電池入口段溫度梯度仍然存在一個迅速下降區域,顯然“入口效應”仍然制約入口處的溫度梯度.為分析其原因,分析了入口延長后陰極入口處溫度分布與原模型溫度分布的對比圖(如圖9 所示).由圖9 可知: 由于氣體入口溫度為均一值,在靠近氣體入口處存在非常大溫差,導致入口處存在較大的對流換熱熱流密度,從而在陰極內形成很大的溫度梯度;而入口段延長后,入口處氣體溫度只是略微降低,陰極入口邊緣處氣體與電池本體間溫差仍然很大,此處必然存在較大的對流換熱量,進而在陰極內形成很大溫度梯度.

圖8 600 s 時,入口延長與原模型得到的SOFC 內溫度梯度沿流道方向分布對比圖Fig.8.Comparison of the distribution of temperature gradient along gas channel in SOFC between the original model and extended entrance model at 600 s.

圖9 600 s 時,入口段延長后陰極入口處溫度分布與原模型溫度分布的對比圖Fig.9.Comparison of temperature distribution near cathode inlet between the original model and extended entrance model at 600 s.

綜上分析結果可知,當氣體入口速度和溫度設置為均一值時,由于入口效應的影響,電池入口處會形成非常大的溫度梯度,此時進行電池安全性分析時,并不能僅僅取電池內最大溫度梯度作為唯一標準;采用入口延長的方式雖然能夠降低入口處的最大溫度梯度,但入口處仍然存在較高的溫度梯度及溫度梯度變化,這是因為入口溫度仍為均一值;最后,當實際過程中采用熱氣體預熱電池時,電池入口處應該是熱應力最大處,因為實際過程中入口效應始終存在.

符號說明

E——內能,J·kg—1

hi——組分i的焓值,J·kg—1

keff——有效導熱系數,W·m—1·K—1

p——壓力,Pa

T——絕對溫度,K

t——時間,s

u——速度矢量,m·s—1

α——多孔介質滲透系數,m2

ε——多孔介質孔隙率;

μ——動力黏度,Pa·s

ρ ——密度,kg·m—3

下角標

f——熱氣體

s——固體

4 結論

本文建立了采用熱氣體預熱的固體氧化物燃料電池動態升溫模型,該模型結果可與文獻模型結果很好地吻合.采用氣體預熱模型,研究了在熱氣體預熱過程中電池內最大溫度梯度的位置及影響因素,并分析了電池內溫度梯度沿流動方向的分布規律,發現了數值模擬結果會在入口處得到異常高的溫度梯度;最后,研究了入口效應對于入口處高溫度梯度形成的影響,最終本文得到如下結論:

1)最大溫度梯度始終在電池入口邊緣處,電池與氣體接觸表面的第一層網格內;最大溫度梯度的位置不受熱氣體溫升速率的變化影響,但熱氣體入口速度會使其位置在垂直于流道方向變化,這主要是由于流速會改變氣體與電極和脊間的對流換熱;當采用雙流道預熱時,最大溫度梯度總位于陰極入口處,這是由于陰極的導熱系數小于陽極的導熱系數.

2)電池入口處存在非常高的溫度梯度,在靠近入口處很小的一段,溫度梯度沿流動方向迅速下降,隨后緩慢變化;當改變預熱氣體參數和通入方式,異常高溫度梯度區域仍然存在;因此,在實際數值模擬研究電池預熱安全性時,僅僅采用最大溫度梯度值作為電池安全性判據必然會得到不合理的結果,而采用平均溫度梯度也會在一定程度上高估電池內的熱應力.

3)在數值模擬中,一般氣體入口速度和溫度設置為均一值,由于入口效應的影響,電池入口處會形成非常大的溫度梯度,此時進行電池安全性分析時,并不能僅僅取電池內最大溫度梯度作為唯一標準;采用入口延長的方式雖然能夠降低入口處的最大溫度梯度,但入口處仍然存在較高的溫度梯度及溫度梯度變化,這是因為入口溫度仍為均一值.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美日韩理论| 国产18在线| 亚洲第一黄片大全| 亚洲一级毛片| 国产91在线|日本| 精品欧美一区二区三区久久久| 亚洲va视频| 毛片三级在线观看| 欧美在线中文字幕| 99伊人精品| 日本黄网在线观看| 亚洲国产精品无码AV| 亚洲第一区在线| 伊人天堂网| 尤物成AV人片在线观看| 真人高潮娇喘嗯啊在线观看| 91在线播放国产| 高潮毛片免费观看| 台湾AV国片精品女同性| 日本三区视频| 久久91精品牛牛| www.91在线播放| 色婷婷亚洲综合五月| 人妻精品久久久无码区色视| aa级毛片毛片免费观看久| 成人福利在线观看| 色综合久久88| 最新国产你懂的在线网址| 国产精品成人一区二区不卡| 国产成人精品男人的天堂| 曰AV在线无码| 国产激情在线视频| 亚洲黄色片免费看| 国产精品漂亮美女在线观看| 丝袜国产一区| 精品国产毛片| 拍国产真实乱人偷精品| 久久精品人人做人人| 久久不卡精品| 国产人人射| 自拍偷拍一区| 黄片在线永久| 无码日韩人妻精品久久蜜桃| 亚洲日韩精品欧美中文字幕 | 久草青青在线视频| 国产亚洲精品在天天在线麻豆| 永久成人无码激情视频免费| 欧美激情视频一区| 国产99精品视频| 欧美性精品不卡在线观看| 亚洲国产综合精品一区| 国产SUV精品一区二区6| 毛片网站在线看| 亚洲一道AV无码午夜福利| 高清久久精品亚洲日韩Av| 毛片久久久| 国产国拍精品视频免费看| 亚洲AV成人一区二区三区AV| 中文字幕久久亚洲一区| 91av国产在线| 亚洲成a人片| 无码精品国产dvd在线观看9久| 国产欧美另类| 91久久性奴调教国产免费| h视频在线播放| 99re66精品视频在线观看| 在线无码九区| 日本高清免费一本在线观看 | 国产精品成人免费视频99| 国产玖玖玖精品视频| 欧美在线中文字幕| 亚洲高清中文字幕| 人妻精品全国免费视频| 亚洲色图欧美在线| 色悠久久综合| 无码人妻热线精品视频| 亚洲欧美一级一级a| 91视频国产高清| 亚洲成av人无码综合在线观看| 91视频首页| 凹凸精品免费精品视频| AV色爱天堂网|