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

沖壓增程制導炮彈氣動特性研究

2020-04-16 01:38:16張寧史金光馬曄璇
兵工學報 2020年3期

張寧,史金光,馬曄璇

(南京理工大學 能源與動力工程學院,江蘇 南京 210094)

0 引言

固體燃料沖壓增程技術是實現炮彈增程的一種有效技術途徑,其增程率約達70%,近年來國內外發展了旋轉穩定的沖壓增程炮彈。隨著射程的增加,火炮射擊精度逐步下降,采用簡易控制技術可改善炮彈射擊精度,為此國內外逐漸開始對沖壓增程制導炮彈進行研究。

針對固體沖壓增程炮彈的氣動布局與氣動特性研究,目前主要采用風洞試驗以及數值仿真相結合的手段進行。沖壓增程炮彈按照穩定方式可分為旋轉穩定彈與尾翼穩定彈。針對旋轉穩定彈,文獻[1-5]在馬赫數Ma為4.3條件下對75 mm沖壓增程炮彈的氣動特性進行了分析。文獻[6]對40 mm沖壓增程高射炮彈的氣動特性與外形進行了分析。文獻[7]在二者分析的基礎上,發現75 mm和40 mm沖壓增程炮彈在內壁壓力及流動分離方面有顯著不同。文獻[8-10]在Ma數4.0條件下對155 mm沖壓增程炮彈進行了數值模擬和飛行試驗,并在Ma數為2.2~3.0情況下對其1/3比例模型進行了風洞實驗。文獻[11]對155 mm沖壓增程炮彈的激波位置和阻力大小進行了模擬。針對尾翼穩定彈,文獻[12]提出了一種155 mm沖壓增程彈的氣動構型。文獻[13]研究了沖壓增程炮彈的錐形激波與入口邊界層損失。相對于傳統常規炮彈,沖壓增程炮彈的頭部進氣道組件直接影響炮彈增程性能,文獻[14]認為冷態飛行試驗的氣動阻力可以直接作為沖壓發動機推力的設計依據。文獻[15]通過風洞實驗發現進氣道有效流通面積與總壓恢復系數呈反比。文獻[16-17]通過數值模擬和風洞實驗,發現雙錐混壓式進氣道性能較高。

對于固體沖壓增程彈的氣動特性模擬,計算網格的質量、離散格式的優劣將直接影響計算精度。對此,文獻[18]使用偏微分法(TTM)貼體網格、文獻[1]選用拼接網格技術與微分方程法、文獻[19]采用塊外形網格分區技術對沖壓增程彈進行了網格劃分。文獻[1]、文獻[19-20]應用2階隱式總變差不增(TVD)離散格式,在沖壓增程炮彈的復雜湍流流場數值仿真上取得了良好效果。沖壓增程炮彈與普通炮彈相比,差異在于由激波和不連續尖銳幾何形狀引起的轉捩流動區域,這些區域良好的模擬對整個流場計算質量尤為重要。文獻[21]認為γ-Reθt-SA(γ為間歇因子,Reθt為當地轉捩雷諾數,SA表示Spalart-Allmaras模型)轉捩模型能比較準確地預測自然轉捩和分離轉捩。文獻[22]對矩形微通道中的流動進行數值研究并與實驗結果進行了對比,認為γ-Reθt-SA模型能準確地捕捉到臨界雷諾數。由于γ-Reθt-SA轉捩模型中間隙因子輸運方程并不依賴于具體的湍流模型,文獻[23-24]認為在小攻角范圍內,γ-Reθt-SA模型的模擬結果與γ-Reθt-剪切應力輸運(SST)模型準確度相當,但計算耗時少。

