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

基于RBF神經(jīng)網(wǎng)絡(luò)的某型裝甲車輛散熱器優(yōu)化設(shè)計(jì)

2021-03-05 14:03:28張玉飛邢俊文馮鵬飛康洪偉
兵器裝備工程學(xué)報(bào) 2021年2期
關(guān)鍵詞:模型設(shè)計(jì)

張玉飛,邢俊文,馮鵬飛,康洪偉,趙 耀

(1.陸軍裝甲兵學(xué)院 車輛工程系, 北京 100072; 2.中國人民解放軍95816部隊(duì), 武漢 432200)

散熱器是裝甲車輛冷卻系統(tǒng)中的重要部件,其性能對裝甲車輛的動力性、經(jīng)濟(jì)性和可靠性有至關(guān)重要的作用[1]。因此,散熱器的優(yōu)化設(shè)計(jì)是冷卻系統(tǒng)設(shè)計(jì)中的重要環(huán)節(jié)[2]。對于車用散熱器的優(yōu)化設(shè)計(jì),國內(nèi)外很多學(xué)者進(jìn)行了相關(guān)研究。閆玉英等在對散熱器流動特性分析的基礎(chǔ)上,建立了數(shù)學(xué)模型,利用罰函數(shù)法進(jìn)行了設(shè)計(jì)[3]。李曉光利用FLUENT建立了某型商用汽車散熱器的數(shù)值計(jì)算模型,并利用風(fēng)洞試驗(yàn)測試了相關(guān)數(shù)據(jù),討論了翅片結(jié)構(gòu)對散熱器性能的影響[4]。王任遠(yuǎn)利用Laminar法針對散熱器建立了流固耦合模型,研究了不同結(jié)構(gòu)參數(shù)對傳熱系數(shù)的影響,最后得到了一組最優(yōu)解[5]。向高利用風(fēng)洞試驗(yàn)數(shù)據(jù)擬合得到了數(shù)學(xué)模型,在此基礎(chǔ)上利用FLUENT軟件分析了各結(jié)構(gòu)因素對散熱性能的影響,最后得到了綜合性能較佳的散熱器結(jié)構(gòu)參數(shù)[6]。景陶敬等針對當(dāng)前汽車散熱器散熱性能不佳的情況,利用試驗(yàn)與仿真結(jié)合的方法,討論了流速散熱性能的影響[7]。目前,針對車輛散熱器優(yōu)化時(shí)考慮的因素較少,多將優(yōu)化簡化為線性問題來解決,隨著高功率密度柴油機(jī)和機(jī)電混合傳動的采用,熱負(fù)荷大幅度增加,因此裝甲車輛的冷卻系統(tǒng)熱負(fù)荷較大,考慮的因素較多,需要進(jìn)行多目標(biāo)優(yōu)化。

本文利用解析法和基于CFX的流體力學(xué)仿真模型分別對某型散熱器進(jìn)行了計(jì)算,驗(yàn)證了CFD仿真模型的精度;選取了目標(biāo)函數(shù)和設(shè)計(jì)變量,利用最優(yōu)拉丁超立方設(shè)計(jì)法構(gòu)建了設(shè)計(jì)空間,建立了基于徑向基神經(jīng)網(wǎng)絡(luò)的替代模型,利用對照組數(shù)據(jù)檢驗(yàn)了模型的精度;利用螢火蟲算法對替代模型進(jìn)行了多目標(biāo)優(yōu)化,得到了最優(yōu)解集;最后利用蒙特卡洛法對影響因素進(jìn)行了統(tǒng)計(jì)學(xué)分析。

1 散熱器建模及計(jì)算

1.1 幾何模型及邊界條件設(shè)置

1.1.1計(jì)算模型的建立

該型散熱器的芯體主要由水管和散熱片等組成。在散熱片上,按水管的截面形狀和尺寸沖孔,再將它串套在水管上,然后在一定溫度下進(jìn)行焊接構(gòu)成芯體。散熱器的幾何實(shí)體模型如圖1所示。散熱器芯體幾何結(jié)構(gòu)的局部示意圖如圖2。

圖1 散熱器幾何模型示意圖

圖2 散熱器幾何結(jié)構(gòu)局部示意圖

