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

HTPB推進(jìn)劑寬泛應(yīng)變率下粘彈性本構(gòu)模型研究*

2017-06-19 19:09:41周海霞李世鵬
固體火箭技術(shù) 2017年3期
關(guān)鍵詞:實(shí)驗(yàn)模型

周海霞,李世鵬,謝 侃,隋 欣,楊 龍

(北京理工大學(xué) 宇航學(xué)院,北京 100081)

HTPB推進(jìn)劑寬泛應(yīng)變率下粘彈性本構(gòu)模型研究*

周海霞,李世鵬,謝 侃,隋 欣,楊 龍

(北京理工大學(xué) 宇航學(xué)院,北京 100081)

為研究HTPB推進(jìn)劑的率相關(guān)性力學(xué)性能,采用材料萬(wàn)能試驗(yàn)機(jī)、液壓試驗(yàn)機(jī)和分離式霍普金森壓桿(SHPB),分別開(kāi)展了低(1.67×10-4~1.67×10-1s-1)、中(1~100 s-1)、高(700~2 500 s-1)應(yīng)變率的單軸壓縮實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明,HTPB推進(jìn)劑的壓縮力學(xué)性能是率相關(guān)性,隨應(yīng)變率的升高,給定應(yīng)變下的應(yīng)力逐漸增大。采用廣義非線性ZWT本構(gòu)模型描述HTPB推進(jìn)劑寬泛應(yīng)變率下的壓縮力學(xué)行為,模型預(yù)測(cè)與實(shí)驗(yàn)數(shù)據(jù)對(duì)比表明,模型中至少需要4個(gè)麥克斯韋元件。

HTPB推進(jìn)劑;壓縮;應(yīng)變率;應(yīng)變率相關(guān);本構(gòu)模型

0 引言

HTPB推進(jìn)劑是一種高含能的聚合物材料[1],其力學(xué)性質(zhì)具有應(yīng)變率相關(guān)性[2]。在其生命周期內(nèi),不僅會(huì)受到藥柱制造、裝配和存儲(chǔ)過(guò)程中的自身重量等低應(yīng)變率載荷作用;而且會(huì)受到高應(yīng)變率載荷作用,如運(yùn)輸過(guò)程中的震動(dòng)沖擊、發(fā)射過(guò)載以及點(diǎn)火瞬間壓力沖擊等沖擊載荷。因此,在設(shè)計(jì)過(guò)程中,考慮和保證推進(jìn)劑裝藥結(jié)構(gòu)的完整性成為火箭武器發(fā)展的關(guān)鍵技術(shù)。從裝藥結(jié)構(gòu)完整性分析中,需要建立推進(jìn)劑相對(duì)應(yīng)的粘彈性本構(gòu)方程,分析應(yīng)力、應(yīng)變及材料特性參數(shù)間的關(guān)系,為數(shù)值模擬提供準(zhǔn)確的本構(gòu)模型[3]。因此,建立準(zhǔn)確的推進(jìn)劑本構(gòu)模型是優(yōu)化火箭發(fā)動(dòng)機(jī)裝藥結(jié)構(gòu)設(shè)計(jì)的基礎(chǔ)。國(guó)內(nèi)外提出的相關(guān)模型有非線性ZWT模型[4-7]、廣義Maxwell模型[8]、基于Yeoh函數(shù)超彈模型建立的一種粘彈本構(gòu)模型[9]等。目前,HTPB的力學(xué)特性的研究主要集中在低應(yīng)變率及高應(yīng)變率方面,而對(duì)固體推進(jìn)劑在中應(yīng)變率下的力學(xué)特性研究報(bào)道較少。

本文基于實(shí)驗(yàn)結(jié)果建立了包括中應(yīng)變率在內(nèi)的寬泛應(yīng)變率下的粘彈性本構(gòu)模型,來(lái)描述HTPB推進(jìn)劑在寬應(yīng)變率范圍內(nèi)的力學(xué)特性,以便為發(fā)動(dòng)機(jī)設(shè)計(jì)及藥柱完整性研究提供理論支持。

1 實(shí)驗(yàn)

1.1 材料和試樣