目前國內外對固體燃料沖壓增程技術開展了大量研究,但是基本都以旋轉穩定炮彈為研究對象,尚未針對沖壓增程制導炮彈的氣動外形開展流場與氣動特性研究。本文依據沖壓增程制導炮彈的工作原理與飛行特點,將其彈道分為3段,即沖壓助推段、爬升飛行段、滑翔控制段[25],并依據不同彈道段特性建立3種氣動外形;通過拼接網格技術進行分區,采用2階迎風矢通量分裂格式與γ-Reθt-SA湍流模型,對3種氣動外形進行數值計算,分析對應的流場與升阻力特性。

1 氣動外形

本文所設計的沖壓增程制導炮彈結構參數為:彈丸總長約為11倍彈徑,采用4對尾翼穩定、2對鴨舵控制。按Ma數為2.5設計進氣道參數,采用雙錐混壓式布局,按照Oswatisch最優波系原理進行配置,為適應更寬的Ma數,設計第1級半錐角為21°,第2級半錐角為16°,唇罩前緣內壁與自由來流夾角為15°.進氣道進口最大捕獲面積約為彈丸最大橫截面積的0.23倍,喉道處橫截面積約為彈丸最大橫截面積的0.33倍。進氣道長度(錐尖點到進氣道亞聲速段出口處)約為彈丸總長的0.29倍。為抵御氣動加熱,進氣道錐尖處及唇罩前緣進行了局部鈍化,錐尖鈍化半徑取為0.15 mm[26-27].按照不同彈道段的技戰術要求,設計3種氣動外形分別如下:

1)沖壓助推外形。炮彈發射出炮口后,尾翼展開,前部膜片破碎,進氣道打開,沖壓空氣流入,與燃料藥柱熱解蒸發產物混合點火并開始燃燒,通過噴管加速、產生推力,使炮彈迅速加速,外形如圖1所示。

圖1 沖壓助推氣動外形Fig.1 Ramjet assisted structure

2)爬升飛行外形。沖壓發動機工作結束,炮彈頭部彈簧鎖打開,中心錐組件向前移動,進氣道關閉,在爬升過程中搜索衛星信號進行定位,外形如圖2所示。

圖2 爬升飛行氣動外形Fig.2 Aerodynamic configuration in climbing flight

3)滑翔控制外形。炮彈在彈道頂點附近,鴨舵展開,控制系統根據炮彈的實際彈道與預設方案彈道偏差,控制炮彈按方案彈道飛行,實現滑翔增程與精度控制,外形如圖3所示。

圖3 滑翔控制氣動外形Fig.3 Sliding control flight structure

2 數值模擬方法

2.1 物理簡化模型、邊界條件與網格劃分

由于不涉及燃燒問題,簡化了燃燒室與補燃室結構,通過改變噴管喉徑來調節燃燒室反壓[1]。彈體外表面與內壁皆采用絕熱無滑移邊界條件,計算域外邊界設置為壓力遠場,來流總溫300 K,總壓101 325 Pa.

對于沖壓助推外形,因為其內部與外部幾何形狀耦合,所以為模擬其黏性邊界層,需要在物面附近采用壓縮比很大的貼體網格。而這些網格又需要與外場的網格相對接,因此本文使用拼接網格技術將炮彈分為內流道與彈體外表面區域,在每個區域中使用貼體網格生成法,并通過設置數據交換面將兩區域結合。彈體內外網格劃分如圖4所示。因為進氣道與中心錐壁面附近氣流速度梯度較大,所以對其進行了加密處理,以提高計算準確性,如圖5所示。

圖4 沖壓助推外形組合網格Fig.4 Combined grid of ramjet assisted structure

圖5 沖壓助推外形頭部網格圖Fig.5 Head grid of ramjet assisted structure

由于爬升飛行、滑翔控制外形的內流道壁仍處于大氣環境中,模擬它們的流場狀態時仍需劃分網格。彈體內外網格如圖6、圖7所示。

圖6 爬升飛行外形組合網格Fig.6 Combined grid of climbing flight structure combined grid

圖7 滑翔控制外形組合網格Fig.7 Combined grid of sliding control flight structure

