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

艦載機著艦的動力學建模

2014-06-23 07:46:14夏桂華董然孟雪朱齊丹
哈爾濱工程大學學報 2014年4期
關鍵詞:模型

夏桂華,董然,孟雪,朱齊丹

(哈爾濱工程大學自動化學院,黑龍江哈爾濱150001)

建立艦載機在著艦狀態下的動力學模型,是完成著艦控制系統綜合仿真的根本前提,該模型總體上分為兩類——描述6自由度運動的非線性模型和線性擾動模型;前者能全面描述艦載機的運動,是對飛行控制系統進行試飛前評定時必然被采用的,但是運動方程的強耦合、非線性的特點導致直接進行非線性系統分析難以得出具有普遍意義的結論;相比之下,線性模型是非線性運動方程經過擾動線性化后導出的,足以描述只圍繞標稱狀態做小幅度機動的著艦飛機對駕駛員輸入以及外部流場作用的響應特性[1],利用完備的線性分析手段,可有效地研究艦載機的穩定性和操縱性,作為飛行控制系統的設計指南。美國海軍的研究經驗表明:在著艦控制系統的分析和完善階段,將能描述艦載機主要特點的線性擾動模型作為地面物理仿真時使用的飛機數學模型是合適的[2-4]。

線性模型的精度很大程度上依賴于線性化方法和側重點的選取。國內已公開的關于艦載戰斗機線性建模的文獻(較有代表性的如文獻[5-8]),主要基于F/A-18飛機的結構及性能參數,在實現非線性建模的基礎上,推導出其在風場中著艦的線性擾動模型;根據“大氣擾動對著艦的主要影響是由擾流的垂直分量引起的艦載機在垂直方向上的軌跡誤差”[9-10],對于考慮紊流因素的線性化過程做了簡化—只考慮擾流垂直分量的作用效果,導致得出的線性模型中不含水平擾流項。

本文驗證了該簡化方式不合理。為此,首先研究了艦尾流場的特點,通過定性理論分析得出其水平分量對艦載機的飛行穩定性有顯著影響;進而,將尾流在飛行平面內的垂直和水平分量引起的附加氣動力和力矩作為干擾項引入已經推導出的無風情況下的線性模型,得到完整的艦載機縱向線性動力學模型;最后,通過艦尾流場下的艦載機線性/非線性模型的對比仿真證明結論。

1 艦載機本體特性的非線性建模

運動方程的小擾動線性化一般是在飛行器基準運動狀態下完成的[11]。假設艦載機著艦階段的基準狀態為:無大氣擾動情況下,艦載機在鉛垂平面內以航跡角-3.5°、迎角 8.1°、航速 70 m/s 的勻速直線降落。建立艦載機非線性模型并求解基準狀態是運動方程線性化的前提。為簡化所研究問題的復雜性,僅針對艦載機的本體特性完成非線性建模,包括發動機推力特性和由氣動布局決定的空氣動力性能,不考慮飛行控制系統,并結合艦載機的運動特點和建模需要,給出下列假設[3,12]:

1)視大地為平面,重力加速度不隨艦載機的位移而變化;

2)將艦載機視為在幾何外形、質量分布和氣動特性方面均對稱的剛體,且不考慮機體彈性變形和旋轉部件的影響;

3)不考慮著艦過程中由燃料消耗造成的艦載機質量和質心的變化,將發動機產生的推力當作外力;

4)處理紊流對飛行的影響時,將飛機近似為在流場中運動的一個質點。

為描述艦載機的動力學方程,本文按照通常的歐美坐標體制規范[12-13],定義的直角坐標系為:大地系Sg,牽連大地系S'g,機體坐標系Sb,氣流坐標系Sa,航跡系Sk。

1.1 艦載機受力分析

1.1.1 發動機推力

載機的總推力位于機體對稱平面內,在仿真條件下,合力作用線通過飛機質心,與機體縱軸夾角為0°。影響發動機推力大小的因素是飛行高度h、馬赫數Ma和油門開度δp。油門桿指令δpl控制油門開度,δp對δpl的響應特性可近似為一階慣性環節1/(1+s)。

1.1.2 空氣動力

