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

改性雙基推進劑剪切壓延過程的數值模擬①

2020-09-05 01:28:04薛曉俊潘浩東尹艷華
固體火箭技術 2020年4期
關鍵詞:改性

英 飛,薛曉俊,潘浩東,尹艷華

(北京理工大學 化學與化工學院,北京 102488)

0 引言

近代推進劑的制造始于19世紀80年代,主要使用擠壓成型工藝。擠壓成型工藝制造的推進劑致密性好,質量均勻、重現性好,生產周期短;但生產過程間斷,只能制造小尺寸的藥柱[1-3]。第二次世界大戰期間,澆鑄工藝逐漸發展起來。澆鑄工藝可以制造大型及復雜形狀的藥柱,但該工藝存在藥柱機械性能較差、燃燒穩定性較差、固化時間長等問題[4-5]。20世紀70年代,美國雷德福兵工廠使用單螺桿擠出工藝完成了單、雙、三基藥的大規模批量生產,但是該工藝不能適用于高填料的推進劑配方,并且加工過程的能耗較高[6-7]。1982年,德國代拿買諾貝爾公司使用雙螺桿擠出機完成了雙基推進劑的制造。雙螺桿擠出工藝可以實現自動化、連續化的供料與出料,提高了生產效率以及生產過程中的安全性,并且具有廢液、廢氣排放量少,混合效率高的優點[8-10]。壓延工藝借助于壓延機兩個輥筒旋轉產生的剪切擠壓力,配以相應的加工溫度,使粘彈性物料多次受到剪切擠壓和延展作用,完成塑化過程,提高材料的可塑性[11]。

剪切壓延機基于壓延機和雙螺桿擠出機,摒除它們的缺點,結合并保留了它們的優點。剪切壓延機包括兩個反向旋轉的輥筒,工作輥(U型輥)和自由輥(V型輥),輥筒表面置有不同形狀、不同數量、不同螺旋角的凹槽,兩輥在不同溫度、不同轉速下運行。改性雙基推進劑吸收藥料通過加料系統均勻、連續地進入剪切壓延機兩個輥筒之間,藥料包裹在工作輥上,沿輥筒軸向向前輸送,藥料在兩個輥筒之間受到加熱作用和剪切擠壓作用,進行驅水干燥、混合和塑化,并得以壓實、平展,由疏松態變為結構致密的塑化態,在出料口經造粒機切割成藥料豆粒。

剪切壓延機加工改性雙基推進劑的生產路線較短,可以實現自動化、連續化生產,并且大大提高了混合效率和塑化性能。但是目前對于剪切壓延機結構參數和工藝參數的選擇仍然憑借加工經驗來確定,具有較大的盲目性,并且耗費大量的時間和資金。因此,本文將計算流體動力學方法(CFD)方法引入到對改性雙基推進劑剪切壓延過程的研究中,利用合適的CFD軟件包對該過程進行數值模擬,分析該過程的流動規律及其影響因素,為剪切壓延機的優化提供參考。

自1970 年以來,人們開始應用CFD軟件對高分子聚合物加工過程的流動規律進行數值模擬研究,得到流動過程中各物理量場(如溫度場、壓力場、剪切速率場、粘度場等)的數值及圖像結果,幫助分析研究系統中存在的問題并提出優化。Zhong Tingting等使用POLYFLOW軟件模擬GR-35雙基推進劑的螺旋擠壓過程,結果表明,溫度、壓力、剪切速率在螺旋推進器處具有最大值,在增大擠出出口產量的同時必須注意安全性[12-13]。馬忠亮等使用POLYFLOW軟件分別對溶劑與硝化棉的體積比為1.35∶1和1.50∶1的單基推進劑在不同工藝條件下的全流道流動情況進行數值模擬。結果表明,變燃速單基推進劑交界面半徑的波動值隨螺桿轉速波動的增加和入口壓力波動的增加而增大,皮料選擇溶劑/硝化棉體積比較大的藥料、芯料選擇溶劑/硝化棉體積比較小的藥料有利于提高變燃速單基推進劑連續加工過程中的尺寸一致性[14-15]。He Zhongqi等使用POLYFLOW軟件模擬研究固體推進劑代料單螺桿擠出過程中的流變參數的分布,結果表明,螺桿轉速的增加會導致流體溫度和壓力的增加[16]。Zhou Ke等使用POLYFLOW軟件對推進劑的單螺桿擠出過程進行建模和分析,結果表明,推進劑藥料的粘度隨螺桿的幾何特性出現周期性波動,呈現"剪切變稀"的特性,在加工過程中提高螺桿轉速,有利于提高推進劑藥料的混合效率和塑化質量,但是會導致壓力和溫度的升高[17]。彭昭宇等應用POLYFLOW軟件對典型含能材料的壓延塑化過程進行數值模擬,探究速度場、壓力場、剪切速率場、粘度場和混合指數場等特征量與塑化效果之間的關系[18]。趙玉蓮等使用POLYFLOW軟件對嚙合同向雙螺桿擠出機進行非等溫數值模擬,對比分析不同轉速對加工效果的影響[19]。