散熱器中流體的流動與傳熱分析是散熱器設(shè)計(jì)和冷卻系統(tǒng)匹配設(shè)計(jì)的重要環(huán)節(jié)[8]。要進(jìn)行數(shù)值計(jì)算,必須首先獲得待求解空間區(qū)域的幾何模型,由于車輛散熱器中,流體在散熱器內(nèi)外流道中流動,內(nèi)外流道之間通過對流、導(dǎo)熱完成熱量的交換,在利用CFX建立模型時(shí)利用完整的幾何模型會造成網(wǎng)格數(shù)量大,收斂慢等問題,因此對散熱器整體進(jìn)行研究時(shí),必須對散熱器模型進(jìn)行簡化處理[9]。為了提高計(jì)算效率,根據(jù)散熱器的結(jié)構(gòu)特點(diǎn),對計(jì)算區(qū)域確定如下:計(jì)算區(qū)域的長度定為芯體的長度,高度定為兩個(gè)散熱片之間的距離,寬度設(shè)定為恰好能涵蓋水管錯(cuò)列位置為一個(gè)周期的距離,選取的區(qū)域計(jì)算模型如圖3所示。

圖3 區(qū)域計(jì)算模型示意圖

1.1.2網(wǎng)格無關(guān)性檢驗(yàn)

由于冷卻空氣和冷卻液之間的水管厚度很薄且面積較大,因此,這里將其作為壁面進(jìn)行處理,以減少網(wǎng)格數(shù)量。為了保證仿真精度,減少網(wǎng)格數(shù)量對計(jì)算的影響,因此需要進(jìn)行網(wǎng)格無關(guān)性檢驗(yàn)[10]。在進(jìn)行網(wǎng)格劃分時(shí),網(wǎng)格相關(guān)度共取5組,范圍由0~100之間。在入口速度為11 m/s情況下進(jìn)行計(jì)算。圖4為空氣流過散熱器內(nèi)部的壓力損失隨網(wǎng)格相關(guān)度的變化情況。從計(jì)算結(jié)果可以看出,隨著網(wǎng)格相關(guān)度的增加,壓力損失出現(xiàn)下降趨勢,網(wǎng)格相關(guān)度由60~100的變化中,相對誤差為1.83%,均值為602.643,標(biāo)準(zhǔn)差為8.264,變異系數(shù)為0.013,大致趨于穩(wěn)定,從提高計(jì)算效率角度可以看作網(wǎng)格的變化對計(jì)算結(jié)果不再產(chǎn)生影響,因此網(wǎng)格無關(guān)度取為60,網(wǎng)格數(shù)量為298 112。

圖4 散熱器壓力損失隨網(wǎng)格相關(guān)度的變化

1.1.3邊界條件設(shè)置

已知空氣的流量為0.81 kg/s,壓力為24.5×104Pa,溫度為160 ℃,冷卻水的流量為14.4×103 kg/h,冷卻水進(jìn)口溫度為25 ℃。

根據(jù)散熱器芯體結(jié)構(gòu)以及分析的目的,確定了管片式散熱器芯體空氣側(cè)CFD分析區(qū)域及邊界條件設(shè)置如圖5所示,在入口處設(shè)置流量入口邊界,出口處設(shè)置自由出口邊界,水管處設(shè)置熱邊界條件,管片設(shè)置為壁面條件,橫向壁面設(shè)置為周期性邊界條件[11]。邊界條件設(shè)置如圖5所示。

圖5 邊界條件設(shè)置示意圖

1.2 解析法

裝甲車輛上使用的散熱器一般為緊湊式散熱器,冷卻空氣流經(jīng)緊湊式散熱器時(shí)的壓力損失可表示為

Δp=Δpi+Δpcf-Δp0

(1)

式中:Δpi為散熱器芯體進(jìn)口處的壓力損失;Δp0為芯體出口處的壓力回升;Δpcf為芯體內(nèi)的壓力損失。

由理論分析得出,流體流經(jīng)散熱器芯體入口和出口處的壓力損失為

(2)

(3)

式中:vm為質(zhì)量流速;v1為芯體進(jìn)口截面處的比體積;v2為芯體出口截面處的比體積;σ為散熱器的流通面積與迎風(fēng)面積之比;Kc為入口壓力損失系數(shù);Ke為出口壓力損失系數(shù)。

散熱器芯體內(nèi)的壓力損失計(jì)算公式為

(4)

式中:f為空氣側(cè)摩擦因子;L為流道長度;vL為流體平均比容;De為流道當(dāng)量直徑。