把艦載機所受空氣動力在氣流坐標系Sa中分解為:沿-x方向的氣動阻力D,沿y方向的氣動側力C,沿-z方向的氣動升力L。把空氣動力對飛行器質心的力矩在機體系Sb中分解為:繞x軸的滾轉力矩Lr,繞y軸的俯仰力矩M,繞z軸的偏航力矩N。計算氣動力、氣動力矩公式分別為

式中:S為機翼面積,m2;ρ為空氣密度,kg·m-3;U為艦載機空速大小,m·s-1;b為翼展,m;CA為機翼平均氣動弦長,m;CL、CD、CC分別為升力系數、阻力系數、側力系數;Cl、Cn、Cm分別為滾轉力矩系數、偏航力矩系數、俯仰力矩系數。

艦載機的氣動力、氣動力矩系數與迎角、側滑角、操縱面偏角等許多因素相關,被表示為多個部分的組成,每一部分以表格形式給出,表格數據主要通過風洞試驗和飛行測試獲得。由于基準飛行狀態下,不考慮橫側向運動和大氣擾動的影響,同時將艦載機襟翼和起落架放出對氣動系數的影響考慮在內。下面僅給出縱向氣動系數的表示:

式中:δe為全動平尾偏角為鴨翼偏角q為俯仰角速度,rad/s;為迎角變化率,rad/s。

1.2 艦載機6自由度運動方程

完整的6自由度運動方程包含質心平動方程和繞質心轉動方程兩部分,每個部分又由動力學方程和運動學方程所組成。本文將方程都建立在機體坐標系下,形成“機體-機體體系”下的閉合方程組。為表達簡明,定義符號(K)n表示矢量K在坐標系Sn中的坐標列陣,定義符號Lξη表示從坐標系Sη到坐標系Sξ的變換矩陣。

若將機體系的角速度矢量表示為 ωb,則有(ωb)b=[p q r]T,p和r分別代表滾轉角速度和偏航角速度。艦載機質心運動動力學方程可寫為

式中:vk為航速, vk( )b=[vxbvybvzb]T;P、A、G分別為艦載機所受發動機推力、空氣動力和重力。

質心運動學方程把艦載機質心的位置變化率與速度相聯系:航速矢量的坐標列陣由本體系變換到大地系的轉換式為式(3),由航跡系變換到地面系的轉換式為式(4),且將式(4)展開后可得到計算航跡角γ與航向角χ的算式。設風速vw在大地系坐標為(vw)g=[wxgwygwzg]T,可得空速U在機體系中的表達式(5),進而可得空速的大小、迎角α和側滑角β的算式(6)。

繞質心轉動的動力學方程式(7)描述了剛體轉動的角速度與所受合外力矩的關系,其中:H為剛體動量矩,I為剛體慣性張量,M為作用于剛體上的外力對質心力矩的合矢量。根據艦載機結構對稱性可將 I()b化簡為對稱矩陣,元素為慣量矩Ix、Iy、Iz和慣量積Izx,kg·m2。

繞質心轉動的運動學方程式(8)描述了機體系角速度和艦載機姿態角之間的關系。其中,θ為俯仰角,φ為滾轉角,ψ為偏航角。

式(1)~(8)聯立起來是封閉的,如圖1所示,圖中δr和δa分別表示方向舵偏角和副翼偏角。

圖1 機體坐標系下的艦載機6自由度運動方程Fig.1 6-DOF motion equations of the carrier-based aircraft in body coordinate system

1.3 艦載機的基準狀態配平

設置艦載機的操縱面偏角和油門使其運動處于平衡狀態,艦載機的升力系數、阻力系數和俯仰力矩系數在配平條件下滿足:式中:下標“*”表示變量在基準狀態的取值,下同。

由于非線性仿真時,氣動系數是在對若干個數據表格差值后再線性組合得到的,因此理論求解準確的配平點難度較大,采取迭代計算的方法,求出不同程度滿足式(9)的初始輸入值。表1給出3組實現近似配平數據,精度依次提高。

第3組初始輸入使艦載機維持基準狀態附近穩定飛行的時間最長,各項性能指標在仿真開始后15 s內的誤差限制在2%以內,隨后響應迅速發散,如圖2所示。造成艦載機難以持久穩定飛行的根源是其本體特性的俯仰靜不穩定性,在仿真中表現為與平衡狀態存在微弱偏離的初始輸入會使誤差隨著時間的積累,最終導致響應發散。有關艦載機著艦狀態的縱向靜不穩定性的詳細驗證和分析可參閱文獻[13-14]。此外,3組近似配平點均在小擾動線性化范圍內,因此由它們導出的線性模型將相差甚微。