1 剪切壓延機幾何模型的建立

本文使用Pro/Engineer軟件建立剪切壓延機的幾何模型。

1.1 U型輥幾何模型的建立

根據剪切壓延機U型輥的直徑、長度、槽寬、槽深和槽數等參數值,建立輥筒的幾何模型(圖1)。從圖1中可看出,U型輥的凹槽形狀為矩形,凹槽方向遵循右手螺旋定則。

圖1 U型輥表面的凹槽結構Fig.1 Grooves on the surface of U-shaped roll

1.2 V型輥幾何模型的建立

V型輥幾何模型的建立方法和U型輥的基本一致,區別在于V型輥表面凹槽的數量和形狀(圖2)。從圖中2可看出,V型輥的凹槽形狀為直角三角形,凹槽方向遵循左手螺旋定則。

圖2 V型輥表面的凹槽結構Fig.2 Grooves on the surface of V-shaped roll

2 剪切壓延機幾何模型的有限元網格劃分

剪切壓延機體積較大,其輥筒長度范圍為500~900 mm,輥筒半徑范圍為350~400 mm;同時兩個輥筒表面開鑿出的大量螺旋上升的不同形狀的凹槽大大增加了剪切壓延機幾何結構的復雜性,使得許多常用的網格劃分軟件與方法無法適用于剪切壓延機的幾何模型。本文通過大量研究,找到了適用于該模型的網格劃分方法,并在保證模型準確性的前提下,大大減少了網格數量、提高了網格精度,以降低模型計算成本,提高模擬計算分析的質量和效率。

本文使用HyperMesh軟件對剪切壓延機幾何模型進行有限元網格劃分。

2.1 U型輥幾何模型的網格劃分

2.1.1 幾何優化

將在Pro/Engineer軟件中建立的U型輥幾何模型導入到HyperMesh軟件中。剪切壓延機的兩個輥筒在改性雙基推進劑剪切壓延過程中作為剛性整體,進行旋轉運動對推進劑藥料產生剪切擠壓作用力并與推進劑藥料進行熱量傳遞使藥料溫度升高。因此,在保證模型準確性的前提下,可適當切割一部分輥筒的實心結構,以減少不必要的網格數量,從而降低模型計算成本,提高模擬計算分析的質量和效率。使用HyperMesh軟件中幾何處理功能中的boolean運算進行幾何體的切割(圖3)。

2.1.2 2D網格的劃分

從圖3中可看出,U型輥的上表面并不是一個規則的平面圖形,如果直接進行2D網格劃分,會導致網格不均勻,進而降低3D網格的質量,不利于后續的求解計算。因此,需要先將該平面分割成若干規則的小四邊形,然后再對各個小四邊形進行2D網格劃分(圖4),從圖中可看出,U型輥的2D網格均為四邊形元素。

圖3 U型輥幾何切割圖Fig.3 Geometric cutting of U-shaped roll

圖4 U型輥的2D網格Fig.4 2D mesh of U-shaped roll

2.1.3 3D網格的劃分

在進行3D網格劃分之前需要先沿垂直軸向將輥筒進行分段,然后再分別對各段輥筒進行3D網格劃分。

在第一段輥筒的外表面上任選一條凹槽棱邊,生成若干個實心圓點(圖5(a));將最下方的點復制平移到輥筒上表面(圖5(b));按照綠點→藍點→紅點的順序測量上表面上的三個點形成的銳角角度值(5(c))。

(a)Creating nodes (b)Projecting the nodes (c)Selecting the nodes圖5 測定旋轉角度Fig.5 Mearsuing the rotation angle

