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

    結構邊界受擾動下翼柱與環(huán)槽型固體火箭發(fā)動機燃燒室壓力振蕩特性研究①

    2023-04-26 01:56:00陳林君侯凱宇
    固體火箭技術 2023年1期
    關鍵詞:燃燒室裝藥振幅

    陳林君,曹 軍,程 書,石 欽,明 爽,郜 冶,侯凱宇

    (1.中國航天科技集團有限公司四院四十一所,西安 710025;2.中國航天科技集團有限公司第四研究院,西安 710025;3.哈爾濱工程大學 航天與建筑工程學院,哈爾濱 150000;4.上海航天技術研究院,上海 201100)

    0 引言

    固體火箭發(fā)動機在航空航天領域和導彈武器的發(fā)射任務中被廣泛使用[1]。為了不斷提高固體火箭發(fā)動機的性能,一直以來研究人員在設計和優(yōu)化上付出了許多努力[2-7]。近年來,為了最大限度地提高發(fā)動機性能,高能推進劑逐漸成為主流,推進劑的裝填系數和發(fā)動機的長徑比不斷提高,這導致固體火箭發(fā)動機不穩(wěn)定燃燒現象頻發(fā)。而這一類型的發(fā)動機具有一些共同特點:發(fā)動機的穩(wěn)定性在整個工作過程會不斷變化,往往在外部載荷或內部壓強脈沖激發(fā)下會出現不穩(wěn)定燃燒現象。

    目前,出現不穩(wěn)定燃燒的發(fā)動機主要有兩種:一種是大推力運載火箭采用的固體火箭發(fā)動機[8-13];另一種是大長徑比的戰(zhàn)術固體火箭發(fā)動機[14-16]。對于大尺寸的固體火箭發(fā)動機而言,產生不穩(wěn)定燃燒的主要為裝藥結構復雜的發(fā)動機,且多為分段式的發(fā)動機。由于裝藥段之間有絕熱擋板,從而易引發(fā)流場中出現渦脫落現象,繼而與燃燒室聲場進行耦合,出現渦-聲耦合導致的不穩(wěn)定燃燒。而對于大長徑比的固體火箭發(fā)動機不穩(wěn)定燃燒的形成機理,學術界并未形成統(tǒng)一的定論。有的學者認為不穩(wěn)定燃燒現象僅由燃燒和聲場相互作用引起的[17-18]。近年來,國內大長徑比的固體火箭發(fā)動機發(fā)生的不穩(wěn)定燃燒現象主要特征是燃燒室內的壓力振蕩可以達到幾百千帕,而導彈結構振動的振幅從發(fā)動機開始向全彈逐步減弱,且結構振動的頻率為燃燒室聲腔基頻或者倍頻,這一現象可以認為燃燒室內的流動同結構發(fā)生了相互耦合振動加強的現象。根據學者們的研究,結構振動的參與可以使燃燒室內的壓力振蕩幅值達到幾百千帕[19-20]。因此,在分析不穩(wěn)定燃燒中需要考慮結構振動所起到的作用。

    從已有的研究成果來看,很少有學者從結構振動和燃燒室內壓力振蕩相結合的角度對不穩(wěn)定燃燒開展研究。因此,本文結合發(fā)動機結構瞬態(tài)響應、推進劑燃速以及燃燒室流場流動,對固體火箭發(fā)動機在受結構振動激勵擾動下燃燒室內壓力振蕩特性進行研究,目的在于從結構振動和燃燒室壓力振蕩相互作用的角度去理解不穩(wěn)定燃燒形成的機理和過程。從而為設計人員在分析不穩(wěn)定燃燒現象產生原因以及尋找抑制不穩(wěn)定燃燒現象的方法時提供理論指導。這對解決我國近年來大長徑比的戰(zhàn)術導彈固體火箭發(fā)動機的不穩(wěn)定燃燒頻發(fā)的現象[14],加快研究進程具有重要的工程指導價值。

    1 數值模型

    1.1 推進劑燃速模型

    根據GREATRIX等的研究[21],推進劑表面的徑向加速度會對推進劑的燃速產生顯著的影響,因而GREATRIX根據實驗數據提出了一種用于描述推進劑燃速受加速度影響的模型。該模型的基本原理是:認為加速度提供了壓縮燃燒區(qū)域的作用,特別是在高密度、低流速接近凝聚相分解層的區(qū)域,這是法向加速度能夠增大推進劑燃速的主要作用機理。下面給出燃速計算模型:

    (1)

    式中cp為燃氣比定壓熱容,J/(kg·K);Tf為火焰溫度,K;Ts為推進劑表面溫度,K;cs為推進劑比熱容,J/(kg·K);ΔHs為推進劑表面熱釋放量,J/kg;Ti為推進劑初始溫度,K;ρs為推進劑密度,kg/m3;rb為推進劑總燃速,m/s;λ為燃氣熱導率,W/(m·K)。

    能量參考層厚度δ0和加速度質量流Ga的表達式如下:

    (2)

    (3)

    其中,R為特定氣體常數,J/(kg·K);an燃面為法向加速度,m/s2;r0為推進劑基礎燃速,m/s;φd為方向角,其表達式如下:

    (4)

    式中al為切向加速度,m/s2;K為比例系數。

    需要注意的是,本文的燃速模型中假設:當加速度的方向指向推進劑表面時,an為負值。根據實驗可知[22],當加速度的方向背離推進劑燃面時,對燃速不起作用,因此最終可使推進劑燃速在徑向振蕩加速度作用下具有凈增加值。為了求解受加速度影響的燃速,根據燃面附近的流體參數和燃面的結構振動加速度對式(1)~式(4)進行迭代求解。上述計算模型中的各參數值保持同GREATRIX等的研究中相一致[21]。

    1.2 結構邊界受擾動下固體火箭發(fā)動機燃燒室振蕩壓強計算模型

    燃燒室在實際的工作過程中會由于邊界受到擾動產生壓力振蕩現象,而結構振動激勵對燃燒室內流場的邊界擾動起到了重要的作用。為了計算燃燒室壓力在結構激勵作用下的振蕩特性,需要將推進劑的燃速和結構的振動相耦合,即計算域中的質量流率入口受結構振動的影響。結構振動對推進劑的燃速影響表現為燃速受推進劑內表面結構振動加速度的作用,這一點已經在1.1節(jié)中給出了半經驗計算公式。因此,在實際的計算中需要知道推進劑內表面的結構振動加速度,同時燃燒室內表面的壓力波動載荷對結構振動的計算來說同樣重要。為了節(jié)省計算資源,且保證計算網格具有足夠的密度,本文的流體計算和瞬態(tài)動力學計算中均使用了二維軸對稱模型。所以,本文提出了一套計算模型用以支持二維模型的流固耦合計算,圖1所示為本文設計的計算方法的流程示意圖。

    圖1 計算模型流程圖

    本文中的結構載荷作用到發(fā)動機結構后通過瞬態(tài)動力學計算模塊計算得到燃面表面的加速度,然后通過推進劑燃速計算模型得到對應的燃速,并以此計算質量流率,從而實現將過載激勵引入計算中以及流場壓力振蕩和結構振動相耦合的過程。

    為了簡化計算,本文設計的流固耦合計算模型為弱耦合,即流固交界面之間的位移沒有聯系,只在交界面上進行加速度和壓強這兩組數據的相互傳遞。同時,本文的研究思路為狀態(tài)凝固化,即取發(fā)動機在整個工作過程中特定的幾個時刻的模型進行瞬態(tài)數值模擬,且在整個計算過程中計算域的網格不發(fā)生變化。除此之外,發(fā)動機的推進劑在整個過程中受壓強和結構振蕩而產生的變形對發(fā)動機內流場的影響可以忽略。在每個求解時間步,每個計算模塊均執(zhí)行一遍,流固耦合數據交換模塊在每個求解時間步提取必要的信息用以流場計算和瞬態(tài)動力學模塊計算時邊界條件的設置。

    流固耦合數據交換模塊主要負責對流場計算輸出的壓強和結構計算輸出的加速度進行線性插值,這一點對于本文非共節(jié)點的兩套網格之間的數據交換來說是十分必要的。假設流體網格和結構網格在發(fā)動機軸向方向上兩節(jié)點之間的間距分別為ΔxFI和ΔxSI。設流體網格的節(jié)點為i-1、i、i+1,結構網格的節(jié)點為j-1、j、j+1,如圖2所示,顯示了流固耦合交界面網格示意圖。則根據流場計算模塊得到的燃面表面的壓強p可得到結構網格中節(jié)點j的壓強pj計算公式如式(5)所示。同樣,對于結構網格輸出的加速度值也是采用和上述相類似的方法進行插值,計算公式如式(6)所示。

    圖2 流固耦合交界面網格示意圖

    (5)

    (6)

    2 數值計算與結果分析

    不同厚度的推進劑(即固體火箭發(fā)動機處于不同的工作時刻)對結構振動的傳遞有不同的特性,所以研究不同工作時刻下發(fā)動機燃燒室受激勵作用下壓力振蕩特性是十分必要的。本文選取了翼柱和環(huán)向開槽兩種燃燒室的三個典型工作時刻的結構和流場模型進行計算,從而對翼柱型和環(huán)向開槽型裝藥燃燒室在結構激勵作用下的壓力振蕩特性進行分析研究。

    2.1 計算模型及邊界條件

    為了后文表述方便,將發(fā)動機的工作時間做無量綱化處理:

    Tn=ta/td

    (7)

    式中Tn為無量綱工作時間;ta為發(fā)動機實際工作時刻,s;td為發(fā)動機總工作時長,s。

    本文中選取了翼柱型和環(huán)向開槽型裝藥燃燒室三個典型無量綱工作時刻Tn=0.4、0.6、0.8的內流場和結構模型進行計算,其中內流場的幾何模型和邊界條件如圖3所示。可以看出,從Tn=0.4時刻開始,翼柱裝藥的尾部翼槽特征已經消失,燃燒室尾部變?yōu)橥粩U結構,燃燒室?guī)缀螢樾D軸對稱幾何,因此可以將燃燒室的三維結構簡化為二維軸對稱模型,而環(huán)向開槽裝藥燃燒室本身就是旋轉軸對稱幾何,故同樣可將三維燃燒室簡化為二維進行計算。

    圖3 燃燒室內流場幾何與邊界條件示意圖(m)

    圖4為翼柱和環(huán)向開槽裝藥發(fā)動機的結構及部分計算邊界條件二維示意圖。圖中綠色部分為發(fā)動機裝藥,發(fā)動機殼體的厚度為5 mm,為圖中的灰色部分。圖中紅色線表示燃燒室的內壓載荷加載邊界,為燃燒室的內表面。在每個計算時間步輸出流場計算中得到的燃燒室內表面的壓力,用于瞬態(tài)結構動力學計算。對發(fā)動機頭部的軸向自由度進行了約束,這一約束同其在實際的工作過程的力學邊界條件相同。殼體和推進劑的材料力學參數見表1,燃燒室內流場計算中燃氣的物性參數見表2。時間步長為1×10-5s,燃燒室中注入的燃氣總溫為3158 K。

    圖4 翼柱和環(huán)向開槽裝藥發(fā)動機二維結構和邊界條件

    表1 不同結構組件材料屬性

    表2 燃氣物性參數

    流場計算中,壁面邊界為無滑移的絕熱壁面,這是因為固體推進劑的導熱系數較低,加之固體火箭發(fā)動機工作時間較短,在整個工作過程中可以將燃燒室內表面視作絕熱壁面。噴管出口為壓力出口,只要出口處壓力與燃燒室內壓力滿足臨界條件,出口處的流動就能夠達到超音速流,出口的壓力并不會對燃燒室內的壓力產生任何影響,故出口壓力取10 kPa。燃燒室的入口為質量流率入口,單位面積上的質量注入量的計算式為

    根據前文中的燃速計算模型對受推進劑表面結構振動加速度作用下的燃速進行計算,從而實現結構激勵對燃燒室內壓力振蕩的影響,具體計算流程如圖1所示。本文中流場計算使用的是k-εRNG湍流模型,結構計算中使用的是ANSYS APDL中的二維單元SOLID273,并設置單元KEYOPT(2)=9,即環(huán)向有9個節(jié)點平面。

    本文中施加的結構激勵為脈沖作用的加速度載荷。加速度載荷大小為80g,0.01 s后停止施加。各個時刻的流場先使用穩(wěn)態(tài)求解器進行計算以快速得到穩(wěn)定的流場,穩(wěn)態(tài)計算的收斂判斷條件為:(1)壓力監(jiān)測點的壓力波動小于0.1%;(2)計算域中入口和出口的質量流量差值波動量小于0.1%。在得到穩(wěn)態(tài)計算收斂的流場后進行瞬態(tài)計算,計算0.1 s后得到用于結構激勵作用下燃燒室壓力振蕩特性計算的初始流場。結構計算的時間步長同流場計算中保持一致。為了監(jiān)測燃燒室中的壓力振蕩特性,從頭部到噴管入口處等距離取31個監(jiān)測點,所有監(jiān)測點在一條直線上,離旋轉對稱軸的徑向距離為0.017 m。為了監(jiān)測發(fā)動機的結構振動特性,在燃燒室殼體外表面上從頭到尾等間距取20個監(jiān)測點。最終,將各測點的數據進行FFT,得到對應的頻譜數據進行分析。

    2.2 流場計算部分網格敏感性驗證

    為了選擇合適的網格尺寸進行計算,必須進行網格尺寸敏感性的驗證。使用翼柱裝藥發(fā)動機在無量綱時間Tn=0.4的燃燒室流場模型,如圖3所示。計算了三種不同尺寸的網格(表3)在經歷脈沖壓力后,燃燒室內壓力振蕩衰減數據。并將計算結果進行處理,通過對比振蕩衰減壓力的相位、幅值衰減特性和幅值頻域特性等三方面,來驗證計算網格的敏感性。

    表3 網格驗證工況列表

    計算中使用的燃氣參數如表2所示,計算域入口的質量流率為13 kg/(s·m2),總溫T0=3158 K,噴管出口的邊界條件為壓力出口,模型的最下方為旋轉軸對稱邊界條件,模型的其他邊界均為無滑移絕熱壁面邊界條件,時間步長為1×10-5s。計算10 000個時間步后,對燃燒室頭部位置施加壓力脈沖。將燃燒室頭部壁面的無滑移壁面邊界條件改為壓力入口,phead=20 MPa。計算50個時間步長之后,停止施加壓力脈沖,重新將燃燒室頭部壁面改回無滑移壁面的邊界條件繼續(xù)求解。

    燃燒室在受到脈沖作用之后,計算中設立的兩個監(jiān)測點的壓力振蕩情況如圖5所示。可以看出,三種不同尺寸的網格對壓力脈沖衰減特性的反應能保證較好的一致性。三種網格對于point-1和point-2的壓力振蕩的相位相差180°左右,這是由這兩個點的軸向位置決定的。且由于脈沖從燃燒室頭部到噴管喉部的傳遞過程中,流體的流動會受到壁面以及流體自身的阻尼作用,在振蕩幅值上會有所衰減,所以point-1的壓強振蕩幅值會比point-2高。三種網格反應出的上述脈沖后燃燒室流力振蕩的特征均符合物理規(guī)律,但在振蕩壓力衰減率以及壓力振幅頻譜圖上有所差異。

    (a)M1 (b)M2 (c)M3

    圖6所示為三種網格point-1和point-2的壓力振蕩曲線的上包絡線??梢钥闯?在脈沖施加后的壓力衰減初期,兩個監(jiān)測點的壓力振幅衰減率較大,這是由于燃燒室內的壓力振蕩為非線性過程,其流場中存在的突變幾何結構的非線性阻尼與壓力振蕩幅值相關,壓力振幅越大阻尼也就越大[23]。三種工況的壓力振幅衰減曲線均具有較好的一致性,但在衰減的初始階段,可以看出不同網格之間的衰減率相同,但其幅值有所差異,隨著網格數量的增多,幅值有所增大,這是由于粗網格計算得到的湍流耗散率相對于實際情況會偏大,因此提高網格的密度,能夠更精確地計算出壓力衰減過程中的振幅振蕩衰減特性。

    (a)point-1 (b)point-2

    為了對比三種尺寸的網格在于壓力振蕩頻域上的計算表現,對兩個監(jiān)測點的振蕩壓力做快速傅里葉變換,得到壓強的均方根(RMS)振幅頻譜圖,如圖7所示。從結果可知,三種尺寸的網格在壓力振蕩峰值頻率上具有良好的一致性,前三階峰值頻率分別為172.5、365、586.7 Hz。這與燃燒室前三階聲頻相接近,第一階壓力振蕩幅值最大,且尖峰表現地最為陡峭,這表明這一階的壓力振蕩能量是最集中的??梢?這三種不同尺度的網格均能在壓力振蕩頻域上有良好表現,計算結果能很好地符合實際物理規(guī)律。

    (a)point-1 (b)point-2

    綜合三種網格的計算結果在燃燒室受脈沖作用后在振蕩衰減壓力的相位、幅值衰減特性和幅值頻域特性上的表現,三種尺度的網格均能得到符合物理規(guī)律的有效結果。在壓力衰減規(guī)律上表現較為一致;在壓力振蕩幅值的計算上,M1~M3依次表現為幅值增大,但三者之間的差別不是很大。網格敏感性的驗證在于探究網格數量對于計算結果的影響,由于網格數量越大,計算量也就越大,因此不斷提高網格密度是不現實的。故在綜合考慮計算量和計算結果精度的情況下,本文后續(xù)的計算中將使用M3工況網格的尺度用于計算,即最大網格尺寸為3 mm。

    2.3 結構計算部分網格敏感性驗證

    為了驗證網格尺寸對結構瞬態(tài)響應的影響,需要開展網格無關性驗證。用于驗證的模型為翼柱型裝藥發(fā)動機在Tn=0.4的模型,對發(fā)動機的頭部進行軸向位移約束,在計算初始時刻在發(fā)動機的內表面施加階躍壓強,大小為10 MPa。本文對三種不同尺寸的網格的計算結果進行了對比。三種網格分別命名為粗網格、中等網格、細網格,具體網格尺寸如表4所示。

    表4 網格無關性驗證各工況參數

    本驗證中取發(fā)動機頭部的推進劑內表面作為監(jiān)測位置,選取從加載階躍壓強載荷(即計算初時時刻)開始到0.4×10-2s這一段時間內監(jiān)測點位移隨時間變化的數據作為對比,計算結果如圖8所示。從圖8中可以看出,三種網格的計算結果具有很好的重合性,這說明此時網格已經收斂。因此,考慮到計算量和結果的精確性,本文結構瞬態(tài)響應計算中的網格尺寸將和中等網格保持相一致。

    圖8 發(fā)動機頭部推進劑內表面徑向位移

    2.4 結構邊界受擾動下燃燒室振蕩壓強計算模型驗證

    為了驗證本文計算模型的有效性,設計了相關實驗進行驗證,根據固體火箭發(fā)動機的實際情況,將燃燒室簡化為圓柱管腔,為模擬真實的固體火箭發(fā)動機,從外到內分別布置殼體(結構鋼)、類橡膠材料(用于模擬推進劑),最終搭建出實驗裝置。

    固體火箭發(fā)動機燃燒室實際工作過程中為高溫高壓的環(huán)境,對于實施本文的實驗是不可行的,學術界對燃燒室聲壓振蕩的研究都是建立在冷流的基礎上,而對發(fā)動機結構振動的研究也是以未工作時的發(fā)動機為對象進行的。因此,本實驗對固體火箭發(fā)動機進行了一定的簡化,為實驗的開展進行了一定的假設。首先,實驗中假設燃燒室內的流動是靜止的。其次,本實驗將發(fā)動機燃燒室簡化為兩端封堵的圓柱管腔進行,且在一端使用揚聲器進行激勵,用以模擬燃燒室中自發(fā)產生的壓力波動。實驗裝置如圖9所示。選擇10 mm厚度橡膠的定頻聲激勵實驗數據同計算結果進行對比驗證。如圖10為10 mm圓柱管腔的計算域示意圖,整個計算模型為旋轉軸對稱。如圖11所示,左側一小塊紅色區(qū)域的網格為激勵區(qū),動量源項的變化頻率保持和實驗中揚聲器激勵頻率一致,為102 Hz。計算中流體的材料為空氣,其中密度由理想氣體狀態(tài)方程計算得到。除了旋轉對稱邊界條件以外,其他邊界均為無滑移絕熱壁面。離散格式為二階迎風,時間步長為1×10-5s。首先使用瞬態(tài)求解器進行流場部分的計算,直到管腔中的壓強出現周期性的振蕩,然后啟動瞬態(tài)結構動力學計算,在每個時間步將流場中管腔內表面的壓力分布輸出作為結構計算下一個時間步的載荷輸入,如此反復迭代計算1 s后停止計算,并對得到的數據進行FFT變化,得到振動頻域上的信息,同實驗值進行對比。

    圖9 實驗系統(tǒng)圖

    圖10 驗證計算模型示意圖

    圖11 圓柱管腔內流場二維旋轉軸對稱網格

    圖12為圓柱管腔外殼測點3和測點5處的實驗和數值計算結果的結構振動加速度振幅頻譜圖。從圖中可以看出,數值計算結果和實驗結果在102、203、306 Hz處均有一致的峰值存在,數值計算值的振幅相對于實驗值來說要高一些,同時實驗中高階頻率的峰值更不明顯,這可能是由于實驗中存在的系統(tǒng)阻尼以及數值計算中并不能保證系統(tǒng)邊界條件同實驗中完全一致造成的。但從整體結果上來看,數值計算的結果同實驗結果吻合性較高,這說明本文提出的用于二維旋轉軸對稱模型的弱流固耦合計算方法是合理且有效的。

    (a)point-3 (b)point-5

    2.5 計算結果及討論

    為了便于對固體火箭發(fā)動機燃燒室在結構邊界受擾動下的壓力振蕩特性開展分析,本文對翼柱和環(huán)向開槽裝藥的整個工作過程中的結構固有頻率、聲腔頻率進行了計算。由于本文的研究目的之一就是燃燒室壓力振蕩頻率在整個工作過程中是否會同結構固有頻率相接近,從而導致兩種振動相互加強。而燃燒室內流場計算中壓力振蕩的各階頻率與燃燒室聲學有限元計算得到的各階聲模態(tài)頻率相接近。因此,對于結構固有頻率的關注重點就是在發(fā)動機燃燒室一階、二階聲頻附近的振型。對于燃燒室來說,其聲腔自發(fā)振蕩能量主要集中在燃燒室一階和二階聲頻,尤其是一階聲頻[24]。對于結構振動來說,彎曲振型不會影響燃燒室內聲腔的壓力振蕩,因為其改變的是軸向尺寸;而呼吸模態(tài)能夠改變燃燒徑向和周向的尺寸,因此對推進劑燃速起到影響,從而能夠影響到燃燒室的壓力波動。

    綜上所述,本部分主要關心在燃燒室一階聲頻附近的呼吸模態(tài),故這里給出的是發(fā)動機燃燒室不同工作時刻的呼吸模態(tài)頻率和一、二階聲腔頻率的變化情況。翼柱裝藥、環(huán)向開槽裝藥發(fā)動機的呼吸模態(tài)頻率和燃燒室前兩階聲頻隨工作時間的變化曲線如圖13所示。

    (a)Finocyl configuration motor (b)Axil configuration motor

    由圖13(a)可以看出,隨著工作時間的增長,發(fā)動機的裝藥量不斷減小,其呼吸模態(tài)出現頻率也在不斷地升高,可見發(fā)動機結構的質量對其振型頻率起到了決定性作用。呼吸模態(tài)的振動頻率的變化范圍為工作初期的140 Hz到工作末期的320 Hz左右。在無量綱工作時刻Tn=0.4,該呼吸模態(tài)頻率與一階聲模態(tài)頻率相交,整個工作過程中這一振型的頻率曲線穿越了一階聲頻。由圖13(b)可知,呼吸模態(tài)頻率先減小后增大,在工作的前中期,其頻率隨裝藥質量減少而減小的現象更明顯,從初始時刻的255 Hz到工作末期的259 Hz,整個工作過程中沒有同發(fā)動機燃燒室的一階和二階聲頻相交。從圖13中可以看出,隨著工作時間的增長,呼吸模態(tài)的頻率不再隨發(fā)動機裝藥質量的減少而增大,這一點同翼柱裝藥發(fā)動機有著明顯的差別。

    翼柱裝藥燃燒室Tn=0.4、Tn=0.6、Tn=0.8時刻的不同軸向位置測點的壓強頻譜圖如圖14所示。從圖14(a)中可以明顯看出,在受結構脈沖加速度激勵后,燃燒室中產生了持續(xù)的壓強振蕩現象,且以一階聲頻下的振蕩為主,頻率為180.4 Hz;二階聲頻下的振蕩強度相對來說就要小很多,其頻率為360.7 Hz。從圖中壓強振蕩幅值隨軸向位置的變化規(guī)律可以看出,燃燒室內的壓強波動表現為明顯的駐波特性,同燃燒室聲模態(tài)振型相符合。具體表現為其一階振型兩端為波峰,中間為波節(jié),故燃燒室兩端的振幅高,中間波節(jié)處的振幅基本為0;其二階振型中存在3個波峰,2個波節(jié)。由于燃燒室內的壓強振蕩會在噴管的喉部發(fā)生反射,因此對于壓強振蕩來說實際波動傳播的長度還需要加上噴管的擴張段,這也就導致圖中一階聲頻下的振型的波節(jié)位置位于燃燒室中間略靠后處。同時由于燃燒室靠近噴管附近的流體流動速度比燃燒室前端要大,導致此處非定常特征更加突出,這就使得尾部的壓強振幅要低于頭部。此外,翼柱裝藥燃燒室的尾部存在突擴的空腔,這一幾何結構對壓強的振蕩具有阻尼作用,同樣降低了尾部的壓強振幅。一階振型中頭部測點的壓強振幅為71.7 kPa,尾部測點的壓強振幅為27.5 kPa,頭部測點壓強振幅是尾部的261%。

    由圖14(b)可見,同Tn=0.4時刻的計算結果相似,Tn=0.6時刻燃燒室中一階聲模態(tài)下的振蕩是最為強烈的,二階聲模態(tài)下的振蕩幅度則要小得多,壓強振幅在軸向位置上分布同燃燒室的聲模態(tài)振型相一致。一階聲模態(tài)的壓強振蕩頻率為182.2 Hz,二階聲模態(tài)的振蕩頻率為364.4 Hz。從圖中可知,燃燒室的壓強振幅較Tn=0.4時刻出現了降低。Tn=0.6時刻燃燒室頭部壓強振幅為49.4 kPa,尾部壓強振幅為28.7 kPa,頭部測點壓強振幅是尾部的172%,頭部壓強的振幅強度同尾部之間的差距進一步減小。由此可見,與Tn=0.4時刻相比,主要是燃燒室頭部的壓強振蕩強度出現了降低。根據圖13的翼柱裝藥發(fā)動機不同工作時刻的呼吸模態(tài)頻率變化規(guī)律可知,在Tn=0.4時刻其呼吸模態(tài)的頻率為179.4 Hz,同燃燒室的一階聲頻相接近,這就使得壓強振蕩同結構振動之間出現了相互耦合加強的現象,兩者的振動得到了相互增強,而燃燒室尾部的壓強振蕩由于流速較大導致的強烈的非定常特性以及燃燒室突擴空腔的聲阻尼作用使得其壓強振幅并沒有出現顯著變化。而Tn=0.6時刻的呼吸模態(tài)頻率為202.7 Hz,遠離了燃燒室的一階聲模態(tài)頻率,這也就使得無法出現耦合振動的現象,壓強的振蕩強度出現了降低,盡管Tn=0.6時刻的燃面面積大于Tn=0.4時刻,即燃燒室內單位時間注入的總擾動能量Tn=0.6時刻更高??梢?燃燒室聲頻同結構的呼吸模態(tài)振動頻率是否相接近是決定燃燒室發(fā)生不穩(wěn)定燃燒現象時壓強振蕩強度的重要因素。

    (a)Tn=0.4

    由圖14(c)可見,同前兩個時刻的計算結果相似,燃燒室中一階聲模態(tài)下的振蕩是最為強烈的,二階聲模態(tài)的振動基本已經不可見,壓強振幅在軸向位置上分布同燃燒室的聲模態(tài)振型相一致。此時,燃燒室一階聲模態(tài)頻率為197 Hz,發(fā)動機的結構呼吸模態(tài)頻率為286 Hz。Tn=0.8時刻燃燒室頭部壓強振幅為41 kPa,尾部壓強振幅為32.9 kPa,頭部測點壓強振幅是尾部的128%,燃燒室頭部后尾部之間壓力振幅之間的差別進一步減小。盡管Tn=0.8時刻燃燒室一階聲頻同結構呼吸模態(tài)頻率相差相較于Tn=0.6時刻要大得多,但是壓強振幅的下降程度并沒有那么大??梢娭灰紵衣曨l同呼吸模態(tài)的頻率一定距離就能夠使得流場的壓力振蕩同結構振動相互作用引發(fā)的耦合振動加強的現象得到較大程度的緩解,從設計角度而言,只需要做有限度的結構調整就能夠獲得較大程度抑制效應,從而削弱燃燒室壓力振蕩同結構耦合振動加強的強度。

    環(huán)向開槽裝藥燃燒室Tn=0.4、Tn=0.6、Tn=0.8時刻不同軸向位置測點的壓強頻譜圖如圖15所示。一階聲頻的壓力振幅最為明顯,其頻率為177.2 Hz,其他階數的壓力振幅在這一時刻并沒有明顯體現,壓強振幅在軸向位置上分布同燃燒室的聲模態(tài)振型相一致。

    (a)Tn=0.4

    由圖15(a)可知,此時的發(fā)動機結構呼吸模態(tài)頻率為212 Hz,同燃燒室的一階聲模態(tài)頻率有較大的差距。因此,壓力振蕩幅值相比于Tn=0.4時刻的翼柱裝藥燃燒室要低得多。該時刻燃燒室頭部壓強振幅為21.9 kPa,尾部壓強振幅為15.8 kPa,頭部測點壓強振幅是尾部的138.6%。

    如圖15(b)所示,一階聲頻的壓力振幅最為明顯,其頻率為172.2 Hz,其他階數聲頻下的壓力振幅在這一時刻并沒有明顯體現,壓力振幅在軸向位置上分布同燃燒室的聲模態(tài)振型相一致。

    如圖15(c)所示,除了一階聲頻下的壓力振蕩,第5階聲頻下的振動也表現得十分明顯,這兩階壓力振蕩的頻率分別為179.8 Hz和913.9 Hz。由此可知對于環(huán)向開槽裝藥燃燒室而言,在工作的后期,其部分振動能量開始向高階聲模態(tài)轉移。一階聲頻下的燃燒室頭部壓強振幅為23.3 kPa,尾部壓強振幅為19.2 kPa,頭部測點的壓強振幅是尾部的121.4%。

    為了對比兩種不同裝藥類型的發(fā)動機在受結構脈沖加速度載荷激勵作用下,燃燒室壓強振蕩特性,將燃燒室各時刻一階聲模態(tài)下頭部和尾部測點的壓強振幅進行對比,如表5所示。表中相對值為環(huán)槽型裝藥的壓強振幅比上翼柱型裝藥。從表5可知,翼柱裝藥燃燒室的壓強振幅始終要高于環(huán)向開槽型裝藥。尤其是頭部測點的壓強振幅,環(huán)向開槽裝藥燃燒室最低只有翼柱型的30.5%。然而,在本文計算的3個時刻,環(huán)向開槽裝藥燃燒室的燃面面積始終大于翼柱型,從線性聲不穩(wěn)定的角度看,燃面面積越大意味著燃燒室的響應增益越大,在受結構激勵時能夠為燃燒室的壓強振蕩和結構的振動提供更多的能量。由此可見,環(huán)向開槽裝藥發(fā)動機的穩(wěn)定性要明顯優(yōu)于翼柱型。這主要有兩個原因:第一,翼柱裝藥燃燒室在Tn=0.4時刻出現了呼吸模態(tài)頻率同燃燒室一階聲頻相接近的情況,造成了燃燒室壓強振蕩同結構振動之間出現了相互耦合振動加強的現象,使得兩者的振動幅度得到了加強,而環(huán)向開槽型裝藥發(fā)動機在整個工作過程中均未出現兩種振動的頻率相接近的情況,從而避免了耦合振動現象的發(fā)生;第二,環(huán)向開槽型裝藥燃燒室的頭部和尾部區(qū)域存在有環(huán)槽,這些環(huán)槽對燃燒室內的壓強振蕩起到了較明顯的阻尼作用,從而使得翼柱裝藥發(fā)動機頭部位置測點的壓強振幅即使在呼吸模態(tài)頻率同一階聲頻不再那么接近的情況下仍要更大一些。而由于兩種燃燒室尾部位置均存在特殊的幾何結構可以對壓強振蕩起到阻滯作用,加之這里的流速遠大于頭部,流動的非定常特性更為明顯。從而使得兩測點的壓強振幅基本相接近。由上述分析可知,對于固體火箭發(fā)動機燃燒室來說,需要避免在整個工作中結構的呼吸模態(tài)頻率同燃燒室聲頻相接近,以及設計推進劑藥型時,盡量使燃燒室的聲腔中可以存在對壓強振蕩起阻尼作用的結構,如本文的環(huán)槽結構。同時也說明,對于這種不斷有能量釋放的燃燒室而言,結構對振動的吸收作用并不是主要的,即推進劑厚度對不穩(wěn)定燃燒的影響并不是因為對振動的吸收能力改變而引起的,更多的是由于結構呼吸模態(tài)的變化導致與聲頻的接近程度產生的。

    表5 燃燒室頭部和尾部測點一階聲頻壓強振幅

    3 結論

    (1)本文以軸向長度相同的翼柱裝藥和環(huán)向開槽裝藥燃燒室為研究對象,使用本文構建的計算方法,研究了兩種燃燒室在三個典型工作時刻受到結構脈沖加速度載荷作用下燃燒室壓力振蕩特性。通過對燃燒室軸向不同點的壓力振幅進行分析可知,發(fā)動機結構呼吸模態(tài)頻率與燃燒室的聲模態(tài)頻率的接近程度以及燃燒室的不同特征結構,將導致在結構激勵作用下燃燒室的壓強振蕩特性表現出明顯的差異。

    (2)在結構脈沖加速度載荷的激勵下,燃燒室中產生的壓強振蕩以一階聲頻為主,發(fā)動機結構振動也表現出相似的振蕩特性??赏浦?當脈沖激勵結束,燃燒室發(fā)生持續(xù)振蕩后,激發(fā)的壓力振蕩帶動了結構一同振動,這一點同實際中觀測到的固體火箭發(fā)動機發(fā)生不穩(wěn)定燃燒現象時彈體的振動特性相類似。因此,可以從避免燃燒室聲頻同結構呼吸模態(tài)頻率相接近方向入手,抑制燃燒不穩(wěn)定現象。此外,燃燒室的特殊幾何結構可以對壓力振蕩起到明顯的阻尼作用。

    (3)本文中的環(huán)向開槽裝藥燃燒室的穩(wěn)定性要明顯優(yōu)于翼柱裝藥,在受結構脈沖激勵作用下產生的壓力振幅要明顯低于翼柱型。主要原因有兩個:第一,環(huán)向開槽型裝藥燃燒室在整個工作過程中其呼吸模態(tài)頻率一直遠離燃燒室的一階聲頻,而翼柱型則在Tn=0.4時刻出現了兩種頻率相交的情況,這就使得壓強振蕩和結構振動出現耦合振動,從而使振幅較大;第二,環(huán)向開槽裝藥的頭部和尾部的環(huán)槽結構在整個工作過程中一直存在,對壓力振蕩起到了阻尼作用。因此,環(huán)向開槽結構替代翼柱結構,可以明顯提高發(fā)動機的穩(wěn)定性。

    猜你喜歡
    燃燒室裝藥振幅
    燃燒室形狀對國六柴油機性能的影響
    火炸藥學報(2022年3期)2022-07-04 07:31:00
    孔內爆炸填塞效應的數值模擬與分析
    某發(fā)射裝藥結構改進設計的新思路
    一種熱電偶在燃燒室出口溫度場的測量應用
    電子制作(2019年19期)2019-11-23 08:41:54
    十大漲跌幅、換手、振幅、資金流向
    十大漲跌幅、換手、振幅、資金流向
    十大漲跌幅、換手、振幅、資金流向
    滬市十大振幅
    深孔不耦合裝藥爆破技術卸壓效果驗證
    日日摸夜夜添夜夜添av毛片| 国产无遮挡羞羞视频在线观看| 国产精品无大码| 国产av一区二区精品久久| 国产伦精品一区二区三区视频9| 亚洲国产日韩一区二区| 乱码一卡2卡4卡精品| 黄色视频在线播放观看不卡| 多毛熟女@视频| 欧美日韩精品成人综合77777| 国产熟女欧美一区二区| 两个人免费观看高清视频 | 丝袜在线中文字幕| 大码成人一级视频| 国产 一区精品| 久久精品国产鲁丝片午夜精品| 婷婷色av中文字幕| 久久久国产一区二区| 爱豆传媒免费全集在线观看| 久久久久久久久久久免费av| 国产亚洲一区二区精品| 国语对白做爰xxxⅹ性视频网站| 亚洲真实伦在线观看| 国产白丝娇喘喷水9色精品| 老女人水多毛片| 久久精品久久久久久久性| 97超视频在线观看视频| 菩萨蛮人人尽说江南好唐韦庄| www.色视频.com| 国产精品福利在线免费观看| 美女主播在线视频| 欧美激情国产日韩精品一区| 十八禁网站网址无遮挡 | 久久国产乱子免费精品| 国产一区二区三区av在线| 看非洲黑人一级黄片| 高清黄色对白视频在线免费看 | 黄色毛片三级朝国网站 | 亚洲图色成人| 国产中年淑女户外野战色| 欧美亚洲 丝袜 人妻 在线| 久久久久国产网址| 亚洲va在线va天堂va国产| 日日啪夜夜撸| 久久婷婷青草| av在线app专区| 精品一区在线观看国产| 亚洲一区二区三区欧美精品| 中文字幕人妻丝袜制服| 在线观看www视频免费| 天堂8中文在线网| 夫妻午夜视频| 熟女电影av网| 日本黄大片高清| 国产女主播在线喷水免费视频网站| 人人妻人人看人人澡| 新久久久久国产一级毛片| 欧美精品亚洲一区二区| 久久久久网色| 免费在线观看成人毛片| 国产乱来视频区| 九九爱精品视频在线观看| 国产成人精品一,二区| 欧美97在线视频| 亚洲av在线观看美女高潮| 国产精品蜜桃在线观看| 欧美 亚洲 国产 日韩一| 久久综合国产亚洲精品| 日日摸夜夜添夜夜爱| 国产免费视频播放在线视频| 九色成人免费人妻av| 久久99蜜桃精品久久| av国产久精品久网站免费入址| 亚洲va在线va天堂va国产| av卡一久久| 欧美区成人在线视频| 国产成人一区二区在线| 伊人久久精品亚洲午夜| 精品国产一区二区三区久久久樱花| 人妻夜夜爽99麻豆av| 国产精品福利在线免费观看| 国产精品人妻久久久影院| 国产高清国产精品国产三级| 国产视频首页在线观看| 啦啦啦视频在线资源免费观看| 亚洲精品日本国产第一区| 久久国产亚洲av麻豆专区| 国模一区二区三区四区视频| 另类精品久久| av在线老鸭窝| 人妻制服诱惑在线中文字幕| 纵有疾风起免费观看全集完整版| 久久99蜜桃精品久久| 草草在线视频免费看| 一边亲一边摸免费视频| 亚洲色图综合在线观看| 丝袜脚勾引网站| 久久久久久久久久久丰满| 久久婷婷青草| 国产在线免费精品| 中文天堂在线官网| 日韩制服骚丝袜av| 91在线精品国自产拍蜜月| 国产无遮挡羞羞视频在线观看| 精品一品国产午夜福利视频| 国产av一区二区精品久久| 黄色一级大片看看| 亚洲激情五月婷婷啪啪| 菩萨蛮人人尽说江南好唐韦庄| 亚洲欧洲日产国产| av一本久久久久| 亚洲综合精品二区| 婷婷色av中文字幕| 国产淫语在线视频| 国产男女超爽视频在线观看| 久久99精品国语久久久| 久久久国产一区二区| 能在线免费看毛片的网站| 中文字幕制服av| 日韩视频在线欧美| 成人美女网站在线观看视频| 国语对白做爰xxxⅹ性视频网站| 一区二区三区四区激情视频| 亚洲婷婷狠狠爱综合网| 国产精品久久久久久精品电影小说| 色婷婷久久久亚洲欧美| 色婷婷久久久亚洲欧美| 能在线免费看毛片的网站| 少妇猛男粗大的猛烈进出视频| 欧美日韩视频高清一区二区三区二| 日本av手机在线免费观看| 卡戴珊不雅视频在线播放| 精品人妻熟女av久视频| 欧美3d第一页| 老女人水多毛片| 久久这里有精品视频免费| 欧美 日韩 精品 国产| 国内精品宾馆在线| 午夜激情久久久久久久| 亚洲av成人精品一二三区| 中文资源天堂在线| 久久久久国产精品人妻一区二区| 亚洲高清免费不卡视频| 91在线精品国自产拍蜜月| 一级a做视频免费观看| 久久精品久久久久久噜噜老黄| 十八禁高潮呻吟视频 | 日本黄大片高清| 春色校园在线视频观看| 久久久久人妻精品一区果冻| 国产精品偷伦视频观看了| 最近中文字幕2019免费版| 国产免费又黄又爽又色| 一级毛片我不卡| 久久久精品94久久精品| 国模一区二区三区四区视频| 嘟嘟电影网在线观看| 午夜影院在线不卡| 国产爽快片一区二区三区| 日韩电影二区| 在线观看人妻少妇| 人妻人人澡人人爽人人| 大片免费播放器 马上看| 国产成人精品久久久久久| 中文字幕人妻丝袜制服| 制服丝袜香蕉在线| 欧美最新免费一区二区三区| 大陆偷拍与自拍| 五月开心婷婷网| 国产白丝娇喘喷水9色精品| a级毛片在线看网站| freevideosex欧美| 国产精品国产三级专区第一集| av网站免费在线观看视频| 国产女主播在线喷水免费视频网站| 亚洲av中文av极速乱| freevideosex欧美| 日产精品乱码卡一卡2卡三| 一区二区三区精品91| 国产亚洲5aaaaa淫片| 日本av免费视频播放| 只有这里有精品99| 2022亚洲国产成人精品| 91久久精品电影网| 一区在线观看完整版| 男女边吃奶边做爰视频| 一本大道久久a久久精品| 嘟嘟电影网在线观看| 午夜免费观看性视频| 午夜日本视频在线| 免费不卡的大黄色大毛片视频在线观看| 一级片'在线观看视频| 男女边吃奶边做爰视频| 成人影院久久| 91在线精品国自产拍蜜月| 性色avwww在线观看| 欧美xxxx性猛交bbbb| 国产精品伦人一区二区| 国产色爽女视频免费观看| 大又大粗又爽又黄少妇毛片口| 久久av网站| 亚洲电影在线观看av| 日本黄大片高清| 自拍偷自拍亚洲精品老妇| 亚洲av男天堂| 国产免费视频播放在线视频| 丝袜在线中文字幕| 国产伦理片在线播放av一区| 亚洲精品aⅴ在线观看| 菩萨蛮人人尽说江南好唐韦庄| 自拍欧美九色日韩亚洲蝌蚪91 | 久久久久网色| 天堂俺去俺来也www色官网| 三上悠亚av全集在线观看 | 日韩av不卡免费在线播放| av在线观看视频网站免费| 我的女老师完整版在线观看| 欧美精品国产亚洲| 色哟哟·www| kizo精华| 国产精品蜜桃在线观看| 国产 精品1| 好男人视频免费观看在线| 狂野欧美激情性xxxx在线观看| 国产伦精品一区二区三区视频9| 亚洲av国产av综合av卡| 熟女电影av网| 黄色日韩在线| 久久精品国产亚洲av涩爱| 少妇的逼水好多| 国产精品99久久久久久久久| 女人久久www免费人成看片| 夜夜看夜夜爽夜夜摸| 中文字幕人妻丝袜制服| 日本vs欧美在线观看视频 | 久久久久国产网址| 制服丝袜香蕉在线| 中文精品一卡2卡3卡4更新| 国产成人a∨麻豆精品| 国产亚洲一区二区精品| 亚洲精华国产精华液的使用体验| 男女啪啪激烈高潮av片| 最新中文字幕久久久久| 欧美日韩精品成人综合77777| 久久人人爽人人爽人人片va| 亚洲一级一片aⅴ在线观看| 亚洲伊人久久精品综合| 亚洲电影在线观看av| 99热这里只有是精品50| 国产黄色视频一区二区在线观看| 啦啦啦视频在线资源免费观看| av播播在线观看一区| 伦精品一区二区三区| 免费黄色在线免费观看| 成人二区视频| 欧美一级a爱片免费观看看| 伊人久久精品亚洲午夜| 成年人免费黄色播放视频 | 国产伦在线观看视频一区| 少妇的逼水好多| 内射极品少妇av片p| 欧美激情国产日韩精品一区| 亚洲精品乱码久久久v下载方式| 青春草视频在线免费观看| 国产欧美亚洲国产| 深夜a级毛片| 精品久久久久久久久av| 亚洲va在线va天堂va国产| 亚洲av免费高清在线观看| 少妇人妻久久综合中文| 国产男人的电影天堂91| 国产精品一区www在线观看| a级一级毛片免费在线观看| 免费少妇av软件| 不卡视频在线观看欧美| 寂寞人妻少妇视频99o| 日韩伦理黄色片| 免费观看无遮挡的男女| 国产黄色视频一区二区在线观看| 久久婷婷青草| 人人妻人人澡人人爽人人夜夜| 在线观看免费日韩欧美大片 | tube8黄色片| 九九在线视频观看精品| 美女主播在线视频| 在线观看免费高清a一片| 人妻 亚洲 视频| 99久久精品热视频| 夜夜骑夜夜射夜夜干| 亚洲av成人精品一区久久| av福利片在线观看| 国产亚洲5aaaaa淫片| 中文精品一卡2卡3卡4更新| 啦啦啦视频在线资源免费观看| 欧美 亚洲 国产 日韩一| 少妇精品久久久久久久| 波野结衣二区三区在线| 成人毛片60女人毛片免费| av免费在线看不卡| av国产精品久久久久影院| 少妇猛男粗大的猛烈进出视频| av黄色大香蕉| 精品酒店卫生间| 99热国产这里只有精品6| 国产精品女同一区二区软件| 99久久精品国产国产毛片| 高清在线视频一区二区三区| 99国产精品免费福利视频| 99re6热这里在线精品视频| 五月伊人婷婷丁香| 看十八女毛片水多多多| 免费观看av网站的网址| 国产成人91sexporn| 午夜福利视频精品| 嘟嘟电影网在线观看| 另类亚洲欧美激情| 日韩欧美 国产精品| 在线观看三级黄色| 一区二区三区乱码不卡18| 高清不卡的av网站| 视频中文字幕在线观看| 天堂俺去俺来也www色官网| 久热这里只有精品99| 男人舔奶头视频| 日韩中文字幕视频在线看片| 人人妻人人澡人人看| 亚洲欧洲日产国产| 国产亚洲最大av| 深夜a级毛片| 少妇被粗大猛烈的视频| 精品熟女少妇av免费看| 日本91视频免费播放| 丁香六月天网| 九九爱精品视频在线观看| 免费人成在线观看视频色| 中国美白少妇内射xxxbb| 大话2 男鬼变身卡| 2022亚洲国产成人精品| 国产极品粉嫩免费观看在线 | 在线精品无人区一区二区三| 一级,二级,三级黄色视频| 少妇被粗大猛烈的视频| 激情五月婷婷亚洲| 最后的刺客免费高清国语| 国产老妇伦熟女老妇高清| 老女人水多毛片| 欧美日本中文国产一区发布| 男女免费视频国产| 亚洲av男天堂| 国精品久久久久久国模美| 欧美日韩av久久| 国产亚洲av片在线观看秒播厂| 国产成人freesex在线| 亚洲精品日韩在线中文字幕| 久久精品久久久久久噜噜老黄| 一级二级三级毛片免费看| 妹子高潮喷水视频| 亚洲性久久影院| 亚洲国产av新网站| 99久久精品国产国产毛片| 国产深夜福利视频在线观看| 中文资源天堂在线| 一区二区三区精品91| 亚洲欧美日韩卡通动漫| 人妻系列 视频| 高清欧美精品videossex| 日韩三级伦理在线观看| 亚洲婷婷狠狠爱综合网| 18禁动态无遮挡网站| 新久久久久国产一级毛片| 国内少妇人妻偷人精品xxx网站| 十分钟在线观看高清视频www | 少妇裸体淫交视频免费看高清| 嘟嘟电影网在线观看| 欧美高清成人免费视频www| 国内少妇人妻偷人精品xxx网站| 亚洲精品色激情综合| 午夜91福利影院| 大又大粗又爽又黄少妇毛片口| 在线天堂最新版资源| 免费播放大片免费观看视频在线观看| 高清不卡的av网站| 在线观看www视频免费| 97超碰精品成人国产| 黄片无遮挡物在线观看| 久久久精品94久久精品| 99热全是精品| 久热这里只有精品99| 你懂的网址亚洲精品在线观看| 91午夜精品亚洲一区二区三区| 亚洲精品乱久久久久久| 麻豆成人午夜福利视频| 观看免费一级毛片| av天堂中文字幕网| 国内揄拍国产精品人妻在线| 在线观看美女被高潮喷水网站| 国产高清三级在线| 久久亚洲国产成人精品v| 交换朋友夫妻互换小说| 九色成人免费人妻av| 国产淫片久久久久久久久| 国产日韩欧美亚洲二区| 日本黄大片高清| 大片电影免费在线观看免费| 亚洲欧美日韩另类电影网站| 久久精品国产亚洲av涩爱| av天堂中文字幕网| 日韩不卡一区二区三区视频在线| 99久久精品国产国产毛片| 亚洲精品乱码久久久v下载方式| 国产在线免费精品| av福利片在线| 中文字幕av电影在线播放| 久久狼人影院| 欧美精品人与动牲交sv欧美| 中文精品一卡2卡3卡4更新| 国产美女午夜福利| 国产淫语在线视频| 国产精品一区二区三区四区免费观看| 美女主播在线视频| 午夜91福利影院| av专区在线播放| 久久99一区二区三区| 国产有黄有色有爽视频| 国产精品人妻久久久影院| 国产在视频线精品| 一级毛片 在线播放| 婷婷色综合www| 欧美高清成人免费视频www| 婷婷色麻豆天堂久久| 日韩强制内射视频| 少妇人妻精品综合一区二区| 欧美激情国产日韩精品一区| xxx大片免费视频| 亚洲人成网站在线播| 日本黄大片高清| 精品久久久久久久久av| 国产黄色视频一区二区在线观看| 午夜福利,免费看| 免费人妻精品一区二区三区视频| 超碰97精品在线观看| 看免费成人av毛片| 精品一区二区免费观看| 丰满少妇做爰视频| 久久精品国产鲁丝片午夜精品| 黄色日韩在线| 婷婷色麻豆天堂久久| 国产精品秋霞免费鲁丝片| 青春草国产在线视频| 欧美+日韩+精品| 99热国产这里只有精品6| 卡戴珊不雅视频在线播放| h视频一区二区三区| av在线观看视频网站免费| 国产免费一区二区三区四区乱码| 国产免费又黄又爽又色| 边亲边吃奶的免费视频| 国产男女内射视频| 久久97久久精品| 韩国高清视频一区二区三区| 岛国毛片在线播放| 久久精品国产自在天天线| 乱系列少妇在线播放| 久久久久久久久久成人| 国产精品秋霞免费鲁丝片| 亚洲va在线va天堂va国产| 亚洲av成人精品一区久久| 少妇被粗大猛烈的视频| 国产黄片视频在线免费观看| 亚洲国产成人一精品久久久| 最近中文字幕2019免费版| 99九九线精品视频在线观看视频| 色视频www国产| 麻豆乱淫一区二区| 国产又色又爽无遮挡免| 人妻夜夜爽99麻豆av| 亚洲,欧美,日韩| 观看av在线不卡| 亚洲,一卡二卡三卡| 中文字幕免费在线视频6| 欧美激情国产日韩精品一区| 99久久中文字幕三级久久日本| 一本一本综合久久| 午夜视频国产福利| 日本色播在线视频| 一级毛片我不卡| 久久久久国产网址| 丰满人妻一区二区三区视频av| 麻豆成人午夜福利视频| 欧美激情国产日韩精品一区| 中国国产av一级| 国产一区有黄有色的免费视频| 国产精品成人在线| 国产av国产精品国产| 亚洲成人手机| 男人添女人高潮全过程视频| 亚洲美女黄色视频免费看| 性色avwww在线观看| 国产色爽女视频免费观看| 亚洲av欧美aⅴ国产| 纵有疾风起免费观看全集完整版| 2018国产大陆天天弄谢| 在线观看www视频免费| 全区人妻精品视频| 有码 亚洲区| 欧美日韩亚洲高清精品| 制服丝袜香蕉在线| 久久久久久久精品精品| 亚洲人成网站在线观看播放| 精品久久久久久久久亚洲| 美女cb高潮喷水在线观看| 丁香六月天网| 国产一区二区在线观看日韩| 大香蕉97超碰在线| 伊人久久精品亚洲午夜| 午夜久久久在线观看| 一级av片app| 日本黄大片高清| 国产精品欧美亚洲77777| 亚洲三级黄色毛片| 亚洲美女黄色视频免费看| 大片免费播放器 马上看| 亚洲人成网站在线观看播放| 久久久久视频综合| 男人狂女人下面高潮的视频| 亚洲欧美一区二区三区黑人 | 人人妻人人爽人人添夜夜欢视频 | 国产欧美日韩精品一区二区| 久久99蜜桃精品久久| 人妻人人澡人人爽人人| 午夜免费男女啪啪视频观看| 国产一区二区在线观看日韩| 最近2019中文字幕mv第一页| 日韩熟女老妇一区二区性免费视频| 少妇的逼好多水| 国产深夜福利视频在线观看| 国产午夜精品久久久久久一区二区三区| 精品卡一卡二卡四卡免费| 亚洲精华国产精华液的使用体验| 久久久久久久久大av| 蜜桃久久精品国产亚洲av| 久久久亚洲精品成人影院| 我要看日韩黄色一级片| 少妇熟女欧美另类| 人人妻人人爽人人添夜夜欢视频 | 国产 一区精品| 男男h啪啪无遮挡| 亚洲第一av免费看| 日本vs欧美在线观看视频 | 国产精品99久久久久久久久| 最新中文字幕久久久久| 热re99久久国产66热| 乱系列少妇在线播放| 色94色欧美一区二区| 成人影院久久| 欧美高清成人免费视频www| av在线观看视频网站免费| 成人漫画全彩无遮挡| 精品一区二区三卡| 欧美精品人与动牲交sv欧美| 多毛熟女@视频| 建设人人有责人人尽责人人享有的| 国产亚洲5aaaaa淫片| 精品少妇内射三级| 尾随美女入室| 久久影院123| 日韩视频在线欧美| 视频区图区小说| 久久久久久久大尺度免费视频| 国产一区二区三区av在线| 国产成人aa在线观看| 美女大奶头黄色视频| 午夜91福利影院| 波野结衣二区三区在线| 一级二级三级毛片免费看| 国产精品国产av在线观看| 女人精品久久久久毛片| 超碰97精品在线观看| 自拍偷自拍亚洲精品老妇| 大香蕉97超碰在线| 三级经典国产精品| 免费不卡的大黄色大毛片视频在线观看| 天堂俺去俺来也www色官网| 欧美精品一区二区大全| 免费观看无遮挡的男女| 国产在线视频一区二区| √禁漫天堂资源中文www| 日本91视频免费播放| 99久久中文字幕三级久久日本| 热re99久久精品国产66热6| 肉色欧美久久久久久久蜜桃| 国产精品久久久久久精品电影小说| 国产精品一区二区性色av| 一级黄片播放器| 日韩一本色道免费dvd| 国产免费一区二区三区四区乱码| 亚洲成色77777| 观看美女的网站| 免费少妇av软件| 各种免费的搞黄视频| 女的被弄到高潮叫床怎么办| 欧美亚洲 丝袜 人妻 在线| 欧美成人精品欧美一级黄| 亚洲精品自拍成人| 热99国产精品久久久久久7| 我的老师免费观看完整版| 成人无遮挡网站| 国产一区二区三区综合在线观看 | 日韩中文字幕视频在线看片| 午夜影院在线不卡| 日日撸夜夜添| 久久久久久人妻| 日韩制服骚丝袜av| 中文字幕亚洲精品专区| 只有这里有精品99|