表1 艦載機的近似配平數據Table 1 Approximate trim data of the carrier-based aircraft

圖2 艦載機在第3組近似配平輸入下的時域響應Fig.2 Time-history response of the carrier-based aircraft under the third approximate trim values

2 艦載機著艦縱向線性建模

2.1 平靜大氣下縱向運動方程的線性化

2.1.1 縱向力和力矩的線性化

基于小擾動線性化的基本原理和方法[12,15]及研究對象的發動機推力特性,可將推力的擾動表示為

由于著艦階段艦載機的高度變化不大且速度基本恒定,因此空氣密度ρ變化很小,音速變化很小,這樣馬赫數基本不變(約為0.2Ma),因此空速(或航速)v對推力的影響可以忽略;此外,雖然艦載機的高度變化對其發動機推力確有影響,但經過后續的線性化處理,該影響將變得很微弱,亦可以忽略。這樣,在艦載機著艦過程中,發動機推力的變化僅與油門開度相關。

高度的變化通過改變空氣密度間接影響氣動力和氣動力矩,它們對高度的偏導數表示為

由標準大氣計算公式可得[11]

指定基準高度和空氣密度后,可得(?ρ/?h)*;此外,由式(10)可知 (?M/?h)*為 0。

根據前面提到的氣動力、力矩計算公式和式(1),可寫出氣動阻力的偏量表達式為

將式(14)中的“D”換為“L”,即為氣動升力的擾動表達式形式:

同理,得到俯仰力矩的擾動形式:

2.1.2 縱向運動方程的線性化

基準運動狀態下,艦載機橫側向通道的物理量均為零,據此可將線性化后的非線性運動方程中只包含縱向通道變量的等式組成縱向小擾動方程組,它又分為2類:一類為耦合方程組(19),另一類為解耦后的方程(20)[12]:

式(19)的前兩式分別由質心動力學方程在航跡系中的表達式(21)的第1、3式線性化而來;后兩式分別由式(7)展開后的第2式和式(3)線性化而來。式(20)的前兩式分別由式(3)、式(8)的第3式線性化而來,第3個等式則根據無風情況下鉛垂平面內對稱飛行的恒等關系[13]。線性方程中各物理量的單位:力為N,舵面偏轉角及油門開度為(°),其他角度和角速度的單位均為弧度制。

將式(20)代入式(19)中進行化簡;為使結果便于表示,引入微分算子符號“K”并定義如下所示的穩定性和控制導數,則推得式(22);再將式(22)的第2式代入第3式,化簡整理得到以狀態空間方程(23)、(24)描述的艦載機縱向小擾動模型,式中Δnz為過載擾動,表示升力擾動與重力之比。

2.1.3 確定狀態空間方程的參數

確定狀態矩陣中的元素需要先求出氣動力和力矩的偏導數,包括氣動系數的基準值和在基準點處的導數值。氣動系數的基準值可通過氣動數據表格在基準狀態處插值獲得,結果列在表2中。

對式(1),求氣動系數偏導數的步驟為:先將系數分量的表格數據做多項式擬合;然后求擬合多項式在基準點處的導數值;最后根據式(14)、(16)和(18)得到結果。要求擬合曲線在基準點及其附近的值與表格數據的誤差盡可能小且過渡平滑,詳細的處理結果見表3。將求出的穩定性和控制導數代入式(23)、(24),即得描述艦載機在平靜大氣下著艦的縱向模型(25)、(26)。

表2 氣動系數在基準狀態的值Table 2 Values of aerodynamic force and moment coefficients at trim condition

表3 氣動系數分量在基準狀態的導數值Table 3 Derivatives of aerodynamic coefficient components at trim condition

2.1.4 線性模型對控制指令響應的精度驗證