表1列出了實(shí)驗(yàn)常用的HTPB推進(jìn)劑成分,在低中應(yīng)變率壓縮實(shí)驗(yàn)中,為消除尺寸變化給實(shí)驗(yàn)帶來(lái)的誤差,將試件尺寸保持一致,實(shí)驗(yàn)采用圓柱形試件,其基本尺寸為φ16 mm×20 mm;在高應(yīng)變率壓縮實(shí)驗(yàn)中,試件尺寸為φ10 mm×5 mm,長(zhǎng)徑比為1∶2,這種較短的長(zhǎng)度可使所記錄的應(yīng)變信號(hào)中波的彌散及試樣的徑向和軸向慣性效應(yīng)均達(dá)到最小。所有試樣由同一批次加工,實(shí)驗(yàn)前儲(chǔ)存在干燥罐中以消除濕度變化對(duì)材料壓縮性能的影響。

表1 HTPB推進(jìn)劑基礎(chǔ)配方

1.2 低中應(yīng)變率壓縮實(shí)驗(yàn)

HTPB推進(jìn)劑低中應(yīng)變速率壓縮實(shí)驗(yàn)條件見(jiàn)表2。每組試驗(yàn)至少重復(fù)5次,取重復(fù)性較好的實(shí)驗(yàn)數(shù)據(jù)的平均值作為實(shí)驗(yàn)結(jié)果。用載荷傳感器和延伸計(jì)分別測(cè)力和形變量,再計(jì)算得到相應(yīng)的應(yīng)力-應(yīng)變曲線。

表2 低中應(yīng)變率壓縮下的實(shí)驗(yàn)條件

1.3 高應(yīng)變率壓縮實(shí)驗(yàn)

HTPB推進(jìn)劑的高應(yīng)變率壓縮實(shí)驗(yàn)在SHPB裝置上完成,實(shí)驗(yàn)工況見(jiàn)表3,每組實(shí)驗(yàn)至少重復(fù)5次。圖1是SHPB裝置的示意圖。實(shí)驗(yàn)時(shí),將試樣放在入射桿和透射桿之間,氣槍發(fā)射子彈,在給定速度下撞擊入射桿前端面,形成入射應(yīng)變脈沖。當(dāng)入射應(yīng)變脈沖傳播至入射桿和試樣的接觸界面時(shí),一部分被反射回入射桿形成反射波,另一部分穿過(guò)試件透射入透射桿。固定在桿上的應(yīng)變片在10 MHz的采樣頻率和500倍的放大倍數(shù)下記錄入射、反射和透射應(yīng)變脈沖。根據(jù)一維應(yīng)力波理論,可計(jì)算出相應(yīng)的應(yīng)力和應(yīng)變。由于實(shí)驗(yàn)時(shí)試樣可能意外點(diǎn)燃,保護(hù)盒可起保護(hù)作用,保證實(shí)驗(yàn)人員的安全性。為保證低阻抗的HTPB推進(jìn)劑能獲得有效的測(cè)試數(shù)據(jù),本文所用的SHPB裝置相對(duì)傳統(tǒng)的裝置有3項(xiàng)改進(jìn):

(1)在桿和試樣的交界面處涂上凡士林,以使摩擦達(dá)到最小;

(2)用高強(qiáng)度的鋁合金代替鋼作為桿材料,且透射桿采用空心結(jié)構(gòu),增加了透射應(yīng)變脈沖的幅值;

(3)用銅片做脈沖整形器,確保試樣加載的恒應(yīng)變率和動(dòng)態(tài)應(yīng)力平衡。

表3 高應(yīng)變率壓縮實(shí)驗(yàn)的實(shí)驗(yàn)條件

圖2給出了HTPB推進(jìn)劑典型的入射波、反射波和透射波信號(hào)圖。由于使用了脈沖整形器,過(guò)濾掉了應(yīng)力波在傳播過(guò)程中由于波形彌散產(chǎn)生的高頻振蕩,促使試件內(nèi)應(yīng)力更快均勻,波形更加平滑。反射信號(hào)中近似平臺(tái)的區(qū)域表明,試樣在恒定的工程應(yīng)變率下發(fā)生變形。

