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打印中的模型分割與打包
主站蜘蛛池模板: 国产91精品最新在线播放| 99久久国产综合精品女同| 在线看片国产| 伊人查蕉在线观看国产精品| 久久一级电影| 国产黄色免费看| 无码中文字幕精品推荐| 91午夜福利在线观看| 亚洲无码视频一区二区三区 | 国产成人精品综合| 2019年国产精品自拍不卡| 欧美综合成人| 欧美国产中文| 亚洲国产中文综合专区在| 久久国产精品电影| 国产人人乐人人爱| 亚洲成a人片在线观看88| 亚洲国产综合第一精品小说| 国产精品思思热在线| 亚洲日产2021三区在线| 国产成人高清精品免费软件| 精品乱码久久久久久久| 国产全黄a一级毛片| 日韩精品中文字幕一区三区| 成人va亚洲va欧美天堂| 欧美激情第一欧美在线| 风韵丰满熟妇啪啪区老熟熟女| 欧美国产在线精品17p| 国产精品成人第一区| 四虎永久免费地址在线网站| 欧美午夜在线播放| 一级毛片在线直接观看| 亚洲国产午夜精华无码福利| 欧洲av毛片| 免费看黄片一区二区三区| 国产视频你懂得| 国产精品私拍在线爆乳| 免费国产小视频在线观看| 亚洲人成网站观看在线观看| 久久精品这里只有国产中文精品| 第一区免费在线观看| 97影院午夜在线观看视频| 亚洲国产精品不卡在线| a毛片在线| 亚洲视频四区| 91午夜福利在线观看精品| 国产精品入口麻豆| 91麻豆国产视频| 亚洲综合专区| 无码专区第一页| 综合色区亚洲熟妇在线| 欧美亚洲激情| 四虎影视8848永久精品| 久久国产精品电影| 日韩在线观看网站| 国产乱人伦AV在线A| 暴力调教一区二区三区| 色妞www精品视频一级下载| 中文天堂在线视频| 久久中文字幕不卡一二区| 手机在线免费不卡一区二| 欧美精品成人一区二区视频一| 久草热视频在线| 国产原创第一页在线观看| 久久人人妻人人爽人人卡片av| 亚洲精品另类| 网友自拍视频精品区| 婷婷六月天激情| 国产凹凸视频在线观看 | 国产91高清视频| 午夜综合网| 色综合久久88| 人妻丰满熟妇AV无码区| 国产男女XX00免费观看| 婷婷中文在线| 亚洲成人精品在线| 国产福利拍拍拍| 无码福利视频| 9久久伊人精品综合| 亚洲六月丁香六月婷婷蜜芽| 看av免费毛片手机播放| 免费高清a毛片|