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

基于耦合傳熱的渦輪增壓器渦輪箱有限元分析

2015-10-28 09:58:48龔金科田應華黃張偉賈國海
中國機械工程 2015年10期

龔金科 田應華 黃張偉 賈國海

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

基于耦合傳熱的渦輪增壓器渦輪箱有限元分析

龔金科田應華黃張偉賈國海

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

基于渦輪增壓器渦輪箱傳熱機理,采用專業(yè)CFD軟件和FEM軟件分別建立了渦輪箱流體區(qū)域和固體區(qū)域網格仿真模型。在流體域建立多重旋轉坐標系,精確計算出渦輪箱流場、壁面?zhèn)鳠嵯禂导皽囟取昧鞴恬詈系姆抡娣椒▽u輪箱進行耦合傳熱分析,得到渦輪箱固體域的溫度場并對其進行熱應力分析。與實驗結果對比發(fā)現,仿真模型的溫度場符合實際渦輪箱溫度分布,最大誤差僅為3.3%。該渦輪箱耦合傳熱模型具有較高的精度,為渦輪增壓器的設計優(yōu)化提供了依據。

渦輪箱;耦合傳熱;溫度場;熱應力

0 引言

渦輪增壓器安裝在發(fā)動機的進排氣歧管上,在高溫、高壓和高速運轉的惡劣環(huán)境下工作,其各部件的溫度分布不均勻[1-3]。隨著渦輪增壓技術的發(fā)展,渦輪增壓器各零部件的工作環(huán)境越來越惡劣,渦輪箱等零部件的熱負荷也逐漸增大,而渦輪箱的溫度場和熱應力狀況很大程度上決定了渦輪增壓器的可靠性[4]。長期以來,國內外渦輪增壓器的研究多集中在渦輪葉片以及軸承體上,缺乏對渦輪箱流場、溫度場和熱疲勞的整體研究[5-7]。

本文基于渦輪箱傳熱機理,利用專業(yè)CFD軟件和FEM軟件對某發(fā)動機渦輪增壓器渦輪箱進行耦合傳熱數值仿真研究,通過在流體區(qū)域建立多重旋轉坐標系,模擬渦輪箱內腔廢氣的流動,得到渦輪高速旋轉時渦輪箱流體域的流場、壁面?zhèn)鳠嵯禂岛蜏囟葓龇植迹约捌涔腆w域的溫度場,分析了其熱應力,并將仿真計算值與實驗值進行了對比驗證。

1 渦輪增壓器渦輪箱傳熱原理

采用有限單元法對流體的流動與傳熱進行數值仿真,首先需要建立反映流體工程本質的數學模型。

高溫廢氣在渦輪箱中的流動與傳熱過程都遵從質量守恒、動量守恒和能量守恒定律。

廢氣在渦流道內的流動屬于湍流運動。簡化的標準k-ε模型方程如下[8]。

湍流動能k方程:

(1)

式中,k為湍流動能;ui為湍流速度;μ為流體動力黏度;ρ為流體密度;μt為湍流動黏度;Gk為由平均速度梯度引起的湍流動能k的產生項;Sk為湍流動能源項;σk為湍流動能k對應的Prandtl數,σk=1.0;ε為湍流耗散率。

湍流耗散率ε方程:

(2)

式中,Sε為湍流耗散源項;σε為湍流動能耗散率ε對應的Prandtl數,σε=1.3;C1ε和C2ε為經驗常數,C1ε=1.44,C2ε=1.92。

渦輪增壓器渦輪箱的傳熱是一個非常復雜的傳熱過程,其內表面以對流和熱傳導的方式與高溫廢氣進行換熱的同時,又以熱傳導的方式向渦輪箱體進行傳熱。

在渦輪箱和高溫氣體接觸的流固邊界上,傳熱過程是一個涉及固體、流體和溫度場等多個物理場同時作用的耦合傳熱過程。流固耦合是典型的弱耦合,只在邊界上存在熱量交換,其邊界上的溫度、傳熱系數都應看成是計算結果的一部分,而不是已知條件,在渦輪箱流固耦合傳熱邊界上有[9]