為考察試樣動(dòng)態(tài)應(yīng)力平衡性,常用的一種方法是對(duì)比其前端面的力F1和后端面的力F2,再基于一維應(yīng)力波理論,計(jì)算F1和F2。這種方法對(duì)本實(shí)驗(yàn)不適用,因?yàn)镠TPB推進(jìn)劑的阻抗較低,使得入射信號(hào)和反射信號(hào)的幅值比較接近,所以會(huì)產(chǎn)生較大的誤差[10]。在本次研究中,F(xiàn)I、FR和FT分別由入射、反射和透射信號(hào)計(jì)算得出,通過(guò)對(duì)比FR和FT-FI的重合度來(lái)評(píng)價(jià)試樣的動(dòng)態(tài)應(yīng)力平衡性。如圖3所示,與FR和FT-FI相關(guān)的2個(gè)力的時(shí)間曲線幾乎重合,表明試樣在恒定的工程應(yīng)變率及動(dòng)態(tài)應(yīng)力平衡下發(fā)生變形。所以,高應(yīng)變率壓縮實(shí)驗(yàn)數(shù)據(jù)是有效的。

1.4 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)所得數(shù)據(jù)經(jīng)過(guò)處理后,得到如圖4所示的應(yīng)力-應(yīng)變曲線。

從圖4可看出,相同應(yīng)變下,應(yīng)力隨應(yīng)變率的升高而增大。另外,隨應(yīng)變率的升高,其初始模量也在不斷增大,因此HTPB推進(jìn)劑具有明顯的率相關(guān)性;推進(jìn)劑在變形過(guò)程中沒(méi)有明顯的屈服現(xiàn)象,可發(fā)生的變形量較大;隨應(yīng)變?cè)龃螅瑧?yīng)力梯度逐漸減小,表明松弛現(xiàn)象越來(lái)越明顯。

圖5給出了不同應(yīng)變下應(yīng)力與應(yīng)變率對(duì)數(shù)的關(guān)系,其中數(shù)據(jù)點(diǎn)由圖4所示的曲線直接獲得。由圖5可知,低應(yīng)變率下,應(yīng)力隨應(yīng)變率對(duì)數(shù)線性增加;中應(yīng)變率下,應(yīng)力隨應(yīng)變率對(duì)數(shù)以指數(shù)方式增加;高應(yīng)變率下,應(yīng)力隨應(yīng)變率對(duì)數(shù)近似線性增加。圖5可定性地看出,HTPB在中應(yīng)變率下的力學(xué)表述與低高應(yīng)變率下的力學(xué)表述是不同的

如果只用中高應(yīng)變率實(shí)驗(yàn)得出的非線性ZWT本構(gòu)模型去預(yù)測(cè)中應(yīng)變率下的力學(xué)特性,有可能是不正確的。因此,下文中會(huì)用得到的寬應(yīng)變率范圍的實(shí)驗(yàn)數(shù)據(jù),去檢驗(yàn)典型非線性ZWT本構(gòu)模型對(duì)中應(yīng)變率力學(xué)性能的預(yù)測(cè)結(jié)果是否可靠。

2 粘彈性本構(gòu)模型的建立

2.1 本構(gòu)模型建立

典型的非線性ZWT模型是由1個(gè)非線性彈簧和2個(gè)分別用于描述低、高應(yīng)變率粘彈特性的Maxwell元件并聯(lián)的本構(gòu)模型。粘彈性力學(xué)中的廣義Maxwell模型也是較常見(jiàn)的模型,是由n個(gè)Maxwell元件組成,具有普適性,可很好地描述松弛時(shí)間相差很大的高聚物。

結(jié)合上面2個(gè)模型,本文采用了圖6所示的模型。該模型由1個(gè)代表準(zhǔn)靜態(tài)力學(xué)特性的非線性彈簧和n個(gè)Mawell元件(n≥2)并聯(lián)組成,將其稱(chēng)為廣義非線性ZWT模型。

其本構(gòu)關(guān)系的積分形式為

(1)

恒應(yīng)變率時(shí),式(1)積分得

(2)

式(1)、式(2)中,等號(hào)右側(cè)前3項(xiàng)是非線性彈性平衡響應(yīng),即與應(yīng)變率無(wú)關(guān)的準(zhǔn)靜態(tài)力學(xué)特性,E0、α、β為非線性彈性常數(shù);后n項(xiàng)分別是用于描述低應(yīng)變率、中應(yīng)變率和高應(yīng)變率的粘彈性響應(yīng)。其中,當(dāng)n≥3時(shí),E1、θ1、Ei、θi(i=2,3,…,n-1)和En、θn分別為描述低、中、高應(yīng)變率粘彈響應(yīng)項(xiàng)的彈性常數(shù)和松弛時(shí)間。當(dāng)n=2時(shí),該模型就是典型的非線性ZWT模型,E1、θ1和E2、θ2分別為描述低、高應(yīng)變率粘彈響應(yīng)項(xiàng)的彈性常數(shù)和松弛時(shí)間。

