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

固定鴨舵雙旋火箭彈超聲速側向氣動特性

2021-03-29 01:02:38黃智康陳少松呂代龍
彈道學報 2021年1期
關鍵詞:模型

黃智康,陳少松,魏 愷,呂代龍

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

為了適應現代戰爭對武器精度的要求,固定鴨舵彈道修正技術憑借其控制方式簡單、制導精度較高和易于集成在現有常規彈藥上而得到廣泛研究。固定鴨舵修正組件(precise guidance kit,PGK)是一種分離式修正組件,即修正組件的滾轉速度與后體不同,因此簡稱雙旋彈。修正組件上固定一對差動舵和一對同向舵,舵偏角是固定的,因此操縱力的大小不能隨意改變,通過改變鴨舵的滾轉角對彈箭進行彈道修正。差動舵在無控階段提供導轉力矩使修正組件相對后體反向旋轉,同向舵在控制階段提供控制力和力矩對彈箭進行彈道修正。

固定鴨舵彈道修正技術的大量研究工作集中在高速旋轉的無尾翼炮彈上,Sahu等[1]對安裝不同修正組件的炮彈的馬格努斯效應進行了數值分析;吳萍[2]通過風洞試驗得到雙旋彈氣動特性的變化規律;張蛟龍[3]建立了雙旋炮彈七自由度動力學模型,結合數值仿真得到的氣動系數,對無控和有控彈道特性進行分析;徐輝雯等[4]對比了雙旋炮彈修正組件反旋與不旋的氣動特性。

在火箭彈上安裝固定鴨舵修正組件的研究集中在動力學與控制方向,在氣動方面的研究比較少。王少飛[5]建立了雙旋火箭彈的六自由度剛體有控彈道模型,并在此基礎上分析了固定鴨舵修正組件用于火箭彈的可行性;王志剛和李偉等[6-7]用凱恩方法對雙旋火箭彈進行多體動力學建模,并對其運動特性進行分析;黨建濤等[8]通過數值模擬分析了雙旋火箭彈在定常狀態下鴨舵對尾翼的氣動干擾;Yin等[9-10]對無鴨舵和有鴨舵的旋轉尾翼導彈的流場進行了數值模擬,研究了有鴨舵彈箭的馬格努斯效應。

鴨式布局的彈箭有著較高的控制效率,但是在安裝鴨舵修正組件后,火箭彈的氣動特性尤其是側向特性會發生較大改變。首先,修正組件本身會在不旋和反旋時產生不同的側向力;其次,由于鴨舵修正組件和后體轉速不一致,鴨舵尾流會干擾彈體周圍的流場,導致彈體的馬格努斯效應發生改變;最后,鴨舵的下洗流會影響尾翼的氣動特性,產生額外的側向力和偏航力矩。因此,本文采用CFD軟件對某雙旋火箭彈的氣動特性進行數值模擬,重點研究該火箭彈后體在低速旋轉狀態下,鴨舵修正組件不旋和反旋狀態對側向力特性的影響。

1 計算方法與研究對象

1.1 計算方法

本文使用商業流體力學計算軟件ANSYS Fluent來計算彈體周圍流場的三維可壓縮非定常Navier-Stokes方程:

式中:V為控制體,?V為包圍控制體的封閉表面,S為控制體表面的面積,W為守恒變量矢量,Fc為對流通量矢量,Fv為黏性通量矢量。采用雷諾平均法對N-S方程進行分解。空間離散采用二階迎風格式,對流通量應用Roe通量差分分裂。非定常計算采用雙時間步法,其中物理時間描述實際運動,內迭代時間使計算收斂,二者皆采用隱式迭代求解。將定常無旋轉狀態的計算結果作為非定常計算的初始條件以快速穩定地得到收斂結果。

為模擬彈箭的旋轉運動,本文采用滑移網格方法。如圖1所示,將流場域劃分為動、靜2個部分,動域根據設置的方式進行平動或旋轉,動、靜域之間通過交界面進行數據傳輸。該方法相對于多重坐標系法,精度更高;而相對于動網格法,因為網格拓撲結果不需要改變而可以獲得較高的網格質量。

圖1 滑移網格區域劃分示意圖

1.2 湍流模型的選擇和計算方法的驗證