當(dāng)在散熱器中流動時(shí),由于流體流經(jīng)每一管排時(shí)都有一次收縮和膨脹,芯體的第一管排和最后一個(gè)管排與中間各管排的阻力情況并無顯著差別,因此,進(jìn)出口處的壓力損失實(shí)際上已經(jīng)包括芯體內(nèi)的黏性摩擦損失中,因此壓力損失的計(jì)算公式可化簡為

(5)

由于散熱器芯體的流體流動加速所引起的壓力損失較小一般可以忽略,上式可以簡化為

(6)

1.3 結(jié)果對比及分析

計(jì)算采用CFX進(jìn)行計(jì)算,設(shè)定參考壓力為101 325 Pa,為了方便與解析解進(jìn)行對比,在散熱器的質(zhì)量流量分別設(shè)為0.41 kg/s、0.55 kg/s、0.69 kg/s、0.77 kg/s、0.81 kg/s,其中根據(jù)文獻(xiàn)[5],該型散熱器在流量為0.81 kg/s時(shí)的試驗(yàn)值為685.2 Pa,通過計(jì)算得到了散熱器在不同的入口質(zhì)量流量下CFD仿真和解析解以及試驗(yàn)值(見圖6)。

圖6 仿真值與解析解曲線

通過對比,在空氣質(zhì)量流量為0.81 kg/s處,仿真值和解析解分別為673.167 Pa和675 Pa,殘差為1.833 Pa,試驗(yàn)值比二者較高,較仿真值,其相對誤差為1.76%。仿真值與解析解之間的最大相對誤差為3.9%,均方誤差為0.000 52,決定系數(shù)為0.999 137 951,相關(guān)系數(shù)為0.999 324 7,從統(tǒng)計(jì)學(xué)角度看,兩組解較為接近,通過以上分析,可以看出仿真值的結(jié)果比較可靠。圖7為局部散熱片速度流線圖,可以看出,最大速度達(dá)到12.29 m/s,且沿空氣流線方向,速度大小分布也不均勻,大致呈階梯狀分布。圖8為局部散熱器片壓力分布圖,從圖8中可以看出出口處壓力較小,最大值位于速度入口處,達(dá)到了223.26 Pa。

圖7 局部散熱器片速度流線圖

圖8 局部散熱器片壓力分布圖

2 基于徑向基神經(jīng)網(wǎng)絡(luò)的替代模型建立

2.1 設(shè)計(jì)變量的確立

本研究的優(yōu)化對象為某型裝甲車輛管帶式散熱器,該散熱器的幾何模型如圖9所示。該種散熱器的散熱片為百葉窗形翅片,冷卻空氣從散熱片中流過,冷卻液從扁平管內(nèi)流過。冷卻液的熱量經(jīng)過管壁表面的對流換熱和翅片的導(dǎo)熱,最后由冷卻空氣將熱量帶走,由于該型散熱器的結(jié)構(gòu)較大,特征重復(fù)較多,因此本文的計(jì)算區(qū)域取兩根水管之間的寬度,高度選為能夠涵蓋散熱片變化恰好為一個(gè)周期的距離[12]。

圖9 散熱器幾何模型示意圖

CFD計(jì)算得到了當(dāng)迎風(fēng)速度為4 m/s、8 m/s、12 m/s情況下的速度云圖如圖10、圖11和圖12所示。從圖10~圖12中可以看出,隨著迎面風(fēng)速的增加,管帶式散熱器內(nèi)部流動特性變化明顯;空氣從入口百葉窗進(jìn)入,沿著百葉窗流動,在百葉窗中部變向,最后從出口百葉窗流出。氣流由一個(gè)百葉窗流入到另一個(gè)百葉窗相當(dāng)于一個(gè)收縮管道,氣流加速運(yùn)動,因此百葉窗翅片間流動速度比較大;流速最大的地方出現(xiàn)在模型中部偏后位置,也就是流動改變方向的第一個(gè)翅片通道,主要原因是氣流在流經(jīng)該區(qū)域之前有一個(gè)氣流加速,經(jīng)過該處氣流換向時(shí)通流截面突然減小所致。迎面速度為 4 m/s 時(shí),空氣以軸向流動為主,而迎面速度為8 m/s和 12 m/s 時(shí),空氣以穿越百葉窗區(qū)域?yàn)橹鳎饕蚴抢字Z數(shù)較低時(shí),在翅片近壁面處,空氣形成很厚的邊界層,阻礙空氣向窗翅區(qū)流動,隨著雷諾數(shù)的增加,空氣流動邊界層變薄,窗翅區(qū)空氣流動阻力降低,空氣的流動方向也逐漸向窗翅區(qū)過渡。