2.2 本構(gòu)模型參數(shù)獲取

由于實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)應(yīng)變率小于1.67×10-4s-1的力學(xué)曲線與1.67×10-4s-1曲線基本重合,因此可將1.67×10-4s-1應(yīng)變率下的力學(xué)曲線看成是準(zhǔn)靜態(tài)曲線,該曲線可用方程式σ=E0ε+αε2+βε3擬合,從而得到E0、α、β3個(gè)參數(shù)。

在低應(yīng)變率和中應(yīng)變率條件下,由于加載時(shí)間遠(yuǎn)大于高應(yīng)變率的松弛時(shí)間,表示高應(yīng)變率的Maxwell流體模型會(huì)很快松弛,使其應(yīng)力值近似為0。所以,此時(shí)的應(yīng)力-應(yīng)變公式為

(3)

當(dāng)n=2時(shí),將低應(yīng)變率下的實(shí)驗(yàn)數(shù)據(jù)代入遺傳算法程序中,可求出低應(yīng)變率的彈性常數(shù)和松弛時(shí)間E1、θ1;當(dāng)n≥3時(shí),將低中應(yīng)變率下的實(shí)驗(yàn)數(shù)據(jù)代入遺傳算法程序中,可將低應(yīng)變率和中應(yīng)變率的彈性常數(shù)和松弛時(shí)間E1、θ1、Ei、θi(i=2,3,…,n-1)求出來(lái)。然后,利用最小二乘法擬合一條2 500 s-1的高應(yīng)變率曲線,即可求出高應(yīng)變率的彈性常數(shù)和松弛時(shí)間En、θn。麥克斯韋元件個(gè)數(shù)為2、3、4時(shí)的模型參數(shù)值如表4所示。2.3 本構(gòu)模型的驗(yàn)證

通過(guò)前面建立的粘彈性本構(gòu)模型以及獲取的模型參數(shù),得到了常溫下HTPB推進(jìn)劑的粘彈性本構(gòu)模型。將低、中、高應(yīng)變率實(shí)驗(yàn)數(shù)據(jù)代入推進(jìn)劑的本構(gòu)模型中,

獲得模型擬合結(jié)果。圖7為實(shí)驗(yàn)結(jié)果與不同麥克斯韋元件個(gè)數(shù)的模型結(jié)果的對(duì)比圖。

表4 粘彈性本構(gòu)模型參數(shù)值

由圖7(a)、(b)可知,典型的非線性ZWT本構(gòu)模型能較好地?cái)M合低、高應(yīng)變率的曲線,但對(duì)中應(yīng)變率力學(xué)性能的預(yù)測(cè)結(jié)果非常差。

由圖7(c)、(d)可知,在三單元本構(gòu)模型(即n=3)中,低中應(yīng)變率過(guò)渡段(即1.67×10-2s-1和1.67×10-1s-1)的曲線擬合效果較差,其他應(yīng)變率下的理論模型與實(shí)驗(yàn)結(jié)果符合程度都很好。說(shuō)明該模型的離散松弛時(shí)間譜較少,不能較好擬合中低應(yīng)變率過(guò)渡段下的力學(xué)特性。

由圖7(e)、(f)可知,實(shí)驗(yàn)數(shù)據(jù)與擬合曲線吻合得很好,說(shuō)明四單元本構(gòu)模型(即n=4)可很好描述出應(yīng)變?cè)?~0.3范圍內(nèi)的力學(xué)特性。因此,采用廣義非線性ZWT模型建立的粘彈性本構(gòu)模型能夠很好地預(yù)測(cè)HTPB推進(jìn)劑的率相關(guān)性。