qw|solid=qw|f luid

(3)

當黏性流體在貼近壁面附近流動且流速很小時,其相對運動可忽略不計。在渦輪箱廢氣與壁面交界處的固體一側,由傅里葉導熱定律可知[10]:

qw|solid=-λgradt

(4)

式中,λ為熱導率;gradt為貼近壁面法線方向上流體溫度變化率。

對流傳熱的牛頓冷卻公式為[10]

qw|f luid=h(tw-tf)

(5)

式中,tw和tf分別為交界面和附近冷卻液的溫度;h為渦輪箱對流傳熱表面?zhèn)鳠嵯禂怠?/p>

由式(4)、式(5)可以得到對流傳熱表面?zhèn)鳠嵯禂蹬c流體溫度場的關系式:

(6)

在采用數值法進行研究時,若只應用單一綜合性有限元軟件對正常工況下渦輪增壓器渦輪箱進行數值模擬仿真計算,很難在減小計算量的同時保證較高的計算精度[9]。

本文采用ABAQUS軟件提供的流固耦合傳熱交界面模型,如圖1所示。在進行流固耦合傳熱仿真計算時,固體域和流體域之間可通過此交界面進行數據傳遞。當利用STAR-CCM+軟件仿真計算出流體區(qū)域的壁面溫度和傳熱系數后,將數據映射到圖1所示的耦合傳熱交界面上,這將成為ABAQUS軟件計算的第三類熱邊界條件。計算得到渦輪箱固體區(qū)域的溫度場,再通過耦合傳熱交接面將固體區(qū)域內腔壁面溫度映射到流體區(qū)域,這將成為流體軟件STAR-CCM+計算的壁面溫度條件。如此反復計算,直到渦輪箱體節(jié)點上的溫度不再發(fā)生變化為止。得到精確的渦輪箱固體區(qū)域溫度場分布后,將此溫度場導入ABAQUS預定義場,并利用其力學模塊計算渦輪箱熱應力,其過程如圖2所示。

圖1 耦合傳熱交界面

圖2 計算過程

2 渦輪增壓器渦輪箱仿真模型

2.1渦輪增壓器渦輪箱網格仿真模型

圖3 渦輪箱固體區(qū)域網格模型

在進行數值仿真計算前,對渦輪箱外表面倒角和細小結構進行簡化,可避免計算出錯并使計算收斂更快。通過幾何處理后采用四面體網格對渦輪箱劃分體網格,圖3為渦輪增壓器渦輪箱固體區(qū)域網格模型,其節(jié)點數為69 230,網格數為312 358。建立渦輪箱流固耦合模型,還需對流體區(qū)域進行網格劃分。為了精確模擬流體流動需要采用邊界層網格,在進出口區(qū)域設置拉伸層網格,以保證流體流動的穩(wěn)定性。運用幾何編輯軟件將導入的渦輪箱的外表面刪除,提取內腔壁面,對進出口處進行封面處理,將其轉換為實體并分割為三部分,即進口段、出口段以及渦輪旋轉區(qū)域,如圖4所示。

圖4 流體區(qū)域網格模型

在進行流體計算時,為了更精確地模擬高溫廢氣的流動,需建立交界面INTERFACE模型和多重旋轉坐標系MRF模型:分別在渦輪旋轉區(qū)域與流體進口段和出口段設置接觸面INTERFACE1和INTERFACE2,用于不同區(qū)域之間的物質及能量傳遞;在渦輪旋轉區(qū)域,建立MRF模型,并定義其旋轉方向和轉速,實現動網格的高速旋轉。圖4為渦輪箱流體區(qū)域網格模型,共有354 380個節(jié)點,1 327 725個單元。

2.2物理模型材料

準確的材料物性參數是仿真分析獲得準確結果的重要前提。在ABAQUS軟件中,對于固體域通常需要提供固體材料的熱導率、質量熱容、密度等參數。渦輪增壓器渦輪箱的材料選用灰鑄鐵,其密度為7010 kg/m3,質量熱容為510 J/(kg·K),泊松比為0.274,熱導率為47W/(m·K),線膨脹系數為1.5×10-5K-1,彈性模量為160 GPa。