計算方法的驗證模型為Basic Finner。遠場邊界為壓力遠場邊界條件,模型邊界設置為無滑移邊界條件。實驗數據參考文獻[11],實驗數據的風洞條件見表1,其中,Ma∞表示來流馬赫數,p0和T0分別表示總壓和總溫,ReD表示基于參考長度的雷諾數,ω表示轉速。

表1 實驗風洞條件

在此基礎上采用2種不同湍流渦黏模型,分別是S-A(Spalart-Allmaras)模型[12]和Transition SST模型[13-15],計算結果見圖2,圖中,Cn表示法向力系數,Cl表示滾轉力矩系數。可以看到,在20°攻角以下Transition SST湍流模型計算的滾轉力矩系數結果與實驗值更加貼合,因此本文選擇Transition SST湍流模型進行后續計算。

圖2 氣動系數計算值與實驗值結果

1.3 研究對象

本文研究對象如圖3、圖4所示。圖中,φ表示后體滾轉角;φcf表示修正組件滾轉角,下標cf表示修正組件。該模型彈頭為固定鴨舵修正組件,彈體和尾翼組合體簡稱為后體。最大彈徑D=106 mm,全彈長L=11.41D,其中修正組件安裝一對差動舵(CANARD1和CANARD3)和一對同向舵(CANARD2和CANARD4),舵偏角δc=3°,彈尾部安裝6片尾翼(FIN1,…,FIN6),尾翼斜置角δw=0.5°。參考長度為全彈長L,參考面積S=πD2/4,力矩參考點距離彈頭0.45L,參考坐標系取彈體系。遠場靜壓為101.325 kPa,靜溫為288.15 K。

圖3 雙旋火箭彈模型

圖4 鴨舵和尾翼布局示意圖

三維六面體網格如圖5所示。計算域邊界參考彈箭尺寸,軸向長度為12L,徑向長度為30D。為了模擬鴨舵修正組件和后體以不同轉速和方向旋轉的過程,使用滑移網格方法,將計算域分為3部分。外部為靜域,內部分為前后2個區域,前部包含鴨舵修正組件,后部包含彈體和尾翼。第1層網格高度設置為6×10-7m,以保證網格在壁面附近的y+<1,并使前20層網格的伸展比保持在1.1。

圖5 模型網格

1.4 無關性驗證

網格無關性驗證的結果如表2所示,來流工況為Ma=2.0,α=6°。時間步長的無關性驗證的結果如表3所示,修正組件旋轉方向為X軸負向,轉速為408.4 rad/s;后體旋轉方向為X軸正向,轉速為119.4 rad/s。表中,N表示網格數量,Δt表示時間步長,CZ表示側向力系數,|ηCn|和|ηCZ|分別表示Cn和CZ的相對誤差絕對值。

表2 網格無關性驗證

表3 時間步長無關性驗證

根據計算結果,450萬網格數和650萬網格數的側向力系數最大相差4.30%;時間步長為1×10-5s和5×10-6s的結果最大相差6.15%。因此,選擇450萬的網格及1×10-5s的時間步長,可以滿足計算對網格密度和時間精度的要求。內迭代步數設置為20可滿足收斂要求。

圖6 雙旋火箭彈氣動力系數隨攻角變化曲線

2 計算結果

本文給出了無鴨舵模型(NCA)、修正組件不旋(NRO)和修正組件反旋(ROT)3種狀態在馬赫數Ma=1.6,2.0,攻角α=0°,2°,6°,10°的平均氣動力,如圖6所示,圖中,Ca,Cn,CZ分別為軸向力系數、法向力系數和側向力系數。

圖6中,鴨舵修正組件不旋設置為俯仰姿態修正的狀態,同向舵提供正向法向力,鴨舵修正組件反旋的轉速和后體轉速設置為相應馬赫數下α=0°對應的平衡轉速(通過計算后體旋轉和不旋轉時的滾轉力矩確定),見表4。

表4 雙旋火箭彈的轉速設置

加裝鴨舵修正組件后軸向力系數增加,其中,Ma=1.6和Ma=2.0時修正組件在不旋狀態下的軸向力系數相對于無鴨舵模型分別增加了約1.3%和1.5%,反旋狀態下的軸向力系數相對于無鴨舵模型分別增加了約0.7%和0.9%。