2.2 計算方法

在三維笛卡爾坐標系下,微分守恒形式Navier-Stokes方程為

(1)

空間離散采用2階迎風型矢通量分裂格式,求解器選用基于密度基的耦合隱式算法。因為沖壓增程彈的內部流場包含由激波和不連續尖銳幾何形狀引起的轉捩流動區域[28],例如在進氣道后擴張段中邊界層在逆壓梯度下容易出現分離,進而造成流動失穩引起分離誘導轉捩,所以湍流模型可選擇適用于大速度梯度與近壁面氣體轉捩流動的γ-Reθt-SA模型,其表達式如下:

(2)

2.3 收斂性分析

為驗證網格數量收斂性,本文設計了4套網格,網格總數分別為410萬、580萬、730萬、820萬,在Ma數為2.5、攻角為0°時計算沖壓助推、爬升飛行、滑翔控制外形的軸向力系數CAr、CAm、CAs,計算結果如表1所示。將410萬、580萬、730萬網格的計算結果分別與820萬網格的計算結果相比,發現對于沖壓助推、滑翔控制外形,網格量為580萬時偏差較小,分別為1.70%、2.13%;對于爬升飛行外形,網格量為410萬時偏差較小,為1.39%.因此,本文中沖壓助推、爬升飛行、滑翔控制3種外形采用的網格總數量分別為580萬、410萬、580萬。

表1 網格收斂性分析Tab.1 Grid convergence analysis

3 數值模擬結果與分析

3.1 數值方法驗證分析

沖壓增程制導炮彈的推力增程性能由其內流道流動特性決定,內壁壓力的計算精度直接影響著整個模型的計算精度。因此,為驗證本文所用網格劃分技術與數值模擬方法的可靠性,需要對內流道的內壁壓力進行驗證。在來流氣壓1 atm、氣溫300 K、Ma數為4.03的工況下,對進氣道喉部直徑48.26 mm、噴管直徑27.94 mm的75 mm沖壓增程炮彈模型進行數值模擬,得到該仿真算例的內壁壓力分布與風洞實驗結果[1]對比如圖8所示。圖9給出了仿真壓力值與實驗值[1]的誤差分布。由圖9可知,計算數據與實驗數據的絕對誤差平均值為1.41%,吻合度較高,表明本文所用的γ-Reθt湍流模型及拼接網格技術對沖壓增程彈丸的仿真模擬合理可行。

圖8 75 mm沖壓增程炮彈的內壁壓力分布Fig.8 Pressure distribution on the inner surface of a 75 mm solid fuel ramjet projectile

圖9 仿真數據與文獻[1]風洞實驗數據的絕對誤差分布Fig.9 Absolute errors of experimental and simulated data

3.2 彈丸流場結構分析

圖10、圖11給出了0°攻角時沖壓助推外形流場的速度、壓力分布情況。由圖10可知:超音速來流在經過中心錐頭錐時產生交于唇口外側的兩道斜激波,波后在唇口外側形成局部高壓區,Ma數減小,但仍大于1;進入喉道后壓力先減小再增加,同時在喉部前收縮段內氣流Ma數不斷減小,喉部后擴張段內Ma數不斷增大;增大到一定程度后超聲速流經一道正激波變為亞音速流,繼續擴張減速后進入后流道,同時在進氣道結束附近存在激波與邊界層的相互作用,迫使一部分氣流向上回流,在彎度較大的后擴張段形成旋流分離區。因為省略了燃燒室與補燃室結構,所以后部速度基本不變,流過尾部拉瓦爾噴管,加速為超音速。

圖10 沖壓助推外形內外流場Ma數圖(0°攻角)Fig.10 Internal and external flow field Mach numbers of the ramjet-assisted structure (0° angle of attack)

圖11 沖壓助推外形內外流場壓力圖(0°攻角)Fig.11 Pressure diagram of internal and external flow fields of the ramjet-assisted structure (0° angle of attack)

