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

    基于誤差估計(jì)的發(fā)動(dòng)機(jī)排氣聲源特性間接識別模型求解方法研究

    2016-09-18 02:45:40劉海濤連小珉劉林芽
    振動(dòng)與沖擊 2016年16期
    關(guān)鍵詞:測管聲壓聲源

    劉海濤, 連小珉, 劉林芽

    (1.華東交通大學(xué) 機(jī)電與車輛工程學(xué)院,南昌 330013; 2.清華大學(xué) 汽車安全與節(jié)能國家重點(diǎn)實(shí)驗(yàn)室,北京 100084)

    ?

    基于誤差估計(jì)的發(fā)動(dòng)機(jī)排氣聲源特性間接識別模型求解方法研究

    劉海濤1,2, 連小珉2, 劉林芽1

    (1.華東交通大學(xué) 機(jī)電與車輛工程學(xué)院,南昌330013; 2.清華大學(xué) 汽車安全與節(jié)能國家重點(diǎn)實(shí)驗(yàn)室,北京100084)

    用于發(fā)動(dòng)機(jī)聲源特性識別的線性時(shí)不變頻域模型仍然是較為實(shí)用的方法,但尚存在求解誤差大、對輸入誤差敏感、識別及求解過程復(fù)雜等問題。為了提高聲源識別的準(zhǔn)確性,建立了間接法聲源特性識別的誤差估計(jì)方法,即聲源聲壓離散度估計(jì)和聲源阻抗值偏差估計(jì)。將識別模型中超定非線性方程組中的方程式進(jìn)行兩兩組合,通過解析幾何的方法直接求解方程組,獲得聲源特性參數(shù)的多個(gè)解,再利用聲源阻抗值偏差估計(jì)選取最優(yōu)解。經(jīng)分析計(jì)算,獲取的聲源特性參數(shù)的誤差估計(jì)值要明顯低于傳統(tǒng)四負(fù)載法。最后,采用三維聲流耦合仿真方法對獲取的聲源特性參數(shù)進(jìn)行驗(yàn)證,仿真預(yù)測的遠(yuǎn)場響應(yīng)點(diǎn)聲壓值與試驗(yàn)結(jié)果吻合良好,表明基于誤差估計(jì)的聲源特性參數(shù)最優(yōu)選取法可以獲取更加準(zhǔn)確的發(fā)動(dòng)機(jī)聲源特性參數(shù),減少識別誤差和識別的復(fù)雜程度。

    發(fā)動(dòng)機(jī);聲源聲壓;聲源阻抗;誤差估計(jì);聲流耦合仿真

    排氣系統(tǒng)的性能與發(fā)動(dòng)機(jī)的特性密切相關(guān),二者之間存在流體耦合和聲學(xué)耦合。流體耦合決定了發(fā)動(dòng)機(jī)的排氣功率損失,而聲學(xué)耦合決定了尾管階次噪聲輻射能量的大小。為了更加有效的控制階次噪聲,需要準(zhǔn)確提取出發(fā)動(dòng)機(jī)聲源特性參數(shù),用于發(fā)動(dòng)機(jī)和排氣系統(tǒng)的匹配分析[1-2]。

    國內(nèi)外學(xué)者在發(fā)動(dòng)機(jī)聲源特性識別模型上做了大量的研究工作。有些學(xué)者將發(fā)動(dòng)機(jī)定義為時(shí)變模型,在時(shí)域內(nèi)計(jì)算發(fā)動(dòng)機(jī)的聲源特性參數(shù)[3-4]。時(shí)變模型能夠考慮發(fā)動(dòng)機(jī)排氣產(chǎn)生的高聲壓級脈動(dòng)的非線性效應(yīng),能夠更好地描述發(fā)動(dòng)機(jī)氣缸和排氣系統(tǒng)之間復(fù)雜的相互作用。時(shí)域計(jì)算模型需要獲取發(fā)動(dòng)機(jī)氣缸及相關(guān)附件的準(zhǔn)確幾何參數(shù),以及發(fā)動(dòng)機(jī)的性能參數(shù),如缸內(nèi)壓力、氣門升程、點(diǎn)火提前角等,但這些關(guān)鍵性能參數(shù)在實(shí)際排氣系統(tǒng)設(shè)計(jì)匹配工作中較難獲取[5]。

    發(fā)動(dòng)機(jī)聲源特性的線性頻域識別模型仍然在排氣系統(tǒng)的匹配設(shè)計(jì)中發(fā)揮重要作用[1,5]。線性頻域識別模型將發(fā)動(dòng)機(jī)簡化成線性時(shí)不變系統(tǒng),通過聲電類比模型,可以將發(fā)動(dòng)機(jī)的聲源特性用聲源聲壓和聲源阻抗兩個(gè)參數(shù)來描述[6]。國內(nèi)外學(xué)者提出了許多方法來確定兩個(gè)聲源特性參數(shù),主要可分為直接法和間接法兩種。直接法[6-8]測試時(shí)需要外部聲源激勵(lì)來獲取源特性參數(shù),但難以找到能夠耐受發(fā)動(dòng)機(jī)排氣管內(nèi)極其惡劣環(huán)境的外部聲源,因而直接法的應(yīng)用受到限制。間接法不需要外部聲源,只需要更換源特性識別負(fù)載測管,并獲取尾管口遠(yuǎn)場輻射噪聲即可,因而間接法在發(fā)動(dòng)機(jī)聲源特性識別中得到廣泛應(yīng)用。經(jīng)過多年發(fā)展,間接法發(fā)展出了多種方法,如雙負(fù)載法[9],三負(fù)載法[10],四負(fù)載法[11]以及多負(fù)載法[12-13],這些方法都有各自的優(yōu)缺點(diǎn)。雙負(fù)載法需要一個(gè)與噪聲產(chǎn)生過程相關(guān)的參考信號來計(jì)算復(fù)聲壓,但是這種參考信號不太容易選取。三負(fù)載法需要求解含有兩個(gè)二元二次方程的方程組,由于測試誤差以及模型簡化的影響,可能在某些頻率點(diǎn)下方程組的解不存在,從而限制了這種方法的應(yīng)用。四負(fù)載法成功地避免了三負(fù)載法沒有解的情況,也是國內(nèi)外學(xué)者研究最多的一種方法。四負(fù)載法可以得到含有三個(gè)二元二次方程的非線性方程組,PRASAD[11]將方程組中的二次項(xiàng)消去,使得三個(gè)二元二次方程簡化成兩個(gè)線性方程,再聯(lián)立求解獲得發(fā)動(dòng)機(jī)聲源阻抗結(jié)果。但是經(jīng)過國外學(xué)者SRIDHARA等[14]的研究,這種求解方法對于輸入量引入的誤差極其敏感,計(jì)算結(jié)果容易出現(xiàn)較大偏差。多負(fù)載方法通常測取十多個(gè)負(fù)載的尾管口輻射噪聲,再采用最小二乘法[12]或者數(shù)值迭代[13]的方法對源特性參數(shù)進(jìn)行數(shù)值求解,從一定程度上克服了傳統(tǒng)四負(fù)載法計(jì)算結(jié)果不穩(wěn)定的缺點(diǎn),提高了識別精度。但是這種方法需要大量重復(fù)的發(fā)動(dòng)機(jī)臺架試驗(yàn)測試,并進(jìn)行復(fù)雜的代數(shù)運(yùn)算去求得源特性參數(shù),會耗費(fèi)大量的時(shí)間和精力,而且缺少對結(jié)果的定量誤差分析。因而有必要研究更合理的聲源特性參數(shù)求解方式,并進(jìn)行定量的誤差估計(jì),以保證求解結(jié)果的準(zhǔn)確性。

    為了解決現(xiàn)有聲源特性識別求解方法存在的問題,本文從聲源特性識別誤差估計(jì)出發(fā),提出一種新的求解方法,相對于傳統(tǒng)四負(fù)載法,在不增加測試負(fù)載數(shù)量的基礎(chǔ)上提高聲源特性參數(shù)的識別精度。

    1 聲源特性間接識別理論模型

    線性頻域求解模型將穩(wěn)態(tài)工況下發(fā)動(dòng)機(jī)排氣噪聲的產(chǎn)生假定成線性時(shí)不變系統(tǒng),并將發(fā)動(dòng)機(jī)簡化成黑箱,用發(fā)動(dòng)機(jī)聲源聲壓和聲源阻抗兩個(gè)集總參數(shù)來表示發(fā)動(dòng)機(jī)的聲源特性。

    1.1集總參數(shù)類比模型

    如圖1(a)表示典型發(fā)動(dòng)機(jī)排氣系統(tǒng)簡化示意圖。通過聲電類比,發(fā)動(dòng)機(jī)排氣系統(tǒng)可以用恒壓源電路來表示,如圖1(b)所示。

    圖1 典型發(fā)動(dòng)機(jī)排氣系統(tǒng)及其聲電類比示意圖Fig.1 Typical engine exhaust system and its electroacoustic analogies

    圖1中,PE是源強(qiáng)度,ZE是源阻抗,PL是源與負(fù)載接口處的聲壓,ZL是負(fù)載阻抗,PR是距尾管口R處響應(yīng)點(diǎn)的輻射聲壓,ZR是尾管口處的輻射阻抗,M表示排氣系統(tǒng)入口處的氣流馬赫數(shù),T表示排氣系統(tǒng)入口處的溫度。在頻域識別模型中,源與負(fù)載之間的相互耦合作用可用式(1)表示。

    PL=PEZL/(ZE+ZL)

    (1)

    對于內(nèi)燃機(jī),只有當(dāng)聲源特性測試中各個(gè)負(fù)載產(chǎn)生的背壓相同或接近的時(shí)候,發(fā)動(dòng)機(jī)的聲源特性參數(shù)才不會隨著測試負(fù)載的變化而產(chǎn)生較大的差異[5],因而采用不同長度的直管作為聲學(xué)負(fù)載會比較合理。另外直管負(fù)載方便使用和加工,而且標(biāo)準(zhǔn)直管的負(fù)載阻抗比較容易精確獲取。對于不同長度的聲源特性測管,式(1)可以改寫成式(2)。

    (i=1,2,3,…,N)

    (2)

    式中:i表示序數(shù),N表示總共使用的測管數(shù)量。

    對于直管聲學(xué)負(fù)載,其進(jìn)出口聲學(xué)參數(shù)可使用四端極子參數(shù)聯(lián)系起來,如式(3)所示。

    (3)

    式中:p是聲壓,U是體積速度;AL,BL,CL,DL分別是直管的四端極子參數(shù)。根據(jù)MUNJAL[15]推導(dǎo)的有流情況下管道內(nèi)平面波的傳播理論,四端極子參數(shù)可以由式(4)表示。

    (4)

    式中:M是氣流馬赫數(shù),M=v/c,v是氣流流速,c是管內(nèi)聲速;k是波數(shù),k=w/c(1-M2),Z0是管內(nèi)氣體的特征阻抗;le=L+Δl,L是測管長度,Δl是聲學(xué)修正長度,采用熱噴射流下的管道末端修正求得[16-17]。

    發(fā)動(dòng)機(jī)排氣系統(tǒng)帶有四端極子參數(shù)的線性頻域類比模型圖,如圖2所示。

    圖2 帶四端極子參數(shù)的發(fā)動(dòng)機(jī)排氣系統(tǒng)聲電類比示意圖Fig.2 Electrical analogue for source-load system with four pole parameters

    圖2可由傳遞矩陣方程式來表示,如式(5)所示。

    (i=1,2,…,N)

    (5)

    式中:UE是發(fā)動(dòng)機(jī)聲源體積振速,Uo是測管末端的體積振速,ZR是有流情況下測管末端的輻射阻抗,可由式(6)[15]計(jì)算得到。

    (6)

    式中Rc是有流情況下尾管口的反射系數(shù),可以通過MUNT[18]的理論計(jì)算得到,如式(7)所示。

    (7)

    式中:θ是反射系數(shù)的相位,θ=π-2kΔl;F+(u)是渦流層位置變換,式(7)的詳細(xì)推導(dǎo)公式可參考文獻(xiàn)[18]。

    1.2聲源阻抗的求解

    在不同的負(fù)載情況下時(shí),將式(2)表示各個(gè)負(fù)載的等式進(jìn)行比值處理,消去聲源聲壓|PE|,得到式(8)。

    (m=1,2,…,N-1)

    (8)

    將源與負(fù)載接口處管內(nèi)聲壓的比值換算到尾管外響應(yīng)點(diǎn)輻射噪聲聲壓的比值,再代入式(8)中,即可獲取αm的值,如式(9)所示。

    (m=1,2,…,N-1)

    (9)

    聯(lián)立式(8)和式(9),根據(jù)PRASAD的公式推導(dǎo),可以得到用于聲源阻抗求解的二元二次非線性方程組,如式(10)[11]所示。

    (m=1,2,…,N-1)

    (10)

    式中:

    ZE=RE+jXE

    通過求解式(10)中的非線性方程組,可以獲得發(fā)動(dòng)機(jī)的源阻抗值。

    1.3聲源強(qiáng)度的求解

    發(fā)動(dòng)機(jī)聲源特性的試驗(yàn)測量通常是在消聲室中進(jìn)行,以避免外界噪聲環(huán)境的干擾。由于排氣氣流屬于低馬赫數(shù)流動(dòng),而且分析的頻率范圍位于低頻段,忽略測管尾端噴射流中氣流以及溫度梯度對末端輻射噪聲聲傳播的影響。測管管口的半徑遠(yuǎn)小于傳聲器到尾管口端的距離,末端噪聲輻射過程可以當(dāng)作點(diǎn)聲源來處理,末端遠(yuǎn)場輻射響應(yīng)點(diǎn)的聲強(qiáng)值可以由式(11)表示。

    (k0r<1,i=1,2,…,N)

    (11)

    式中:ρ0和c0分別為消聲室中空氣的密度和聲速,IR為響應(yīng)點(diǎn)處的聲強(qiáng),可由式(12)來計(jì)算。

    (12)

    式中:PR是尾管末端遠(yuǎn)場輻射響應(yīng)點(diǎn)傳聲器測量的聲壓。

    將式(11)代入式(12)中,可以得到測管末端尾管口處的體積振速,如式(13)所示。

    (13)

    聯(lián)立式(5)和式(13)可以獲得聲源聲壓PE和末端遠(yuǎn)場輻射響應(yīng)點(diǎn)的聲壓PR之間的關(guān)系,如式(14)所示。

    (i=1,2,…,N)

    (14)

    發(fā)動(dòng)機(jī)的聲源聲壓值可由式(14)計(jì)算得到。

    2 間接識別模型的誤差估計(jì)方法

    2.1聲源阻抗值偏差估計(jì)

    由式(10)可知,當(dāng)聲源阻抗結(jié)果沒有誤差時(shí),將其回代入聲源阻抗求解方程組中,可以滿足所有的方程。但當(dāng)聲源阻抗求解結(jié)果中帶有誤差時(shí),式(10)中的方程等式右邊將不再是零,而會出現(xiàn)一些偏差。而這個(gè)偏差可以用來估計(jì)聲源阻抗求解結(jié)果的誤差大小,并可以用來分析方程組中的參量,如測管的負(fù)載阻抗值ZL對聲源阻抗識別結(jié)果的影響。

    (m=1,2,…,N-1)

    (15)

    式(15)中可以看出,對于N個(gè)負(fù)載,將會產(chǎn)生N-1個(gè)偏差。為了方便統(tǒng)計(jì)偏差的大小,總的偏差可以由式(16)表示。

    (16)

    式中:E′t代表聲源阻抗求解方程的偏差估計(jì)函數(shù)值,可以用來量化源阻抗求解結(jié)果的誤差大小。

    2.2聲源聲壓值偏差估計(jì)

    由式(14)可以看出,聲源聲壓PE的結(jié)果主要由測管四端極子參數(shù)、測管末端輻射阻抗、聲源阻抗以及末端遠(yuǎn)場輻射響應(yīng)點(diǎn)聲壓決定。末端遠(yuǎn)場輻射響應(yīng)點(diǎn)聲壓通過試驗(yàn)測試得到,其精度主要由試驗(yàn)環(huán)境和測試系統(tǒng)來保證。而測管四端極子參數(shù)、測管末端輻射阻抗、聲源阻抗的誤差將直接影響聲源聲壓計(jì)算結(jié)果的精度。根據(jù)式(14)可知,每一個(gè)負(fù)載都可以計(jì)算出一個(gè)聲源聲壓值,如果沒有誤差,則理論上所有負(fù)載對應(yīng)的聲源聲壓值都相同。但是當(dāng)參與求解的參數(shù)存在誤差時(shí),每個(gè)負(fù)載計(jì)算出的聲源聲壓值就會有偏差,而且誤差越大,偏差值也會越大。因而每個(gè)負(fù)載計(jì)算出的聲源聲壓值之間的離散度可以用來估計(jì)由測管四端極子參數(shù)、測管末端輻射阻抗、聲源阻抗引入的誤差大小,從而用于聲源特性識別誤差的分析。聲源聲壓值的離散度可由式(17)表示。

    (17)

    式中:DPE表示聲源聲壓值的離散度,而EPE表示各個(gè)負(fù)載求出的聲源聲壓值的平均值,如式(18)所示。

    (18)

    聲源聲壓值的離散度DPE可以用來估計(jì)聲源特性識別模型中引入的誤差大小,而EPE作為最終求解得出的聲源聲壓值。

    3 聲源阻抗值最優(yōu)選取求解方法

    聲源阻抗值的獲取在于求解式(10)中的非線性超定方程組。傳統(tǒng)四負(fù)載法是將非線性方程組經(jīng)過替代轉(zhuǎn)化成線性方程組,從而求解出聲源阻抗值。但這種傳統(tǒng)求解方法,會使得求解結(jié)果對于輸入?yún)?shù)的誤差極為敏感,容易產(chǎn)生極大的誤差[14]。傳統(tǒng)四負(fù)載法產(chǎn)生誤差的核心問題是消除了式(10)中的二次項(xiàng),因而本小節(jié)從這點(diǎn)出發(fā),將超定非線性方程組中的方程式進(jìn)行兩兩組合,通過解析幾何的方法分別直接求解方程組,獲得多個(gè)解,然后結(jié)合第2節(jié)中的偏差估計(jì)方法選取最優(yōu)解,從而有效控制誤差。

    選用四種測管負(fù)載進(jìn)行計(jì)算分析,則式(10)可寫成方程組的形式,如式(19)所示。

    (19)

    從式(19)可知,非線性方程組從代數(shù)幾何的角度來看,其實(shí)質(zhì)是一系列圓方程。非線性方程組的解,就是這一系列圓的交點(diǎn)。當(dāng)沒有誤差干擾時(shí),三個(gè)圓將相交于一點(diǎn),即聲源阻抗值的精確解;而當(dāng)有輸入誤差干擾時(shí),三個(gè)圓方程將會有多個(gè)交點(diǎn)。選取式(19)中任意兩個(gè)方程,通過解析幾何的方法推導(dǎo)出交點(diǎn)的求解公式,兩圓相交的幾何示意圖如圖3所示。

    圖3 兩圓相交的幾何示意圖Fig.3 Geometric diagram for the intersection of two circles

    如圖3所示,O1(ai/2,bi/2),O2(aj/2,bj/2)分別表示兩個(gè)圓的圓心,其中i,j∈{1,2,3}。Ic(RE0,XE0)點(diǎn)是兩圓心連線與兩圓公共弦的交點(diǎn),Is1(RE1,XE1),Is2(RE2,XE2)分別是兩個(gè)圓的交點(diǎn)。

    由兩圓之間的平面幾何關(guān)系可得到Ic的坐標(biāo),如式(20)所示。

    (20)

    式中:lo12為兩個(gè)圓的圓心間距離,ro1和ro2分別為兩個(gè)圓的半徑。其中:

    Is1和Ic之間的距離lcs可以由式(21)表示。

    (21)

    根據(jù)相交兩圓的幾何關(guān)系可以求得兩圓交點(diǎn)的表達(dá)式,如式(22)和式(23)所示,式中Re表示取實(shí)部值。

    (22)

    (23)

    聯(lián)立以上等式,即可求解得到交點(diǎn)的坐標(biāo)值。當(dāng)兩圓相交時(shí),有兩個(gè)交點(diǎn);當(dāng)兩圓相切時(shí),求出的兩個(gè)交點(diǎn)的坐標(biāo)值相同,即為一個(gè)交點(diǎn);而當(dāng)兩圓相離時(shí),求出的兩個(gè)交點(diǎn)的坐標(biāo)值是復(fù)數(shù),這時(shí)沒有幾何意義上的交點(diǎn)。相離的情況下,應(yīng)該取距離兩圓最近的點(diǎn)作為圓方程的解,從而有效減小偏差。相離時(shí)式(22)和式(23)復(fù)數(shù)坐標(biāo)點(diǎn)的實(shí)部值相同,構(gòu)成的點(diǎn)位于圓心的連線上,處于兩圓的中間,距離兩圓是最近的點(diǎn),因而將復(fù)數(shù)坐標(biāo)點(diǎn)的實(shí)部值作為兩圓的虛交點(diǎn)。

    按對圓方程進(jìn)行兩兩組合的求解方法,獲得了多個(gè)交點(diǎn),而選取其中誤差最小的點(diǎn)作為發(fā)動(dòng)機(jī)聲源阻抗值是問題的關(guān)鍵。2.1小節(jié)中推導(dǎo)出了聲源阻抗的誤差估計(jì)方法,可以量化聲源阻抗求解結(jié)果的誤差大小。因而本文中以式(16)作為衡量標(biāo)準(zhǔn),選取使偏差估計(jì)函數(shù)值最小的交點(diǎn),可以用式(24)表示。

    (24)

    4 實(shí)際發(fā)動(dòng)機(jī)聲源特性參數(shù)求解分析

    4.1試驗(yàn)裝置

    選用某款四缸汽油發(fā)動(dòng)機(jī)進(jìn)行聲源特性識別試驗(yàn)。發(fā)動(dòng)機(jī)臺架采用電渦流測功機(jī),負(fù)載測管以及響應(yīng)點(diǎn)的傳聲器放置于半消聲室中,同時(shí)在尾管口的地面上鋪滿吸聲劈尖以減少地面的聲反射,從而使測試環(huán)境接近于自由場。試驗(yàn)裝置的示意圖如圖4所示。

    圖4 臺架試驗(yàn)裝置示意圖Fig.4 Experimental set-up

    響應(yīng)點(diǎn)距離測管的尾管口500 mm,與尾管口的軸心線成45°夾角。為了防止尾管口噴射的熱氣流對傳聲器產(chǎn)生影響,在傳聲器的頭部安裝了風(fēng)罩。試驗(yàn)現(xiàn)場的傳感器布置照片如圖5所示。

    圖5 測試裝置現(xiàn)場實(shí)物圖Fig.5 The photograph of test equipment

    聲壓測量采用的是PCB公司生產(chǎn)的377C01型自由場傳聲器,可以用來測量高聲壓。而數(shù)據(jù)采集采用LMS公司生產(chǎn)的SCM02型數(shù)據(jù)分析儀。

    采用了四個(gè)不同長度的直管當(dāng)作聲學(xué)負(fù)載用于發(fā)動(dòng)機(jī)聲源特性的測試。通過對壓力傳感器的信號進(jìn)行分析發(fā)現(xiàn),同一工況下,各種長度的直管負(fù)載所產(chǎn)的排氣壓力損失小于2 kPa,遠(yuǎn)小于排氣系統(tǒng)中三元催化器的排氣壓力損失,因而不同直管負(fù)載所產(chǎn)生的背壓對發(fā)動(dòng)機(jī)源特性的影響可以忽略。四種不同長度的直管負(fù)載的幾何參數(shù)如表1所示。

    表1 試驗(yàn)測管幾何參數(shù)

    4.2聲源特性參數(shù)計(jì)算分析

    選取四個(gè)測管作為聲學(xué)負(fù)載,根據(jù)線性頻域源特性求解模型,其構(gòu)成的三個(gè)圓方程的交點(diǎn)求解結(jié)果如圖6所示。圖6中,三個(gè)圓是由四個(gè)負(fù)載根據(jù)式(19)構(gòu)成的,橫軸是聲源阻抗的實(shí)部,而縱軸是聲源阻抗的虛部,并都進(jìn)行了歸一化處理。由于誤差的存在,圖中三個(gè)圓沒有相交于一點(diǎn),兩個(gè)圓相交,另外一圓相離。圖中的符號“*”標(biāo)識出了三個(gè)圓所有的交點(diǎn),包括相離情況下構(gòu)成的虛交點(diǎn)。

    圖6 圓方程交點(diǎn)求解結(jié)果圖Fig.6 The results of intersection points of circle equations

    如圖6所示,按對圓方程進(jìn)行兩兩組合的求解方法,獲得了多個(gè)交點(diǎn),而選取其中誤差最小的點(diǎn)作為發(fā)動(dòng)機(jī)聲源阻抗值是問題的關(guān)鍵。本文以式(24)作為衡量標(biāo)準(zhǔn),選取使偏差估計(jì)函數(shù)值最小的交點(diǎn),選取結(jié)果如圖7(a)所示。圖中的符號“*”標(biāo)識出了求解結(jié)果。

    圖7 發(fā)動(dòng)機(jī)聲源特性參數(shù)求解結(jié)果Fig.7 The results of engine source impedance

    圖7(b)顯示的是采用傳統(tǒng)四負(fù)載法的求解結(jié)果。如圖所示,按滿足式(24)的阻抗值最優(yōu)選取方法的計(jì)算結(jié)果基本位于三個(gè)圓接近相交的地方,而傳統(tǒng)四負(fù)載方法的求解結(jié)果則偏離三個(gè)圓很遠(yuǎn)。圖7充分說明,本小節(jié)提出的圓方程交點(diǎn)最優(yōu)選取的求解方法能夠有效控制輸入誤差對聲源特性識別結(jié)果的影響。

    按本文中提出的圓方程交點(diǎn)最優(yōu)選取法和傳統(tǒng)四負(fù)載法計(jì)算得到的發(fā)動(dòng)機(jī)聲源特性參數(shù)如圖8所示。

    圖8 聲源阻抗計(jì)算結(jié)果對比圖Fig.8 The comparison of engine source impedance solved by the two methods

    圖9 聲源聲壓計(jì)算結(jié)果對比圖Fig.9 The comparison of engine source strength solved by the two methods

    從圖8看出,兩種計(jì)算方式的聲源阻抗結(jié)果存在一定的差異,傳統(tǒng)四負(fù)載法的聲源阻抗值波動(dòng)明顯較大。從圖9可以看出,傳統(tǒng)四負(fù)載法計(jì)算的聲源聲壓值在階次峰值處略要高于本文圓交點(diǎn)最優(yōu)選取法的計(jì)算結(jié)果,另外在非階次峰值處還存在較多的波動(dòng)峰值。

    將傳統(tǒng)四負(fù)載法和本文提出的方法獲取的聲源阻抗結(jié)果分別代入式(16)和式(17)之中,可以獲取兩種方法求解結(jié)果的誤差估計(jì),如圖10所示。

    圖10 兩種求解方法結(jié)果的誤差估計(jì)Fig.10 The error estimation of the results of two methods

    求得的聲源阻抗值偏差估計(jì)如圖10(a)所示,聲源聲壓離散度估計(jì)如圖10(b)所示。從圖10(a)可以看出,傳統(tǒng)四負(fù)載法的阻抗值偏差估計(jì)結(jié)果在整個(gè)頻帶上都遠(yuǎn)大于圓交點(diǎn)最優(yōu)選取法。圖10(b)可以看出,傳統(tǒng)四負(fù)載法的聲源聲壓離散度估計(jì)值在全頻帶內(nèi)也高于圓交點(diǎn)最優(yōu)選取法,平均高出4 dB左右。

    圖7和圖10充分說明,傳統(tǒng)四負(fù)載法的求解結(jié)果偏差太大,準(zhǔn)確性難以保證,而本小節(jié)提出的圓交點(diǎn)最優(yōu)選取法,通過誤差估計(jì)使聲源特性求解結(jié)果更加穩(wěn)定,從而能夠有效對發(fā)動(dòng)機(jī)聲源識別誤差進(jìn)行控制。

    5 聲源特性識別結(jié)果驗(yàn)證

    本小節(jié)對傳統(tǒng)四負(fù)載法和本文提出的方法獲取的聲源特性結(jié)果進(jìn)行驗(yàn)證。驗(yàn)證方法為選取一個(gè)帶內(nèi)插管的膨脹腔進(jìn)行發(fā)動(dòng)機(jī)臺架試驗(yàn),測量得到尾管輻射噪聲的聲壓級。然后對帶內(nèi)插管的膨脹腔進(jìn)行三維聲流耦合仿真,將獲得的聲源特性結(jié)果施加到聲學(xué)仿真模型的入口邊界上,再計(jì)算尾管遠(yuǎn)場輻射聲壓級。最后將試驗(yàn)結(jié)果與仿真結(jié)果進(jìn)行對比即可驗(yàn)證兩種計(jì)算方式得到的發(fā)動(dòng)機(jī)聲源特性結(jié)果的準(zhǔn)確性。

    內(nèi)插管膨脹腔的幾何尺寸以及試驗(yàn)測試中的實(shí)物照片如圖11所示。

    圖11 內(nèi)插管膨脹腔的幾何尺寸及實(shí)物照片F(xiàn)ig.11 The physical dimension and photograph of the expansion chamber with insert tube

    采用CFD方法[19]對內(nèi)插管膨脹腔內(nèi)部的流場進(jìn)行計(jì)算,計(jì)算結(jié)果如圖12所示。再將流場計(jì)算結(jié)果通過積分插值運(yùn)算映射到聲場網(wǎng)格上,從而使得聲場的計(jì)算能夠充分考慮非均一流場場量的影響。從圖12可以看出,內(nèi)插管膨脹腔內(nèi)部呈現(xiàn)復(fù)雜的流場以及溫度場分布,必定會對管內(nèi)的聲傳播產(chǎn)生影響。因而將流場信息耦合到聲場計(jì)算中,可以大大減小仿真誤差。在完成流場信息插值計(jì)算以后,對聲場仿真模型施加激勵(lì),然后進(jìn)行耦合場情況下的聲傳播計(jì)算。

    圖12 內(nèi)插管膨脹腔內(nèi)部流場仿真結(jié)果Fig.12 The internal flow field simulation results of the expansion chamber with insert tube

    要將聲源特性結(jié)果施加到聲場仿真模型的入口邊界上,先要將聲源阻抗換算成聲源導(dǎo)納,將聲源聲壓換算成加速度激勵(lì)。換算方法如式(25)和式(26)所示。式中YE是入口邊界處的聲源導(dǎo)納,ɑE是入口處的加速度激勵(lì)。

    (25)

    (26)

    內(nèi)插管膨脹腔最終的聲場仿真模型如圖13所示。

    將仿真模型中遠(yuǎn)場響應(yīng)點(diǎn)在各階次處的輻射聲壓值提取出來,與試驗(yàn)測試結(jié)果進(jìn)行對比如圖14所示。

    圖13 聲源特性參數(shù)驗(yàn)證的聲學(xué)仿真模型Fig.13 The acoustic simulation model for the verification of engine source characteristics results

    圖14 內(nèi)插管膨脹腔遠(yuǎn)場輻射噪聲仿真結(jié)果與試驗(yàn)對比Fig.14 The comparison of the radiated sound pressure level between simulation and experiment for the expansion chamber with insert tube

    從圖14中可以看出,在各階次處傳統(tǒng)四負(fù)載法的遠(yuǎn)場輻射噪聲仿真結(jié)果與試驗(yàn)結(jié)果偏差較大,尤其在8階以下的結(jié)果明顯高于試驗(yàn)測試的結(jié)果。而本文中提出的基于誤差估計(jì)的圓交點(diǎn)最優(yōu)選取法的遠(yuǎn)場輻射噪聲仿真結(jié)果與試驗(yàn)結(jié)果吻合良好,充分說明本文中提出發(fā)動(dòng)機(jī)聲源特性求解方法可以有效地減小識別誤差,從而獲取更為準(zhǔn)確的發(fā)動(dòng)機(jī)聲源特性參數(shù)。

    6 結(jié) 論

    本文從誤差估計(jì)的角度出發(fā),建立了間接法聲源特性識別的聲源聲壓離散度估計(jì)和聲源阻抗值偏差估計(jì),將識別模型中的超定非線性方程組中的方程式進(jìn)行兩兩組合,通過解析幾何的方法直接求解方程組,然后再利用聲源阻抗值偏差估計(jì)選取最優(yōu)解。本文提出的方法克服了傳統(tǒng)四負(fù)載法通過替代消除方程組中的二次項(xiàng),將超定非線性方程組轉(zhuǎn)化成線性方程組所帶來的對輸入誤差敏感、容易產(chǎn)生極大誤差的缺點(diǎn)。最后通過仿真和試驗(yàn)相結(jié)合的方法驗(yàn)證了本文中提出的方法的有效性和準(zhǔn)確性。本文研究主要得到以下結(jié)論:

    (1) 建立的間接法聲源特性識別的誤差估計(jì)方法可對識別結(jié)果的誤差大小進(jìn)行估計(jì),為聲源特性參數(shù)求解結(jié)果的選取提供了判別依據(jù)。

    (2) 提出的基于誤差估計(jì)的圓交點(diǎn)最優(yōu)選取法可以使聲源特性求解結(jié)果更加穩(wěn)定,從而能夠有效控制輸入誤差對聲源特性識別結(jié)果的影響。

    (3) 采用仿真和試驗(yàn)相結(jié)合的方法可以對源特性識別結(jié)果進(jìn)行驗(yàn)證。驗(yàn)證結(jié)果對比表明本文提出的方法克服了傳統(tǒng)四負(fù)載法的缺陷,提高了發(fā)動(dòng)機(jī)聲源特性識別精度,為源特性參數(shù)的準(zhǔn)確識別提供了新方法。

    [1] DOKUMACI E. Prediction of source characteristics of engine exhaust manifolds[J]. Journal of Sound and Vibration, 2005, 280(3/4/5): 925-942.

    [2] 劉海濤,鄭四發(fā),康鐘緒,等.基于四負(fù)載方法的汽車發(fā)動(dòng)機(jī)排氣源特性研究[J].振動(dòng)工程學(xué)報(bào),2011,24(5):573-577.

    LIU Haitao, ZHENG Sifa, KANG Zhongxu, et al. Acoustical source characteristics of vehicle exhuast system based on the four-load method[J]. Journal of Vibration Engineering, 2011, 24(5):573-577.

    [3] BODEN H. The multiple load method for measuring the source characteristics of time-variant sources[J]. Journal of Sound and Vibration, 1991, 148(3): 437-453.

    [4] RAMMAL H, BODEN H. Modified multi-load method for nonlinear source characterisation[J]. Journal of Sound and Vibration, 2007, 299(4/5): 1094-1113.

    [5] MACIAN V, TORREGROSA A J, BROATCH A. A view on the internal consistency of linear source identification for I.C. engine exhaust noise prediction[J]. Mathematical and Computer Modelling, 2013, 57(7/8): 1867-1875.

    [6] GALAITISIS G A, BENDER K E. Measurement of the acoustic impedance of an internal combustion engine[J]. Journal of the Acoustical Society of America, 1975, 58(1): S8.

    [7] PARSAD G M, CROCKER J M. Acoustical source characteriztion studies on a multi-cylinder engine exhaust system[J]. Journal of Sound and Vibration, 1983, 90(4): 490-497.

    [8] ROSS F D, CROCKER J M. Measurement of the acoustic internal source impedance of an internal combustion engine[J]. Journal of the Acoutical Society of America,1983,74(1): 18-27.

    [9] DOIGE G A, ALVES S H. Experimental characterization of noise source for duct acoustics[J]. Journal of Vibration and Acoustics, 1989, 111(1): 108-114.

    [10] ALVES H S. Characterization of noise sources in ducts[D].Calgary, Canada: The University of Calgary, 1986.

    [11] PRASAD G M. A four load method for evaluation of acoustic source impedance in a duct[J]. Journal of Sound and Vibration, 1987, 114(2): 347-356.

    [12] BODEN H. On multi-load methods for determination of the source data of acoustic one-port sources[J]. Journal of Sound and Vibration, 1995, 180(5): 725-743.

    [13] JANG H S, IH G J. Refined multiload method for measuring acoustical source characteristics of an intake or exhaust system[J]. Journal of the Acoustical Society of America, 2000, 107(6): 3217-3225.

    [14] SRIDHARA S B, CROCKER J M. Error analysis for the four-load method used to measure the source impedance inducts[J]. The Journal of the Acoustical Society of America,1992,92(5):2924-2931.

    [15] MUNJAL L M. Acoustic of ducts and mufflers[M]. 2nd ed. New York: John Wiley, 2014.

    [16] RIENSTRA W S. A small strouhal number analysis for acoustic wave-jet flow-pipe interaction[J]. Journal of Sound and Vibration, 1983, 86(4): 539-556.

    [17] TIIKOJA H, LAVRENTJEV J, RAMMAL H. Experimental investigations of sound reflection from hot and subsonic flow duct termination[J]. Journal of Sound and Vibration, 2014, 333(3): 788-800.

    [18] MUNT R M. Acoustic transmission properties of a jetpipe with subsonic jet flow: I. the cold jet reflection coefficient[J]. Journal of Sound and Vibration, 1990, 142(3): 413-436.

    [19] 徐航手,季振林,康鐘緒. 抗性消聲器傳遞損失預(yù)測的三維時(shí)域計(jì)算方法[J]. 振動(dòng)與沖擊,2010, 29(4): 107-110.

    XU Hangshou, JI Zhenlin, KANG Zhongxu.Three-dimensional time-domain computational approach for predicting transmission loss of reactive silencers [J]. Journal of Vibration and Shock, 2010, 29(4): 107-110.

    On the solving method of the indirect source characteristics identification model for an engine exhaust system based on error estimation

    LIU Haitao1,2, LIAN Xiaomin2, LIU Linya1

    (1. School of Mechatronics & Vehicle Engineering, East China Jiaotong University, Nanchang 330013, China;2. State Key Laboratory of Automotive Safety and Energy, Tsinghua University, Beijing 100084, China)

    Linear time-invariant hypothetical model solved in the frequency domain for the determination of engine source characteristics plays a practical role in the engine exhaust system design. However, the error of existing methods for solving the model is large, the results are sensitive to input error, and the identification and solving processes are complicated. In order to improve the accuracy of engine source characteristics identification results, an error estimation method was established for the indirect identification model, namely source strength level dispersion estimation and source impedance deviation estimation. The equation from the over-determined nonlinear equations in the identification model was pairwise combined, which were directly solved by an analytic geometry method. Multiple source characteristic parameter solutions were obtained. And then the optimal solution was selected by the source impedance deviation estimation. According to analysis and calculation, the error value of the source characteristic parameters acquired by the method above, is significantly lower than the traditional four-load method. Finally, the three-dimensional coupled simulation of acoustic and flow method was used to verify the results of source characteristic identification. The sound pressure level of far-field response point predicted by simulations agrees well with experimental results, which indicates that the optimal source characteristic parameters selection method based on error estimation can achieve more accurate source identification results, and can efficiently reduce the identification error and complexity.

    engine; source strength; source impedance; error estimation; coupled simulation of acoustic and flow

    國家自然科學(xué)基金(51268014);華東交通大學(xué)科研基礎(chǔ)經(jīng)費(fèi)(26541022)

    2015-05-26修改稿收到日期:2015-08-01

    劉海濤 男,博士,講師,1986年1月生

    U464.134.4

    A

    10.13465/j.cnki.jvs.2016.16.016

    猜你喜歡
    測管聲壓聲源
    基于嘴唇處的聲壓數(shù)據(jù)確定人體聲道半徑
    虛擬聲源定位的等效源近場聲全息算法
    聲測管對聲波透射法檢測樁身完整性的影響
    基于GCC-nearest時(shí)延估計(jì)的室內(nèi)聲源定位
    電子制作(2019年23期)2019-02-23 13:21:12
    車輛結(jié)構(gòu)噪聲傳遞特性及其峰值噪聲成因的分析
    汽車工程(2018年12期)2019-01-29 06:46:36
    螺旋式聲測管在基樁工程的應(yīng)用
    運(yùn)用內(nèi)積相關(guān)性結(jié)合迭代相減識別兩點(diǎn)聲源
    基于GIS內(nèi)部放電聲壓特性進(jìn)行閃絡(luò)定位的研究
    電測與儀表(2016年9期)2016-04-12 00:30:02
    高速鐵路橋梁樁基聲測管防堵的控制措施
    力-聲互易在水下聲源強(qiáng)度測量中的應(yīng)用
    色播亚洲综合网| 久久欧美精品欧美久久欧美| 免费看a级黄色片| 搡女人真爽免费视频火全软件| 97超视频在线观看视频| 欧美一区二区国产精品久久精品| 国产精华一区二区三区| 亚洲精品乱码久久久v下载方式| 国产精品野战在线观看| 在线天堂最新版资源| 黄片无遮挡物在线观看| 麻豆国产97在线/欧美| 欧美xxxx黑人xx丫x性爽| 干丝袜人妻中文字幕| 日韩一区二区视频免费看| 黄色日韩在线| 国产一区二区在线观看日韩| 69人妻影院| 久久久久性生活片| 国产在线一区二区三区精 | 国产午夜福利久久久久久| 黄片wwwwww| 久久久精品欧美日韩精品| 日韩亚洲欧美综合| 观看免费一级毛片| 伦理电影大哥的女人| 亚洲乱码一区二区免费版| 国产高清视频在线观看网站| 大香蕉97超碰在线| 亚洲国产精品成人久久小说| 欧美成人一区二区免费高清观看| 极品教师在线视频| 国产黄片视频在线免费观看| 一个人观看的视频www高清免费观看| 卡戴珊不雅视频在线播放| 久99久视频精品免费| 一卡2卡三卡四卡精品乱码亚洲| 99久久成人亚洲精品观看| 少妇的逼水好多| 欧美高清成人免费视频www| www.av在线官网国产| 亚洲色图av天堂| 国产熟女欧美一区二区| 亚洲成人av在线免费| 国产精品国产三级专区第一集| 中文欧美无线码| 国产真实伦视频高清在线观看| 美女大奶头视频| 最新中文字幕久久久久| 五月玫瑰六月丁香| 国产视频首页在线观看| av视频在线观看入口| 日韩高清综合在线| 久久久国产成人免费| 麻豆精品久久久久久蜜桃| 男人狂女人下面高潮的视频| 中文字幕制服av| 一级爰片在线观看| 国产精品一区二区在线观看99 | 亚洲电影在线观看av| 日韩精品有码人妻一区| 69人妻影院| 99热6这里只有精品| 国产精品电影一区二区三区| 国产高潮美女av| 国产成人免费观看mmmm| 99国产精品一区二区蜜桃av| 精品久久久久久成人av| 岛国毛片在线播放| 日本五十路高清| 国产一区二区在线av高清观看| 一边摸一边抽搐一进一小说| 99久国产av精品| 秋霞伦理黄片| 在线播放无遮挡| 一边亲一边摸免费视频| 亚洲综合色惰| or卡值多少钱| 免费大片18禁| 国产精品久久视频播放| 熟女人妻精品中文字幕| 国产精品一二三区在线看| 亚洲精品一区蜜桃| 一级爰片在线观看| 免费av观看视频| 99久久九九国产精品国产免费| 人体艺术视频欧美日本| 大话2 男鬼变身卡| 亚洲自拍偷在线| 亚洲av男天堂| 午夜福利在线观看吧| 午夜精品国产一区二区电影 | 有码 亚洲区| 国产片特级美女逼逼视频| 久久久a久久爽久久v久久| 九九在线视频观看精品| www.色视频.com| 69av精品久久久久久| 18禁动态无遮挡网站| 国产成人一区二区在线| 欧美日韩在线观看h| 听说在线观看完整版免费高清| 在线免费观看不下载黄p国产| 亚洲丝袜综合中文字幕| 只有这里有精品99| 国产精品三级大全| 色视频www国产| 日韩成人伦理影院| 中文欧美无线码| 嫩草影院新地址| 免费看av在线观看网站| 黄色配什么色好看| 久久精品久久久久久噜噜老黄 | 国产白丝娇喘喷水9色精品| 亚洲精品乱码久久久久久按摩| АⅤ资源中文在线天堂| 亚洲人与动物交配视频| 欧美+日韩+精品| 国产午夜精品久久久久久一区二区三区| 国产视频首页在线观看| 免费看av在线观看网站| 免费看光身美女| 精品久久久久久久久亚洲| 99国产精品一区二区蜜桃av| 国产 一区 欧美 日韩| 白带黄色成豆腐渣| 国产成年人精品一区二区| 老司机影院成人| 精品久久久久久久久久久久久| 人体艺术视频欧美日本| 有码 亚洲区| 久久精品夜色国产| 内地一区二区视频在线| 三级男女做爰猛烈吃奶摸视频| 丰满人妻一区二区三区视频av| 精品欧美国产一区二区三| 十八禁国产超污无遮挡网站| 日韩精品青青久久久久久| 久久综合国产亚洲精品| 亚洲精品一区蜜桃| 久久精品人妻少妇| av卡一久久| 丝袜喷水一区| 亚洲精品亚洲一区二区| 噜噜噜噜噜久久久久久91| 婷婷色av中文字幕| 中文亚洲av片在线观看爽| 亚洲熟妇中文字幕五十中出| 淫秽高清视频在线观看| 如何舔出高潮| 在线观看66精品国产| 日韩大片免费观看网站 | 熟妇人妻久久中文字幕3abv| 精品一区二区免费观看| 日韩 亚洲 欧美在线| 成人欧美大片| 成人毛片a级毛片在线播放| 成人特级av手机在线观看| 熟女人妻精品中文字幕| 国产成人午夜福利电影在线观看| 欧美日韩在线观看h| 日韩成人av中文字幕在线观看| 亚洲精品456在线播放app| www.av在线官网国产| 亚洲欧美日韩无卡精品| 一二三四中文在线观看免费高清| 国产精品永久免费网站| 日本免费一区二区三区高清不卡| 一个人免费在线观看电影| 一卡2卡三卡四卡精品乱码亚洲| 国产av码专区亚洲av| 亚洲美女视频黄频| 久久久国产成人免费| 国产又色又爽无遮挡免| 看非洲黑人一级黄片| 91精品伊人久久大香线蕉| 尤物成人国产欧美一区二区三区| 国产乱人视频| 性插视频无遮挡在线免费观看| 老司机福利观看| 国产精品1区2区在线观看.| 丰满人妻一区二区三区视频av| 在线免费十八禁| 国产在线男女| av免费观看日本| 久久午夜福利片| 久久久精品大字幕| 亚洲精品乱码久久久久久按摩| 国产精品爽爽va在线观看网站| 成人亚洲精品av一区二区| 亚洲欧美日韩东京热| 国产精品av视频在线免费观看| 欧美成人一区二区免费高清观看| 亚洲丝袜综合中文字幕| 日产精品乱码卡一卡2卡三| 最近中文字幕高清免费大全6| 日韩欧美国产在线观看| 日韩精品青青久久久久久| 日本爱情动作片www.在线观看| 亚洲国产精品国产精品| 亚洲精品aⅴ在线观看| 男女边吃奶边做爰视频| 欧美成人午夜免费资源| 日韩av在线免费看完整版不卡| 亚洲成人精品中文字幕电影| 村上凉子中文字幕在线| 久久99精品国语久久久| 欧美bdsm另类| 视频中文字幕在线观看| 中文在线观看免费www的网站| 女人久久www免费人成看片 | 亚洲精品乱码久久久v下载方式| 日本黄大片高清| 最近视频中文字幕2019在线8| 91精品国产九色| 国产老妇女一区| 国国产精品蜜臀av免费| 成年av动漫网址| 寂寞人妻少妇视频99o| 精品久久久久久久久av| 晚上一个人看的免费电影| 在线观看66精品国产| 日本猛色少妇xxxxx猛交久久| 麻豆乱淫一区二区| .国产精品久久| 一边亲一边摸免费视频| 晚上一个人看的免费电影| 亚洲中文字幕日韩| 麻豆国产97在线/欧美| 久久精品国产99精品国产亚洲性色| 亚洲成色77777| 亚洲av福利一区| 建设人人有责人人尽责人人享有的 | 免费看av在线观看网站| 三级男女做爰猛烈吃奶摸视频| 床上黄色一级片| 久久精品国产99精品国产亚洲性色| 小说图片视频综合网站| 深爱激情五月婷婷| 精品99又大又爽又粗少妇毛片| 熟女电影av网| av在线天堂中文字幕| 亚洲精品乱码久久久v下载方式| 国产成人精品久久久久久| 乱人视频在线观看| 国产日韩欧美在线精品| 精华霜和精华液先用哪个| 亚洲av免费高清在线观看| 亚洲在久久综合| 亚洲精品乱久久久久久| 亚洲人成网站在线播| 国产午夜福利久久久久久| 欧美一区二区亚洲| 亚洲美女搞黄在线观看| 亚洲精品亚洲一区二区| 午夜精品在线福利| 日韩一区二区视频免费看| 日本色播在线视频| 亚洲欧美清纯卡通| 国产成人一区二区在线| 久久久久久国产a免费观看| 免费黄网站久久成人精品| 久久欧美精品欧美久久欧美| 18禁裸乳无遮挡免费网站照片| 级片在线观看| 日韩中字成人| 国产精品一二三区在线看| 夜夜看夜夜爽夜夜摸| 免费观看性生交大片5| 成人欧美大片| 中文精品一卡2卡3卡4更新| 午夜福利成人在线免费观看| av播播在线观看一区| 春色校园在线视频观看| 看免费成人av毛片| 视频中文字幕在线观看| 中文资源天堂在线| 淫秽高清视频在线观看| 99久久精品一区二区三区| 日日撸夜夜添| 亚洲精品一区蜜桃| 久久精品夜夜夜夜夜久久蜜豆| 国产精品国产三级国产专区5o | 麻豆精品久久久久久蜜桃| 寂寞人妻少妇视频99o| 一个人观看的视频www高清免费观看| 91精品一卡2卡3卡4卡| 国产精品一区二区三区四区久久| 中文字幕人妻熟人妻熟丝袜美| 亚洲欧美日韩无卡精品| 成年免费大片在线观看| 午夜福利在线在线| 久久久国产成人免费| 成人亚洲欧美一区二区av| 免费观看精品视频网站| 国产精品国产三级国产av玫瑰| 国产精品蜜桃在线观看| 人人妻人人澡欧美一区二区| 最近2019中文字幕mv第一页| 我的老师免费观看完整版| 嫩草影院入口| 国产精品国产高清国产av| 舔av片在线| 老女人水多毛片| 久久欧美精品欧美久久欧美| 我的女老师完整版在线观看| 欧美精品一区二区大全| 久久久成人免费电影| 一级毛片电影观看 | 欧美xxxx性猛交bbbb| av福利片在线观看| 亚洲熟妇中文字幕五十中出| 亚洲欧美一区二区三区国产| 小蜜桃在线观看免费完整版高清| 国内揄拍国产精品人妻在线| 国产黄色视频一区二区在线观看 | 丰满乱子伦码专区| 一级av片app| 美女脱内裤让男人舔精品视频| 亚州av有码| 99热这里只有是精品在线观看| 免费看光身美女| av天堂中文字幕网| www日本黄色视频网| av视频在线观看入口| 建设人人有责人人尽责人人享有的 | 成人综合一区亚洲| 欧美激情久久久久久爽电影| 精品人妻偷拍中文字幕| 亚洲av成人精品一二三区| 国产精品熟女久久久久浪| 国产一级毛片在线| 欧美高清成人免费视频www| 五月伊人婷婷丁香| 99九九线精品视频在线观看视频| 在线免费十八禁| 国产午夜精品论理片| 日本免费一区二区三区高清不卡| 2021天堂中文幕一二区在线观| 长腿黑丝高跟| 久久99热6这里只有精品| 久久久久久久久久久免费av| 亚洲av电影在线观看一区二区三区 | 成人午夜高清在线视频| 亚洲aⅴ乱码一区二区在线播放| 久久久久性生活片| 国产91av在线免费观看| 丝袜喷水一区| 国产精品一区二区在线观看99 | 日本爱情动作片www.在线观看| 亚洲精品乱码久久久v下载方式| 国产午夜精品论理片| 天天躁夜夜躁狠狠久久av| 亚洲国产欧洲综合997久久,| 综合色丁香网| 国产精品一区www在线观看| 国产精品电影一区二区三区| 欧美bdsm另类| 日韩 亚洲 欧美在线| 搞女人的毛片| 亚洲精品国产成人久久av| 亚洲av电影不卡..在线观看| 欧美一区二区精品小视频在线| av女优亚洲男人天堂| 久久久欧美国产精品| 成年女人永久免费观看视频| 99久久人妻综合| 国产女主播在线喷水免费视频网站 | 国产爱豆传媒在线观看| 国产免费视频播放在线视频 | 91狼人影院| 亚洲欧美成人综合另类久久久 | 午夜福利在线观看免费完整高清在| 国产国拍精品亚洲av在线观看| 一级爰片在线观看| 精品久久久久久久久亚洲| 国内精品宾馆在线| 看十八女毛片水多多多| 舔av片在线| 日本三级黄在线观看| 男人舔奶头视频| 亚洲aⅴ乱码一区二区在线播放| 国产毛片a区久久久久| 国产av不卡久久| 老司机福利观看| 最新中文字幕久久久久| 狂野欧美白嫩少妇大欣赏| 国产精品久久久久久久电影| 日本av手机在线免费观看| av在线老鸭窝| 麻豆国产97在线/欧美| 中文字幕制服av| 免费播放大片免费观看视频在线观看 | 综合色av麻豆| videossex国产| 亚洲美女搞黄在线观看| 日本免费在线观看一区| 国产亚洲最大av| 国产精品不卡视频一区二区| 国产精品一及| av线在线观看网站| 免费看a级黄色片| av卡一久久| 国产精品人妻久久久久久| 亚洲最大成人中文| 午夜a级毛片| 国产精品一区www在线观看| 麻豆成人av视频| 久久精品综合一区二区三区| 亚洲一区高清亚洲精品| 国产一级毛片七仙女欲春2| 亚洲av.av天堂| 国产精品蜜桃在线观看| 国产极品精品免费视频能看的| 久久午夜福利片| 国产精品国产三级国产专区5o | 高清av免费在线| 日本黄色片子视频| 精品一区二区三区人妻视频| 中文在线观看免费www的网站| 午夜福利在线在线| 国产精品久久久久久av不卡| 久久久久久久久中文| 午夜福利在线在线| 国产成人福利小说| av黄色大香蕉| 国产免费一级a男人的天堂| 免费搜索国产男女视频| 天堂影院成人在线观看| 日本-黄色视频高清免费观看| 国产高清有码在线观看视频| 免费看av在线观看网站| 乱码一卡2卡4卡精品| 亚洲电影在线观看av| 少妇的逼好多水| 亚洲欧美中文字幕日韩二区| 黄色日韩在线| 国产黄a三级三级三级人| 久久久久久久久久久丰满| 日本色播在线视频| 精品无人区乱码1区二区| 大话2 男鬼变身卡| 欧美xxxx黑人xx丫x性爽| 亚洲国产精品合色在线| 国产一区二区亚洲精品在线观看| 久久人妻av系列| 久久亚洲精品不卡| 久久精品国产99精品国产亚洲性色| 国产乱人视频| 国产淫片久久久久久久久| 国产熟女欧美一区二区| 嫩草影院入口| 欧美激情久久久久久爽电影| 超碰av人人做人人爽久久| 亚洲人成网站在线观看播放| 天堂√8在线中文| 美女高潮的动态| 中文资源天堂在线| 国产男人的电影天堂91| 成人亚洲精品av一区二区| 中文字幕av成人在线电影| 91久久精品国产一区二区三区| 国产毛片a区久久久久| 国产高清不卡午夜福利| 久久人人爽人人爽人人片va| 色播亚洲综合网| 国产av一区在线观看免费| 亚洲第一区二区三区不卡| 嫩草影院精品99| 国产精品国产高清国产av| 在线观看美女被高潮喷水网站| 国产精品1区2区在线观看.| 国产亚洲最大av| 内射极品少妇av片p| 国产精品日韩av在线免费观看| 日本三级黄在线观看| 爱豆传媒免费全集在线观看| 老司机福利观看| 爱豆传媒免费全集在线观看| 国产大屁股一区二区在线视频| 国产av不卡久久| 国产av一区在线观看免费| 有码 亚洲区| 成人综合一区亚洲| 久久精品综合一区二区三区| 国产一区有黄有色的免费视频 | 国产精品久久久久久精品电影小说 | 国产大屁股一区二区在线视频| 成人国产麻豆网| 黄色欧美视频在线观看| 国产高清视频在线观看网站| 国产伦一二天堂av在线观看| 变态另类丝袜制服| 91精品一卡2卡3卡4卡| 五月玫瑰六月丁香| 国产一区有黄有色的免费视频 | 波野结衣二区三区在线| 亚洲欧美日韩卡通动漫| 亚洲第一区二区三区不卡| 欧美成人a在线观看| 日韩欧美 国产精品| 国产成人午夜福利电影在线观看| 精品国产三级普通话版| 国产成人a区在线观看| 哪个播放器可以免费观看大片| 精品人妻一区二区三区麻豆| 亚洲色图av天堂| 午夜精品在线福利| 天天一区二区日本电影三级| 亚洲精品国产av成人精品| 成年版毛片免费区| 床上黄色一级片| 老师上课跳d突然被开到最大视频| 日本猛色少妇xxxxx猛交久久| 精品久久久久久久久久久久久| 黄片无遮挡物在线观看| 在线天堂最新版资源| 91aial.com中文字幕在线观看| 一区二区三区乱码不卡18| 国产伦理片在线播放av一区| 一级毛片久久久久久久久女| 十八禁国产超污无遮挡网站| 久久精品人妻少妇| 成人鲁丝片一二三区免费| 波多野结衣高清无吗| 成人av在线播放网站| 亚洲无线观看免费| 特级一级黄色大片| 日韩一本色道免费dvd| 2021少妇久久久久久久久久久| 日本与韩国留学比较| 久久精品国产亚洲网站| 日本wwww免费看| 国产人妻一区二区三区在| 村上凉子中文字幕在线| 精品人妻熟女av久视频| 少妇被粗大猛烈的视频| 精品欧美国产一区二区三| 免费观看a级毛片全部| 欧美一级a爱片免费观看看| 久久人人爽人人爽人人片va| 午夜福利成人在线免费观看| 国产精品三级大全| 黄色一级大片看看| 春色校园在线视频观看| 一区二区三区高清视频在线| 亚洲欧美精品综合久久99| av在线蜜桃| 午夜福利成人在线免费观看| 色综合色国产| 最后的刺客免费高清国语| 亚洲精品亚洲一区二区| 最近视频中文字幕2019在线8| 看十八女毛片水多多多| 日韩欧美精品免费久久| 成人亚洲欧美一区二区av| 人妻系列 视频| 国产精品国产三级专区第一集| 午夜久久久久精精品| 床上黄色一级片| 97超碰精品成人国产| 亚洲综合精品二区| 网址你懂的国产日韩在线| 亚洲国产精品成人久久小说| 日韩欧美 国产精品| 亚洲av电影不卡..在线观看| 久久久久久久午夜电影| 国产色婷婷99| 亚洲av福利一区| 免费黄网站久久成人精品| 菩萨蛮人人尽说江南好唐韦庄 | 欧美成人精品欧美一级黄| 在线观看av片永久免费下载| 免费av毛片视频| a级毛色黄片| 亚洲最大成人中文| 亚洲精品影视一区二区三区av| 久久久久久久亚洲中文字幕| 久久99热这里只有精品18| 精品无人区乱码1区二区| 国产精品国产三级专区第一集| 亚洲国产欧洲综合997久久,| 中文字幕亚洲精品专区| 只有这里有精品99| 久久草成人影院| 日韩av在线免费看完整版不卡| 美女高潮的动态| 精品熟女少妇av免费看| 亚洲欧美中文字幕日韩二区| 精品国产露脸久久av麻豆 | 舔av片在线| 日韩,欧美,国产一区二区三区 | 少妇人妻精品综合一区二区| 亚洲成人精品中文字幕电影| 国内精品宾馆在线| 一个人看视频在线观看www免费| 中文字幕av成人在线电影| 三级毛片av免费| 欧美三级亚洲精品| 少妇猛男粗大的猛烈进出视频 | 好男人视频免费观看在线| av黄色大香蕉| 内射极品少妇av片p| av在线观看视频网站免费| 国产精品熟女久久久久浪| 国产精品国产高清国产av| 日本一二三区视频观看| 精品国内亚洲2022精品成人| 国产一区二区在线av高清观看| 国产免费又黄又爽又色| 亚洲怡红院男人天堂| 一本久久精品| 成人综合一区亚洲| 久久久精品欧美日韩精品| av黄色大香蕉|