2.3邊界條件

渦輪箱承受的熱載荷與發(fā)動機氣缸數、廢氣溫度等密切相關,是一個隨時間變化的交變載荷。本文在研究渦輪箱溫度及熱應力時,視高溫廢氣的流動為三維不可壓縮的黏性湍流流動,湍流模型采用k-ε湍流模型。入口采用質量流量邊界條件,根據提供資料,選取渦輪轉速為220 kr/min,使用GT-Power/Endenvironment模塊建模,計算得到正常工況下廢氣溫度T及質量流量qm隨時間變化情況,并導入STAR-CCM+,成為流體計算的入口邊界條件,如圖5所示。

圖5 增壓器入口端廢氣質量流量及溫度隨時間變化

渦輪箱采用壓力出口邊界條件,由于渦輪葉片高速旋轉導致氣流在出口端產生湍動,故其出口壓力p同樣隨時間變化,具體情況如圖6所示。

圖6 渦輪增壓器出口端廢氣壓力隨時間變化

3 計算結果分析

3.1渦輪箱流場分析

通過流體計算得到渦輪箱內部廢氣壓力以及速度分布。圖7顯示出了渦流道的壓降過程,越靠近渦輪處,壓降越明顯。在渦輪葉片處,由于廢氣高速沖擊渦輪,導致葉片與廢氣直接接觸的一側壓力較高,而葉片的另一側,壓力最低,為0.116 40 MPa。圖8為廢氣在渦流道內的速度矢量圖,在入口端,速度最小且分布均勻。當廢氣進入渦流道后,由于流道的引流作用,廢氣速度增大,對比其壓降云圖,由于能量守恒,流道內廢氣的勢能轉換為動能,在達到渦輪葉片時呈現最高速度,為594.28 m/s,保證了渦輪獲得動能最大。

圖7 渦輪箱渦流道壓力場分布

圖8 渦輪箱渦流道速度場分布

3.2渦輪箱內壁面?zhèn)鳠嵯禂蹬c溫度分布

發(fā)動機廢氣具有很高的動能和勢能,渦輪箱與之直接接觸,會接收從廢氣中傳出的大量能量。渦輪增壓器渦輪箱內壁面的傳熱系數直接反映了高溫廢氣和渦輪箱傳熱情況。

圖9為流體計算得到的渦輪箱內壁面?zhèn)鳠嵯禂档姆植荚茍D,其傳熱系數處于569.07~4083.9 W/(m2·K)之間。從渦流道至出口端面,傳熱系數先減小后增大,顯示出明顯的傳熱系數梯度。最高傳熱系數分布在進口端渦流道與渦輪旋轉區(qū)域交界面處,在此區(qū)域高溫廢氣高速沖擊渦輪葉片,并伴隨著大量勢能和動能之間的轉換。整個渦流道由于在能量轉換之前與高溫廢氣充分接觸,故此區(qū)域整體傳熱系數較大。

圖9 渦輪箱內壁面?zhèn)鳠嵯禂捣植紙D

傳熱系數直接影響溫度的分布,如圖10所示,可以看出渦輪箱內腔溫度分布明顯不均勻。總體從渦輪箱入口端至出口端呈現明顯的下降趨勢,高溫區(qū)域分布在渦輪旋轉區(qū)域,此區(qū)域由于劇烈傳熱導致局部溫度較高,最高達到1185.7 K。進口渦流道總體溫度較高,在1063.1 K以上,出口區(qū)域溫度最低,僅為572.72 K。

圖10 渦輪箱內壁面溫度分布圖

3.3渦輪增壓器渦輪箱溫度及熱應力分析

渦輪箱在正常工作時,由于高溫廢氣傳熱,故將其流體區(qū)域傳熱系數及溫度映射到渦輪箱固體網格內腔壁面上,作為計算固體域溫度場的第三類熱邊界條件。