圖10 迎面風(fēng)速為4 m/s時(shí)速度云圖

圖11 迎面風(fēng)速為8 m/s時(shí)速度云圖

圖12 迎面風(fēng)速為12 m/s時(shí)速度云圖

圖13為風(fēng)速為12 m/s時(shí)百葉窗散熱器流動方向的壓力變化曲線,壓力隨著流動方向長度的增加而降低。在流動單元中部附近變化曲率較小,而在其他地方基本呈線性變化。流動單元前半部分壓降略小于后半部分壓降,這也是由于后半部分流速變大,雷諾數(shù)增加,流體穿越百葉窗也越多,流速增加導(dǎo)致摩擦阻力損失增大,從而使得壓力損失也變大。

圖13 管帶式散熱器流動方向壓力分布曲線

2.2 設(shè)計(jì)變量和目標(biāo)函數(shù)

車輛散熱器在設(shè)計(jì)時(shí),應(yīng)根據(jù)車輛冷卻系統(tǒng)的要求來設(shè)計(jì)散熱器的結(jié)構(gòu)及其參數(shù)[13],車輛冷卻系統(tǒng)對其系統(tǒng)內(nèi)散熱器的要求,可以是在給定的空間容積和允許的壓力損失條件下,使得散熱器具有盡可能大的散熱能力,散熱器熱力設(shè)計(jì)優(yōu)化模型應(yīng)該包括目標(biāo)函數(shù)、設(shè)計(jì)變量和約束方程。該散熱器的翅片結(jié)構(gòu)為百葉窗形翅片,散熱器幾何結(jié)構(gòu)參數(shù)如圖14所示,其中l(wèi)、a、Fh、Fp和δf為百葉窗間距、百葉窗開度、翅片高、翅片間距和翅片厚度,本研究選取這5個(gè)參數(shù)作為設(shè)計(jì)變量。目標(biāo)函數(shù)為散熱量Φ和空氣側(cè)壓力損失Δp。設(shè)計(jì)變量的約束條件為

(7)

圖14 散熱器幾何結(jié)構(gòu)參數(shù)示意圖

2.3 樣本點(diǎn)抽取

最優(yōu)拉丁超立方設(shè)計(jì)改進(jìn)了隨機(jī)拉丁超立方設(shè)計(jì)的均勻性,使因子和響應(yīng)的擬合更加精確真實(shí),具有非常好的空間填充性和均衡性[14]。本次共抽取了15組樣本點(diǎn),最后抽取樣本點(diǎn)如表1所示。

不同種類的車輛,由于其發(fā)動機(jī)經(jīng)常所處的工況不相同,故其散熱器的設(shè)計(jì)工況點(diǎn)和校核工況點(diǎn)也不相同。本文為裝甲車輛用的水散熱器,發(fā)動機(jī)較多地處于最大扭矩工況點(diǎn)附近工作,風(fēng)扇和水泵的工作轉(zhuǎn)速都比較低,使得散熱器處在比較惡劣的條件下工作[15]。所以,應(yīng)該把最大扭矩點(diǎn)工況作為散熱器的設(shè)計(jì)工況,而把發(fā)動機(jī)的額定工況作為校核工況。

表1 樣本點(diǎn)數(shù)據(jù)

2.4 替代模型的建立

以最大扭矩工況點(diǎn)為例,利用CFD模型對抽取的15個(gè)樣本點(diǎn)進(jìn)行了計(jì)算,利用徑向基神經(jīng)網(wǎng)絡(luò)建立替代模型。根據(jù)實(shí)驗(yàn)設(shè)計(jì)點(diǎn)的值及其響應(yīng),得到了這5個(gè)參數(shù)的部分替代模型,如圖15所示。

圖15 部分替代模型示意圖

為了驗(yàn)證替代模型的精度,選用決定系數(shù)來評價(jià)擬合效果,決定系數(shù)R2的計(jì)算公式為

圖16 決定系數(shù)計(jì)算曲線

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