使用軟件的“translate”功能,將2.1.2中劃分的2D網格復制平移到第一段輥筒的下表面,使用“rotate”功能對下表面的2D網格進行旋轉,旋轉角度為圖5中測定的角度值。然后使用軟件的“solid map”中的“line drag”方法,利用生成的2D網格和點軌跡,進行第一段輥筒的3D網格劃分(圖6)。

圖6 U型輥第一段的3D網格Fig.6 3D mesh of the first segment of U-shaped roll

重復上述步驟,直至完成整個U型輥的3D網格劃分(圖7)。從圖中7可看出,U型輥的3D網格均為六面體,六面體網格可以大大減少模型的網格數量,并且提高求解過程的穩定性和準確性。

圖7 U型輥的3D網格Fig.7 3D mesh of U-shaped roll

2.2 V型輥幾何模型的網格劃分

V型輥的網格劃分方法與U型輥的基本一致。區別是V型輥表面的凹槽形狀為直角三角形,因此其2D網格除了四邊形元素以外,還存在三角形元素(圖8)。

圖8 V型輥的2D網格Fig.8 2D mesh of V-shaped roll

為對剪切壓延機的結構進行優化,需要進行大量不同結構參數的剪切壓延機加工改性雙基推進劑過程的模擬計算,通過求解結果的分析與比較,選擇較優的結構,為實際設備優化提供參考。1.1節和2.1節中概述的幾何建模和網格劃分方法可以應用到各種不同結構的剪切壓延機中,并可為其他形狀復雜的設備的有限元分析提供參考。

3 改性雙基推進劑剪切壓延過程模擬計算

3.1 改性雙基推進劑剪切壓延過程的數值模型

使用1.1節和2.1節中概述的方法,建立剪切壓延機和推進劑的幾何模型,并對其進行有限元網格劃分(圖9)。然后,將其導入到POLYFLOW軟件中進行裝配(圖10)。剪切壓延機的U型輥和V型輥關于Y軸對稱,輥筒直徑為350 mm,輥筒長度為500 mm,出口邊界位于XZ(Y=0)平面。

(a)U-shaped roll (b)V-shaped roll (c)Fluid圖9 網格模型Fig.9 3D mesh

圖10 U型輥、V型輥、物料層模型裝配Fig.10 Assembly of U-shaped roll,V-shaped roll and fluid

改性雙基推進劑藥料的物性參數和剪切壓延機的設備參數、工藝參數分別如表1~表3所示。

表1 藥料物性參數Table 1 Properties of propellant

表2 剪切壓延機設備參數Table 2 Parameters of shear calender

表3 剪切壓延機工藝參數Table 3 Parameters of shear calender

物料層的邊界設置如圖11所示。物料層的邊界條件設置如表4所示。剪切壓延機為開放式結構,因此物料層的外層邊界和出口邊界的流動邊界條件均設置為大氣壓的壓力條件,外層邊界的熱邊界條件設置為與空氣的對流傳熱條件。

圖11 物料層邊界設置Fig.11 Boundaries of fluid

表4 物料層邊界條件設置Table 4 Boundary conditions of fluid

3.2 改性雙基推進劑剪切壓延過程的模擬結果及分析

使用POLYFLOW軟件對數值模型進行求解計算,根據得到的計算結果分析改性雙基推進劑藥料剪切壓延過程的流動規律及其影響因素。

3.2.1 改性雙基推進劑剪切壓延過程的速度分析

速度分布對于推進劑藥料的成型加工質量具有直接的影響,如果速度分布差異較大,藥料容易發生局部變形,嚴重影響產品的使用性能,高質量成型要求的藥料必須保證在加工過程中具有較好的速度均勻性。改性雙基推進劑剪切壓延過程的XZ平面和XY平面速度分布云圖分別如圖12(a)、(b)所示,輥隙區域速度分布矢量圖如圖12(c)所示,圖12中的灰色區域表示剪切壓延機的兩個輥筒。

(a)Contour of velocities on XZ plane (b)Contour of velocities on XY plane (c)Vector ofvelocities in nip region圖12 速度分布示意圖Fig.12 Velocity distribution