在ABAQUS軟件中進行計算,得到渦輪箱溫度場分布情況,如圖11所示,渦輪箱溫度分布趨勢總體上向出口端遞減,這是由于渦輪箱主要的熱源來自發(fā)動機排放的廢氣,高溫廢氣通過渦輪箱流道時與之充分接觸,熱量迅速傳遞到渦流道壁面上,且廢氣大量勢能轉化為渦輪葉片的動能,越靠近出口端溫度越低。渦輪箱整體溫度較高,在624.2~1176 K之間。最高溫度在進氣渦道、進氣法蘭及V形圈邊,都超過1130 K,最低溫度分布在渦輪箱尾部出口端下緣,局部甚至溫度不足700 K。渦輪箱入口端和出口端溫差較大,從1176 K迅速降低到624 K。渦輪箱進氣法蘭端面溫度分布不均勻,在900~1100 K之間,越靠近渦流道的區(qū)域溫度越高。

圖11 渦輪箱整體溫度分布

圖12為渦輪箱熱應力分布情況,渦流道入口、V形圈邊等部位,熱應力遠大于其他部分熱應力,渦輪箱有兩邊向中間擠壓的變形趨勢,這主要由于渦輪箱進出口端被安裝固定,溫度變化使得渦輪箱高溫的部位產生膨脹,而這種膨脹受到法蘭、閥門圈座等固定部位的限制,從而產生局部熱應力。

圖12 渦輪箱整體熱應力分布

圖13顯示了渦輪箱流道截面溫度分布情況,在正常工作狀態(tài)下,渦流道長時間承受高溫,其內壁面的溫度達到最高,接近排氣溫度,渦輪箱流道整體溫度在1100 K以上,在渦流道導向噴嘴處,其內外側溫差較大。

圖13 渦輪箱渦流道截面溫度分布

圖14所示為渦輪箱渦流道截面熱應力分布情況,當發(fā)動機冷卻60 s后,內外側受熱狀況不一樣,渦流道外壁面迅速降溫,內表面由于殘氣繼續(xù)加熱,溫度下降較慢,導致溫度梯度增大,故渦流道局部溫度不均而產生熱應力,越靠近固定法蘭處熱應力越大。渦流道導向噴嘴處由于壁厚較薄,在冷卻過程中其內外側溫差較大,易產生較大熱應力。設計時在保證能量轉換效率下適當增加導向噴嘴厚度,可有效減小熱應力。

圖14 渦輪箱渦流道截面熱應力分布

圖15為渦輪箱整體截面的溫度分布云圖,從渦流道至出口端,溫度下降較快,形成明顯溫度梯度。針對本模型的單流道渦輪箱,其入口端的渦流道溫度較高,都處于1100 K以上,因此區(qū)域在高溫、高壓、高速的惡劣工況下工作并伴隨大量能量的轉換。渦輪箱整體渦流道部分溫度分布不均勻,越靠近廢氣出口端的渦流道部位溫度越低。

圖15 渦輪箱截面溫度分布

圖16顯示了渦輪箱截面熱應力分布情況。對于入口端,停機冷卻時,由于內部受到殘余廢氣的繼續(xù)加熱,其溫度較為穩(wěn)定,而法蘭的邊緣溫度下降較快,由于螺栓固定住了法蘭面,導致法蘭面高溫區(qū)域膨脹受擠壓,故產生熱應力及變形。在法蘭端螺栓孔處設置凸臺,可減小螺栓孔熱應力,同時在保證密封前提下采用彈性墊圈可緩解法蘭端部分熱應力。在V形圈邊處可適當增加其厚度,并在保證密封前提下采用彈性墊圈,可緩解由于限制熱膨脹而產生的熱應力。

圖16 渦輪箱截面熱應力分布

4 實驗驗證