圖17(a)和圖17(b)為百葉窗間距l(xiāng)按照正態(tài)分布進(jìn)行1 000次蒙特卡洛抽樣計(jì)算得到的散熱量Φ和空氣側(cè)壓力損失Δp的關(guān)系曲線,另外4個(gè)設(shè)計(jì)變量均取取值范圍的中間值。從圖17中可以看出:散熱量隨著百葉窗間距的增加而變大,壓力損失隨著百葉窗間距的增加先增大后減小。圖18(a)和圖18(b)為百葉窗開度a按照正態(tài)分布進(jìn)行1 000次蒙特卡洛抽樣計(jì)算得到的散熱量Φ和空氣側(cè)壓力損失Δp的關(guān)系曲線。從圖18中可以看出:散熱量隨著百葉窗開度的增加而變大,壓力損失隨著百葉窗開度的增加先增大后減小。圖19(a)和圖19(b)為翅片高Fh按照正態(tài)分布進(jìn)行 1 000 次蒙特卡洛抽樣計(jì)算得到的散熱量Φ和空氣側(cè)壓力損失Δp的關(guān)系曲線。從圖19中可以看出:散熱量隨著百葉窗開度的增加而降低,壓力損失隨著百葉窗開度的增加而減小。

圖17 百葉窗間距與散熱量和壓力損失的關(guān)系曲線

圖18 百葉窗開與散熱量和壓力損失的關(guān)系曲線

圖19 翅片高與散熱量和壓力損失的關(guān)系曲線

圖20(a)和圖20(b)為翅片間距Fp按照正態(tài)分布進(jìn)行1 000次蒙特卡洛抽樣計(jì)算得到的散熱量Φ和空氣側(cè)壓力損失Δp的關(guān)系曲線。從圖20中可以看出:散熱量隨著百葉窗開度的增加而降低,壓力損失隨著百葉窗開度的增加先變小后增大。圖21(a)和圖21(b)為翅片厚度δf按照正態(tài)分布進(jìn)行1 000次蒙特卡洛抽樣計(jì)算得到的散熱量Φ和空氣側(cè)壓力損失Δp的關(guān)系曲線。

圖20 對翅片間距抽樣得到的分布

圖21 對翅片厚度抽樣得到的分布

從圖21中可以看出:散熱量隨著百葉窗開度的增加先減小后增加,壓力損失隨著百葉窗開度的增加先變小后增大。從上述這些規(guī)律中可以看出,每個(gè)設(shè)計(jì)變量相對于目標(biāo)函數(shù)的變化趨勢是不同甚至是矛盾的,因此不可能通過調(diào)整設(shè)計(jì)變量讓兩個(gè)目標(biāo)函數(shù)同時(shí)達(dá)到最優(yōu),只能在各個(gè)目標(biāo)之間進(jìn)行協(xié)調(diào)和權(quán)衡,盡可能使每個(gè)目標(biāo)盡可能均達(dá)到最優(yōu),而這種優(yōu)化解也不可能是單一解,而是一個(gè)解集,一般稱為帕累托前沿。

3 基于螢火蟲算法的多目標(biāo)優(yōu)化

自然界中的群居昆蟲,它們雖然個(gè)體結(jié)構(gòu)簡單,但是通過個(gè)體間的合作卻能夠表現(xiàn)極其復(fù)雜的行為能力。螢火蟲算法是一種模擬螢火蟲行為的群集只能優(yōu)化算法,由于它具有控制參數(shù)少且易于實(shí)現(xiàn)等特點(diǎn),為解決存在于科學(xué)領(lǐng)域的多目標(biāo)優(yōu)化問題提供了一種新的方法。基于螢火蟲算法的多目標(biāo)優(yōu)化流程如圖22所示。

利用螢火蟲算法進(jìn)行計(jì)算時(shí),需要設(shè)置相關(guān)參數(shù)。螢火蟲算法參數(shù)設(shè)置如表2所示。

利用螢火蟲算法求得的帕累托前沿如圖23所示。從圖23中可以看出,散熱量的取值范圍為797~807 kW,壓力損失的取值范圍為400~460 Pa,從散點(diǎn)的分布情況來看,不可能使二者同時(shí)達(dá)到最優(yōu),如果決策者追求散熱量盡可能地大,應(yīng)該往右上角取值,如果對壓力損失比較重視,則應(yīng)該往左下角取值。本文對壓力損失和散熱量綜合考慮,選擇了一組使各個(gè)目標(biāo)盡可能最優(yōu)的折衷方案,該組的位置如圖23所示,設(shè)計(jì)變量的尺寸分別為l=8.22,a=1.43,F(xiàn)h=6.5,F(xiàn)p=5,δf=0.12,將優(yōu)化后的尺寸利用CFD仿真可以得到散熱量為803.5 kW,壓力損失為417.559 Pa,與優(yōu)化前相比,散熱量增加了6%,壓力損失降低了13.8%。