在外彈體頭部與圓柱部連接處,由于存在向外折轉角,形成了膨脹波扇區,氣流在此區不斷加速,壓力減小。彈身部分近壁面Ma數由于遠場影響與來流Ma數相同。在尾翼末端,超音速來流向外折轉,與彈體尾部截面一道形成拱形膨脹波,氣流穿過時壓力、密度連續減小、速度增加。在翼面后端與彈底截面處形成擾流,但是因為內流道氣流的加入,所以空氣密度稀薄區域向后移動些許。圖12給出了6°攻角時沖壓助推外形進氣道的壓力分布情況。由圖12可見,外部流場迎風面與背風面流場出現不對稱,迎風面壓力高于背風面壓力。內部流場呈不對稱性。迎風側兩道斜激波交于唇口內側,背風側激波與外罩的斜激波相互干擾,在唇口外側形成局部高壓區,迎風側的波后壓強略高于背風側,流場結構不均勻性增加。同時由于頭部激波對氣流的壓縮不夠充分,流道內的Ma數增加,如圖13所示。

圖12 沖壓助推外形進氣道壓力圖(6°攻角)Fig.12 Pressure diagram of the ramjet-assisted structure (6° angle of attack)

圖13 沖壓助推外形進氣道Ma數圖(6°攻角)Fig.13 Mach number of the ramjet-assisted structure (0° angle of attack)

圖14、圖15給出了0°攻角時爬升飛行外形流場的速度、壓力分布情況。由圖14、圖15可以看出,爬升飛行與沖壓助推外形的最大區別是堵住了進氣道,無氣流流入,流道中氣體稀薄,壓強較低,造成此現象的原因主要有兩個:一是內流道開口背向來流方向;二是炮彈底部的渦流區空氣稀薄。如圖16所示,有攻角時,迎風面與背風面流場出現不對稱,迎風面壓力高于背風面壓力。

圖14 爬升外形內外流場速度圖(0°攻角)Fig.14 Velocity diagram of the internal and external flow fields of the climbing structure (0° angle of attack)

圖15 爬升外形內外流場壓力圖(0°攻角)Fig.15 Pressure diagram of the internal and external flow fields of the climbing structure (0° angle of attack)

圖16 爬升結構內外流場壓力圖(6°攻角)Fig.16 Pressure diagram of the internal and external flow fields of the climbing structure (6° angle of attack)

圖17、圖18給出了0°攻角時滑翔控制外形流場的速度、壓力分布情況。由圖17、圖18可見,滑翔控制與爬升飛行外形的最大區別是展開了鴨舵,通過對比圖15、圖18可以發現兩種外形內流場相似。

圖17 滑翔控制外形內外流場速度圖(0°攻角)Fig.17 Velocity diagram of the internal and external flow fields of the sliding control flight structure (0° angle of attack)

圖18 滑翔控制外形內外流場壓力圖(0°攻角)Fig.18 Pressure diagram of the internal and external flow fields of the sliding control flight structure (0° angle of attack)

滑翔控制與爬升飛行流場差異主要表現在彈頭與鴨舵激波耦合,在遠處合并為一道激波,在鴨舵后形成低壓區。在鴨舵末端,超音速來流向外折轉,形成膨脹波系,氣流穿過膨脹波系時壓力、密度連續減小,速度增加,到達尾翼處形成新激波,氣流壓力升高、速度減小,此間形成了較長的渦流區。如圖19所示,有攻角時迎風面與背風面流場不對稱。

圖19 滑翔控制結構內外流場壓力圖(6°攻角)Fig.19 Pressure diagram of the internal and external flow fields of sliding control flight structure (6° angle of attack)

3.3 彈丸阻力特性分析