對于機械零部件而言,熱應力產生的主要原因有兩個方面:一方面是不同材料由于不同的熱膨脹系數導致變形不均;另一方面是相同材料本身冷熱不均致使熱脹冷縮產生變形以及變形位移受到某些限制而導致。對于渦輪箱而言,箱體由相同材料構成,熱膨脹系數的不同由溫度差異造成,而這些差異不足以對熱膨脹系數產生太大的影響,可以忽略,因而,渦輪箱熱應力的主要來源是冷熱不均。渦輪箱處于冷熱交替變化的工作環(huán)境中,因而對其熱應力的研究實際在很大程度上可以等價為對其溫度場的研究[11]。基于本文中的渦輪增壓器實驗臺架,對渦輪增壓器渦輪箱的溫度分布進行實驗研究,能夠為流固耦合數值模擬計算提供更準確的邊界條件,同時也可以用于實驗結果與數值模擬結果的對比分析。

在渦輪增壓器渦輪箱的傳熱實驗中,采用WRTK-112工業(yè)鎧裝熱電偶分別對渦輪箱內外表面進行溫度測量并驗證模擬仿真的準確性[12]。由于測試條件的限制,部分測點的溫度無法直接測量,需要進行間接測量,如本文的渦輪箱溫度實驗,常用的方法是在渦輪箱壁面上,從外部向內部開測量孔,測量孔直徑為1 mm,測點底部安裝熱電偶球頭,位置距離渦輪箱內壁面為0.5 mm,溫度測量結果作為估算渦輪箱內壁面溫度值參考。考慮零件工藝要求和安裝的便捷性,溫度測點布置如圖17所示,其中測點1~6、10~12主要用于進行渦輪箱溫度測量實驗及結果分析;測點7~9、13~16采用間接測量的方式,在渦輪箱壁面上開測量孔,測量渦輪箱內壁面的溫度,主要用于進行傳熱邊界測量實驗。

圖17 渦輪箱測溫點分布

流固耦合傳熱邊界實驗在實驗臺架穩(wěn)態(tài)運行工況下進行,使用燃燒時加熱后排出的氣體模擬發(fā)動機廢氣驅動渦輪增壓器,在渦輪增壓器穩(wěn)定運行20 min后讀取各測點溫度值(圖18)。測量實驗測點溫度與數值模擬實驗結果對比如表1所示,各測點計算溫度值與實測值最大誤差僅為3.3%。即經過實驗驗證,利用流固耦合數值模擬的方法對渦輪增壓器渦輪箱內氣體流動及傳熱分析是準確及可靠的。

圖18 實驗設備

測點編號實測值(K)計算值(K)誤差(%)測點編號實測值(K)計算值(K)誤差(%)11152.11137.31.39822.4819.10.421119.91126.7-0.610812.8800.41.531162.51130.82.811895.9914.22.041104.61127.1-2.012979.8967.21.351143.41126.51.513929.8960.7-3.361083.41066.31.6141025.21056.9-3.071014.4991.62.3151051.31078.7-2.68888.9898.8-1.8161113.41138.7-2.2

5 結論

(1)基于MRF模型對渦輪箱進行流體計算,得到流體域的流場、壁面?zhèn)鳠嵯禂导皽囟确植记闆r。計算結果表明,廢氣的勢能和動能在渦輪處轉化最劇烈;從渦流道至出口端,傳熱系數先減小后增大,顯示出明顯的傳熱系數梯度,最大傳熱系數分布在流固耦合交界面處,其次是渦流道區(qū)域;傳熱系數直接影響渦輪箱的溫度分布,溫度總體上呈現從進口端向出口端遞減趨勢,在渦輪箱出口端溫度最低。

(2)渦輪箱整體溫度較高使得高溫的部位產生膨脹,各部位由于溫度不均使膨脹程度不一致,且受到法蘭等固定部位的限制,從而產生局部熱應力。其中進氣法蘭端面、渦流道以及導向噴嘴等處熱應力最大,應在設計時進行優(yōu)化。

(3)通過測量實驗測點溫度并與數值模擬實驗結果進行對比,可知數值模擬結果非常接近實驗值,驗證了利用流固耦合數值模擬對渦輪增壓器渦輪箱內廢氣流動及傳熱分析的方法準確性和可靠性。