由圖12(a)、(b)可知,XZ平面的速度關于輥隙區域呈近似對稱分布;藥料在進入輥隙區域時,由于受到兩輥筒旋轉產生的強大剪切力和擠壓力,運動阻力驟然增大,速度有所減小;進入輥隙區域后,在兩輥筒的剪切擠壓作用下,進行強制輸送,流速明顯增大。輥隙區域內靠近輥筒表面的藥料速度較小,中心部分速度較大,這是因為與輥筒表面接觸的藥料的速度幾乎與輥筒旋轉產生的速度相同,該速度與藥料強制輸送的速度相比較小。由圖12(c)可知,輥隙區域的速度方向與輥筒旋轉方向一致,并且具有較好的均勻性,有利于提高產品成型質量。

從物料進口邊界處每隔20 mm建立一系列與Y軸垂直的平面,計算各平面的平均速度,圖13為藥料平均速度沿Y軸的變化。由圖13可知,由于進出口效應的影響,進出口區域的速度變化相對較大,忽略進出口效應之后,流動過程中的平均速度變化較平穩;由于藥料在兩輥筒的剪切擠壓作用下,進行強制輸送,流動過程中的平均速度隨Y軸距離的減小略有增大,即藥料速度沿輥筒輸送方向略有增大。

圖13 平均速度沿Y軸的變化Fig.13 The change of average velocity along Y-axis

3.2.2 改性雙基推進劑剪切壓延過程的壓力分析

改性雙基推進劑剪切壓延過程的XZ平面和XY平面壓力分布云圖如圖14所示,灰色區域表示剪切壓延機的兩個輥筒。從物料進口邊界處每隔20 mm建立一系列與Y軸垂直的平面,計算各平面的平均壓力,圖15為藥料平均壓力沿Y軸變化。

(a)On XZ plane

(b)On XY plane圖14 壓力分布云圖Fig.14 Contour of pressure

圖15 平均壓力沿Y軸的變化Fig.15 The change of average pressure along Y-axis

由于剪切壓延機為開放式設備,藥料在加工過程中所受到的壓力主要來自于兩輥筒旋轉產生的剪切擠壓作用,因此輥隙區域的壓力明顯高于其他區域。藥料在進入剪切壓延機時,壓力梯度較大,沿藥料的輸送方向,壓力梯度逐漸減小,加工過程中無壓力突變行為,藥料所受壓力變化較平穩,說明剪切壓延過程的安全性較高。

3.2.3 改性雙基推進劑剪切壓延過程剪切速率分析

藥料的剪切速率越大則粘性生熱越多,一方面使藥料溫度升高,另一方面導致熱量堆積,形成“熱點”,容易引起藥料的著火、燃爆,降低生產過程的安全性;同時,如果剪切速率分布差異較大,會影響藥料加工成型的質量。改性雙基推進劑剪切壓延過程的XZ平面和XY平面剪切速率分布云圖如圖16所示,圖中的灰色區域表示剪切壓延機的兩個輥筒。從物料進口邊界處每隔20 mm建立一系列與Y軸垂直的平面,計算各平面的最大剪切速率,圖17為藥料最大剪切速率沿Y軸的變化。

(a)On XZ plane

(b)On XY plane圖16 剪切速率分布云圖Fig.16 Contour of shear rate

圖17 最大剪切速率沿Y軸的變化Fig.17 The change of maximum shear rate along Y-axis

由圖16可知,由于藥料在輥隙區域受到兩輥筒旋轉產生的剪切擠壓作用,該區域的剪切速率明顯高于其他區域。由圖17可知,由于進出口效應的影響,進出口區域剪切速率波動相對較大,忽略進出口效應之后,流動過程中的最大剪切速率在25~30 s-1范圍內波動,變化較平穩。

3.2.4 改性雙基推進劑剪切壓延過程的粘性熱分析

改性雙基推進劑剪切壓延過程的XZ平面和XY平面粘性熱分布云圖如圖18所示,圖中的灰色區域表示剪切壓延機的兩個輥筒。

由圖18(a)可知,粘性生熱的分布與剪切速率的分布基本一致,輥隙區域的粘性生熱明顯高于其他區域。

(a)On XZ plane

(b)On XY plane圖18 粘性熱分布云圖Fig.18 Contour of viscous heating

3.2.5 改性雙基推進劑剪切壓延過程的溫度分析