圖20給出了0°攻角時3種外形與參考彈(除舵翼與頭部母線外形參數外,其他外形參數相同,不采用沖壓增程形式的鴨式布局制導炮彈)在不同馬赫數下的阻力系數變化曲線[24]。由圖20可見,沖壓助推、滑翔控制、爬升飛行外形在相同條件下對應的阻力系數依次遞減,分別較參考彈阻力系數增大約50.5%、42.9%、33%.圖21給出了在Ma數為2.5下沖壓助推外形的內壁壓力分布曲線。由圖21可以發現:由于激波在內壁面形成一系列激波反射區,內壁表面壓力系數出現周期性波動,隨著激波反射強度逐漸變弱,壓力系數波動幅度逐漸減小,最大壓力值出現在進氣道內部擴展段。因此,通過優化進氣道內部擴展段彎度,可以減小沖壓助推段的阻力;由于噴管入口處為亞音速來流,根據變截面一維流理論,來流在噴管收斂段加速減壓,在噴管喉部壓力系數曲線減小,故可以通過調整噴管喉徑得到適宜的燃燒室壓強。

圖20 4種氣動外形的阻力系數隨Ma數變化曲線Fig.20 Drag coefficients of four structures as a function of Mach number

圖21 沖壓助推外形內壁壓力分布曲線(0°攻角)Fig.21 Wall pressure distribution curve of the ramjet boost structure (0°angle of attack)

表2給出了Ma數為2.5下沖壓助推外形各組件阻力系數。由表2可見:壓差阻力系數占主導地位,約為總阻力系數的87.16%;在各組件阻力系數中內流道阻力系數占59.73%,這是因為氣流經過進氣道后速度降低,在進氣道內壁形成了極大的壓力;且由于在進氣道后擴張段內的氣流靜壓遠高于來流,使得中心錐的前部壓力小于后部壓力,壓差阻力為負值。

表2 沖壓助推外形各部分阻力系數(Ma=2.5)Tab.2 Drag coefficient of each part of ramjet extended range structure (Ma=2.5)

表3給出了Ma數為2.5下爬升飛行外形各部分阻力系數。由表3可以發現,壓差阻力系數仍占主導地位,約為總阻力系數的91.25%.摩擦阻力系數減小,是因為進氣道不再進氣,彈體內空氣密度下降;阻力貢獻最大的組件變為復合前罩(中心錐前部與外彈體前罩組合)與彈翼,占全彈阻力系數的77.67%;因為進氣道關閉使內流道壓差阻力迅速下降,因此貢獻阻力最大的組件發生變化。

表3 爬升飛行外形各部分阻力系數(Ma=2.5)Tab.3 Drag coefficient of each part of the climbing flight structure (Ma=2.5)

表4給出了Ma數為2.5下的滑翔控制外形各組件阻力系數。由表4可見:壓差阻力系數仍占據主導地位,約為總阻力系數的91.26%;阻力系數占比最大的組件變成復合前罩、彈翼與鴨舵,約為80.98%.

表4 滑翔制導外形各部分阻力系數表(Ma=2.5)Tab.4 Drag coefficient of each part of gliding guidance structure (Ma=2.5)

3.4 彈丸升力特性分析

圖22給出了6°攻角時3種外形與參考彈在不同Ma數下的升力系數變化曲線。由圖22可見:4種外形的升力系數規律一致,在超音速范圍內都隨Ma數增大而減小,其中沖壓助推與爬升飛行外形升力系數相近,相差4%;同一工況下,滑翔控制外形較參考彈的升力系數約小11.9%,這是因為其進氣道外形限制了鴨舵的存放空間,從而其較參考彈的鴨舵面積約小50%.圖23給出了Ma數為2.5時4種外形的升力系數隨攻角的變化曲線,從中可以發現4種外形的升力都隨攻角增大而增大。

圖22 4種外形的升力系數隨Ma數變化(6°攻角)Fig.22 Change of lift coefficients of four structures with Mach number (6° angle of attack)