[1]張俊紅,李志剛,王鐵寧,等.車用渦輪增壓技術的發(fā)展回顧、現狀及展望[J].小型內燃機與摩托車,2007,36(1):66-69.Zhang Junhong,Li Zhigang,Wang Tiening,et al.The Past,Present and Prospects of the Development on Vehicle Turbo Technology[J].Small Internal Combustion Engines and Motorcycles,2007(1):66-69.[2]周虹偉.渦輪增壓器主要部件結構特性及改進研究[J].中國鐵道科學,2004,25(2):72-77.

Zhou Hongwei.The Structure Characteristics and Improvement of Main Components of Turbo-supercharger[J].China Railway Science,2004,25(2):72-77.

[3]李志江.柴油機廢氣渦輪增壓器的常見故障與使用保養(yǎng)[J].工程機械與維修,2008(4):198-199.

Li Zhijiang.The Common Fault and User Maintenance on Diesel Exhaust Turbocharger[J].Construction Machinery & Maintenance,2008(4):198-199.

[4]廖愛華.增壓器的非線性力學分析[D].大連:大連理工大學,2007.

[5]胡友安,李曉東,陳圖鈞,等.渦殼溫度場和熱應力的有限元分析[J].機械強度,2007,29(1):130-134.Hu Youan,Li Xiaodong,Chen Tujun,et al.Analysis on Temperature Field and Thermal Stress of Turbocharger Housing by Finite Element Method[J].Journal of Mechanical Strength,2007,29(1):130-134.[6]Errera M P, Chemin S. Optimal Solutions of NumericalInterface Conditions in Fluid-structure Thermal Analysis[J].Journal of Computational Physics,2013,245:431-455.

[7]Pesiridis A,Martinez-botas R F.Experimental Evaluation of Active Flow Control Mixed-flow Turbine for Automotive Turbocharger Application[J].Journal of Turbo Machinery,2007,129(1):44-52.

[8]徐行軍.柴油機冷卻系統(tǒng)結構優(yōu)化及缸蓋熱應力分析[D].天津:天津工業(yè)大學,2010.

[9]傅松,胡玉平,李新才,等.柴油機缸蓋水腔流動與沸騰傳熱的流固耦合數值模擬[J].農業(yè)機械學報,2010,41(4):26-30.

Fu Song,Hu Ping,Li Xincai,et al.Chamber Cylinder Head Water Flow and Heat Transfer Simulation Boiling Fluid-structure Interaction[J].Agricultural Machinery,2010,41(4):26-30.

[10]賈力,方肇洪. 高等傳熱學[M]. 北京:高等教育出版社,2008.

[11]胡友安,李曉東.渦輪增壓器蝸殼熱裂紋的試驗研究[J].河海大學學報,2008,36(6):846-849.

Hu Youan,Li Xiaodong.Experimental Study on Heat Cracks of Turbine Housing on Turbochargers[J].Journal of Hehai University,2008,36(6):846-849.

[12]徐思友,吳新濤,閆瑞乾,等.增壓器軸承和密封環(huán)溫度試驗研究[J].車用發(fā)動機,2010(2):35-37.

Xu Siyou,Wu Xintao,Yan Ruiqian,et al.Experimental Study on Temperature for Turbocharger Bearing and Sealing Ring[J].Vehicle Engine,2010(2):35-37.

(編輯王艷麗)

Finite Element Analysis on Turbocharger Turbine Box Based on Coupled Heat Transfer

Gong JinkeTian YinghuaHuang ZhangweiJia Guohai

State Key Laboratory of Advanced Design and Manufacture for Vehicle Body of Hunan University,Changsha,410082

Based on the heat transfer mechanism of the turbocharger turbine box,a fluid mesh model and a solid one for simulation were built by using the professional software CFD and FEM respectively.The multiple coordinate system of rotation in the fluid domain was established to calculate accurately the flow field of turbine boxes,heat transfer coefficient and temperature.Coupled heat transfer analysis of the turbocharger turbine box was carried out by the fluid-solid coupling method,and then the temperature field for solid areas was obtained.Based on the temperature,its thermal stress analysis was carried out. In comparison with the experiments,the temperature distribution of the simulation model accords with that of the realistic turbocharger turbine box, the maximum deviation is only 3.3%.The coupled heat transfer model for the turbocharger turbine box has a high accuracy,which provides the basis for better design of turbocharger turbine.