HTPB的結(jié)構(gòu)復(fù)雜,運(yùn)動(dòng)單元大小不等,不同尺度結(jié)構(gòu)的運(yùn)動(dòng)單元有不同的粘性系數(shù)和彈性系數(shù),它們的松弛時(shí)間就不止一個(gè),而是形成一個(gè)范圍寬廣的連續(xù)譜,所以在很寬的應(yīng)變率范圍內(nèi),高聚物均呈現(xiàn)有粘彈性。但在實(shí)際工程應(yīng)用中,較難求得連續(xù)的松弛時(shí)間譜,所以可考慮用不連續(xù)的離散譜來(lái)近似計(jì)算,且Maxwell元件越多,預(yù)測(cè)越精確。

結(jié)合圖5可看出,與低、高應(yīng)變率不同,中應(yīng)變率下的應(yīng)力隨應(yīng)變率對(duì)數(shù)以指數(shù)方式增加。低、高應(yīng)變率粘彈特性分別用1個(gè)Maxwell元件表示,則中應(yīng)變率粘彈特性應(yīng)至少用2個(gè)Maxwell元件表示,由圖7(e)、(f)知,可很好地描述類(lèi)似的率相關(guān)性材料的力學(xué)特性的n的最小值為n=4。

3 結(jié)論

常溫下HTPB推進(jìn)劑的實(shí)驗(yàn)曲線表明,HTPB推進(jìn)劑具有明顯的率相關(guān)性和粘彈性變形特性。通過(guò)理論模型與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比發(fā)現(xiàn),典型非線性ZWT本構(gòu)模型(即n=2)預(yù)測(cè)的中應(yīng)變率力學(xué)特性與實(shí)驗(yàn)結(jié)果不符。廣義非線性ZWT本構(gòu)模型能夠很好地描述HTPB推進(jìn)劑在實(shí)驗(yàn)應(yīng)變范圍內(nèi)的單抽壓縮力學(xué)性能的n的最小值為4。

[1] Kakavas P A. Mechanical properties of propellant composite materials reinforced with ammonium perchlorate particles[J]. International Journal of Solids and Structures,2014,51(10): 2019-2026.

[2] Zalewski R,Wolszakiewicz T. Analysis of uniaxial tensile tests for homogeneous solid propellants under various loading conditions[J]. Central European Journal of Energetic Materials,2011,8(4): 223-231.

[3] Kalaycioglu B,Dirikolu M H,Celik V. An elasto-viscoplastic analysis of direct extrusion of a double base solid propellant[J]. Advances in Engineering Software,2010,41(9): 1110-1114.

[4] Chou S C,Robertson K D,Rainey J H. The effect of strain rate and heat developed during deformation on the stress-strain curve of plastics[J]. Experimental Mechanics,1973,13(10): 422-432.

[5] Siviour C R,Walley S M,Proud W G,et al. The high strain rate compressive behaviour of polycarbonate and polyvinylidene difluoride[J]. Polymer,2005,46(26): 12546-12555.

[6] 鄧凱,陽(yáng)建紅,陳飛,等. HTPB復(fù)合固體推進(jìn)劑本構(gòu)方程[J]. 宇航學(xué)報(bào),2010,31(7): 1815-1818.

[7] Walley S M,Field J E. Strain rate sensitivity of polymers in compression from low to high rates[J]. DYMAT J.,1994,1(3): 211-227.

[8] 胡少青,鞠玉濤,孟紅磊,等. 雙基推進(jìn)劑壓縮力學(xué)性能的應(yīng)變率相關(guān)性研究[J]. 彈道學(xué)報(bào),2011,23(4):75-79.

[9] 周相榮,王強(qiáng),王寶珍. 一種基于Yeoh函數(shù)的非線性粘超彈本構(gòu)模型及其在沖擊仿真中的應(yīng)用[J]. 振動(dòng)與沖擊,2007,26(5): 33-37.

[10] Kolsky H. An investigation of the mechanical properties of materials at very high rates of loading[J]. Proceedings of the Physical Society,Section B,1949,62(11): 676.

(編輯:劉紅利)

Research on the viscoelastic constitutive model of HTPB propellant over a wide range of strain rates

ZHOU Hai-xia,LI Shi-peng,XIE Kan,SUI Xin,YANG Long

(School of Aerospace Engineering,Beijing Institute of Technology,Beijing 100081,China)