圖23 4種外形的升力系數隨攻角變化曲線(Ma=2.5)Fig.23 Change of lift coefficients of three structures with angle of attack (Ma=2.5)

表5給出了Ma數為2.5時3種外形各組件的升力系數,從中可以看出:3種外形的主要升力都由外彈體貢獻,但滑翔控制外形因鴨舵展開,升力系數較其他兩種外形約大12%;沖壓助推外形中內流道升力貢獻約為6.1%,后兩種外形與之相比,由于進氣道關閉,在復合流道內產生微弱負升力約1.8%.

圖24給出了參考彈與滑翔控制外形在不同Ma數下的壓心位置(彈丸壓力中心距彈頂的距離)變化曲線。由圖24可見,在超音速范圍內,兩種外形的壓心位置都隨Ma數增大而減小,但同一工況下,滑翔控制外形較參考彈的壓心位置約靠后3.12%.圖25給出了參考彈與滑翔控制外形在不同Ma數下的穩定力矩系數導數變化曲線。由圖25可見,在超音速范圍內,兩種外形的穩定力矩系數導數都隨Ma數增大而減小,但同一工況下滑翔控制外形較參考彈的穩定力矩系數導數約大5.2%.對沖壓增程制導炮彈滑翔控制外形,經計算滑翔控制飛行段的全彈靜穩定性儲備量約6%,操縱性關系為0.4,即舵面偏轉10°約產生4°的平衡攻角,因此該外形在滑翔控制飛行時的氣動性能良好,穩定性與操縱性、舵偏角與平衡攻角匹配較好。

表5 3種外形各組件升力系數表(Ma=2.5)Tab.5 Lift coefficient of each component of the three structures (Ma=2.5)

圖24 參考彈與滑翔控制外形的壓心位置(口徑倍數)隨Ma數變化曲線Fig.24 Change of pressure center positions (aperture multiple) of reference projectile and sliding control profile with Mach number

圖25 參考彈與滑翔控制外形的穩定力矩系數導數隨Ma數變化曲線Fig.25 Change of stability moment coefficient derivatives of reference projectile and sliding control profile with Mach number

3.5 彈體擺動對沖壓發動機進氣特性的影響

圖26給出了沖壓助推段進氣道流量系數隨攻角的變化曲線。由圖26可知,在相同Ma數下,流量系數隨攻角的增大而減小。隨著攻角增大,進氣口的有效面積減小,迎風側激波角增大,波系交于唇口外,進氣道出現溢流,使流量系數減小。

圖26 進氣道的流量系數隨攻角變化曲線(Ma=2.5)Fig.26 Change of flow coefficient of inlet with angle of attack (Ma=2.5)

圖27給出了沖壓助推段進氣道總壓恢復系數隨攻角的變化曲線。由圖27可知,在相同Ma數下,總壓恢復系數隨攻角的增大而減小。隨著攻角增大,進氣道中的結尾正激波在迎風側和背風側的位置差異增大,增強了激波的扭曲程度,降低了進氣道的總壓恢復系數。

圖27 進氣道的總壓恢復系數隨攻角變化曲線(Ma=2.5)Fig.27 Change of total pressure recovery coefficient of inlet with angle of attack (Ma=2.5)

4 結論

本文針對沖壓增程制導炮彈在不同彈道階段的氣動特性,依據其工作原理與飛行特點,設計沖壓助推、爬升飛行、滑翔控制狀態所對應的氣動外形,采用拼接網格技術與γ-Reθt-SA湍流模型進行了三維流場模擬與數值仿真,獲得了炮彈流場與升阻力特性。通過分析得到以下結論:

1)75 mm沖壓增程炮彈的內壁壓力仿真結果與風洞實驗數據的絕對誤差平均值為1.41%,吻合度較高,表明采用拼接網格技術與γ-Reθt-SA模型進行沖壓增程制導炮彈三維流場數值模擬,具有良好的精度。