驗證線性擾動模型對非線性模型的近似精度,具體方法是:在相同的偏差指令輸入下,將非線性模型的輸出相對其基準狀態的偏離量同線性擾動模型的對應項相比較。鑒于艦載機本體特性的靜不穩定性和小擾動方法的局限性,設定的2組仿真條件為:仿真時長6 s,在2 s時刻向升降舵引入半階躍信號指令;仿真時長8 s,在2 s時刻向油門加入單位階躍信號。對比結果如圖3、4所示。

仿真結果表明,線性模型反映了非線性模型的不穩定性,且兩類模型的對應項輸出在變化趨勢和幅度上基本吻合,特別是在指令加入后約4 s內,之后線性系統的誤差不斷變大,主要原因是線性化的精度只在基準狀態附近有保證。

圖3 非線性和線性模型對升降舵階躍指令的仿真結果對比Fig.3 Comparisons of the simulated linear and nonlinear model responses to the elevator step signal command

圖4 非線性和線性模型對油門階躍指令的仿真結果對比Fig.4 Comparisons of the simulated linear and nonlinear model responses to the throttle step signal command

2.2 非平靜大氣下艦載機縱向運動方程的線性化

2.2.1 理論分析

艦尾氣流場是造成著艦偏差的主要因素之一,集中于著艦階段的最后0.5英里(約846 m)對飛行產生影響[10];軍標 MIL-HDBK-1797將其分成 4個部分:海面自由大氣紊流分量、艦尾流穩態分量、艦尾流周期性分量和隨機分量[16]。若設甲板風為利于著艦的情況—大小27 kn(約14 m/s),方向大致沿航母斜角甲板中線[10],則結合軍標對艦尾流模型的規定和本文的仿真假設,得到作用在艦載機上的尾流場的時域仿真結果,如圖5所示。

通過仿真艦尾流的過程發現:甲板風越大,尾流場強度越大,其位于艦載機垂直對稱面內的分量隨艦機距離的減小呈現變大的趨勢,且水平分量的強度明顯大于垂直分量。據此,可將影響艦載機縱向運動的艦尾紊流視為垂直和水平方向上連續分布的強度各異的突風,嘗試把由突風帶來的附加力和力矩作為干擾項引入已經導出的平靜大氣下的艦載機線性動力學模型。

圖5 艦尾流場時域仿真Fig.5 Time-history simulation of ship burble

設垂直方向存在下洗流(wzg>0),它在短時間內對飛行的影響如圖6所示。艦載機的航跡傾角變大且有式(27)成立,式中慣性迎角αd表示航速矢量與機體系軸xb的夾角,迎角和空速的變化過程難以直觀判斷。由式(27)可得式(28)。

圖6 垂直突風對飛行的影響Fig.6 Influence of vertical gust on flight

若紊流方向與進場方向平行(wxg>0),則在氣動效果上相當于艦載機空速降低,有式(29)滿足,式中的Δv代表空速擾動,Δvk代表航速擾動。聯系圖5的仿真結果和艦載機氣動力、力矩計算公式,可以估算出艦載機升力和阻力的變化幅度有幾千甚至上萬牛頓,對于控制精度要求很高的著艦過程,這是值得重視的要素。因此在線性動力學建模時只考慮流場的垂直分量的作用效果或者認為水平分量和垂直分量一起只對(慣性)迎角產生突出影響的做法都欠妥。

綜上分析,將平靜大氣下的線性模型“改造”為風場中動力學模型的方法是:將式(28)第1式代入式(19),除了速度微分項不變外,式(19)的Δv均由式(29)取代;用式(28)第2式取代式(20)第3式;采用與前面推導線性模型完全類似的做法,即得描述艦載機在艦尾流場中著艦的縱向動力學模型(30)~(33)。

2.2.2 線性模型對紊流響應的精度驗證

驗證風場下的線性擾動模型對非線性模型的近似程度。根據前面設置仿真條件的經驗,并考慮艦尾流場的特點,截取圖5在12~16 s的艦尾流前向分量和垂直分量的仿真數據作為模型擾動輸入,得到3組對比結果,如圖7~9所示。

圖7 非線性和線性模型對艦尾流垂直分量的響應Fig.7 Simulated linear and nonlinear model responses to the vertical component of ship burble