法向力系數隨攻角和馬赫數的增大而增大。圖7給出了修正組件反旋時修正組件法向力系數隨滾轉角周期性變化曲線。圖7及后文圖中,下標b表示彈體,f表示尾翼,aft表示后體,cf表示修正組件。由圖7可知,修正組件反旋時的法向力較不旋時提供的平均法向力小,但由于修正組件提供的法向力始終為正,因此反旋狀態的平均法向力比無鴨舵模型大。

圖7 Ma=1.6,α=2°,6°時鴨舵修正組件反旋時修正組件法向力系數隨滾轉角φcf變化曲線

在加裝鴨舵修正組件后側向力系數隨攻角非線性增加,修正組件反旋的狀態下其側向力系數增加量有所減少,即反旋的狀態有助于抑制側向力隨攻角的增加。

圖8給出了Ma=1.6時修正組件彈體和尾翼的部件側向力隨攻角變化曲線。無鴨舵模型的彈體和尾翼提供的側向力大小大致相同,方向相反。安裝鴨舵修正組件后,彈體和尾翼的側向力顯著增加,其中彈體提供的側向力比尾翼大很多。小攻角時,修正組件反旋的彈體和尾翼側向力隨攻角變化的增量小于不旋狀態;隨著攻角的進一步增加,反旋狀態的彈體和尾翼側向力逐漸超過不旋狀態。同時α=10°時,修正組件不旋時尾翼的側向力并不會隨著攻角增加而持續增加,而是保持與無鴨舵模型相近的大小。

圖8 Ma=1.6時部件側向力系數隨攻角變化曲線

彈頭的鴨舵修正組件本身提供的側向力在不旋和反旋時方向相反,相對于尾翼和彈體而言是一小量,但是由于其離質心距離較遠,所以仍會提供不容忽視的偏航力矩。鴨舵修正組件和后體偏航力矩曲線見圖9。由圖可見,反旋狀態下修正組件提供的偏航力矩最高能達到全彈偏航力矩的40%。

2.1 無鴨舵模型側向氣動特性分析

對于無鴨舵模型,小攻角時無論是彈體還是尾翼的側向力都很小,僅在α=10°時側向力才有明顯的增加。圖10為無鴨舵模型在Ma=1.6,α=10°時,彈體段4個截面上的無量綱渦量Ω*(Ω*=ΩD/v∞,其中Ω表示渦量大小,v∞表示來流速度)云圖和空間流線圖。如圖所示,彈體前段x/L=0.22處并未發生流動分離。在彈體的中后段,來流繞過彈體形成分離渦附著在彈體上表面,并且體渦區隨著軸向距離增加而擴大。盡管彈體低速旋轉會使彈體渦隨軸向距離增加,略微沿旋轉方向偏轉,但是整體流場仍保持左右對稱,產生的側向力系數也相應較小。

圖10 Ma=1.6,α=10°,無鴨舵模型x/L=0.22,0.42,0.62,0.82處無量綱渦量云圖和空間流線圖

圖11為無鴨舵模型彈體4個截面處的壓力系數Cp沿周向分布曲線。彈體旋轉產生馬格努斯力的主要原因在于φ=90°和φ=270°附近的壓力系數。在x/L=0.22處,彈體繞流并未分離,壓力系數曲線分布左右對稱。隨著邊界層沿軸向發展,流動發生分離,導致左右兩側的壓力系數不對稱加劇。同時彈體分離渦處的壓力系數分布也受旋轉影響,產生的正向側向力與上面所描述的側向力方向正好相反,但由于這部分影響作用在彈體上表面,因此產生的側向力并不大。

圖11 Ma=1.6,α=10°,無鴨舵模型不同x/L處的壓力系數沿周向分布曲線