To investigate the rate-dependent mechanical properties of hydroxyl-terminated polybutadiene (HTPB) propellant,low(1.67×10-4~1.67×10-1s-1), intermediate(1~100 s-1),and high strain rate (700~2 500 s-1) compressive testing was performed,by using a universal testing machine,a hydraulic testing machine,and a split Hopkinson pressure bar (SHPB) system,respectively. Results indicate that the compressive mechanical properties of HTPB propellant are rate-dependent and the stress at a given strain increases with an increase in strain rate.Generalized nonlinear ZWT model was employed to describe the compressive mechanical behaviours of HTPB propellant over a wide range of strain rates.Compasrison of the predicted values from the model and the experimental data shows that the number of Maxwell elements in the model should be at least 4.

HTPB propellant;compression;strain rate;strain-rate correlation;constitutive model

2016-01-19;

2016-04-12。

周海霞(1992—),女,博士生,主要從事高能推進(jìn)劑研究。E-mail:1354587242@qq.com

謝侃(1982—),男,副教授。E-mail:xiekan@bit.edu.cn

V512

A

1006-2793(2017)03-0325-05

10.7673/j.issn.1006-2793.2017.03.010

猜你喜歡
實(shí)驗(yàn)模型
一半模型
記一次有趣的實(shí)驗(yàn)
微型實(shí)驗(yàn)里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
做個(gè)怪怪長(zhǎng)實(shí)驗(yàn)
3D打印中的模型分割與打包
NO與NO2相互轉(zhuǎn)化實(shí)驗(yàn)的改進(jìn)
實(shí)踐十號(hào)上的19項(xiàng)實(shí)驗(yàn)
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 精品国产免费人成在线观看| 国内精品一区二区在线观看| 亚洲青涩在线| 色天天综合| 成人免费视频一区二区三区 | 色亚洲激情综合精品无码视频 | www.精品国产| 国产成人一区| 欧美久久网| 亚洲成在线观看 | 四虎在线观看视频高清无码 | 亚洲午夜天堂| 91精品国产麻豆国产自产在线| 久久亚洲欧美综合| 国产亚洲精品97在线观看| 国产成人永久免费视频| 成年A级毛片| 色综合中文字幕| 久久窝窝国产精品午夜看片| 无码内射中文字幕岛国片| 国产制服丝袜91在线| 激情在线网| 成年人福利视频| 国内精品一区二区在线观看| 亚洲日本中文字幕天堂网| 国产亚洲视频在线观看| 亚洲另类国产欧美一区二区| 亚洲国产精品人久久电影| 刘亦菲一区二区在线观看| 久久久久亚洲AV成人网站软件| 高清欧美性猛交XXXX黑人猛交| 亚洲视频无码| 亚洲日韩精品伊甸| 国产精品不卡永久免费| 国产在线观看一区精品| 国产成人高清亚洲一区久久| 亚洲bt欧美bt精品| 99精品免费欧美成人小视频 | 亚洲av日韩综合一区尤物| 青青青国产精品国产精品美女| 啪啪永久免费av| 日韩麻豆小视频| 成AV人片一区二区三区久久| 久久五月天综合| 国产真实乱了在线播放| 国产AV毛片| 精品国产免费观看一区| 国产一级特黄aa级特黄裸毛片| 日韩美毛片| 色有码无码视频| 欧美日韩午夜视频在线观看| 国内精自视频品线一二区| 在线日韩一区二区| 欧美曰批视频免费播放免费| 美女内射视频WWW网站午夜| 免费视频在线2021入口| 久久亚洲日本不卡一区二区| 欧美精品成人一区二区视频一| 免费人成视频在线观看网站| 美女扒开下面流白浆在线试听 | 日韩a在线观看免费观看| 91色在线视频| 婷婷六月在线| 国产日韩欧美黄色片免费观看| 亚洲国产午夜精华无码福利| 国产黑丝一区| 成人福利在线免费观看| 成年人国产视频| 国产精品自拍露脸视频| 国产喷水视频| 亚洲欧美h| 欧美日韩在线国产| 99精品免费欧美成人小视频 | 国产成人久久综合一区| 婷婷激情五月网| 四虎影视国产精品| 无码网站免费观看| 欧美一区二区三区欧美日韩亚洲| 国产av剧情无码精品色午夜| 国产在线视频欧美亚综合| 少妇露出福利视频| 欧美日韩中文字幕二区三区|