• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    脈沖燃燒風(fēng)洞測(cè)力系統(tǒng)動(dòng)態(tài)標(biāo)定方法

    2017-09-15 09:09:42樂嘉陵
    實(shí)驗(yàn)流體力學(xué) 2017年4期
    關(guān)鍵詞:測(cè)力天平標(biāo)定

    武 龍, 王 鋒, 樂嘉陵

    (中國(guó)空氣動(dòng)力研究與發(fā)展中心 高超聲速?zèng)_壓發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000)

    脈沖燃燒風(fēng)洞測(cè)力系統(tǒng)動(dòng)態(tài)標(biāo)定方法

    武 龍, 王 鋒*, 樂嘉陵

    (中國(guó)空氣動(dòng)力研究與發(fā)展中心 高超聲速?zèng)_壓發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000)

    針對(duì)脈沖燃燒風(fēng)洞中的測(cè)力系統(tǒng),提出了一種動(dòng)態(tài)標(biāo)定方法。利用力錘在模型表面上不同位置,沿不同方向施加一系列集中載荷,由輸入載荷和天平輸出辨識(shí)出該表面對(duì)應(yīng)的單位脈沖響應(yīng)函數(shù)(UIRF),再將各表面對(duì)應(yīng)的UIRF加權(quán)得到系統(tǒng)的UIRF,加權(quán)系數(shù)由試驗(yàn)狀態(tài)下各表面的壓力分布確定。辨識(shí)某表面對(duì)應(yīng)的UIRF時(shí),通過將其參數(shù)化使反卷積問題轉(zhuǎn)化為參數(shù)優(yōu)化問題以回避問題的病態(tài)特性。求解參數(shù)優(yōu)化問題時(shí),先用遺傳算法搜索到參數(shù)全局最優(yōu)解的近似值,再以此作為單純形方法的初值繼續(xù)優(yōu)化得到參數(shù)最優(yōu)值。在ANSYS中模擬了動(dòng)態(tài)標(biāo)定過程,考慮了實(shí)際試驗(yàn)中輸出應(yīng)變含有較大噪聲的情況,驗(yàn)證了這種動(dòng)態(tài)標(biāo)定方法的準(zhǔn)確性。

    脈沖燃燒風(fēng)洞;動(dòng)態(tài)標(biāo)定;參數(shù)辨識(shí);遺傳算法;單純形方法

    0 引 言

    脈沖燃燒風(fēng)洞的建造和運(yùn)行成本較低,能夠以不同尺度、方式進(jìn)行基本流動(dòng)和工程問題的研究,是目前國(guó)內(nèi)開展大尺度超燃沖壓發(fā)動(dòng)機(jī)和機(jī)體推進(jìn)一體化高超聲速飛行器試驗(yàn)研究的主要設(shè)備之一[1]。然而,目前脈沖燃燒風(fēng)洞的工作時(shí)間僅能達(dá)到300ms左右[2],模型振動(dòng)在試驗(yàn)時(shí)間內(nèi)來(lái)不及衰減,測(cè)力天平得到的信號(hào)含有較大的振動(dòng)成分,傳統(tǒng)方法難以獲得良好的結(jié)果。如果能先獲得試驗(yàn)系統(tǒng)的動(dòng)態(tài)輸入輸出特性,就可以用天平的輸出信號(hào)反算出模型所受的載荷歷程而不需要模型達(dá)到穩(wěn)定,這種載荷辨識(shí)方法已經(jīng)在激波風(fēng)洞內(nèi)得到應(yīng)用[3-4]。該方法能夠得到載荷的時(shí)變過程,這對(duì)于飛行器帶動(dòng)力試驗(yàn)[5-6]尤為重要。

    載荷辨識(shí)的前提是已知系統(tǒng)的動(dòng)態(tài)輸入輸出特性,動(dòng)態(tài)標(biāo)定試驗(yàn)是獲得這一特性的可靠途徑。昆士蘭大學(xué)的Abdel-Jawad等對(duì)超高速膨脹管中的多分量應(yīng)力波天平進(jìn)行了動(dòng)態(tài)標(biāo)定,通過加載一系列標(biāo)定載荷,得到了模型對(duì)分布載荷的脈沖響應(yīng)矩陣[7]。牛津大學(xué)的L.J.Doherty等在Stalker管風(fēng)洞中對(duì)一帶動(dòng)力飛行器進(jìn)行測(cè)力試驗(yàn)時(shí),用力錘敲擊方法對(duì)一臺(tái)三分量應(yīng)力波天平進(jìn)行了動(dòng)態(tài)標(biāo)定。中國(guó)空氣動(dòng)力研究與發(fā)展中心的王鋒等采用瞬時(shí)卸載方式對(duì)系統(tǒng)進(jìn)行激勵(lì),分別用階躍函數(shù)和斜坡階躍函數(shù)描述卸載過程,對(duì)某單分量天平進(jìn)行了動(dòng)態(tài)標(biāo)定。本文針對(duì)一個(gè)六分量測(cè)力試驗(yàn)系統(tǒng),給出了一種動(dòng)態(tài)標(biāo)定方法并對(duì)這種方法的準(zhǔn)確性進(jìn)行了驗(yàn)證,為后續(xù)開展載荷辨識(shí)工作提供必要條件[8]。

    1 測(cè)力試驗(yàn)系統(tǒng)介紹

    試驗(yàn)系統(tǒng)包括蒙皮框架結(jié)構(gòu)的試驗(yàn)?zāi)P?、天平以及固定底?部分,如圖1所示。模型包含5個(gè)面,受載荷沖擊時(shí)通過框架將振動(dòng)傳遞至天平浮動(dòng)框,天平固定框與底座相連,底座固定于地面。

    六分量天平由沿軸向布置的前后2個(gè)環(huán)形結(jié)構(gòu)組合而成,如圖2所示,上部為浮動(dòng)框,下部為固定框。浮動(dòng)框與固定框之間是測(cè)量元件,其上貼有應(yīng)變片,可測(cè)量應(yīng)變并通過后續(xù)電路將應(yīng)變信號(hào)輸出。

    2 測(cè)力試驗(yàn)系統(tǒng)模態(tài)分析

    模態(tài)分析是結(jié)構(gòu)動(dòng)力學(xué)分析的基礎(chǔ),系統(tǒng)在載荷作用下發(fā)生振動(dòng),可視為各階模態(tài)振動(dòng)的疊加[9]。通過模態(tài)分析可近似了解系統(tǒng)的基本動(dòng)力學(xué)特征,為開展動(dòng)態(tài)標(biāo)定提供重要的參考信息。利用有限元軟件計(jì)算得到系統(tǒng)的前6階模態(tài),如表1所示。

    表1 試驗(yàn)系統(tǒng)的前6階模態(tài)Table 1 The first six order modes of the testing system

    3 測(cè)力試驗(yàn)系統(tǒng)動(dòng)態(tài)標(biāo)定方法

    3.1 動(dòng)態(tài)標(biāo)定目的

    動(dòng)態(tài)標(biāo)定的目的是獲得單位脈沖載荷作用下系統(tǒng)的響應(yīng)函數(shù)(UIRF)。這個(gè)單位脈沖載荷是指由模型表面等效到天平中心點(diǎn)(或模型質(zhì)心)后的單位脈沖載荷。試驗(yàn)?zāi)P捅旧肀M管整體剛度較大,但也是一個(gè)彈性體,在模型不同位置處加載對(duì)應(yīng)的UIRF不同。這里近似認(rèn)為模型同一個(gè)面上各位置對(duì)應(yīng)的UIRF相同,因此可以在1個(gè)面上施加一系列集中載荷,通過輸入載荷和輸出應(yīng)變求得該面對(duì)應(yīng)的UIRF。在5個(gè)面上分別進(jìn)行上述操作,然后將得到的5個(gè)UIRF加權(quán)即可得到整個(gè)試驗(yàn)系統(tǒng)的UIRF。加權(quán)系數(shù)與試驗(yàn)狀態(tài)下對(duì)應(yīng)面上的平均壓力大小成正比[10]。各面上的平均壓力變化,加權(quán)系數(shù)要相應(yīng)調(diào)整。仿真計(jì)算中,平均壓力由事先的CFD計(jì)算近似得到;風(fēng)洞試驗(yàn)時(shí),在各面上布置測(cè)壓傳感器,每個(gè)時(shí)刻都可以測(cè)出1組壓力值進(jìn)而算得1組加權(quán)系數(shù),實(shí)際上同一車試驗(yàn)中各面上平均壓力的比例關(guān)系隨時(shí)間變化不大,近似計(jì)算也可將不同時(shí)刻的加權(quán)系數(shù)取平均用于后續(xù)計(jì)算。

    定義模型的5個(gè)表面為e1~e5,在模型的面e(e=1,2,3,4,5)上施加一時(shí)變載荷,其相對(duì)于天平中心有6個(gè)分量,即x、y、z方向的力和這3個(gè)方向的力矩,分別記為Fx(t),F(xiàn)y(t),F(xiàn)z(t),Mx(t),My(t),Mz(t)。天平輸出Y有6個(gè)通道,將第k(k=1,2,…,6)個(gè)通道的輸出記為yk(t)。測(cè)力系統(tǒng)的輸入輸出關(guān)系為

    *

    其中,“*”表示卷積。第k行表示第k通道的輸入輸出關(guān)系:在面e上施加載荷,載荷相對(duì)天平中心分別為x、y、z方向單位脈沖力或力矩時(shí),第k通道輸出分別為Ixek、Iyek、Izek、Imxek、Imyek、Imzek。令

    則(1)式可寫為

    Y=Ie*

    Ie即為面e對(duì)應(yīng)的單位脈沖響應(yīng)函數(shù),是一個(gè)6×6的矩陣,整個(gè)系統(tǒng)的單位脈沖響應(yīng)函數(shù)I為,

    其中,αe為面e對(duì)應(yīng)的加權(quán)系數(shù)。對(duì)于某個(gè)時(shí)間點(diǎn),I是一個(gè)6×6的矩陣,動(dòng)態(tài)標(biāo)定的目的即是獲得I。

    3.2 動(dòng)態(tài)標(biāo)定過程

    本文僅針對(duì)模型的1個(gè)表面,以天平的1個(gè)通道為例介紹動(dòng)態(tài)標(biāo)定方法,即只求出Ie的第k行:Iek。其他各面、各通道的動(dòng)態(tài)標(biāo)定方法是相同的。因此只需要考慮(1)式中的第k行

    yk=Iek*

    本文中僅針對(duì)e5面(e=5),第1通道(k=1)求解對(duì)應(yīng)的I51,為了簡(jiǎn)潔,下文中略去下標(biāo)e和k,將(5)式、(8)式分別簡(jiǎn)寫為:

    y=I*

    下文中I均指代I51,對(duì)于某個(gè)時(shí)間點(diǎn)為一1×6向量,下面介紹I的求法。

    動(dòng)態(tài)標(biāo)定采用力錘敲擊法,在e5面上安裝加載部件(圖1中有放大顯示),為力錘敲擊提供5個(gè)不同方向的加載面。加載部件的加載面面積遠(yuǎn)大于力錘錘頭的面積,這樣便于操作并有利于保證敲擊方向的準(zhǔn)確性。定義加載方向?yàn)閐1~d5,如圖3所示。

    在點(diǎn)p(p=1,2,3,…)處沿方向d(d=1,2,3,4,5)敲擊,敲擊力要足夠大以保證天平有明顯輸出但又不能超出量程,力錘可記錄其輸入載荷大小的時(shí)間歷程f(t),結(jié)合加載位置和方向計(jì)算出其相對(duì)于天平中心的載荷向量F

    其中

    cdp為單位載荷相對(duì)于天平中心的載荷矢量,其6個(gè)分量不獨(dú)立,由加載位置和方向決定。給定d、p后,可由簡(jiǎn)單的幾何關(guān)系求得cdp,是已知的。

    設(shè)輸出信號(hào)為ydp(t),由(10)式

    ydp=I*F=I*

    定義

    ydp=Idp*

    這里,Idp即為模型表面加載點(diǎn)對(duì)應(yīng)的脈沖響應(yīng)函數(shù),將在3.3節(jié)中給出其確定算法,此處暫認(rèn)為已知。于是Idp、cdp已知,只有I待求。

    由(9)式,I有6個(gè)未知分量,需要6個(gè)系數(shù)不相關(guān)的方程組成非齊次線性方程組才可解。因此,需要改變d、p進(jìn)行6次不相關(guān)的敲擊。記敲擊次數(shù)為h(h=1,2,3,4,5,6),并令:

    由(14),(16)和(17)式得到

    因此,只需C可逆,即cdp_1,…,cdp_6線性無(wú)關(guān),便可解得I

    3.3Idp的求法

    3.3.1 參數(shù)化Idp

    設(shè)有n個(gè)時(shí)間點(diǎn)t1…tj…tn,那么就存在Idp(t1)…Idp(tj)…Idp(tn)共n個(gè)未知數(shù),若直接將f和ydp(t)代入(15)式解卷積,問題的病態(tài)特征可能使結(jié)果誤差很大。

    考慮結(jié)構(gòu)是線彈性的,可將脈沖響應(yīng)函數(shù)Idp(t)參數(shù)化

    參數(shù)化將Idp寫成m個(gè)不同頻率振動(dòng)疊加的形式。其中,m為選取的模態(tài)數(shù)量,Ai是第i個(gè)模態(tài)的幅值,ωi、γi分別是第i個(gè)模態(tài)的頻率和阻尼比。參數(shù)化后,未知數(shù)變?yōu)?m個(gè),通常,n?m,未知數(shù)的個(gè)數(shù)大大減少。

    設(shè)實(shí)際測(cè)量得到的輸出信號(hào)為yodp(t),令R(t)=yodp(t)-ydp(t),定義殘差r

    求Idp(t)的問題轉(zhuǎn)化為用Idp*f擬合yodp(t)的問題,即確定能使r取到最小值的m及其對(duì)應(yīng)的Ai、ωi、γi這一優(yōu)化問題。

    3.3.2 優(yōu)化算法

    分析3.3.1節(jié)中的這一優(yōu)化問題,待定參數(shù)是m階模態(tài)參數(shù)(ω1,A1,γ1)…(ωm,Am,γm),優(yōu)化目標(biāo)是使殘差r取最小值。優(yōu)化過程中,r逐漸減小,給定一個(gè)值rl,以r

    阻尼比γ一般在10-3量級(jí)或更小,在振動(dòng)頻率不太高、測(cè)量時(shí)長(zhǎng)較短的情況下對(duì)優(yōu)化結(jié)果影響很小,故先不考慮阻尼比,對(duì)每個(gè)模態(tài)僅考慮幅值和頻率2個(gè)參數(shù),Idp(t)簡(jiǎn)化為

    求解前m未知,可預(yù)估m(xù)值,然后對(duì)3m個(gè)參數(shù)同時(shí)優(yōu)化。但這對(duì)預(yù)估的要求很高,m取小了導(dǎo)致Idp(t)中必然缺少高階振動(dòng)成分;m取大了,則優(yōu)化參數(shù)過多,優(yōu)化難度較大。

    這里設(shè)計(jì)了一個(gè)算法,將整個(gè)優(yōu)化過程分若干次進(jìn)行,每次應(yīng)用一定的優(yōu)化算法只辨識(shí)出1個(gè)模態(tài)對(duì)應(yīng)的成分。第1次優(yōu)化時(shí),待擬合曲線為yodp(t)。第i次(i=1,2…)優(yōu)化終止時(shí),計(jì)算Ri(t)并輸出辨識(shí)出的第i階頻率ωi和幅值A(chǔ)i。第i+1次優(yōu)化時(shí),以Ri(t)作為待擬合曲線,重復(fù)上一步的優(yōu)化過程,得到ωi+1和Ai+1。因?yàn)槊看蝺?yōu)化的目標(biāo)都是使殘差盡可能減小,被識(shí)別出來(lái)的都是當(dāng)次待擬合曲線中最主要的振動(dòng)成分,所以各振動(dòng)成分是按照幅值由大到小的順序被辨識(shí)出來(lái)的,即Ai+1≥Ai。當(dāng)進(jìn)行到第M次優(yōu)化后,AM?A1,主要的模態(tài)成分已經(jīng)被辨識(shí)出來(lái),若達(dá)到了r

    在第i次優(yōu)化中,如圖4虛線框中部分,未知數(shù)為ωi和Ai,待擬合曲線為Ri-1(t)(i=1時(shí)為yodp(t)),目標(biāo)函數(shù)為ri(ωi,Ai)??紤]r(ωi,Ai)是一個(gè)多峰函數(shù),使用傳統(tǒng)的搜索方法極易陷入局部最優(yōu)解。遺傳算法是一種概率搜索算法,可以從一個(gè)種群開始并行搜索,并且不依賴于目標(biāo)函數(shù)的梯度信息,具有很強(qiáng)的全局優(yōu)化性能[11-12],適合提取當(dāng)前幅值最大的振動(dòng)成分。這里采用Matlab的遺傳算法工具箱(GA)進(jìn)行優(yōu)化[13-14]。收斂準(zhǔn)則為本代殘差r與上一代殘差rlg的相對(duì)差值小于ε,即

    此處ε取10-4。求解前給定初代種群的取值范圍,求解中如果某一代滿足了收斂準(zhǔn)則,就輸出結(jié)果;如果不滿足收斂準(zhǔn)則,就通過選擇性復(fù)制、交叉、變異

    生成下一代種群,直到滿足收斂準(zhǔn)則為止。

    在上述過程中,若rl取得較小則會(huì)出現(xiàn)ωi≈ωj,Ai?Aj,(i

    前面已經(jīng)得到了(ω1,A1)…(ωm,Am)和Idp(t)?,F(xiàn)加入對(duì)阻尼比γ的考慮,若繼續(xù)基于上述的遺傳算法進(jìn)行優(yōu)化,減小ε和rl,是可以得到預(yù)期的Idp(t)的,但是求解效率很低。進(jìn)一步優(yōu)化時(shí)采用Matlab優(yōu)化工具箱的fminsearch函數(shù),它是基于單純形算法的優(yōu)化工具。將前面得到的(ω1,A1)…(ωm,Am)作為頻率和幅值的初值,0作為阻尼比γ的初值。由于這些初值已經(jīng)很接近最優(yōu)解,fminsearch函數(shù)只需在最優(yōu)解附近搜索即可,這就避免了其易于陷入局部最優(yōu)解的缺點(diǎn),而充分發(fā)揮了其收斂速度更快的優(yōu)勢(shì)[15]。經(jīng)過若干次迭代,殘差迅速下降并達(dá)到穩(wěn)定,此時(shí)終止優(yōu)化,輸出(ω1,A1,γ1)…(ωm,Am,γm),并代入(20)式計(jì)算出Idp(t)。

    4 測(cè)力試驗(yàn)系統(tǒng)動(dòng)態(tài)標(biāo)定方法驗(yàn)證

    第2節(jié)中介紹了動(dòng)態(tài)標(biāo)定試驗(yàn)的操作方法和求系統(tǒng)單位脈沖響應(yīng)函數(shù)的算法。本節(jié)將利用有限元軟件模擬試驗(yàn)過程,利用上述算法求得I,并檢驗(yàn)I的準(zhǔn)確性。

    4.1 模擬動(dòng)態(tài)標(biāo)定試驗(yàn)

    在Ansys Workbench中進(jìn)行瞬態(tài)動(dòng)力學(xué)分析,模擬上述的動(dòng)態(tài)標(biāo)定過程。分別在點(diǎn)p1(d1,d2,d4方向)、點(diǎn)p2(d1,d2方向)、點(diǎn)p3(d1方向)模擬力錘加載,各加載點(diǎn)位置如圖3所示。由幾何關(guān)系分別算得這6次加載的載荷矢量:c11、c21、c41、c12、c22、c13,代入(17)式得到矩陣C

    C的秩r(C)=6,載荷矢量滿足線性無(wú)關(guān)的要求。

    載荷歷程為f1(t)

    測(cè)得對(duì)應(yīng)的輸出,記為yo11(t)、yo21(t)、yo41(t)、yo12(t)、yo22(t)、yo13(t)。

    4.2 求I(t)

    4.2.1 求6次敲擊分別對(duì)應(yīng)的Idp(t)

    利用3.3中的算法求這6次敲擊對(duì)應(yīng)的Idp(t),即I11(t)、I21(t)、I41(t)、I12(t)、I22(t)、I13(t)。

    以I11(t)為例,利用遺傳算法得到前10個(gè)振動(dòng)成分,如表2所示。

    表2 前10個(gè)振動(dòng)成分Table 2 The first ten frequencies and amplitudes

    可見,幅值遞減,前5個(gè)幅值明顯大于后續(xù)幅值,符合預(yù)期。前10個(gè)頻率都在5.8Hz、9.4Hz、12.3Hz、16.9Hz、23.0Hz附近,這也與第2節(jié)中模態(tài)分析所得到的系統(tǒng)前幾階固有頻率相符。進(jìn)一步利用fminsearch函數(shù)優(yōu)化時(shí)只需考慮這5個(gè)振動(dòng)成分,按照(24)、(25)式算得優(yōu)化初值,如表3所示。

    表3 優(yōu)化初值Table 3 Initial values of optimization

    由于初值已經(jīng)比較準(zhǔn)確,殘差r在前1000步迭代中便快速下降,經(jīng)過約4000步穩(wěn)定在10-9量級(jí),得到模態(tài)參數(shù),如表4所示。

    表4 優(yōu)化結(jié)果Table 4 Results of optimization

    將表4中的模態(tài)參數(shù)代入(20)式,即可得到I11(t)。同理,可求得I21(t)、I41(t)、I12(t)、I22(t)、I13(t)。再按照(16)中的定義即可得到

    4.2.2 求I(t)

    將(26)、(28)式代入(19)式,即可求得I(t)

    I(t)=S·

    由上文,此處I(t)省略了下標(biāo)e和k,其中

    (29)式即為模型e5面上天平1通道對(duì)應(yīng)的UIRF,其他各面、各通道對(duì)應(yīng)的UIRF求法相同。

    得到其他4個(gè)面上天平1通道對(duì)應(yīng)的UIRF,均與(29)式相近,其中表征頻率和阻尼比的S相同。表征幅值的系數(shù)矩陣略有不同,但其中較大的系數(shù)相對(duì)變化很小,說(shuō)明在不同面上敲擊激發(fā)的起主導(dǎo)作用的模態(tài)是相同的。由CFD計(jì)算算得5個(gè)面上的平均壓力進(jìn)而算出各面對(duì)應(yīng)的加權(quán)系數(shù)分別為:0.3865,0.1725,0.0062,0.1150,0.3197。由(7)式可得整個(gè)模型對(duì)天平1通道的UIRF,

    I(t)=S·

    測(cè)力系統(tǒng)在天平1通道的動(dòng)態(tài)標(biāo)定完成。

    4.3 驗(yàn)證I(t)的準(zhǔn)確性

    在p4點(diǎn)施加d5方向的載荷,載荷歷程為f2(t)

    測(cè)得輸出應(yīng)變yo54(t),并利用4.2中得到的I(t)計(jì)算出理論輸出y54(t)。將yo54(t)、y54(t)繪制在圖5中。兩條曲線幾乎完全重合,說(shuō)明動(dòng)態(tài)標(biāo)定得到的I(t)是準(zhǔn)確的。

    Fig.5 Comparison between computed strainy54(t)andmeasuredstrainyo54(t)

    4.4 驗(yàn)證標(biāo)定方法在測(cè)得應(yīng)變含噪聲時(shí)的準(zhǔn)確性

    實(shí)際試驗(yàn)中,輸出應(yīng)變中含有噪聲,這里驗(yàn)證上述動(dòng)態(tài)標(biāo)定方法在ydp(t)含較大噪聲時(shí)的準(zhǔn)確性。

    以p1點(diǎn)處沿d1方向的敲擊為例,在yo11(t)中加入以該信號(hào)標(biāo)準(zhǔn)差15%為幅值的隨機(jī)信號(hào),記為yn11(t),模擬含有較大噪聲的輸出信號(hào),局部圖如圖6所示。

    利用3.3中的算法由yn11(t)計(jì)算出In11(t),對(duì)比In11(t)和I11(t)。先利用遺傳算法得到In11(t)的前10階振動(dòng)成分,發(fā)現(xiàn)與表2中的結(jié)果十分接近。然后利用fminsearch函數(shù)繼續(xù)優(yōu)化,殘差r穩(wěn)定在約1.77,得到In11(t)的各振動(dòng)成分,如表5所示。

    表5 優(yōu)化結(jié)果Table 5 Results of optimization

    表5與表4結(jié)果非常接近,只有低頻成分的阻尼比稍有差異。這是因?yàn)楸竟?jié)為了驗(yàn)證算法,加入的噪聲很大,低頻阻尼比的影響遠(yuǎn)遠(yuǎn)小于噪聲的影響。實(shí)際試驗(yàn)中噪聲并不會(huì)這么大,低頻阻尼比的誤差會(huì)更小。繪制I11(t)、In11(t),局部圖如圖7所示。兩條曲線幾乎重合,說(shuō)明即使yn11(t)中含有明顯的噪聲,利用3.3中的算法,得到的In11(t)仍然準(zhǔn)確。

    可見,yn11(t)中的噪聲基本沒有被擬合,而是保留在了殘差r中。

    5 結(jié) 論

    本文給出了一種針對(duì)脈沖燃燒風(fēng)洞測(cè)力系統(tǒng)的動(dòng)態(tài)標(biāo)定方法并利用ANSYS仿真對(duì)其進(jìn)行了驗(yàn)證,得到了以下結(jié)論:

    (1) 求解UIRF時(shí)通過參數(shù)化將解卷積問題轉(zhuǎn)化為參數(shù)優(yōu)化問題,可以有效回避問題的病態(tài)特性。

    (2) 求解參數(shù)優(yōu)化問題時(shí)先利用遺傳算法搜索到全局最優(yōu)解的近似值,再以其作為單純形方法的初值繼續(xù)優(yōu)化,可以迅速求得全局最優(yōu)解。

    (3) 即使實(shí)際試驗(yàn)測(cè)得的應(yīng)變信號(hào)含有較大的噪聲,這種動(dòng)態(tài)標(biāo)定方法仍具有很高的精度。

    (4) 本動(dòng)態(tài)標(biāo)定方法在實(shí)際應(yīng)用中可能存在的誤差主要由力錘敲擊方向的偏差引起,試驗(yàn)人員在操作前需要進(jìn)行適當(dāng)?shù)木毩?xí)以盡量減小敲擊方向的偏差。

    [1]樂嘉陵, 劉偉雄, 賀偉, 等. 脈沖燃燒風(fēng)洞及其在火箭和超燃發(fā)動(dòng)機(jī)研究中的應(yīng)用[J]. 實(shí)驗(yàn)流體力學(xué), 2005, 19(1): 1-10.

    Le J L, Liu W X, He W, et al. Impulse combustion wind tunnel and its application in rocket and scramjet research[J]. Journal of Experiments in Fluid Mechanics, 2005, 19(1): 1-10.

    [2]劉偉雄, 譚宇, 毛雄兵, 等. 一種新運(yùn)行方式脈沖燃燒風(fēng)洞研制及初步應(yīng)用[J]. 實(shí)驗(yàn)流體力學(xué), 2007, 21(4): 59-64.

    Liu W X, Tan Y, Mao X B, et al. The development and preliminary application of a pulse combustion wind tunnel with new running way[J]. Journal of Experiments in Fluid Mechanics, 2007, 21(4): 59-64.

    [3]Robinson M J, Mee D J, Tsai C Y, et al. Three-component force measurements on a large scramjet in a shock tunnel[J]. Journal of Spacecraft and Rockets, 2004, 41(3): 416-425.

    [4]Robinson M J, Hannemann K, Schramm J M. Design and implementation of an internal stress wave force balance in a shock tunnel[J]. CEAS Space Journal, 2011, 1(1): 45-57.

    [5]賀偉, 于時(shí)恩, 李宏斌. 高超聲速一體化飛行器推阻特性測(cè)量研究[J]. 實(shí)驗(yàn)流體力學(xué), 2010, 24(2): 65-68.

    He W, Yu S E, Li H B. Experimental investigation on thrust drag performance of hypersonic integrative vehicle[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(2): 65-68.

    [6]賀偉, 童澤潤(rùn), 李宏斌. 單模塊超燃發(fā)動(dòng)機(jī)推力測(cè)量天平研制[J]. 航空動(dòng)力學(xué)報(bào), 2010, 25(10): 2285-2289.

    He W, Tong Z R, Li H B. Investigation of thrust balance for the single module scramjet[J]. Journal of Aerospace Power, 2010, 25(10): 2285-2289.

    [7]Abdel-Jawad M M, Mee D J, Morgan R G. New calibration technique for multiple-component stress wave force balances[J]. Review of Scientific Instruments, 2007, 78(6): 065101-1-065101-7.

    [8]王鋒, 任虎, 周正, 等. 載荷辨識(shí)方法用于脈沖風(fēng)洞模型阻力測(cè)量研究[J]. 振動(dòng)與沖擊, 2015, 34(24): 202-208.

    Wang F, Ren H, Zhou Z, et al. Drag force measurement in impulse facilities by using load identification method[J]. Journal of Vibration and Shock, 2015, 34(24): 202-208.

    [9]李東旭. 高等結(jié)構(gòu)動(dòng)力學(xué)[M]. 北京: 科學(xué)出版社, 2010: 273-298.

    Li D X. Advanced structural dynamics[M]. Beijing: The Science Publishing Company, 2010: 273-298.

    [10]Doherty L J, Smart M K, Mee D J. Measurement of three-components of force on an airframe integrated scramjet at Mach 10[R]. AIAA-2015-3523, 2015.

    [11]劉國(guó)春, 費(fèi)強(qiáng), 趙武云, 等. 基于Matlab 遺傳算法優(yōu)化工具箱的應(yīng)用[J]. 機(jī)械研究與應(yīng)用, 2014, 27(2): 71-73.

    Liu G C, Fei Q, Zhao W Y, et al. Application of genetic algorithm optimization toolbox based on matlab[J]. Mechanical Research and Application, 2014, 27(2): 71-73.

    [12]林鴻彬. 基于遺傳算法的數(shù)據(jù)擬合在 MATLAB 環(huán)境中的實(shí)現(xiàn)[J]. 湖南農(nóng)機(jī), 2010, 37(3): 92-97.

    Lin H B. Data fitting based on genetic algorithm implementation in MATLAB environment[J]. Hunan Agricultural Machinery, 2010, 37(3): 92-97.

    [13]羅述全. 傳統(tǒng)優(yōu)化算法與遺傳算法的比較[J]. 湖北工業(yè)大學(xué)學(xué)報(bào), 2007, 22(3): 32-35.

    Luo S Q. Comparison between traditional optimized algorithm and heredity algorithm[J]. Hubei University of Technology Journal, 2007, 22(3): 32-35.

    [14]郭海雙, 梁佳雯, 張劭昀. MATLAB 遺傳算法工具箱GADS優(yōu)化及應(yīng)用[J]. 電子設(shè)計(jì)工程, 2015, 23(10): 27-30.

    Guo H S, Liang J W, Zhang S Y. Optimization and examples in Matlab GA toolbox GADS[J]. Electronic Design Engineering, 2015, 23(10): 27-30.

    [15]楊改強(qiáng), 霍麗娟, 楊國(guó)義, 等. 利用MATLAB擬合van Genuchten方程參數(shù)的研究[J]. 土壤, 2010, 42(2): 268-274.

    Yang G Q, Huo L J, Yang G Y, et al. Research on fitting van genuchten equation parameter with MATLAB software[J]. Soils, 2010, 42(2): 268-274.

    (編輯:張巧蕓)

    A dynamic calibration method for a dynamometric systemin impulse combustion facilities

    Wu Long, Wang Feng*, Le Jialing

    (Science and Technology on Scramjet Laboratory, China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China)

    A new dynamic calibration method for a dynamometric system in impulse combustion facilities is proposed. The calibration involved uses an instrumented impact hammer to apply individual calibration forces in different directions at different positions on a face of the model and calculates the unit impulse response function (UIRF) of the face from input loads and output strains. UIRFs of different faces are weighted to obtain the UIRF of the dynamometric system and the weighting coefficients are determined by the pressure on each face under the experimental condition. By parameterization, the problem is converted into a parameter optimization problem to solve the UIRF. Using a genetic algorithm to obtain the approximation of the global optimal solution of the parameters and setting it as the initial value of a simplex algorithm, the exact solution is obtained by the simplex algorithm then. ANSYS simulation of the dynamic calibration is presented. Input loads and output strains are recorded and noises are added to the output strains to simulate the actual experimental situation. The simulation validates the accuracy and feasibility of the dynamic calibration method.

    impulse facilities;dynamic calibration;parameter identification;genetic algorithm;simplex algorithm

    1672-9897(2017)04-0051-08

    10.11729/syltlx20160158

    2016-10-21;

    2017-01-09

    國(guó)家自然科學(xué)基金項(xiàng)目(11372339);高超聲速?zèng)_壓發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室基金項(xiàng)目(STSKFKT2012001)

    WuL,WangF,LeJL.Adynamiccalibrationmethodforadynamometricsysteminimpulsecombustionfacilities.JournalofExperimentsinFluidMechanics, 2017, 31(4): 51-58. 武 龍, 王 鋒, 樂嘉陵. 脈沖燃燒風(fēng)洞測(cè)力系統(tǒng)動(dòng)態(tài)標(biāo)定方法. 實(shí)驗(yàn)流體力學(xué), 2017, 31(4): 51-58.

    V211.7

    A

    武 龍(1991-),男,黑龍江省大慶市人,碩士研究生。研究方向:載荷辨識(shí)。通信地址:四川省綿陽(yáng)市涪城區(qū)劍門路西段278號(hào)(621000)。E-mail: 779483196@qq.com。

    *通信作者 E-mail: wfscholar@163.com

    猜你喜歡
    測(cè)力天平標(biāo)定
    說(shuō)說(shuō)天平的使用
    主向力作用下壓電測(cè)力儀內(nèi)部側(cè)向力計(jì)算方法
    天平使用前后的兩次平衡
    使用朗仁H6 Pro標(biāo)定北汽紳寶轉(zhuǎn)向角傳感器
    測(cè)力延度在膠粉改性瀝青低溫性能評(píng)價(jià)中的應(yīng)用
    石油瀝青(2019年1期)2019-03-05 08:25:46
    天平的平衡
    基于勻速率26位置法的iIMU-FSAS光纖陀螺儀標(biāo)定
    船載高精度星敏感器安裝角的標(biāo)定
    基于Harris-張正友平面標(biāo)定法的攝像機(jī)標(biāo)定算法
    剛?cè)峄旌先攘S力傳感器測(cè)力性能分析
    精品国产亚洲在线| 欧美激情高清一区二区三区| 欧美日韩一级在线毛片| 国产日韩一区二区三区精品不卡| 久久久久精品国产欧美久久久| 精品一区二区三区av网在线观看| 色老头精品视频在线观看| 男女午夜视频在线观看| 久久久久久久午夜电影| 久久久精品国产亚洲av高清涩受| 51午夜福利影视在线观看| 在线观看日韩欧美| 久久国产亚洲av麻豆专区| 亚洲三区欧美一区| 久久天堂一区二区三区四区| 国产成人影院久久av| 久久草成人影院| 窝窝影院91人妻| 在线观看www视频免费| 黄色女人牲交| 国产成人影院久久av| 日本三级黄在线观看| 国产成+人综合+亚洲专区| 男人舔女人的私密视频| 90打野战视频偷拍视频| 最新在线观看一区二区三区| 亚洲自偷自拍图片 自拍| 国产精品久久久久久亚洲av鲁大| 色综合站精品国产| 日本一区二区免费在线视频| 热re99久久国产66热| 国产亚洲精品久久久久久毛片| 久久久国产欧美日韩av| 国产av精品麻豆| 久久久精品国产亚洲av高清涩受| 亚洲午夜理论影院| 久久国产乱子伦精品免费另类| 日韩三级视频一区二区三区| 日韩有码中文字幕| 精品国产国语对白av| ponron亚洲| 亚洲熟妇中文字幕五十中出| 可以在线观看的亚洲视频| 国产成人av激情在线播放| 国产精品亚洲一级av第二区| 亚洲久久久国产精品| a在线观看视频网站| 色尼玛亚洲综合影院| 少妇熟女aⅴ在线视频| 亚洲美女黄片视频| 大型av网站在线播放| 欧美黄色淫秽网站| 一本综合久久免费| 99re在线观看精品视频| 18禁裸乳无遮挡免费网站照片 | 午夜成年电影在线免费观看| 女同久久另类99精品国产91| 亚洲熟妇中文字幕五十中出| 一进一出抽搐动态| 欧美日本视频| 国产成人系列免费观看| 中亚洲国语对白在线视频| 日本三级黄在线观看| 九色亚洲精品在线播放| 黄片小视频在线播放| 免费搜索国产男女视频| 久久人妻av系列| 亚洲成人免费电影在线观看| 久久香蕉国产精品| 午夜精品久久久久久毛片777| av天堂久久9| 亚洲精品一卡2卡三卡4卡5卡| 搡老岳熟女国产| 国产成人精品久久二区二区91| 老司机午夜十八禁免费视频| 色婷婷久久久亚洲欧美| 国产精品久久久久久精品电影 | 老司机福利观看| 欧美性长视频在线观看| 午夜福利一区二区在线看| a级毛片在线看网站| 欧美另类亚洲清纯唯美| 国产视频一区二区在线看| 91麻豆精品激情在线观看国产| 男女床上黄色一级片免费看| av福利片在线| 看黄色毛片网站| 日韩精品青青久久久久久| 90打野战视频偷拍视频| 精品日产1卡2卡| 又黄又爽又免费观看的视频| 免费av毛片视频| 这个男人来自地球电影免费观看| 国产人伦9x9x在线观看| 亚洲精品国产区一区二| 亚洲av成人av| 久99久视频精品免费| 国产1区2区3区精品| ponron亚洲| 淫秽高清视频在线观看| 中出人妻视频一区二区| 色婷婷久久久亚洲欧美| 久久中文看片网| 精品欧美国产一区二区三| 啦啦啦观看免费观看视频高清 | 很黄的视频免费| 搡老熟女国产l中国老女人| 久久国产乱子伦精品免费另类| 俄罗斯特黄特色一大片| 亚洲第一电影网av| 丝袜在线中文字幕| 美女高潮到喷水免费观看| 国产精品秋霞免费鲁丝片| 欧美日本视频| 他把我摸到了高潮在线观看| 久久久国产成人免费| 日日爽夜夜爽网站| 亚洲第一av免费看| 在线观看免费日韩欧美大片| 窝窝影院91人妻| 国产av在哪里看| 欧美日韩乱码在线| 中文字幕最新亚洲高清| 亚洲 国产 在线| 国产一区二区激情短视频| 亚洲精品久久国产高清桃花| 女人被狂操c到高潮| 精品久久久久久久久久免费视频| www日本在线高清视频| 久久久久久人人人人人| 久久精品aⅴ一区二区三区四区| ponron亚洲| 亚洲久久久国产精品| 十八禁人妻一区二区| √禁漫天堂资源中文www| 色综合婷婷激情| 成人av一区二区三区在线看| 色综合亚洲欧美另类图片| 成在线人永久免费视频| 中亚洲国语对白在线视频| 午夜精品在线福利| 欧美中文日本在线观看视频| 高潮久久久久久久久久久不卡| av电影中文网址| 麻豆久久精品国产亚洲av| 18禁国产床啪视频网站| 99国产综合亚洲精品| 久久久久国内视频| 亚洲avbb在线观看| 欧美日本视频| 一区二区三区高清视频在线| 精品乱码久久久久久99久播| av网站免费在线观看视频| 自拍欧美九色日韩亚洲蝌蚪91| 丝袜美足系列| 精品国产一区二区三区四区第35| 老鸭窝网址在线观看| 久久午夜综合久久蜜桃| 黑人巨大精品欧美一区二区蜜桃| av视频在线观看入口| 校园春色视频在线观看| 老熟妇仑乱视频hdxx| 久久中文字幕一级| 久久热在线av| 午夜福利18| 精品免费久久久久久久清纯| 亚洲 欧美一区二区三区| 乱人伦中国视频| 欧美成人性av电影在线观看| 精品国内亚洲2022精品成人| 午夜福利欧美成人| 国产亚洲av嫩草精品影院| 看片在线看免费视频| 精品久久久久久,| 一级,二级,三级黄色视频| 午夜免费鲁丝| 久久国产精品男人的天堂亚洲| 日本 欧美在线| 美女被艹到高潮喷水动态| 日本与韩国留学比较| 国产精品久久久久久av不卡| 成人欧美大片| 国产主播在线观看一区二区| 国产欧美日韩一区二区精品| 亚洲经典国产精华液单| 久久久久久国产a免费观看| 久久草成人影院| av.在线天堂| 国产精品爽爽va在线观看网站| 99热这里只有是精品在线观看| 国产亚洲精品综合一区在线观看| 国产午夜精品久久久久久一区二区三区 | 婷婷精品国产亚洲av在线| 少妇猛男粗大的猛烈进出视频 | 夜夜夜夜夜久久久久| 国产精品1区2区在线观看.| 99在线人妻在线中文字幕| 久久久国产成人精品二区| 女同久久另类99精品国产91| 国产黄色小视频在线观看| 国产大屁股一区二区在线视频| 亚洲中文字幕一区二区三区有码在线看| 国产av不卡久久| 性欧美人与动物交配| 国产高清不卡午夜福利| 丰满乱子伦码专区| 亚洲天堂国产精品一区在线| 亚洲精品一区av在线观看| 日本 欧美在线| 真人做人爱边吃奶动态| 久久久久久久久大av| 一级黄片播放器| 网址你懂的国产日韩在线| 欧美一区二区国产精品久久精品| 亚洲国产精品久久男人天堂| 亚洲自拍偷在线| 日韩欧美国产在线观看| 亚洲性久久影院| 91麻豆精品激情在线观看国产| 亚洲国产色片| av福利片在线观看| 狂野欧美激情性xxxx在线观看| 啦啦啦啦在线视频资源| 国产精品一区二区性色av| 久久久精品欧美日韩精品| 亚洲美女搞黄在线观看 | 欧美黑人欧美精品刺激| 中文字幕免费在线视频6| 国产女主播在线喷水免费视频网站 | 88av欧美| 中出人妻视频一区二区| 久久久午夜欧美精品| videossex国产| 欧美潮喷喷水| 美女黄网站色视频| 国产精品久久久久久久久免| 成人国产一区最新在线观看| 男人的好看免费观看在线视频| 中文字幕av在线有码专区| 色综合亚洲欧美另类图片| 亚洲成a人片在线一区二区| 在线a可以看的网站| 五月玫瑰六月丁香| 国内精品久久久久久久电影| 99久久精品一区二区三区| 欧美激情在线99| 国产精品一区二区免费欧美| 99久久久亚洲精品蜜臀av| 99热网站在线观看| 美女cb高潮喷水在线观看| 在线播放国产精品三级| 99久久久亚洲精品蜜臀av| 69人妻影院| 淫妇啪啪啪对白视频| 午夜影院日韩av| 午夜福利视频1000在线观看| 在线播放无遮挡| 国产精品一及| 国产中年淑女户外野战色| 国产乱人伦免费视频| 国产黄a三级三级三级人| 欧美中文日本在线观看视频| 亚洲中文字幕一区二区三区有码在线看| 长腿黑丝高跟| 欧美丝袜亚洲另类 | 88av欧美| 国产精品久久久久久久久免| 在线观看免费视频日本深夜| 男人舔女人下体高潮全视频| 成人午夜高清在线视频| 91麻豆av在线| 亚洲欧美日韩东京热| 日韩精品有码人妻一区| 最后的刺客免费高清国语| 免费观看精品视频网站| 国产精品久久久久久亚洲av鲁大| 欧美性猛交黑人性爽| 国产视频内射| 亚洲av五月六月丁香网| 国产免费av片在线观看野外av| 男女视频在线观看网站免费| 中文在线观看免费www的网站| 日本黄色视频三级网站网址| 女生性感内裤真人,穿戴方法视频| 99精品久久久久人妻精品| 黄色一级大片看看| 国内揄拍国产精品人妻在线| 成人无遮挡网站| 男人狂女人下面高潮的视频| 日韩精品中文字幕看吧| 国产免费av片在线观看野外av| 久久午夜福利片| 日日夜夜操网爽| 此物有八面人人有两片| 我要搜黄色片| 精品人妻一区二区三区麻豆 | 少妇的逼水好多| 日本色播在线视频| 亚洲va在线va天堂va国产| 天堂av国产一区二区熟女人妻| 久久久午夜欧美精品| 丝袜美腿在线中文| 国内精品一区二区在线观看| 亚洲性夜色夜夜综合| 精品久久久噜噜| 久久精品国产亚洲av香蕉五月| 夜夜夜夜夜久久久久| 91久久精品电影网| 真人一进一出gif抽搐免费| 国语自产精品视频在线第100页| 变态另类丝袜制服| 大又大粗又爽又黄少妇毛片口| 一级a爱片免费观看的视频| 亚洲人成伊人成综合网2020| 此物有八面人人有两片| 韩国av在线不卡| 亚洲专区中文字幕在线| 国产成人aa在线观看| 一区二区三区激情视频| 嫁个100分男人电影在线观看| 欧美精品啪啪一区二区三区| 美女xxoo啪啪120秒动态图| 可以在线观看毛片的网站| 久久久久久久久久黄片| 老司机午夜福利在线观看视频| 国产高清视频在线播放一区| 午夜精品在线福利| 亚洲成人久久性| 亚洲va在线va天堂va国产| 婷婷精品国产亚洲av| 69av精品久久久久久| 国产精品99久久久久久久久| 欧美性猛交╳xxx乱大交人| 亚洲中文日韩欧美视频| 色综合站精品国产| 欧美高清成人免费视频www| АⅤ资源中文在线天堂| 日韩一区二区视频免费看| 久久久久久久久久黄片| 国内精品久久久久精免费| 久久精品国产清高在天天线| 最近最新中文字幕大全电影3| 五月伊人婷婷丁香| 精品久久久久久成人av| 国产一区二区在线观看日韩| 亚洲久久久久久中文字幕| 黄片wwwwww| 国产v大片淫在线免费观看| 国产av在哪里看| 亚洲人成伊人成综合网2020| av在线观看视频网站免费| 亚洲人成伊人成综合网2020| 日本与韩国留学比较| 欧美色欧美亚洲另类二区| 日韩一本色道免费dvd| 人妻少妇偷人精品九色| 亚洲av中文av极速乱 | 亚洲自拍偷在线| 夜夜夜夜夜久久久久| 欧美又色又爽又黄视频| 欧美丝袜亚洲另类 | 欧美日韩亚洲国产一区二区在线观看| 一个人免费在线观看电影| 在线观看舔阴道视频| 久久精品综合一区二区三区| 中文字幕高清在线视频| a级毛片a级免费在线| 国产伦精品一区二区三区四那| 久久中文看片网| 国产在视频线在精品| 99热网站在线观看| 亚洲在线自拍视频| 久久精品影院6| 国产精品电影一区二区三区| 久久精品综合一区二区三区| 尾随美女入室| 18禁在线播放成人免费| 国产精品亚洲美女久久久| 极品教师在线免费播放| 亚洲最大成人中文| 亚洲黑人精品在线| av福利片在线观看| 午夜福利视频1000在线观看| 12—13女人毛片做爰片一| 成人性生交大片免费视频hd| 亚洲精品在线观看二区| 啦啦啦啦在线视频资源| 男人舔女人下体高潮全视频| 精品人妻一区二区三区麻豆 | 久久精品国产亚洲av香蕉五月| 国产成人av教育| 久久久久久国产a免费观看| 日韩欧美一区二区三区在线观看| 99九九线精品视频在线观看视频| 婷婷丁香在线五月| 日本黄大片高清| 听说在线观看完整版免费高清| 性欧美人与动物交配| av在线蜜桃| 欧美成人a在线观看| av在线亚洲专区| 哪里可以看免费的av片| 日本三级黄在线观看| 日本与韩国留学比较| 亚洲国产色片| 又紧又爽又黄一区二区| 舔av片在线| 啦啦啦观看免费观看视频高清| av国产免费在线观看| 国产麻豆成人av免费视频| 国产高清不卡午夜福利| 国产黄色小视频在线观看| 国产激情偷乱视频一区二区| 久久久成人免费电影| 性欧美人与动物交配| 美女免费视频网站| 尤物成人国产欧美一区二区三区| 久久久久久伊人网av| 亚洲在线自拍视频| 亚洲无线在线观看| 日韩,欧美,国产一区二区三区 | 欧美日韩国产亚洲二区| 丝袜美腿在线中文| 成人特级黄色片久久久久久久| 亚洲va在线va天堂va国产| 永久网站在线| 国内精品久久久久久久电影| eeuss影院久久| 高清在线国产一区| 长腿黑丝高跟| 亚洲无线观看免费| 亚洲aⅴ乱码一区二区在线播放| 色综合色国产| 久久九九热精品免费| 精品一区二区三区人妻视频| 无人区码免费观看不卡| 91av网一区二区| 欧美另类亚洲清纯唯美| 亚洲精品色激情综合| 啦啦啦韩国在线观看视频| 噜噜噜噜噜久久久久久91| 伊人久久精品亚洲午夜| 欧美成人一区二区免费高清观看| 成人特级黄色片久久久久久久| 赤兔流量卡办理| 热99在线观看视频| 国产色爽女视频免费观看| 国产一区二区在线av高清观看| 成人一区二区视频在线观看| 舔av片在线| 国产毛片a区久久久久| a在线观看视频网站| 12—13女人毛片做爰片一| xxxwww97欧美| 男女做爰动态图高潮gif福利片| 日本三级黄在线观看| 亚洲自偷自拍三级| 人人妻人人看人人澡| 亚洲七黄色美女视频| 在线观看免费视频日本深夜| 熟女人妻精品中文字幕| av国产免费在线观看| 乱码一卡2卡4卡精品| АⅤ资源中文在线天堂| 国产精品99久久久久久久久| 午夜老司机福利剧场| 国产精品一区二区性色av| 日本撒尿小便嘘嘘汇集6| 91av网一区二区| 日韩欧美国产在线观看| 91久久精品国产一区二区三区| 欧美日韩国产亚洲二区| 一边摸一边抽搐一进一小说| 一区二区三区高清视频在线| 国产精品久久久久久久久免| 99热只有精品国产| 亚洲成av人片在线播放无| 观看美女的网站| 亚洲经典国产精华液单| 国产精品一区www在线观看 | 亚洲图色成人| 国产欧美日韩精品亚洲av| 看黄色毛片网站| 欧美在线一区亚洲| 国产极品精品免费视频能看的| 亚洲国产色片| av国产免费在线观看| 最近最新中文字幕大全电影3| 久久精品影院6| 欧美+日韩+精品| 久久久久国内视频| 国产三级中文精品| 别揉我奶头 嗯啊视频| 亚洲国产欧洲综合997久久,| 日日摸夜夜添夜夜添av毛片 | 免费观看在线日韩| 两人在一起打扑克的视频| 男女做爰动态图高潮gif福利片| 亚洲人成网站高清观看| 我要看日韩黄色一级片| 在线a可以看的网站| 亚洲人与动物交配视频| 午夜福利欧美成人| 大又大粗又爽又黄少妇毛片口| 亚洲最大成人av| 成人高潮视频无遮挡免费网站| 久久久久国内视频| 神马国产精品三级电影在线观看| 成熟少妇高潮喷水视频| 长腿黑丝高跟| 欧美日韩瑟瑟在线播放| 精品欧美国产一区二区三| 日韩欧美精品v在线| 真实男女啪啪啪动态图| 超碰av人人做人人爽久久| 又紧又爽又黄一区二区| aaaaa片日本免费| 他把我摸到了高潮在线观看| 日韩在线高清观看一区二区三区 | 午夜免费男女啪啪视频观看 | 蜜桃久久精品国产亚洲av| 国产精品一区二区三区四区久久| 韩国av一区二区三区四区| 中国美女看黄片| 亚洲av中文av极速乱 | 久久精品影院6| 亚洲最大成人中文| 午夜精品久久久久久毛片777| 日本与韩国留学比较| 国产又黄又爽又无遮挡在线| 国产一区二区三区在线臀色熟女| 人妻制服诱惑在线中文字幕| 欧美色视频一区免费| 最好的美女福利视频网| 乱码一卡2卡4卡精品| 99九九线精品视频在线观看视频| 在线观看66精品国产| 美女 人体艺术 gogo| 在线观看一区二区三区| 亚洲av中文字字幕乱码综合| 欧美激情在线99| 午夜福利高清视频| 欧美国产日韩亚洲一区| 国产精品自产拍在线观看55亚洲| 老熟妇乱子伦视频在线观看| 国产美女午夜福利| 久久99热这里只有精品18| 免费av观看视频| 99riav亚洲国产免费| 亚洲乱码一区二区免费版| 国产精品一区二区免费欧美| 日本a在线网址| 色哟哟哟哟哟哟| 成人亚洲精品av一区二区| 成年女人永久免费观看视频| 亚洲专区国产一区二区| xxxwww97欧美| 日本精品一区二区三区蜜桃| 久久久久国产精品人妻aⅴ院| 真人一进一出gif抽搐免费| 久久精品国产自在天天线| 免费一级毛片在线播放高清视频| 成人午夜高清在线视频| 亚洲 国产 在线| 成人av一区二区三区在线看| 亚洲精品乱码久久久v下载方式| 色哟哟·www| 亚洲专区国产一区二区| xxxwww97欧美| 亚洲国产欧洲综合997久久,| 伦理电影大哥的女人| 99热6这里只有精品| 国产高清激情床上av| 日韩在线高清观看一区二区三区 | 少妇熟女aⅴ在线视频| 真人一进一出gif抽搐免费| 少妇丰满av| 人人妻人人看人人澡| 日本一二三区视频观看| 国产一区二区在线观看日韩| 少妇的逼水好多| 一本久久中文字幕| 深夜a级毛片| 日韩一区二区视频免费看| videossex国产| 无遮挡黄片免费观看| 午夜精品一区二区三区免费看| 舔av片在线| 精品一区二区三区人妻视频| 亚洲不卡免费看| 国产一区二区激情短视频| 日本成人三级电影网站| 久久草成人影院| 色综合站精品国产| 在线观看免费视频日本深夜| 99在线视频只有这里精品首页| 国产 一区精品| 啦啦啦韩国在线观看视频| 69人妻影院| 国产精品久久久久久久电影| 可以在线观看的亚洲视频| 久久人人精品亚洲av| 国产单亲对白刺激| 97人妻精品一区二区三区麻豆| 亚洲精品456在线播放app | 非洲黑人性xxxx精品又粗又长| 男人舔奶头视频| 久久精品国产99精品国产亚洲性色| 国产真实伦视频高清在线观看 | 一区福利在线观看| 精品久久久久久久久久免费视频| 精品午夜福利在线看| 日日摸夜夜添夜夜添小说| 国产黄色小视频在线观看| 小说图片视频综合网站|