圖12是無鴨舵模型尾翼FIN1旋轉一周提供的側向力。小攻角時,彈體上表面渦跡較小甚至沒有,對尾翼影響也只限于翼根,此時尾翼產生的側向力很小。如圖12中α=2°的曲線,在平衡轉速下側向力曲線呈現振幅相等的2個振蕩周期。隨著攻角增加,彈體渦影響范圍擴大,對尾翼的側向力產生影響。如圖12中α=10°的曲線,在背風區(φ∈[0°,90°]∪[270°,360°]),側向力曲線左右出現不對稱,同時在φ=0°附近產生了第3個振蕩。圖13給出了α=10°時尾翼在φ=0°附近的流線圖,由圖可見,彈體兩側的繞流在頂部匯聚,在φ=0°附近,當地流動方向會發生偏折,改變了尾翼側向力的方向,由此產生了第3個振蕩周期。

圖12 Ma=1.6,無鴨舵模型FIN1提供的側向力系數隨滾轉角φ變化曲線

圖13 Ma=1.6,α=10°,無鴨舵模型表面壓力系數云圖和空間流線圖

2.2 修正組件不旋時側向氣動特性分析

當鴨舵修正組件不旋時,其本身產生的側向力來自攻角的影響。當同向舵處于橫向平面而提供額外升力時,差動舵處于縱向平面并提供一個負的滾轉力矩,如果此時攻角不為0,那么此時處于迎風面的差動舵提供的正向側向力將大于背風面差動舵提供的負向側向力。圖14為迎風側和背風側鴨舵表面壓力系數云圖,可以清晰地觀察到迎風面鴨舵上下表面壓差大于背風面鴨舵上下表面的壓差。

圖14 Ma=1.6,α=6°,鴨舵修正組件不旋時背風側鴨舵(CANARD1)和迎風側鴨舵(CANARD3)表面壓力系數云圖

圖15給出了無鴨舵模型和修正組件不旋狀態下彈體截面側向力系數軸向分布曲線。加裝鴨舵修正組件后,側向力系數增加明顯。從x/L=0.26至x/L=0.88,截面側向力系數負向增加,且攻角越大,其變化的斜率越大。在x/L=0.26之前和x/L=0.88之后,因截面不連續變化,彈體截面側向力出現震蕩變化。并且彈體尾部區域受到尾翼干擾,截面側向力隨軸向距離增加迅速減小,方向變為Z軸正向。由此可見,彈體提供的側向力主要集中在彈體中后段。

圖15 Ma=1.6,無鴨舵模型和修正組件不旋時彈體截面側向力系數沿軸向分布曲線

鴨舵的尾渦和彈體的相互干擾是彈體側向力增加的主要原因。為了進一步分析鴨舵尾渦對彈體周圍壓力分布的影響,給出彈體不同截面的壓力系數周向分布曲線和渦量云圖,如圖16和圖17所示。

由圖可見,在小攻角時,鴨舵尾渦分布在彈體φ=0°,90°,180°,270° 4個方向上。在尾渦對應的位置,彈體的壓力系數會出現波動,波動的幅度與當地尾渦的強度有關。顯然,軸向距離越小,鴨舵尾渦更強,波動幅度就會更大;同向舵的當地攻角是舵偏角和來流攻角之和,其尾渦強度比處于俯仰平面的差動舵更強,因此,φ=90°和φ=180°方向的波動幅度更大。從壓力系數曲線和渦量云圖中可以觀察到,壓力系數發生波動的位置隨著軸向距離的增加而改變,攻角不為0,使得尾渦呈現向上偏移的趨勢;彈體旋轉使得尾渦隨著軸向距離增加略微向著旋轉方向偏移。

圖16 Ma=1.6,無鴨舵模型和修正組件不旋時x/L=0.22,0.42,0.62,0.82處壓力系數周向分布曲線

圖17 Ma=1.6,修正組件不旋時彈體截面的無量綱渦量云圖和流線圖

α=10°時,受迎風面的影響,180°方向的尾渦很快就消失了,因此彈體底部的壓力系數與無鴨舵模型完全一致。φ=0°方向的尾渦因為向上偏移,對彈體表面的影響變小。φ=90°和φ=270°方向的鴨舵尾渦在隨著軸向距離增加并上移的過程中,不僅會直接影響彈體表面的流動,而且和彈體體渦相互干擾,且旋轉使渦流順時針偏移,從而使得渦流在彈體背風面呈現非對稱的分布,在壓力系數曲線圖上表現為左右不對稱的加劇;同時,在彈體頂部區域,無鴨舵模型彈體背風面存在的正向側向力也因受鴨舵尾渦的影響而大幅減弱。