圖22 基于螢火蟲算法的多目標(biāo)優(yōu)化流程框圖

表2 螢火蟲算法參數(shù)設(shè)置

圖23 帕累托前沿示意圖

4 結(jié)論

1) 本文所建立的基于CFX的裝甲車輛散熱器模型的精度較高,與解析解之間的最大相對誤差僅為3.9%。

2) 通過對散熱器的仿真計(jì)算可以散熱片中的空氣流速分布不均勻,大致呈階梯狀分布。散熱器片出口處壓力較小,最大值位于速度入口處。

3) 所建立的基于RBF神經(jīng)網(wǎng)絡(luò)的替代模型精度較高,目標(biāo)函數(shù)為散熱量Φ和壓力損失Δp的決定系數(shù)計(jì)算結(jié)果為0.989 7和0.964 28,完全滿足計(jì)算精度。

4) 利用螢火蟲算法所得到的優(yōu)化后的結(jié)構(gòu)散熱量增加6%,壓力損失降低13.8%。

猜你喜歡
模型設(shè)計(jì)
一半模型
重要模型『一線三等角』
何為設(shè)計(jì)的守護(hù)之道?
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
《豐收的喜悅展示設(shè)計(jì)》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設(shè)計(jì)萌到家
設(shè)計(jì)秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設(shè)計(jì)叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 久久婷婷五月综合97色| www.国产福利| 香蕉久人久人青草青草| 国产正在播放| 成人午夜精品一级毛片| 国产sm重味一区二区三区| 色偷偷一区| 国产激情第一页| 欧美成人午夜在线全部免费| 麻豆精品在线| 国产无码高清视频不卡| 一区二区三区四区精品视频 | 国产96在线 | 欧美成人午夜视频免看| 日韩无码精品人妻| 欧美激情第一欧美在线| 亚洲成A人V欧美综合| 国产成人亚洲欧美激情| www中文字幕在线观看| 片在线无码观看| 国产成a人片在线播放| 亚洲无码免费黄色网址| 国产精品天干天干在线观看| 97亚洲色综久久精品| 久久成人18免费| 亚洲欧洲综合| 色综合中文综合网| 日韩国产亚洲一区二区在线观看| 国产靠逼视频| a毛片免费观看| 高清国产在线| 欧美第二区| 爱色欧美亚洲综合图区| 久久国产成人精品国产成人亚洲| 欧美亚洲国产精品第一页| 亚洲最大福利网站| 国产福利在线免费| 国产男女免费完整版视频| 波多野结衣中文字幕久久| 亚洲天堂精品在线观看| 亚洲国产看片基地久久1024| 久久无码高潮喷水| 2021最新国产精品网站| 国产成人综合日韩精品无码首页| 久久久久人妻一区精品色奶水| 大陆国产精品视频| 精品国产欧美精品v| 日本道综合一本久久久88| 亚洲国产精品无码久久一线| 婷婷综合缴情亚洲五月伊| 国产区在线观看视频| 久久九九热视频| 天天摸夜夜操| 亚洲人成网站日本片| 日韩国产综合精选| 美女毛片在线| 国产玖玖视频| 草草影院国产第一页| 一区二区在线视频免费观看| 久久性妇女精品免费| 亚洲国产清纯| 欧美无遮挡国产欧美另类| 日本道中文字幕久久一区| 亚洲天堂视频在线观看免费 | 国产一区二区三区日韩精品| 啊嗯不日本网站| 日本午夜在线视频| 亚洲一区二区三区在线视频| 四虎永久在线精品国产免费| 三上悠亚在线精品二区| 婷婷色中文网| 国产乱子精品一区二区在线观看| 国产视频自拍一区| AⅤ色综合久久天堂AV色综合| 欧美在线天堂| 久久无码免费束人妻| 国产XXXX做受性欧美88| 98超碰在线观看| 毛片免费高清免费| 亚洲第一精品福利| A级毛片无码久久精品免费| 国产91熟女高潮一区二区|