圖7表明,加入擾流垂直分量后的2.5 s內,除航速以外,線性模型同非線性模型的各對應輸出響應均很接近。研究線性模型在航速變化方面誤差較大的原因發現,這是本文針對艦載機著艦動力學所采取的線性化方法的固有缺陷。一個簡便的驗證方案是:數值仿真環境下,用脈寬為0.5 s的單位脈沖信號代表穩態垂直擾流(下洗風),在1 s時加入2類模型中;結果顯示,脈沖存在時,2類模型唯獨在航速響應上呈現出完全相反的趨勢,如圖10(a)所示;對于非線性模型,依據圖6,分析航速方向上的受力,仿真結果如圖10(b)所示,并由此可見突風的存在確實會引起航速減小,如圖10(a)中虛線所示;而對于線性模型,由式(30)可得

式中:除 xvΔvk外的各項均為正;雖然xv< 0,但xv數量級太小,因此線性模型的慣性速度為正。

圖8 非線性和線性模型對艦尾流水平分量的響應Fig.8 Simulated linear and nonlinear model responses to the horizontal component of ship burble

即便如此,圖9表明,補充艦尾流水平分量后,該“缺陷”明顯減弱,至少在受擾后的1.5 s,線性模型各項輸出指標的誤差較小。綜合艦載機本體特性的靜不穩定性并結合前面的仿真結果可知,本文的小擾動線性化過程是合理的,線性模型符合在著艦控制系統仿真中使用的飛機動力學模型的要求。此外,從很大程度上,圖7中的實線課被視為僅考慮了艦尾流垂直分量影響的線性動力學模型在完整艦尾流作用下的響應情況,這樣與圖9的虛線對比后可知:對艦載機著艦的縱向線性建模必須重視紊流在水平方向上的分量。

圖9 非線性和線性模型對艦尾流的響應Fig.9 Simulated linear and nonlinear model responses to ship burble

圖10 線性模型的缺陷驗證Fig.10 Verification on the flaw of linear model

3 結束語

本文在建立艦載機6自由度非線性模型的基礎上,應用小擾動理論,完成了對艦載機著艦縱向動力學的線性建模。數值仿真環境下,動力學模型對控制指令和艦尾流場的響應表明,在基準狀態附近,線性模型相對非線性模型具有較高精度,可作為在著艦仿真中使用的飛機動力學模型。

理論分析以及線性、非線性對比仿真都證明,影響著艦過程的紊流水平分量會對艦載機穩定性產生顯著干擾,對線性模型的精度起決定作用,因此是艦載機縱向線性建模所必須考慮的因素。

[1]DURAND T S,WASICKO R J.Factors influencing glide path control in carrier landing[J].Journal of Aircraft,1967,4(2):146-158.

[2]DURAND T S,WASICKO R J.An analysis of carrier landing[C]//Aircraft Design and Technology Meeting.Los Angeles,1965.

[3]DAVIES W D T,NOURY R.AN/SPN-42 automatic carrier landing system[C]//Proceedings of the First Annual Advanced Control Conference.Lafayette,USA,1974:99-110.

[4]RUDOWSKY T,COOK S,HYNES M,et al.Review of the carrier approach criteria for carrier-based aircraft.NAWCADPAX/TR-2002/71[R].[S.l.]:Naval Air Warfare Center Aircraft Division,2002.

[5]彭兢.艦載艦載機進艦著艦的自動引導和控制研究[D].北京:北京航空航天大學,2001:167-170.PENG Jing.Research on the automatic guide and control of a carrier-based airplane approach and landing[D].Beijing:Beihang University,2001:167-170.

[6]鄧娟.艦載艦載機自動著艦縱向控制系統設計的理論與仿真研究[D].上海:復旦大學,2010:30-53.DENG Juan.Theory and simulation for the design of longitudinal automatic carrier landing system[D].Shanghai:Fudan University,2010:30-53.

[7]劉冰.基于LQG/LTR方法的艦載艦載機自動著艦系統設計的理論及仿真研究[D].上海:復旦大學,2010:25-41.LIU Bing.Research on theory and simulation for the design of automatic carrier landing system based on LQG/LTR[D].Shanghai:Fudan University,2010:25-41.

[8]吳鑫.艦載機自動著艦控制系統設計與研究[D].哈爾濱:哈爾濱工程大學,2009:8-20.WU Xin.Design and simulation for the automatic carrier landing control system[D].Harbin:Harbin Engineering U-niversity,2009:8-20.