turbine box;coupled heat transfer;temperature field;thermal stress

2014-07-04

國家高技術研究發(fā)展計劃(863計劃)資助項目(2008AA11A116);汽車車身先進設計制造國家重點實驗室自主課題資助項目(61075002);湖南省自然科學基金資助項目(10JJ6080)

TK411.8DOI:10.3969/j.issn.1004-132X.2015.10.012

龔金科,男,1954年生。湖南大學機械與運載工程學院教授、博士研究生導師。主要研究方向為汽車排放及控制技術,熱動力設備新技術等。田應華,男,1988年生。湖南大學機械與運載工程學院碩士研究生。黃張偉,男,1986年生。湖南大學機械與運載工程學院碩士研究生。賈國海,男,1986年生。湖南大學機械與運載工程學院博士研究生。

主站蜘蛛池模板: 亚洲天堂首页| 亚洲中文字幕在线精品一区| 国产裸舞福利在线视频合集| 波多野结衣第一页| 免费人成视网站在线不卡 | 精品无码一区二区三区在线视频| 日本一本正道综合久久dvd | 亚洲国产中文在线二区三区免| 久久青青草原亚洲av无码| 中文字幕人成人乱码亚洲电影| 国产性精品| 亚洲欧美日本国产综合在线| 精品国产美女福到在线直播| 国产精品美女在线| 国产网友愉拍精品视频| 欧美激情二区三区| 99在线视频免费观看| 精品人妻一区无码视频| 日本免费福利视频| 久久天天躁夜夜躁狠狠| 亚洲三级影院| 亚洲中字无码AV电影在线观看| 无码'专区第一页| 亚洲永久色| 久久精品女人天堂aaa| 国产欧美精品一区二区| 欧洲欧美人成免费全部视频| 婷婷伊人五月| 97视频免费在线观看| 欧美全免费aaaaaa特黄在线| 亚洲性影院| 亚洲伊人电影| 欧美成人午夜在线全部免费| 亚洲国产91人成在线| 国产丰满成熟女性性满足视频| 国产在线观看成人91| 婷婷综合缴情亚洲五月伊| 欧美一区二区精品久久久| 国产天天射| 亚洲国产精品VA在线看黑人| 一级全黄毛片| 日本一区二区不卡视频| Jizz国产色系免费| 国产欧美日韩综合在线第一| 丁香五月亚洲综合在线| 日本免费新一区视频| 亚洲AV永久无码精品古装片| 亚洲熟女中文字幕男人总站| 国产日本欧美亚洲精品视| 中文字幕人成人乱码亚洲电影| 91小视频在线观看免费版高清| 麻豆国产精品| 久久精品亚洲中文字幕乱码| 亚洲精品动漫在线观看| 在线看片中文字幕| 国产精品福利在线观看无码卡| 美女视频黄频a免费高清不卡| 国产女人在线| 国产自无码视频在线观看| 国产丝袜丝视频在线观看| 亚洲综合狠狠| 亚洲AV无码乱码在线观看裸奔 | 91久久偷偷做嫩草影院| 国产呦精品一区二区三区下载| 狠狠做深爱婷婷综合一区| v天堂中文在线| 99久久国产综合精品2020| 区国产精品搜索视频| 欧美中日韩在线| 久草中文网| 蜜臀av性久久久久蜜臀aⅴ麻豆| 日韩少妇激情一区二区| 日本在线亚洲| 国产在线观看91精品亚瑟| 亚洲成人福利网站| 亚洲视频欧美不卡| 日本91视频| 99热国产这里只有精品9九| 无码精品一区二区久久久| 久久熟女AV| 欧美日本在线观看| 毛片基地视频|