改性雙基推進劑配方中含有大量炸藥及金屬粉,能量高、感度高,如果在剪切壓延過程設備和工藝參數選擇不當,造成高溫點的出現,極易發生著火、燃爆,造成人員傷亡與財產損失。改性雙基推進劑藥料在剪切壓延過程中溫度升高主要有兩個源項,一是剪切壓延機輥筒的加熱作用,兩個輥筒的溫度通常設置為90~100 ℃,藥料在室溫溫度下進入剪切壓延機中,輥筒與藥料之間進行熱量傳遞,使藥料溫度升高;二是粘性熱的作用,改性雙基推進劑藥料為高粘性流體,在剪切壓延過程抵抗變形的粘性力做功,把流體運動的機械能不可逆地轉換為熱能,造成藥料溫度升高。

改性雙基推進劑剪切壓延過程的XZ平面和XY平面溫度分布云圖如圖19所示,圖中的灰色區域表示剪切壓延機的兩個輥筒。從物料進口邊界處每隔20 mm建立一系列與Y軸垂直的平面,計算各平面的最高溫度,圖20為藥料最高溫度沿Y軸的變化。

(a)On XZ plane

(b)On XY plane圖19 溫度分布云圖Fig.19 Contour of temperature

圖20 最高溫度沿Y軸的變化Fig.20 The change of maximum temperature along Y -axis

由圖19(a)可知,XZ平面的溫度分布呈圓環狀,圓環狀溫度分布從輥筒表面向外呈現先增大、后減小的趨勢。靠近輥筒表面的藥料與輥筒溫度基本一致,并且由于V型輥的溫度低于U型輥,靠近V型輥表面的藥料的溫度低于靠近U型輥表面的藥料的溫度。中心區域由于粘性熱的作用,藥料溫度逐漸有所升高。物料層外層邊界與空氣直接接觸,進行對流傳熱,熱量有所散失,因此溫度逐漸下降。除進出口區域由于進出口效應,溫度有相對較大的波動之外,剪切壓延過程的最高溫度變化較平穩,在380~385 K范圍內波動。

3.2.6 改性雙基推進劑剪切壓延過程的混合指數分析

為了提高產品質量,改性雙基推進劑藥料要在剪切壓延過程中進行充分混合。改性雙基推進劑剪切壓延過程的XZ平面和XY平面混合指數分布云圖如圖21所示,圖中的灰色區域表示剪切壓延機的兩個輥筒。從物料進口邊界處每隔20 mm建立一系列與Y軸垂直的平面,計算各平面的最大混合指數,圖22為藥料最大混合指數沿Y軸的變化。

(a)On XZ plane

(b)On XY plane圖21 混合指數分布云圖Fig.21 Contour of mixing index

由圖22可知,在進入和離開剪切壓延機輥隙區域時,藥料受到兩輥筒旋轉產生的剪切擠壓作用較弱,藥料的混合程度較低;藥料在剪切壓延機輥隙區域受到兩個輥筒旋轉產生的強大的剪切擠壓作用而進行混合,混合指數明顯高于其他區域;混合指數整體較高且變化平穩,大致穩定在0.88附近,說明藥料在兩輥筒之間受力均勻且穩定,有利于提高產品成型質量。

圖22 最大混合指數沿Y軸的變化Fig.22 Change of maximum mixing index along Y-axis

4 結論

考慮到剪切壓延機幾何結構和推進劑流變性質的復雜性,以及實際試驗修正的時間和資金成本,本文將CFD方法引入到改性雙基推進劑剪切壓延過程的研究中,利用合適的CFD軟件包對該過程進行數值模擬,分析該過程的流動規律及其影響因素。使用Pro/Engineer軟件建立剪切壓延機的三維幾何模型,使用HyperMesh軟件對幾何模型進行有限元網格劃分,使用POLYFLOW軟件建立改性雙基推進劑剪切壓延過程的數值模型并進行模擬計算,分析該過程的速度、壓力、剪切速率、粘性生熱、溫度、混合指數等參數的分布規律:

(1)改性雙基推進劑藥料在剪切壓延機的輥隙區域受到兩個輥筒旋轉產生的強大的剪切擠壓作用,該區域的壓力、剪切速率、粘性熱、混合指數明顯高于其他區域;

(2)藥料在剪切壓延過程中受到兩輥筒旋轉產生的剪切擠壓作用而進行混合,最大混合指數可達0.88左右;