為了確定鴨舵對尾翼的影響,圖18給出了在鴨舵組件不旋狀態下單獨FIN1旋轉一周的側向力,并與無鴨舵模型對比。對比結果表明,鴨舵尾流對迎風面尾翼基本沒有影響,主要影響區域集中在背風區;且隨著攻角增大,影響區域逐漸減小。這種影響主要來自于鴨舵尾渦和彈體體渦的互相干擾,使得尾翼段流場更加復雜并且失去對稱性。當α<6°時,影響區域有側向力正向增加的現象,而在α=10°時尾翼的側向力并沒有明顯增加的趨勢。

圖18 Ma=1.6,無鴨舵模型和修正組件不旋時FIN1提供的側向力系數隨滾轉角φ的變化曲線

圖19給出了不同攻角時鴨舵尾流流線圖和尾翼處截面的壓力云圖。由圖可見,隨著攻角增加,鴨舵尾流對尾翼的影響區域逐漸變小,同時彈體兩側的鴨舵受來流攻角的影響最大,因此尾渦更強。迎風面鴨舵的尾流對尾翼影響不大,因此圖中沒有給出。當α=2°和α=6°時,背風側和彈體兩側鴨舵尾流能夠直接影響尾翼區域,因此處于φ=0°附近的尾翼上下表面壓力系數分布改變非常明顯。當α=10°時,彈身兩側的鴨舵尾渦更強,對背風面的壓力系數改變更加明顯,但是其對尾翼段流場的直接影響區域超過了尾翼的展長,因此并沒有加劇尾翼側向力,這與圖8和圖18中尾翼側向力并沒有隨攻角進一步增加的結果相符。

圖19 Ma=1.6,無鴨舵模型和修正組件不旋時鴨舵尾流流線圖和x/L=0.97處截面壓力云圖

2.3 修正組件反旋時側向氣動特性分析

鴨舵修正組件在反旋時提供的側向力會比不旋狀態大,同時方向相反,因為此時鴨舵不僅產生側向力,也會干擾彈頭旋轉體反旋的馬格努斯效應,二者隨修正組件滾轉角φcf變化的曲線見圖20,圖中,CANARDs表示4片鴨舵之和。

圖20 Ma=1.6,修正組件反旋時鴨舵側向力系數隨修正組件滾轉角φcf變化曲線

單片鴨舵提供的側向力的幅值隨攻角增加而增加,但4片鴨舵的總側向力卻并沒有增加很多。總側向力曲線的極值對應φcf=90°和φcf=270°,此時同向舵正處于縱向平面內,產生最大的瞬時側向力。

3片舵偏角與滾轉方向相同的鴨舵(CANARD1,CANARD2,CANARD3),其側向力曲線幅值相同,隨滾轉角變化只相差90°的相位,而另外的1片鴨舵(CANARD4),由于其舵偏角方向與滾轉方向相反,因而側向力曲線的幅值更大。這就導致當攻角不為0時,CANARD4于迎風面和背風面產生的側向力的差值也要比其余鴨舵大得多。修正組件反旋時,鴨舵會干擾彈頭旋轉體段的馬格努斯力,如圖20所示,最主要的變化集中在φcf=270°附近,且隨著攻角增大,該影響也隨之增大,結果使側向力增加。

圖21給出了修正組件反旋時彈體截面側向力系數沿軸向分布。由圖可見,修正組件反旋相比不旋情況,彈體前半段側向力增加較小。同時,當α=2°和α=6°時,反旋時彈體的側向力增量相比不旋狀態也小很多,僅在α=10°時,側向力系數隨軸向距離增加而大幅增加。

圖21 Ma=1.6,無鴨舵模型和修正組件反旋時彈體截面側向力系數沿軸向分布曲線

圖22 Ma=1.6,彈體瞬時側向力系數隨時間變化曲線

圖23對比了FIN1在組件不旋和反旋狀態下旋轉一周的側向力系數。當α<6°時,側向力系數與修正組件不旋時大致相同,反旋時鴨舵對尾翼的主要影響同樣集中在背風區。主要的區別在于,α=10°時,尾翼側向力系數有明顯增加,這也與圖8中所述一致。原因在于,修正組件反旋后,鴨舵的尾流在背風區周期變化,不僅對尾翼段的影響范圍變大,而且此時鴨舵尾流會直接影響尾翼,導致尾翼側向力增加。