2)沖壓助推、滑翔控制、爬升飛行外形在相同條件下對應的阻力系數依次遞減,在超音速范圍內都隨Ma數增大而減小,分別較參考彈阻力系數約大50.5%、42.9%、33%,表明后兩種外形向前推進中心錐組件可以減小阻力;沖壓助推外形內流場阻力占支配地位,為59.73%,后兩種外形外流場阻力占支配地位,為80.98%和77.67%.

3)沖壓助推、滑翔控制、爬升飛行外形的升力系數規律一致,在超音速范圍內相同攻角時都隨Ma數的增大而減小,相同Ma數時隨攻角的增大而增大。同一工況下,滑翔控制外形因鴨舵展開,升力系數較其他兩種外形約大12%,又因進氣道外形限制了鴨舵的面積,升力系數較參考彈約小11.9%;內流道影響輕微,沖壓助推外形內流道升力貢獻約6.1%,后兩種形態因進氣道關閉,內流道產生微弱負升力約1.8%.

4)彈體擺動對沖壓發動機進氣道的總體性能有不利影響,即在相同Ma數下,隨著攻角的增大,進氣道的流量系數和總壓恢復系數減小,在6°攻角時分別較0°攻角時減小了13%、7.2%.

主站蜘蛛池模板: 国产日韩欧美中文| 99久久精品视香蕉蕉| 欧美一区二区啪啪| 一级不卡毛片| 极品国产一区二区三区| 欧美一级大片在线观看| 午夜人性色福利无码视频在线观看| 韩国自拍偷自拍亚洲精品| 最新国产成人剧情在线播放| 国产精品自在线拍国产电影 | 中文字幕第4页| 高清无码不卡视频| 波多野结衣一区二区三区四区视频| 国产亚洲一区二区三区在线| 露脸真实国语乱在线观看| 国产在线无码一区二区三区| 人妻精品全国免费视频| 美女被操黄色视频网站| 久久久精品久久久久三级| 色婷婷综合在线| 97se亚洲综合在线天天| 欧美国产另类| 91在线无码精品秘九色APP| 九九这里只有精品视频| 久久99热这里只有精品免费看| 一级毛片免费不卡在线视频| 国产精品香蕉在线| 欧美一级爱操视频| 中文无码伦av中文字幕| 国产黑丝一区| 欧美亚洲第一页| 久久综合五月| 2020国产精品视频| 久久99蜜桃精品久久久久小说| 成人福利在线观看| 国产高清自拍视频| 亚洲欧美自拍视频| 新SSS无码手机在线观看| 国产无码精品在线| 国产门事件在线| 亚洲精品在线影院| 日韩无码黄色网站| 欧美日韩亚洲国产| 国产免费人成视频网| 久久不卡精品| 尤物亚洲最大AV无码网站| 丝袜亚洲综合| 最新国产精品鲁鲁免费视频| 欧洲日本亚洲中文字幕| 国产亚洲视频播放9000| 国产凹凸视频在线观看| 国产成人免费| 国产日韩精品欧美一区灰| 久热re国产手机在线观看| 尤物国产在线| 91娇喘视频| 欧美精品1区| 久操中文在线| 国产成人精品男人的天堂| 中国黄色一级视频| 亚洲国产日韩在线观看| 日本国产在线| 在线另类稀缺国产呦| 久久77777| 国产欧美日韩精品第二区| 久久精品国产在热久久2019 | 国产综合在线观看视频| 伦伦影院精品一区| 欧美日韩免费在线视频| 99人体免费视频| 99久久精品无码专区免费| 在线播放国产一区| 草草线在成年免费视频2| 亚洲一级毛片在线播放| 国产美女在线观看| 午夜福利在线观看入口| 91福利国产成人精品导航| 91视频区| 伊伊人成亚洲综合人网7777| 99免费在线观看视频| a毛片免费在线观看| 一级高清毛片免费a级高清毛片|