(3)藥料的速度分布和剪切速率分布的均勻性較好,有利于提高產品成型質量,忽略進出口效應后,平均速度在0.376 m/s附近波動,最大剪切速率在25~30 s-1范圍內波動;

(4)該過程中溫度和壓力的變化較平穩,沒有出現突變,工藝的安全性較高。

后續將建立大量不同結構參數的剪切壓延機加工改性雙基推進劑過程的數值模型,并進行模擬計算,通過求解結果的分析與對比,選擇較優的結構,為實際工藝設備的改進提供參考。

猜你喜歡
改性
尼龍6的改性研究進展
P(3,4HB)/PHBV共混改性及微生物降解研究
中國塑料(2016年12期)2016-06-15 20:30:07
汽車用ABS/PC共混物改性研究進展
中國塑料(2016年2期)2016-06-15 20:30:00
有機磷改性納米SiO2及其在PP中的應用
中國塑料(2016年2期)2016-06-15 20:29:59
我國改性塑料行業“十二·五”回顧與“十三·五”展望
中國塑料(2016年5期)2016-04-16 05:25:36
三聚氰胺甲醛樹脂的改性研究進展
聚乳酸的阻燃改性研究進展
中國塑料(2015年3期)2015-11-27 03:41:38
ABS/改性高嶺土復合材料的制備與表征
中國塑料(2015年11期)2015-10-14 01:14:14
聚甲醛增強改性研究進展
中國塑料(2015年9期)2015-10-14 01:12:17
聚乳酸擴鏈改性及其擠出發泡的研究
中國塑料(2015年4期)2015-10-14 01:09:19
主站蜘蛛池模板: 伊人国产无码高清视频| 极品av一区二区| 亚洲人成人无码www| 国产簧片免费在线播放| аv天堂最新中文在线| 国产精品夜夜嗨视频免费视频| 久久国产亚洲欧美日韩精品| 欧美三级不卡在线观看视频| 国产精品污污在线观看网站| 992tv国产人成在线观看| 国产精品 欧美激情 在线播放| 精品撒尿视频一区二区三区| 热久久国产| 精品无码一区二区三区在线视频| 日本午夜网站| 欧美国产日韩在线| 91在线无码精品秘九色APP| 国产福利影院在线观看| 国产成人久视频免费| 91在线精品免费免费播放| 久久semm亚洲国产| 亚洲国模精品一区| 欧美国产在线精品17p| 看国产毛片| 亚洲中文字幕在线一区播放| 久久男人视频| 亚洲大学生视频在线播放| 91在线日韩在线播放| 国产精品区视频中文字幕| 日本成人不卡视频| 久久青草免费91线频观看不卡| 国产资源站| 狠狠五月天中文字幕| 中国国产一级毛片| 亚洲欧美极品| av一区二区三区在线观看| 久久中文字幕2021精品| 91久久精品国产| 亚洲免费福利视频| 国产永久在线视频| 亚洲综合经典在线一区二区| 美女无遮挡拍拍拍免费视频| 综合网久久| 中文字幕伦视频| 色综合中文| 久久狠狠色噜噜狠狠狠狠97视色 | 亚洲国产精品无码久久一线| 国内熟女少妇一线天| 国产欧美成人不卡视频| 亚洲欧美自拍中文| 亚洲日韩第九十九页| 亚洲人成色在线观看| 国产美女91视频| 女人18一级毛片免费观看| 精品国产一区二区三区在线观看 | 亚洲精品成人片在线观看| 亚洲无码高清视频在线观看| 在线看国产精品| 最新国产你懂的在线网址| 91小视频在线| 超清无码熟妇人妻AV在线绿巨人| 在线99视频| 无码中文字幕精品推荐| AV不卡在线永久免费观看| 国产91丝袜| 欧美五月婷婷| 色综合a怡红院怡红院首页| 91久久大香线蕉| 国产99在线观看| 亚洲码在线中文在线观看| 91久久天天躁狠狠躁夜夜| 国产欧美日韩视频怡春院| 国产国模一区二区三区四区| 永久免费精品视频| 亚洲欧洲日产国产无码AV| 99视频精品全国免费品| www.狠狠| 国产午夜一级毛片| 亚洲精品不卡午夜精品| 六月婷婷精品视频在线观看| 成人午夜久久| 波多野结衣爽到高潮漏水大喷|