[9]張明廉,徐軍.艦載艦載機自動著艦系統的研究[J].北京航空航天大學學報,1994,20(4):386-391.ZHANG Minglian,XU Jun.Research on the automatic carrier landing system[J].Journal of Beijing University of Aeronautics and Astronautics,1994,20(4):386-391.

[10]SCHUST A P,YOUNG P N,SIMPSON W R.Automatic carrier landing system(ACLS)categoryⅢcertification manual.AD-A1181817[R].[S.l.]:ARINC Research Corporation,1982.

[11]肖業倫.飛行器運動方程[M].北京:航空工業出版社,1987:12-35,62-81,90-103.

[12]方振平,陳萬春,張曙光.航空飛行器飛行動力學[M].北京:北京航空航天大學出版社,2012:16-31,174-214.

[13]ROBERT C N.Fligt stability and automatic control[M].2nd ed.New York:McGraw-Hill,1997:42-53.

[14]董然.ACLS縱向內回路控制系統研究[D].哈爾濱:哈爾濱工程大學,2013:17-21.DONG Ran.Research on the longitudinal inner loop configuration in ACLS[D].Harbin:Harbin Engineering University,2013:17-21.

[15]錢杏芳,林瑞雄,趙亞男.導彈飛行力學[M].北京:北京理工大學出版社,2011:158-162,179-184.

[16]MIL-HDBK-1797,Department of defense interface standard,flying qualities of piloted aircraft[S].[S.l.]:Department of Defense,1997.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 日本高清在线看免费观看| 真人免费一级毛片一区二区| 国产农村妇女精品一二区| 欧美综合区自拍亚洲综合天堂| 亚洲综合日韩精品| 一级毛片基地| 无码中文字幕乱码免费2| 亚洲嫩模喷白浆| 爱色欧美亚洲综合图区| 免费可以看的无遮挡av无码| 亚洲日韩日本中文在线| 国产成人高精品免费视频| 四虎国产在线观看| 国产拍在线| 国产欧美在线观看一区| 欧美伊人色综合久久天天| 国产成人麻豆精品| 91视频青青草| 亚洲swag精品自拍一区| 玖玖精品在线| 色婷婷在线播放| 手机在线免费不卡一区二| 欧美日韩午夜| 97国产一区二区精品久久呦| 欧美在线观看不卡| 免费又爽又刺激高潮网址| 热热久久狠狠偷偷色男同| 国产成人精品日本亚洲| 伊人久久影视| 欧美色伊人| 国产亚洲成AⅤ人片在线观看| 成人第一页| 青青草一区| 妇女自拍偷自拍亚洲精品| 久爱午夜精品免费视频| 国产99视频精品免费观看9e| 天天色天天综合| 成人免费午夜视频| 色婷婷啪啪| 亚洲国产理论片在线播放| 日韩欧美中文字幕在线精品| 伊人婷婷色香五月综合缴缴情| 亚洲A∨无码精品午夜在线观看| 特级精品毛片免费观看| 久久久久久久久久国产精品| 青青草a国产免费观看| 免费久久一级欧美特大黄| 亚洲成a人片77777在线播放| 四虎精品黑人视频| 高清不卡一区二区三区香蕉| 黄色福利在线| 久久大香伊蕉在人线观看热2| 久久精品丝袜高跟鞋| 色哟哟国产精品一区二区| 日本五区在线不卡精品| 亚洲综合18p| 亚洲无码日韩一区| 午夜丁香婷婷| 在线国产毛片| 成人国内精品久久久久影院| 欧美高清视频一区二区三区| 五月激情综合网| 精品无码人妻一区二区| 国产老女人精品免费视频| 欧美视频二区| 波多野结衣在线se| 特级毛片8级毛片免费观看| 国产无码高清视频不卡| 国产成人综合网| 色妞永久免费视频| 亚洲最新地址| 久久大香伊蕉在人线观看热2| 精品色综合| 999国内精品视频免费| www.99在线观看| 久久人搡人人玩人妻精品| 老熟妇喷水一区二区三区| 成人一级免费视频| 无码中文字幕乱码免费2| 在线视频97| 亚卅精品无码久久毛片乌克兰| 三级欧美在线|