亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于單純形優(yōu)化算法的火箭彈增程研究

        2021-09-01 09:33:20陳國光
        關(guān)鍵詞:優(yōu)化

        張 斌,陳國光,王 捷

        (1.中北大學(xué),山西太原 030051;2.西南技術(shù)物理研究所,四川成都 610041)

        為滿足當(dāng)前遠(yuǎn)距離、精確化、大威力的現(xiàn)代戰(zhàn)爭(zhēng)作戰(zhàn)要求,火箭彈增程技術(shù)當(dāng)前已成為當(dāng)前研究的熱點(diǎn)之一。目前,火箭彈增程方法有提高發(fā)動(dòng)機(jī)總沖、改良火箭彈結(jié)構(gòu)以及增添舵機(jī)控制的滑翔增程等。單就結(jié)構(gòu)和推進(jìn)劑的改良而言,增程效果并不佳。依據(jù)彈道修正彈彈道特性分析,采用舵機(jī)控制的滑翔增程方法不僅能有效增加射程,還能減少因火箭彈增大射程而帶來的散布問題。

        滑翔增程是通過調(diào)整舵翼偏轉(zhuǎn)來控制火箭彈的升力,使火箭彈緩慢下降,進(jìn)而延長(zhǎng)飛行時(shí)間,達(dá)到增程的目的。這種增程方法是基于彈箭飛行特性提出的,對(duì)制導(dǎo)控制及彈道修正有著重要指導(dǎo)意義。經(jīng)查閱文獻(xiàn)了解,目前增程方案有:分析增程原理和彈道特點(diǎn),控制重力法向方向始終平衡的方案;采用控制火箭彈法向過載最優(yōu)的方案;采用數(shù)值解算最大升阻比的思想建立滑翔段模型;以能量法研究滑翔彈道,最小能量損失的彈道設(shè)計(jì);對(duì)滑翔段進(jìn)行分段控制,保證舵偏角在每1小段都不變,探究分段次數(shù)與增程的關(guān)系;依據(jù)氣動(dòng)力參數(shù),擬合最大升阻比下的攻角隨馬赫數(shù)的變化曲線及攻角隨時(shí)間的函數(shù),實(shí)現(xiàn)滑翔增程;以射程為優(yōu)化變量,采用序列二次規(guī)劃優(yōu)化算法,進(jìn)行數(shù)值仿真;采用現(xiàn)代優(yōu)化算法,對(duì)彈道進(jìn)行全局約束搜索最優(yōu)攻角方案。

        目前,多采用地面雷達(dá)監(jiān)測(cè)及發(fā)射指令的方式修正彈道,對(duì)于高速短時(shí)作戰(zhàn)的彈藥而言,這種方法存在時(shí)間誤差和預(yù)測(cè)偏差等問題,修正效果并不理想。

        基于此,本文從能量法角度分析并提出了1種基于彈載計(jì)算機(jī)實(shí)時(shí)優(yōu)化的改進(jìn)單純形優(yōu)化算法,在SIMULINK 中建立彈道模型,并將優(yōu)化算法模塊化,嵌套在仿真程序中,跟隨彈道解算過程對(duì)攻角進(jìn)行實(shí)時(shí)優(yōu)化。通過該算法,能快速找到當(dāng)前彈道諸元(速度、高度、彈道傾角等)下的最優(yōu)控制量。通過仿真分析驗(yàn)證了算法有效性。

        1 模型建立

        本文的研究對(duì)象是某型制導(dǎo)火箭彈,其采用雙通道舵機(jī),卷弧尾翼。三維彈體模型,如圖1所示。

        圖1 制導(dǎo)火箭彈模型Fig.1 Guided rocket model

        對(duì)于雙通道控制舵機(jī),升降舵和方向舵是分開控制的:升降舵主要控制火箭彈在縱向平面的射程修正;方向舵則控制水平面的偏差量修正。本文主要研究火箭彈在縱向平面的增程問題。

        1.1 滑翔增程彈原理

        制導(dǎo)火箭彈通常在彈道前期無控飛行,飛行到彈道頂點(diǎn)附近打開舵機(jī),通過電機(jī)驅(qū)動(dòng)舵翼偏轉(zhuǎn),調(diào)整彈體的升力,減小或平衡彈丸在速度矢量法向方向的重力分量,從而延長(zhǎng)飛行時(shí)間,增加飛行距離。

        在彈道縱向平面內(nèi),火箭彈沿速度矢量方向的受力分析,如圖2所示。

        圖2 火箭彈在速度矢量方向的受力分析Fig.2 Force analysis of rocket in the direction of velocity vector

        從能量角度分析,因高度降低所產(chǎn)生的勢(shì)能轉(zhuǎn)化成動(dòng)能和克服空氣的阻力,動(dòng)能的增加勢(shì)必加快高度的降低,這對(duì)增程而言是不利的。因此,將勢(shì)能盡可能地用來克服空氣阻力,維持動(dòng)能不變或緩慢減小是1種增加射程的有效方法。

        由圖2受力分析可知,速度是由切向合力控制,高度的降低是由空氣作用力和重力兩者的合力控制。控制速度變化的切向合力最小,以及控制彈道傾角變化的法向合力最小,即可達(dá)到緩慢下降及速度緩慢減小的目的。將上述理論轉(zhuǎn)化為數(shù)值表達(dá)式,為更易實(shí)現(xiàn)這一方法,采用標(biāo)量的合力之和來替換矢量合力。為避免各合力出現(xiàn)負(fù)值的情況,對(duì)各合力做絕對(duì)值運(yùn)算。整合以上理論,如式(1):

        式(1)中:

        X

        為空氣阻力;

        Y

        為空氣升力;

        m

        為火箭彈質(zhì)量;

        θ

        為彈道傾角;

        g

        為當(dāng)?shù)刂亓铀俣?,與高度

        y

        有關(guān)。

        式(2)中:

        ρ

        為空氣密度;

        X

        為空氣阻力;

        Y

        為空氣升力;

        S

        為彈丸特征面積;

        C

        為零升阻力系數(shù);

        C

        為誘導(dǎo)阻力系數(shù);

        C

        ′為升力系數(shù)導(dǎo)數(shù);

        α

        為攻角?;鸺龔検芰﹄S飛行位置不斷發(fā)生變化,受飛行穩(wěn)定性及修正裝置等的限制,此處的合力最小為優(yōu)化問題,須通過優(yōu)化方法獲取。為此,本文提出實(shí)時(shí)優(yōu)化的算法以解決這一問題。將優(yōu)化算法嵌套在彈載計(jì)算機(jī)程序中,每隔

        T

        時(shí)刻,接收并處理火箭彈的位置及速度信息,求解出采樣時(shí)刻的最優(yōu)攻角值,再轉(zhuǎn)化成修正指令,控制舵翼偏轉(zhuǎn),實(shí)現(xiàn)增程目的。

        1.2 彈道模型及方程建模

        通過對(duì)火箭彈增程原理分析,結(jié)合原火箭彈特點(diǎn),總結(jié)滑翔增程火箭彈彈道特點(diǎn)如下。

        1)主動(dòng)段。該段彈丸受重力、空氣動(dòng)力及發(fā)動(dòng)機(jī)推力作用,受推進(jìn)劑燃?xì)馔苿?dòng)作用,彈丸做加速運(yùn)動(dòng)。發(fā)動(dòng)機(jī)停止工作后,彈丸達(dá)到最大速度。

        2)被動(dòng)段。該段彈丸僅受空氣動(dòng)力和自身重力作用,與無控彈的彈道解法一致。

        3)控制段。在彈丸到達(dá)彈道頂點(diǎn)附近時(shí),打開舵機(jī)。在電機(jī)驅(qū)動(dòng)下,舵翼按照舵機(jī)指令進(jìn)行偏轉(zhuǎn),空氣作用力隨之發(fā)生改變,進(jìn)而達(dá)到調(diào)整姿態(tài)、修正彈道的作用。

        根據(jù)以上彈道特點(diǎn)描述,飛行彈道特性模擬圖,如圖3所示。

        圖3 滑翔增程火箭彈彈道模擬圖Fig.3 Ballistic simulation diagram of glide extended range rocket

        依據(jù)前文分析,本文只研究有控彈箭在縱向平面的運(yùn)動(dòng),由于采用卷弧尾翼,彈丸旋轉(zhuǎn)速度會(huì)降低,飛行至控制段時(shí)近似看作無旋運(yùn)動(dòng),即忽略彈箭繞質(zhì)心的轉(zhuǎn)動(dòng),只研究彈箭質(zhì)心的運(yùn)動(dòng)。由此,縱向平面可操縱質(zhì)點(diǎn)方程組如下:

        2 優(yōu)化算法建模

        2.1 基于單純形的攻角優(yōu)化建模

        單純形法(Simplex)是1962 年由Spendley、Hext和Himsworth 提出,再由Nelder 和Mead 于1965 年改進(jìn),故又稱為NM 算法。它是1 種不使用導(dǎo)數(shù)求解無約束極小化問題的直接搜索方法。相比最優(yōu)化方法中的解析法,此法僅需要相關(guān)函數(shù)值信息,就可尋找最優(yōu)解。不同于其他優(yōu)化的算法,它從最初構(gòu)造的單純形,不斷替換最大值點(diǎn),構(gòu)造更適合尋優(yōu)的新單純形?;驹砣缦拢?/p>在

        n

        維空間中,以

        n

        +1 個(gè)頂點(diǎn)構(gòu)成一個(gè)單純形(

        x

        ,

        x

        ,

        x

        ,…,

        x

        )。并將單純形中各頂點(diǎn)代入到目標(biāo)函數(shù)中求解函數(shù)值。記各函數(shù)值中最小值為

        f

        (

        x

        ),最大值為

        f

        (

        x

        )。并求除最大值以外各點(diǎn)的重心,及最大值頂點(diǎn)關(guān)于重心的反射點(diǎn),依據(jù)反射點(diǎn)函數(shù)值的大小選擇延伸、收縮或壓縮,得到新的較好頂點(diǎn)并替換原單純形中最大值頂點(diǎn),構(gòu)造新的單純形。通過反復(fù)求解較小點(diǎn)來替換新單純形中最大值頂點(diǎn),迭代單純形,逼近極小值點(diǎn)。單純形中重心、反射點(diǎn)、延伸點(diǎn)、收縮點(diǎn)、壓縮點(diǎn)以及精度的計(jì)算公式如下:

        式(4)中:

        x

        為單純形頂點(diǎn);

        x

        為最大值頂點(diǎn);

        x

        為最小值頂點(diǎn);

        x

        為重心;

        x

        為反射點(diǎn);

        x

        為延伸點(diǎn);

        x

        為收縮點(diǎn);

        x

        為壓縮點(diǎn);

        h

        為反射系數(shù);

        β

        為延伸系數(shù);

        γ

        為收縮系數(shù);

        δ

        為壓縮系數(shù);

        ε

        為計(jì)算精度。

        2.1.1 選取設(shè)計(jì)變量

        依據(jù)公式(1)分析,為獲取采樣時(shí)刻最優(yōu)攻角,實(shí)時(shí)彈道優(yōu)化以采樣時(shí)刻彈道諸元作為優(yōu)化輸入量,以平衡攻角為優(yōu)化變量。隨采樣時(shí)刻的更新進(jìn)行實(shí)時(shí)優(yōu)化。設(shè)計(jì)變量

        Γ

        見式(5)。

        式(5)中:

        y

        為火箭彈采樣時(shí)刻高度;

        v

        為采樣時(shí)刻速度;

        θ

        為采樣時(shí)刻彈道傾角;

        ρ

        為采樣時(shí)刻空氣密度;

        g

        為采樣時(shí)刻的重力加速度;

        α

        為采樣時(shí)刻待優(yōu)化的攻角。

        2.1.2 建立目標(biāo)函數(shù)

        基于公式(1)分析,將采樣時(shí)刻的位置和速度信息作為優(yōu)化輸入量,以采樣時(shí)刻合力最小值為目標(biāo)函數(shù)。結(jié)合式(1),目標(biāo)函數(shù)

        F

        見式(6):

        2.1.3 添加參數(shù)約束

        攻角的取值,取決于彈箭物理參數(shù)及舵機(jī)的修正能力,更重要的是不能影響飛行穩(wěn)定性,這限制了攻角只能在一定范圍內(nèi)變化,見式(7):

        式(7)中:

        α

        為攻角最小值;

        α

        為攻角最大值。

        2.1.4 算法終止條件

        單純形算法優(yōu)點(diǎn)是尋優(yōu)效率高,能快速迭代到最優(yōu)解附近;缺點(diǎn)是在最優(yōu)點(diǎn)附近收斂速度較慢。而實(shí)時(shí)修正的核心即時(shí)效性,短時(shí)多修,因此對(duì)于精度要求并不是很高。為了更好滿足實(shí)時(shí)解算要求,縮短算法終止時(shí)間,本文在原算法基礎(chǔ)上再增添終止條件,以連續(xù)3 次輸出迭代優(yōu)化值之差小于迭代誤差

        ε

        作為終止條件。雙重終止條件能有效縮短迭代時(shí)間。迭代誤差

        ε

        見公式(8):

        式(8)中:

        ε

        為迭代誤差;

        x

        為第

        i

        次迭代結(jié)果;

        x

        為第

        i

        -1次迭代結(jié)果。

        2.1.5 單純形法優(yōu)化攻角流程

        1)設(shè)定控制段采樣時(shí)刻

        T

        =0.1 s,獲取采樣時(shí)刻的火箭彈位置和速度信息。統(tǒng)一設(shè)定待優(yōu)化的攻角初始值為0。并設(shè)定優(yōu)化函數(shù)中的參數(shù)值:初始步長(zhǎng)、反射系數(shù)、收縮系數(shù)、壓縮系數(shù)及計(jì)算精度。

        2)將初始攻角值及采用初始步長(zhǎng)獲取的新攻角值代入插值函數(shù)中,分別插值計(jì)算出當(dāng)前攻角及馬赫數(shù)下的氣動(dòng)力參數(shù),并將氣動(dòng)力參數(shù)值代入目標(biāo)函數(shù)中,求解當(dāng)前目標(biāo)函數(shù)值。

        3)將2)中得到的目標(biāo)函數(shù)值返回到優(yōu)化算法中,確定最小值頂點(diǎn)

        x

        及最大值頂點(diǎn)

        x

        ,并計(jì)算重心、反射點(diǎn)。4)通過對(duì)反射點(diǎn)的函數(shù)值求解,選擇延伸處理、收縮處理或壓縮處理。使目標(biāo)函數(shù)更小的點(diǎn)為新的單純形頂點(diǎn),返回并替換步驟3)中的

        f

        ( )

        x

        ,重新確定新單純形的最大值頂點(diǎn)及最小值頂點(diǎn),重復(fù)以上步驟,反復(fù)迭代新單純形,直至滿足計(jì)算精度為止。

        5)將優(yōu)化算法得到的攻角值與式(6)比對(duì),將不在攻角設(shè)定范圍的攻角值進(jìn)行輸出限制或替換。最后輸出攻角,作為當(dāng)前火箭彈飛行的控制量。直至得到下1 個(gè)最優(yōu)攻角。通過這種交替優(yōu)化解算方法,獲取控制全段的最優(yōu)攻角值。

        2.2 SIMULINK建模仿真

        工程上常用C語言求解彈道,這對(duì)編程能力要求高。SIMULINK 是將MATLAB 的相關(guān)函數(shù)進(jìn)行模塊化封裝,對(duì)使用者來說要求低、易上手,對(duì)于所有設(shè)定的參數(shù)均可實(shí)時(shí)讀取,通用性高,可視性強(qiáng),已成為學(xué)者和工程師求解彈道的常用軟件。本文采用SIMULINK 仿真平臺(tái)進(jìn)行滑翔增程火箭彈模型建立與仿真分析。

        依據(jù)滑翔增程彈的彈道特性,整個(gè)彈道模型由彈道初始模塊、氣象標(biāo)準(zhǔn)模塊、氣動(dòng)力模塊、彈道方程組模塊、優(yōu)化算法模塊組成。各模塊內(nèi)部以子模塊或函數(shù)組成?;柙龀袒鸺龔棌椀澜馑憧驁D,如圖4所示。

        圖4 滑翔增程火箭彈彈道解算框圖Fig.4 Diagram of ballistic calculation for glide extended range rocket

        2.2.1 彈道初始化模塊

        在SIMULINK中設(shè)置滑翔增程火箭彈的參數(shù),見表1。

        表1 火箭彈參數(shù)Tab.1 Rocket parameters

        2.2.2 NM優(yōu)化算法模塊

        依據(jù)算法,將函數(shù)的所有輸入的參數(shù)封裝到函數(shù)模塊中,該函數(shù)模塊包括3 個(gè)函數(shù):氣動(dòng)力插值函數(shù),目標(biāo)函數(shù)及單純形優(yōu)化函數(shù)。NM子模塊中的單純形優(yōu)化算法函數(shù)模型,見圖5所示。

        圖5 單純形優(yōu)化算法函數(shù)模塊Fig.5 Function module of simplex optimization algorithm

        本文火箭彈模型,依據(jù)舵機(jī)修正能力及飛行穩(wěn)定性,設(shè)定攻角取值范圍為[0 °,8°] 。單純形優(yōu)化算法是一種無約束優(yōu)化,對(duì)于優(yōu)化結(jié)果超出優(yōu)化變量約束范圍的,將以攻角最值輸出。也就是優(yōu)化結(jié)果大于8°則取8°,小于0°則取0°。

        2.2.3 氣動(dòng)力參數(shù)模塊

        本文采用DATCOM 軟件獲取滑翔增程火箭彈的氣動(dòng)參數(shù)。將建好的滑翔火箭彈模型的彈體參數(shù)(見表1)及舵機(jī)參數(shù)編制文件,設(shè)定不同馬赫數(shù)和攻角值,分別計(jì)算原火箭彈氣動(dòng)參數(shù),以及添加舵機(jī)后的滑翔火箭彈的氣動(dòng)參數(shù)。根據(jù)滑翔增程火箭彈彈道特點(diǎn),在主動(dòng)段和被動(dòng)段使用無舵翼氣動(dòng)力參數(shù),打開舵機(jī)后使用有舵翼的氣動(dòng)力參數(shù)。

        3 仿真結(jié)果及分析

        由于滑翔增程彈的舵機(jī)是在頂點(diǎn)附近打開的,因此彈道的主動(dòng)段和被動(dòng)段均與原火箭彈相同。忽略無控段受外部隨機(jī)干擾造成的攻角變化,采用零攻角仿真無控段。本文采用SIMULINK 自帶算法(ode45),四階龍格庫塔解算彈道。

        設(shè)定滑翔增程火箭彈初始條件是:

        θ

        =50°,

        v

        =50 m/s,

        x

        =0,

        y

        =0。將初始條件設(shè)置到SIMULINK模塊中,設(shè)置積分步長(zhǎng)與采樣時(shí)刻同為0.1。設(shè)置程序終止條件為

        y

        <0。優(yōu)化函數(shù)算法中各參數(shù)取值分別為:初始步長(zhǎng)

        h

        =2,延伸系數(shù)

        β

        =2,收縮系數(shù)

        γ

        =0.5,壓縮系數(shù)

        δ

        =0.5,計(jì)算精度

        ε

        =0.05,迭代終止精度

        ε

        =0.001。

        考慮到火箭彈飛行中受到(風(fēng)等)擾動(dòng),以及彈載計(jì)算機(jī)獲取的位置坐標(biāo)誤差及地磁測(cè)量姿態(tài)誤差等因素,轉(zhuǎn)化各種誤差因素為攻角的響應(yīng)白噪,即在實(shí)時(shí)優(yōu)化后的攻角值上隨機(jī)添加[ ]-0.2°,0.2° 的偏差量。

        圖6為零攻角下原火箭彈的彈道曲線。圖7為有噪與無噪的優(yōu)化彈道對(duì)比圖。圖8為有噪和無噪的優(yōu)化攻角對(duì)比圖。圖9 為實(shí)時(shí)攻角下的優(yōu)化目標(biāo)函數(shù)值。圖10是在原單純形算法下,對(duì)不同采集時(shí)刻的彈道信息分別進(jìn)行實(shí)時(shí)優(yōu)化的迭代次數(shù)與優(yōu)化攻角的對(duì)比圖。

        圖6 零攻角火箭彈彈道Fig.6 Ballistic curve of the rocket at zero angle of attack

        圖7 實(shí)時(shí)優(yōu)化后的彈道Fig.7 Ballistic curve of the rocket after real-time optimization

        圖8 實(shí)時(shí)優(yōu)化后的火箭彈攻角Fig.8 Angle of attack of the rocket after real-time optimization

        圖9 實(shí)時(shí)優(yōu)化下的目標(biāo)函數(shù)值Fig.9 Objective function value under real-time optimization

        圖10 不同時(shí)刻下迭代次數(shù)與優(yōu)化結(jié)果Fig.10 Number of iterations and optimization results at different moments

        圖6、7 對(duì)比表明,在射角50°、零攻角下的原無控火箭彈的射程僅22 000 m。而采用本文提出的改進(jìn)單純形優(yōu)化算法解算實(shí)時(shí)攻角,無噪射程為38 344.87 m ,有噪射程為38 205.04 m ,偏差僅139.83 m,無噪下增程率達(dá)到74.3%。說明此優(yōu)化方法能有效增加火箭彈的射程。

        圖7~9表明,在相同初始條件下,無噪和有噪解算結(jié)果相差甚小,實(shí)時(shí)解算具有較高的抗噪性。實(shí)時(shí)解算的最優(yōu)攻角變化平穩(wěn),并無大的突變,對(duì)于舵機(jī)修正而言,易于實(shí)現(xiàn)。此改進(jìn)算法能有效地應(yīng)用于彈載計(jì)算機(jī)上。

        圖10 表明,雖然達(dá)到原算法的計(jì)算精度,但迭代次數(shù)較多,耗時(shí)較長(zhǎng),不利于快速多次的實(shí)時(shí)修飛理念。

        通過表2 的數(shù)據(jù)比對(duì),原算法實(shí)時(shí)求解的時(shí)間較長(zhǎng),從而使舵機(jī)修正的時(shí)間間隔延長(zhǎng)。而隨著修正間隔時(shí)間的增加,精度必然降低。采用改進(jìn)后的算法,迭代次數(shù)更少,而且優(yōu)化結(jié)果與原算法相差較小。改進(jìn)算法能有效地縮短間隔時(shí)間,增加修正次數(shù),更能滿足實(shí)時(shí)解算的要求。

        表2 原算法和改進(jìn)算法尋優(yōu)結(jié)果對(duì)比Tab.2 Comparison of optimization results between the original algorithm and the improved algorithm

        4 結(jié)論

        針對(duì)火箭彈增程問題,在理想滑翔狀態(tài)理論基礎(chǔ)上,分析并提出了基于改進(jìn)單純形算法的實(shí)時(shí)優(yōu)化攻角的增程方法,通過在SIMULINK仿真平臺(tái)上建立彈道模型并仿真分析,得出以下結(jié)論:

        1)與無控火箭彈射程相比,基于改進(jìn)單純形優(yōu)化算法的增程方案使得火箭彈增程率達(dá)到74.3%;

        2)改進(jìn)后的單純形優(yōu)化算法更適合彈載計(jì)算機(jī)實(shí)時(shí)優(yōu)化。

        猜你喜歡
        優(yōu)化
        超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
        PEMFC流道的多目標(biāo)優(yōu)化
        能源工程(2022年1期)2022-03-29 01:06:28
        民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
        關(guān)于優(yōu)化消防安全告知承諾的一些思考
        一道優(yōu)化題的幾何解法
        由“形”啟“數(shù)”優(yōu)化運(yùn)算——以2021年解析幾何高考題為例
        圍繞“地、業(yè)、人”優(yōu)化產(chǎn)業(yè)扶貧
        事業(yè)單位中固定資產(chǎn)會(huì)計(jì)處理的優(yōu)化
        4K HDR性能大幅度優(yōu)化 JVC DLA-X8 18 BC
        幾種常見的負(fù)載均衡算法的優(yōu)化
        電子制作(2017年20期)2017-04-26 06:57:45
        免费国产在线视频自拍白浆| 91精品国产91久久久久久青草 | 亚洲悠悠色综合中文字幕| 日本大片免费观看视频| 四虎国产精品永久在线无码| 亚洲欧美变态另类综合| 久久久噜噜噜久久熟女| 亚洲国产精品无码久久一线 | 亚洲av无码av吞精久久| 午夜精品久视频在线观看| 国产一区二区资源在线观看| 欧美国产激情二区三区| 国产乱人伦av在线a| 91福利精品老师国产自产在线| 极品少妇人妻一区二区三区| 少妇被又大又粗又爽毛片 | 欧美巨鞭大战丰满少妇 | 亚洲中文无码久久精品1| 国产精品一区二区三区黄片视频| 精品国产亚洲av麻豆| 久久人人爽人人爽人人片av东京热 | 国产精品一区二区韩国AV| 小草手机视频在线观看| 色偷偷久久久精品亚洲| 亚洲av无码不卡久久| 亚洲熟女av超清一区二区三区| 中文字幕日韩精品人妻久久久| 天堂中文а√在线| 国产av一区二区三区日韩| 天堂女人av一区二区| 公厕偷拍一区二区三区四区五区| 成人免费看www网址入口| 免费观看久久精品日本视频| 加勒比av在线一区二区| 乱中年女人伦av三区| 波多野结衣亚洲一区二区三区| av免费网站在线免费观看| 国产女人精品视频国产灰线| 亚洲va中文字幕无码| 九月色婷婷免费| 国产av一级黄一区二区三区|