圖23 Ma=1.6,修正組件不旋和修正組件反旋時FIN1提供的側向力系數隨滾轉角φ變化曲線

3 結論

本文通過計算流體力學軟件ANSYS Fluent計算了雙旋火箭彈在超聲速下修正組件不旋和反旋的氣動特性,并與無鴨舵模型進行對比,重點分析了加裝鴨舵修正組件后的側向力特性,得到了以下結論:

①在加裝鴨舵修正組件后,全彈的側向力隨攻角增加而大幅增加,其中彈體側向力增加最多,即受鴨舵組件的影響最大,是全彈側向力的主要組成部分。彈頭修正組件不旋和反旋時的側向力方向相反。尾翼受鴨舵影響,其側向力也會增加。

②當修正組件不旋時,其彈頭鴨舵修正組件提供的側向力來源于迎風面鴨舵與背風面鴨舵提供的側向力的差值。彈體側向力隨攻角大幅增加,主要源于鴨舵尾渦與彈體體渦相互干擾,加劇了彈體兩側壓力分布的不對稱。鴨舵尾渦對尾翼影響區域集中在背風區。

③當修正組件反旋時,彈頭鴨舵修正組件的側向力來源于舵偏角與滾轉方向相反的鴨舵(CANARD4);彈體側向力隨時間變化曲線的周期被打破,振蕩幅值大幅增加;同時,修正組件反旋使鴨舵對尾翼的影響加劇。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 91小视频版在线观看www| 国产呦精品一区二区三区网站| aⅴ免费在线观看| 日韩无码白| 欧美精品另类| 亚洲无码高清一区二区| 国产aⅴ无码专区亚洲av综合网| 国产99视频免费精品是看6| a天堂视频| 欧美日韩在线国产| 亚洲成人精品| 日韩麻豆小视频| AV在线天堂进入| 在线无码九区| 久久精品一卡日本电影| 亚洲人精品亚洲人成在线| 一区二区自拍| 激情亚洲天堂| 亚洲浓毛av| 国产在线观看成人91| 91精品日韩人妻无码久久| 日本一本正道综合久久dvd| 无码免费的亚洲视频| 国产成人三级| 欧美日韩国产在线播放| 精品无码一区二区三区在线视频| 91欧洲国产日韩在线人成| 中文字幕啪啪| 亚洲精品麻豆| 色九九视频| 一区二区三区成人| 欧亚日韩Av| 蜜桃视频一区二区| 久久婷婷六月| 欧美亚洲国产精品久久蜜芽| 黄色网站不卡无码| 中文字幕资源站| 在线观看亚洲天堂| 国产人人干| 欧日韩在线不卡视频| 这里只有精品国产| 亚洲国产系列| 在线欧美一区| 欧美啪啪一区| 欧美一级在线看| 在线观看av永久| 欧美在线综合视频| 91久久国产综合精品| 99久久国产精品无码| 天天爽免费视频| 欧美日韩亚洲国产| 制服丝袜国产精品| 久久美女精品| 丁香婷婷激情综合激情| 欧美日韩成人在线观看| 又大又硬又爽免费视频| 久久精品国产精品国产一区| 成人福利在线视频| 国产综合精品一区二区| 色亚洲成人| 青青操视频在线| 69免费在线视频| 性69交片免费看| 国产乱人视频免费观看| 国内99精品激情视频精品| 亚洲国产综合精品中文第一| 国产在线自乱拍播放| 亚洲精品视频免费观看| 99re在线免费视频| A级全黄试看30分钟小视频| 女人毛片a级大学毛片免费| 久久semm亚洲国产| 亚洲中文精品久久久久久不卡| 久久久久国产精品嫩草影院| 亚洲丝袜第一页| 亚洲成人一区二区| 亚洲无码高清免费视频亚洲| 久久精品亚洲中文字幕乱码| 欧美亚洲第一页| 国产日韩欧美中文| 国产成人无码综合亚洲日韩不卡| 国产99免费视频|