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

    渦輪盤篦齒裂紋擴(kuò)展的有限元數(shù)值模擬

    2016-11-04 09:32:31張成成高靖云李孟光
    工程設(shè)計學(xué)報 2016年2期
    關(guān)鍵詞:裂紋有限元模型

    任 遠(yuǎn),張成成,高靖云,李孟光

    (1.中航商用航空發(fā)動機(jī)有限責(zé)任公司上海商用飛機(jī)發(fā)動機(jī)工程技術(shù)研究中心,上?!?00240;2.上海量維信息科技有限公司,上?!?00233)

    渦輪盤篦齒裂紋擴(kuò)展的有限元數(shù)值模擬

    任 遠(yuǎn)1,張成成1,高靖云1,李孟光2

    (1.中航商用航空發(fā)動機(jī)有限責(zé)任公司上海商用飛機(jī)發(fā)動機(jī)工程技術(shù)研究中心,上海200240;2.上海量維信息科技有限公司,上海200233)

    為了詳細(xì)考察篦齒裂紋的擴(kuò)展規(guī)律,對篦齒裂紋從齒尖一直擴(kuò)展到即將完全穿透篦齒環(huán)的過程進(jìn)行了數(shù)值模擬.含篦齒裂紋的渦輪盤有限元模型采用子模型法建立,使用M積分計算裂紋前沿的應(yīng)力強(qiáng)度因子;在確定篦齒裂紋前沿每一節(jié)點(diǎn)處的局部擴(kuò)展方向及距離后,通過樣條曲線擬合出新裂紋前沿,并依靠自適應(yīng)網(wǎng)格劃分實現(xiàn)裂紋區(qū)有限元網(wǎng)格的更新.數(shù)值模擬結(jié)果表明,篦齒根部過渡圓角頂部可以視為裂紋緩慢擴(kuò)展階段與快速擴(kuò)展階段的分界點(diǎn),在此之前篦齒裂紋以穿透型裂紋的形態(tài)以較低的速度進(jìn)行擴(kuò)展,在此之后篦齒裂紋開始向表面裂紋進(jìn)行轉(zhuǎn)化,并且平均擴(kuò)展速度大大增加,分界點(diǎn)前的裂紋擴(kuò)展壽命是之后的數(shù)倍.此外,由數(shù)值模擬結(jié)果還可以發(fā)現(xiàn),增大篦齒根部過渡圓角半徑以及減小相鄰篦齒的間距,均有助于延緩篦齒裂紋的擴(kuò)展.

    篦齒裂紋;渦輪盤;裂紋擴(kuò)展;應(yīng)力強(qiáng)度因子;疲勞壽命

    篦齒封嚴(yán)結(jié)構(gòu)安裝在高低壓動靜轉(zhuǎn)子之間,利用獨(dú)特的幾何形狀形成一個曲折通道,增大流體流動阻力與能量耗散,從而阻隔高低壓腔間的流體泄露[1].由于結(jié)構(gòu)簡單、維護(hù)方便并且封嚴(yán)特性好,它被廣泛應(yīng)用在航空發(fā)動機(jī)內(nèi)流腔室間的氣體封嚴(yán)系統(tǒng)中.

    渦輪盤的封嚴(yán)篦齒在巨大的熱-機(jī)械載荷交變作用下,很容易在結(jié)構(gòu)強(qiáng)度最為薄弱的篦齒齒尖處出現(xiàn)疲勞裂紋[2].在周向拉伸應(yīng)力的驅(qū)動下,篦齒齒尖萌生的微觀裂紋將持續(xù)擴(kuò)展到篦齒根部并沿篦齒連接段向渦輪盤體方向發(fā)展,形成半穿透或穿透性裂紋.此外,短時燃燒室局部排氣溫度過高、燃?xì)獾构?、封?yán)篦齒與封嚴(yán)環(huán)間可能產(chǎn)生的碰磨,都大大增大了渦輪盤篦齒裂紋故障的危害.據(jù)外場不完全統(tǒng)計,大約3%的某機(jī)I級渦輪盤由于篦齒裂紋故障在達(dá)到使用壽命之前報廢[3];據(jù)某廠故障檢測原始資料,在翻修的918臺次某乙型發(fā)動機(jī)中,有154臺次出現(xiàn)渦輪盤封嚴(yán)篦齒裂紋,故障率為16.78%[4].不僅在經(jīng)濟(jì)上造成巨大損失,渦輪盤篦齒裂紋所引發(fā)的頻繁檢修與盤件更換也嚴(yán)重影響到了飛機(jī)的出勤率.

    國內(nèi)圍繞著渦輪盤篦齒裂紋問題開展了很多研究,通過殘余應(yīng)力測量、疲勞試驗、斷裂力學(xué)估算以及斷口分析等多種技術(shù)手段,對封嚴(yán)篦齒裂紋的故障起因與特征進(jìn)行了歸納與分析,并為防止危害性結(jié)果的發(fā)生給出了許多指導(dǎo)性的建議[3-7].

    限于計算力學(xué)發(fā)展水平,以往對篦齒裂紋故障的研究比較缺乏精確定量的數(shù)值模擬,因而無法詳細(xì)考察裂紋前沿的應(yīng)力強(qiáng)度因子分布與裂紋擴(kuò)展規(guī)律,在確定臨界裂紋長度與結(jié)構(gòu)剩余壽命時缺乏理論依據(jù),并且難以制定出兼顧經(jīng)濟(jì)性與安全性的檢修周期.

    本文使用有限元軟件ANSYS 13.0與三維裂紋擴(kuò)展分析軟件FRANC3D6.0對渦輪盤篦齒裂紋擴(kuò)展進(jìn)行直接數(shù)值模擬.首先采用子模型法建立了含篦齒裂紋的渦輪盤有限元模型,使用M積分[8-9]計算裂紋前沿的應(yīng)力強(qiáng)度因子;在確定篦齒裂紋前沿每一節(jié)點(diǎn)處的局部擴(kuò)展方向與擴(kuò)展距離后,通過樣條曲線擬合出新裂紋前沿,并依靠自適應(yīng)網(wǎng)格劃分[10]實現(xiàn)裂紋區(qū)有限元網(wǎng)格的更新;得出篦齒裂紋前沿的歷史軌跡后,在應(yīng)力強(qiáng)度因子計算結(jié)果的基礎(chǔ)上,對擴(kuò)展速率模型進(jìn)行積分得到篦齒裂紋擴(kuò)展壽命[11].

    1 篦齒渦輪盤的應(yīng)力分析

    某篦齒渦輪盤的實體模型如圖1所示,依次將其上的2個篦齒記為篦齒1與篦齒2.輪盤材料采用變形高溫合金,對應(yīng)的有限元模型如圖2所示,采用SOLID186單元進(jìn)行離散,單元總數(shù)為329 954個.有限元模型包含篦齒渦輪盤0°~360°的完整回轉(zhuǎn)體,而沒有采用只對一個扇區(qū)進(jìn)行建模的周期對稱方法.其原因在于,如果采用周期對稱,在1/N扇區(qū)的輪盤模型上插入的一條篦齒裂紋,等效于渦輪盤上沿周向均布著N條大小、形狀完全相同的篦齒裂紋,渦輪盤由此轉(zhuǎn)化為包含裂紋間相互干涉作用的多裂紋結(jié)構(gòu).為從有限元數(shù)值模擬中突出渦輪盤篦齒裂紋擴(kuò)展的主要特征,本文對裂紋出現(xiàn)的條數(shù)進(jìn)行了簡化,僅考慮含單一篦齒裂紋的渦輪盤結(jié)構(gòu).為了便于后續(xù)的裂紋區(qū)子模型的建立,在網(wǎng)格劃分時對即將插入裂紋的位置進(jìn)行了局部網(wǎng)格加密處理.

    圖1 篦齒渦輪盤的實體模型Fig.1 Solid model of a turbine disc with labyrinth seal

    圖2 篦齒渦輪盤的有限元模型Fig.2 Finite element model of a turbine disc with labyrinth seal

    以高溫起飛為計算工況,在篦齒渦輪盤有限元模型上施加對應(yīng)的轉(zhuǎn)速與溫度場,將葉片質(zhì)量產(chǎn)生的離心載荷轉(zhuǎn)換為等效的徑向面載荷施加在圖1中的截面A上,將渦輪轉(zhuǎn)子軸對稱分析得出的軸向位移與徑向位移施加在圖1中的截面B上,取自網(wǎng)格加密區(qū)域的徑向、周向、軸向應(yīng)力云圖如圖3(a)、圖3(b)、圖3(c)所示,分布范圍分別為-9.6~10.9 MPa,350.3~386.4MPa,-0.7~89.6MPa,可見該區(qū)域以周向拉伸應(yīng)力為主導(dǎo).

    圖3 篦齒區(qū)域的有限元應(yīng)力分析結(jié)果Fig.3 The stress results of finite element analysis in the vicinity of labyrinth seal

    篦齒齒尖厚度通常僅為幾百微米,在周向拉伸應(yīng)力與溫度載荷的周期性作用下,很容易萌生微小裂紋,裂紋面往往與周向拉伸應(yīng)力相垂直形成最為危險的Ⅰ型裂紋,在交變載荷作用下繼續(xù)向齒根方向擴(kuò)展.

    2 含篦齒裂紋有限元模型的建立

    為了適應(yīng)裂紋區(qū)網(wǎng)格動態(tài)更新的需要,并同時對模型整體的計算規(guī)模加以控制,采用子模型法建立含篦齒裂紋的渦輪盤有限元模型[12].將篦齒裂紋在擴(kuò)展過程中可能影響到的輪盤體積取出作為子模型,剩余的絕大部分輪盤體積作為全局模型,如圖4所示.

    圖4 篦齒渦輪盤有限元網(wǎng)格中的子模型與全局模型Fig.4 Submodel and global model for the finite element mesh of a turbine disc with labyrinth seal

    篦齒1的周向拉伸應(yīng)力水平高于篦齒2,因此選擇在篦齒1的齒尖插入初始裂紋.根據(jù)故障檢測經(jīng)驗[3-5],假設(shè)初始裂紋為貫穿齒尖厚度的穿透型裂紋,初始裂紋深度為0.254mm,如圖5所示.在插入初始裂紋后,采用自適應(yīng)網(wǎng)格劃分對子模型網(wǎng)格進(jìn)行更新:①裂紋前沿被一圈楔形奇異單元所包圍,該單元在裂尖相鄰兩邊上的中間節(jié)點(diǎn)被移動到1/4邊長位置,以反映裂紋尖端應(yīng)力的奇異性[13];②裂尖奇異元包裹在2層六面體單元圍成的管道內(nèi),以保證裂紋前沿區(qū)域應(yīng)力計算的精確性;③全局模型與子模型交界面上的節(jié)點(diǎn)保持不變,在子模型交界面內(nèi)側(cè)生成一層五面體單元,實現(xiàn)自由網(wǎng)格向映射網(wǎng)格的過渡;④在子模型其他區(qū)域采用四面體自由網(wǎng)格進(jìn)行填充,以保證裂紋區(qū)網(wǎng)格對復(fù)雜幾何形狀的適應(yīng)性.

    含裂紋有限元模型在建模過程中,裂紋前沿的網(wǎng)格尺寸及密度是一項十分重要的考量因素.經(jīng)驗法則認(rèn)為,為確保有限元計算結(jié)果的精確性,裂紋區(qū)域的單元尺寸至少要小于裂紋面的幾何特征尺寸10倍[14]:①對于穿透性裂紋,裂紋面寬度與長度均應(yīng)保持10個單元以上的網(wǎng)格密度;②對于表面裂紋,裂紋面在半徑方向及周長方向也需要保持10個單元以上的網(wǎng)格密度,如圖6所示.經(jīng)測試發(fā)現(xiàn),對于深度為0.254mm的初始裂紋,裂尖奇異元與2層六面體單元圍成的管道半徑控制在0.06mm,圍繞管道的周向布置8個單元(見圖5中的虛線框出部分)即可滿足上述經(jīng)驗法則.在隨后每一步的擴(kuò)展計算中,都需要根據(jù)當(dāng)前裂紋面的幾何尺寸對管道半徑進(jìn)行調(diào)整,在滿足網(wǎng)格密度需求的同時對單元總數(shù)加以控制.

    圖5 初始篦齒裂紋插入前后的子模型與裂紋區(qū)的自適應(yīng)網(wǎng)格劃分Fig.5 Submodels with and without initial crack and self-adapting mesh generation in the vicinity of the crack

    圖6 裂紋區(qū)域網(wǎng)格尺寸及密度的經(jīng)驗法則Fig.6 Experience rules for element size and density in crack region

    3 篦齒裂紋擴(kuò)展的數(shù)值模擬

    圖7 初始情況下篦齒裂紋前沿的應(yīng)力強(qiáng)度因子分布Fig.7 The distributions of the stress intensity factors along the initial crack of the labyrinth seal disc

    根據(jù)含篦齒裂紋有限元模型的位移計算結(jié)果,使用M積分法計算初始裂紋前沿的應(yīng)力強(qiáng)度因子KⅠ,KⅡ,KⅢ(對應(yīng)Ⅰ型、Ⅱ型、Ⅲ型裂紋),結(jié)果如圖7所示.沿篦齒裂紋前沿的KⅠ遠(yuǎn)大于KⅡ及KⅢ,這與周向拉伸應(yīng)力占主導(dǎo)的應(yīng)力計算結(jié)果吻合;KⅠ沿裂紋前沿變化不大,處于中部的大部分長度上數(shù)值基本保持恒定,約為375MPa;在篦齒裂紋與外部自由表面的2個交點(diǎn)P1,P2附近(見圖5,篦齒1齒尖靠近篦齒2的一側(cè)交點(diǎn)記為P1,對應(yīng)圖7中橫坐標(biāo)為0的位置),KⅠ有一定的降低,交點(diǎn)P1,P2附近的裂紋擴(kuò)展速率將因此低于中部區(qū)域,裂紋前沿將從初始的直線逐步過渡到一條中部下凹、兩端上翹的曲線.

    方案階段選取0-Max-0三角轉(zhuǎn)速譜,將初始的篦齒裂紋前沿記為C0,根據(jù)C0各節(jié)點(diǎn)處的應(yīng)力強(qiáng)度因子可以通過擴(kuò)展操作得出新裂紋前沿C1,過程如下:

    1)根據(jù)C0前沿KⅠ,KⅡ,KⅢ的M積分結(jié)果,使用最大環(huán)向拉伸應(yīng)力準(zhǔn)則[15]計算C0每一節(jié)點(diǎn)處的裂紋局部擴(kuò)展方向;

    2)確定擴(kuò)展步長Δa(m)(Δa(m)為裂紋前沿“中位節(jié)點(diǎn)”處的擴(kuò)展距離,該節(jié)點(diǎn)處的應(yīng)力強(qiáng)度因子為裂紋前沿的中位數(shù),本文C0前沿對應(yīng)的Δa(m)取0.038mm,后續(xù)擴(kuò)展過程將逐步增大Δa(m)的數(shù)值),然后根據(jù)Δa(m)的取值,利用Paris公式推算出C0前沿其他節(jié)點(diǎn)處的擴(kuò)展距離;

    3)根據(jù)C0前沿各節(jié)點(diǎn)處的擴(kuò)展方向與擴(kuò)展距離,得到C1前沿的離散預(yù)測結(jié)果,然后采用三次樣條曲線進(jìn)行擬合,得到經(jīng)過光順化處理的C1前沿;

    4)圍繞C1前沿生成一圈楔形奇異元并在外部裹以2層六面體單元,在子模型交界面內(nèi)側(cè)生成1層五面體單元并在子模型其他區(qū)域填充四面體自由網(wǎng)格,完成裂紋區(qū)網(wǎng)格的更新,同時全局模型保持不變;

    重復(fù)上述步驟依次得到C2,C3,C4,…,直至裂紋即將穿透篦齒環(huán)為止,共經(jīng)歷53步擴(kuò)展,其中C52向C53的擴(kuò)展如圖8所示.C0至C53描述了擴(kuò)展過程中篦齒裂紋前沿變化的歷史軌跡,如圖9所示.

    圖8 篦齒裂紋前沿C52向前沿C53的擴(kuò)展Fig.8 The extension from crack front C52to crack front C53for the crack of the labyrinth seal disc

    圖9 C0至C53所描述的篦齒裂紋前沿變化歷史軌跡Fig.9 The history tracks for the variations of the crack fronts of the labyrinth seal disc represented by C0to C53

    需要指出,上述步驟中所提到的裂紋擴(kuò)展步長并非裂紋每一循環(huán)的擴(kuò)展量.整個分析過程中裂紋擴(kuò)展步的數(shù)量,以及其中每一步的步長值,都需要根據(jù)所研究問題的實際情況具體指定.過大的擴(kuò)展步長一方面可能導(dǎo)致不精確的結(jié)果,另一方面也有可能使預(yù)測出的裂紋前沿出現(xiàn)不穩(wěn)定的變化(因為未能精確捕捉到擴(kuò)展過程中裂紋前沿載荷重新分布的過程);過小的擴(kuò)展步長同樣應(yīng)該避免,因為它一方面會導(dǎo)致過多的擴(kuò)展步數(shù)與過大的計算開銷,另一方面還有可能導(dǎo)致網(wǎng)格劃分困難.

    第i步的擴(kuò)展步長Δa(m),i有2種確定方法.

    第1種是指定第i個擴(kuò)展步內(nèi)所經(jīng)歷的循環(huán)次數(shù),即

    其中:(da/dN)(m),i是中位節(jié)點(diǎn)的當(dāng)前裂紋擴(kuò)展速率;ΔNi是第i步所經(jīng)歷的載荷循環(huán)數(shù),需由人工指定.由于裂紋擴(kuò)展速率與裂紋尺寸間高度的非線性關(guān)系,往往需要通過“trial-and-error”的方式獲得使Δa(m),i大小合適的ΔNi.

    確定Δa(m),i的另一種方法是,以一個裂紋當(dāng)前的幾何特征尺寸為參照,以該尺寸的一定百分比作為Δa(m),i的數(shù)值.本文采用篦齒在渦輪盤半徑方向上的深度作為特征尺寸,初始時取其數(shù)值的15%作為擴(kuò)展步長.初始裂紋前沿C0的深度為0.254mm,則與其對應(yīng)的Δa(m)=0.254mm×15%≈0.038mm.根據(jù)筆者的經(jīng)驗,該比例取在5%~15%時,可以避免裂紋前沿出現(xiàn)不穩(wěn)定的變化,同時將計算開銷控制在合理范圍內(nèi).一般來說,由于初始裂紋尺寸往往十分微小,這時宜采用相對較大的比例;隨著裂紋在擴(kuò)展過程中尺寸的不斷增大,該比例可以視情況酌情采用較小的數(shù)值.影響擴(kuò)展步長比例的并非單純只有裂紋幾何方面的因素,由于裂紋區(qū)域網(wǎng)格劃分的復(fù)雜性,某一擴(kuò)展步長比例可能會導(dǎo)致網(wǎng)格劃分失敗或網(wǎng)格高度畸形,這時也需要對比例數(shù)值作出適當(dāng)調(diào)整.

    篦齒渦輪盤材料為某國產(chǎn)變形高溫合金,根據(jù)文獻(xiàn)[16]中基于標(biāo)準(zhǔn)CT試樣的裂紋擴(kuò)展速率試驗結(jié)果(圖10),在篦齒1高溫起飛時的工作溫度下,其材料的裂紋擴(kuò)展速率da/dN與應(yīng)力強(qiáng)度因子變程ΔK在雙對數(shù)坐標(biāo)系下可以用線性關(guān)系進(jìn)行近似,擬合方程見式(2),擬合優(yōu)度為0.984 2,轉(zhuǎn)化為Paris冪函數(shù)表達(dá)式后得出的材料常數(shù)C=1.445 4×10-17,n=4.333,采用的單位制為:長度,mm;應(yīng)力強(qiáng)度因子,MPa;裂紋擴(kuò)展速率,mm/循環(huán).

    在裂紋擴(kuò)展模擬過程中,Paris公式的意義不僅僅在于通過積分得出裂紋擴(kuò)展壽命,它還影響著每一個擴(kuò)展步中裂紋前沿不同節(jié)點(diǎn)之間擴(kuò)展距離的相對大小.設(shè)j,k是當(dāng)前裂紋前沿2個不同的節(jié)點(diǎn),2處的應(yīng)力強(qiáng)度因子變程分別為ΔK(j),ΔK(k),若施齒1的對稱面與圖9所示徑向平面的交線作為擴(kuò)展壽命的計算路徑(圖9中的虛線),該路徑上的應(yīng)力強(qiáng)度因子KⅠ隨路徑長度的變化曲線如圖12所示,路徑長度的起點(diǎn)取在篦齒1齒尖上端面的O點(diǎn)(圖9),A點(diǎn)的路徑長度等于初始裂紋深度,為0.254mm.

    圖10 根據(jù)裂紋擴(kuò)展速率試驗得出的Paris模型Fig.10 Paris model based on the crack growth rate test

    加在當(dāng)前裂紋擴(kuò)展步的載荷循環(huán)數(shù)為ΔN,則根據(jù)Paris公式可得出該擴(kuò)展步中節(jié)點(diǎn)j,k各自的擴(kuò)展距離為

    式(3)與式(4)左右兩邊分別相除,消去相同的因子,可得

    由此可見,Δa(j)與Δa(k)之比不僅取決于裂紋前沿的應(yīng)力強(qiáng)度因子分布,還取決于Paris公式中的材料常數(shù)n.本文以當(dāng)前裂紋前沿“中位節(jié)點(diǎn)”的擴(kuò)展距離Δa(m)為輸入,根據(jù)式(5)計算出其他各節(jié)點(diǎn)處的擴(kuò)展距離,進(jìn)而預(yù)測擴(kuò)展步結(jié)束后的新裂紋前沿,如圖11所示.

    圖11 以擴(kuò)展步長Δa(m)為輸入預(yù)測新裂紋前沿Fig.11 Prediction of new crack front takingΔa(m)as input

    4 裂紋擴(kuò)展壽命計算與分析

    由于初始裂紋的裂紋面與占主導(dǎo)的周向拉伸應(yīng)力相垂直,后續(xù)的裂紋擴(kuò)展基本保持在一個徑向平面內(nèi)進(jìn)行,該徑向平面經(jīng)過發(fā)動機(jī)的旋轉(zhuǎn)軸線.將篦

    圖12 計算路徑上KⅠ隨路徑長度的變化曲線Fig.12 The variation curve of KⅠversus path length along the path

    根據(jù)圖12中KⅠ曲線的轉(zhuǎn)折變化,從中提取A,B,…,G共7個特征點(diǎn),并在圖9中依次標(biāo)出它們在裂紋前沿軌跡上的對應(yīng)位置,可以看出:①A-B-C是篦齒裂紋從齒尖沿徑向擴(kuò)展直至達(dá)到齒根的過程,在該過程中,KⅠ隨裂紋擴(kuò)展逐步增大;②C-D-E是篦齒裂紋逐步穿透齒根過渡圓角的過程,在該過程中,篦齒裂紋由穿透型裂紋轉(zhuǎn)變?yōu)楸砻媪鸭y,過渡圓角致使裂紋在擴(kuò)展過程中遇到結(jié)構(gòu)厚度的迅速增加,因此C-D-E階段KⅠ隨裂紋擴(kuò)展有逐步減小的趨勢;③E-F-G是篦齒裂紋轉(zhuǎn)變?yōu)楸砻媪鸭y后的擴(kuò)展過程,該階段KⅠ隨裂紋擴(kuò)展逐步增大.

    由圖12擬合出KⅠ與路徑長度間的函數(shù)關(guān)系,然后沿計算路徑對裂紋擴(kuò)展速率模型(Pairs公式)進(jìn)行積分,得出路徑長度隨載荷循環(huán)數(shù)的變化曲線,如圖13所示,圖中標(biāo)出了A,B,…,G各點(diǎn)的對應(yīng)位置以及所經(jīng)歷的載荷循環(huán)數(shù),可以看出:①篦齒裂紋擴(kuò)展的全過程(從A到G)共經(jīng)歷20 953循環(huán);②篦齒裂紋從齒尖沿徑向擴(kuò)展直至達(dá)到齒根的A-B-C階段經(jīng)歷15 414循環(huán),占總壽命的73.6%;③篦齒裂紋穿透齒根過渡圓角的C-D-E過程經(jīng)歷2 687循環(huán),占總壽命的12.8%;④篦齒裂紋轉(zhuǎn)變?yōu)楸砻媪鸭y后的擴(kuò)展過程E-F-G經(jīng)歷2 852循環(huán),占總壽命的13.6%.

    對圖13中的曲線進(jìn)行微分,可以得出路徑長度增長速率隨載荷循環(huán)數(shù)的變化曲線,以此衡量各階段篦齒裂紋擴(kuò)展的快慢,如圖14所示.從中可以看出:①篦齒裂紋從齒尖沿徑向擴(kuò)展直至達(dá)到齒根的A-B-C階段,裂紋擴(kuò)展處于不斷加速的狀態(tài);②由于過渡圓角引起的結(jié)構(gòu)厚度增加,C-D-E階段的裂紋擴(kuò)展速度有所減緩;③在篦齒裂紋完全轉(zhuǎn)化為表面裂紋的E-F-G階段,裂紋擴(kuò)展速度首先很快地增大到最高點(diǎn),然后出現(xiàn)一定的下降,這是因為在E-F-G的后半段,裂紋前沿的擴(kuò)張開始受到來自于相鄰篦齒過渡圓角的約束作用,如圖9所示.

    圖13 路徑長度隨載荷循環(huán)數(shù)的變化曲線Fig.13 The variation of path length versus load cycle numbers

    圖14 路徑長度增長速率隨載荷循環(huán)數(shù)的變化曲線Fig.14 The variation curve of path length growth rate versus load cycle numbers

    前面的裂紋擴(kuò)展計算中采用了如圖10所示的Paris公式作為裂紋擴(kuò)展速率模型,圖10中150個試驗點(diǎn)ΔK的變動范圍是558.7~2 117.6MPa;而結(jié)合圖9與圖12中可以看出,篦齒裂紋擴(kuò)展的全過程(從A到G)計算路徑上KⅠ的變動范圍是376.8~1 054.6MPa,由此可見,圖10中Paris公式的適用范圍并不能保證覆蓋篦齒裂紋擴(kuò)展的全過程,因為初始篦齒裂紋C0之后0.3~0.4mm的擴(kuò)展是在擴(kuò)展速率試驗的ΔK范圍之外進(jìn)行的.

    由于試驗數(shù)據(jù)的缺乏,目前無法確定C0之后0.3~0.4mm的擴(kuò)展是否依然符合Paris公式,因此需對壽命計算結(jié)果的誤差進(jìn)行分析.在雙對數(shù)坐標(biāo)系中,裂紋擴(kuò)展速率da/dN隨ΔK的完整變化關(guān)系曲線近似呈反“S”形[17-18],分為初始穩(wěn)定擴(kuò)展區(qū)、線性穩(wěn)定擴(kuò)展區(qū)(Paris區(qū))以及快速擴(kuò)展區(qū)[13];在初始穩(wěn)定擴(kuò)展區(qū)da/dN-ΔK曲線存在一個垂直漸近線ΔK=ΔKth,ΔKth通常稱裂紋擴(kuò)展門檻值.

    對于所研究的篦齒盤材料,由文獻(xiàn)[16]中的擴(kuò)展速率試驗數(shù)據(jù)只能估算出初始穩(wěn)定擴(kuò)展區(qū)與Paris區(qū)的分界點(diǎn)ΔK*≤558.7MPa;如果ΔK*的實際值小于等于376.8MPa,則Paris公式可以覆蓋篦齒裂紋擴(kuò)展的全過程,包括C0之后0.3~0.4mm的擴(kuò)展;如果ΔK*的實際值在376.8~558.7MPa之間,則ΔK<ΔK*時的裂紋擴(kuò)展將不再服從Paris公式,此時da/dNΔK曲線位于Paris直線的下方,依然使用Paris公式等效于高估裂紋擴(kuò)展速率,在這種情況下對于圖9中A-B段,基于Paris公式將得出偏于保守的壽命計算結(jié)果(也就是說ΔK*>376.8MPamm會致使A-B段實際經(jīng)歷的載荷循環(huán)數(shù)高于圖13中給出的12610循環(huán)).

    5 結(jié) 論

    根據(jù)應(yīng)力強(qiáng)度因子以及裂紋擴(kuò)展速率的變化,可以將篦齒裂紋的擴(kuò)展過程分為3個階段:穿透型裂紋擴(kuò)展階段(A-B-C)、穿透型裂紋向表面裂紋的轉(zhuǎn)化階段(C-D-E)以及表面裂紋擴(kuò)展階段(E-F-G). A-B-C是篦齒裂紋從齒尖沿徑向擴(kuò)展直至達(dá)到齒根的過程,應(yīng)力強(qiáng)度因子KⅠ隨裂紋擴(kuò)展逐步增大,裂紋擴(kuò)展處于不斷加速的狀態(tài);C-D-E是篦齒裂紋逐步穿透齒根過渡圓角的過程,KⅠ隨裂紋擴(kuò)展有逐步減小的趨勢,裂紋擴(kuò)展速度有所減緩;E-F-G是篦齒裂紋轉(zhuǎn)變?yōu)楸砻媪鸭y后的擴(kuò)展過程,KⅠ隨裂紋擴(kuò)展逐步增大,裂紋擴(kuò)展速度首先增大到最高點(diǎn),然后在相鄰篦齒的約束作用下擴(kuò)展速度出現(xiàn)一定下降.

    A-B-C與C-D-E-F-G的路徑長度比為1∶1.5,擴(kuò)展壽命比為3∶1,所以裂紋擴(kuò)展速度之比為1∶4.5,這說明篦齒裂紋前沿在達(dá)到齒根過渡圓角后平均擴(kuò)展速度會大大增加,所以筆者建議以齒根過渡圓角頂部(圖9中的C點(diǎn))為分界,將篦齒裂紋擴(kuò)展分為緩慢擴(kuò)展階段(A-B-C)與快速擴(kuò)展階段(C-D-E-F-G).

    緩慢擴(kuò)展的A-B-C階段,應(yīng)力強(qiáng)度因子與裂紋擴(kuò)展速度均處于較低水平,裂紋擴(kuò)展局限于齒根以上區(qū)域,裂紋在發(fā)動機(jī)旋轉(zhuǎn)時產(chǎn)生的不平衡力較小,并且能夠承受相當(dāng)數(shù)量的載荷循環(huán)(本例中為15 414循環(huán)),所以渦輪盤即使在檢修中發(fā)現(xiàn)了尚未擴(kuò)展至齒根的篦齒裂紋,也不必立即作更換處理,但需在裂紋進(jìn)入齒根以下的快速擴(kuò)展階段之前采取措施處理.

    文獻(xiàn)[16]中裂紋擴(kuò)展速率試驗的ΔK范圍有限,因此對壽命計算結(jié)果可能出現(xiàn)的誤差進(jìn)行了分析:①若ΔK*≤376.8MPa,則可全程使用Paris公式進(jìn)行裂紋擴(kuò)展壽命計算;②若ΔK*>376.8MPa,則使用Paris公式會在篦齒裂紋擴(kuò)展的起始階段高估裂紋擴(kuò)展速率,從而得出偏于保守的壽命計算結(jié)果.為了獲得更精確的篦齒裂紋擴(kuò)展壽命,有必要進(jìn)一步補(bǔ)充ΔKth<ΔK<ΔK*區(qū)域內(nèi)的裂紋擴(kuò)展速率試驗.

    由本文分析結(jié)果還可以得出2條有關(guān)篦齒結(jié)構(gòu)設(shè)計的參考性建議:①適當(dāng)增大篦齒齒根部的過渡圓角半徑有助于延緩篦齒裂紋從穿透型裂紋向表面裂紋的轉(zhuǎn)化;②適當(dāng)減小篦齒的排列間距,可以利用相鄰篦齒齒根對裂紋后期的快速擴(kuò)展而起一定約束作用.

    現(xiàn)階段筆者對于渦輪盤篦齒裂紋擴(kuò)展的研究側(cè)重于理論方法與數(shù)值模擬,隨著研制工作的不斷深入,后續(xù)將安排相關(guān)試驗,其目的一方面是對數(shù)值模擬的結(jié)果進(jìn)行驗證并加以必要的修正,另一方面是結(jié)合數(shù)值模擬與試驗結(jié)果建立一套行之有效的渦輪盤篦齒損傷容限設(shè)計方法.

    [1]黃彬彬.航空發(fā)動機(jī)篦齒封嚴(yán)特性數(shù)值模擬研究[D].北京:北京交通大學(xué)機(jī)械與電子控制工程學(xué)院,2014:2-9.

    HUANG Bin-bin.Numerical simulation of leakage character of aero-engine labyrinth seals[D].Beijing:Beijing Jiaotong University,School of Mechanical,Electronic & Control Engineering,2014:2-9.

    [2]陶春虎,鐘培道,王仁智,等.航空發(fā)動機(jī)轉(zhuǎn)動部件的失效分析與預(yù)防[M].北京:國防工業(yè)出版社,2000:138-147.

    TAO Chun-hu,ZHONG Pei-dao,WANG Ren-zhi,et al. Failure analysis and prevention for rotor in aero-engine[M].Beijing:National Defense Industry Press,2000:138-147.

    [3]龔夢賢,王衛(wèi)國,古遠(yuǎn)興,等.某發(fā)動機(jī)I級渦輪盤篦齒裂紋故障分析[J].燃?xì)鉁u輪試驗與研究,1999,12(4):34-39.

    GONG Meng-xian,WANG Wei-guo,GU Yuan-xing,et al.Failure analysis of sealing labyrinth crack for the first stage turbine disc in an aero-engine[J].Gas Turbine Experiment and Research,1999,12(4):34-39.

    [4]唐光輝.某乙型發(fā)動機(jī)Ⅰ級渦輪盤封嚴(yán)篦齒裂紋故障分析及維修措施[J].航空發(fā)動機(jī),1995,21(4):50-58.

    TANG Guang-hui.Failure analysis and maintenance measures of sealing labyrinth crack for the first stage turbine disc in a type-B areo-engine[J].Aeroengine,1995,21(4):50-58.

    [5]李運(yùn)菊,孫智軍.某發(fā)動機(jī)高壓Ⅱ級渦輪盤封嚴(yán)篦齒裂紋分析[J].失效分析與預(yù)防,2006,1(3):27-30.

    LI Yun-ju,SUN Zhi-jun.Crack analysis of sealing labyrinth on high-pressure turbine discⅡin an aero-engine[J].Failure Analysis and Prevention,2006,1(3):27-30.

    [6]田武剛,潘孟春,羅飛路,等.某型航空發(fā)動機(jī)篦齒盤裂紋的原位渦流檢測[J].測試技術(shù)學(xué)報,2008,22(3):241-244.

    TIAN Wu-gang,PAN Meng-chun,LUO Fei-lu,et al.In situ eddy current testing for the cracks in an aeroengine labyrinth disc[J].Journal of Test and Measurement Technology,2008,22(3):241-244.

    [7]卓正敏,李文闖.渦輪盤篦齒補(bǔ)焊工藝研究[J].金屬加工(熱加工),2013(2):74-75.

    ZHUO Zheng-ming,LI Wen-chuang.Research of repair welding technology for the sealing labyrinth in a turbine disc[J].Machinist Metal Forming,2013(2):74-75.

    [8]YAU J F,WANG S S,CORTEN H T.A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity[J].Journal of Applied Mechanics,1980,47(2):335-341.

    [9]WAWRZYNEK P A,CARTER B J,BANKS-SILLS L. The M-integral for computing stress intensity factors in generally anisotropic materials:NASA/CR-2005-214006[R].Alabama:Marshall Space Flight Center,2005.

    [10]HOMBAL V K,MAHADEVAN S.Surrogate modeling of 3Dcrack growth[J].International Journal of Fatigue,2013,47(1):90-99.

    [11]WAWRZYNEK P A,CARTER B J,INGRAFFEA A R.Advances in simulation of arbitrary 3Dcrack growth using FRANC3D/NG[C]//Proceedings of the 12th International Conference on Fracture.Ottawa,Jul.12-17,2009.

    [12]CARTER B J,SCHENCK E C,WAWRZYNEK P A,et al.Three-dimensional simulation of fretting crack nucleation and growth[J].Engineering Fracture Mechanics,2012,96(1):447-460.

    [13]酈正能,關(guān)志東,張紀(jì)奎,等.應(yīng)用斷裂力學(xué)[M].北京:北京航空航天大學(xué)出版社,2012:64-73.

    LI Zheng-neng,GUAN Zhi-dong,ZHANG Ji-kui,et al. Application fracture mechanics[M].Beijing:Beihang University Press,2012:64-73.

    [14]FAC.FRANC3Dreference manual[EB/OL].(2012-05-08)[2015-07-15].http://www.franc3d.com/wp-content/uploads/2012/05/FRANC3D_V7_Reference.pdf.[15]BUCZEK M B,HERAKOVICH C T.A normal stress criterion for crack extension direction in orthotropic composite materials[J].Journal of Composite Materials,1985,19(6):544-553.

    [16]航空發(fā)動機(jī)設(shè)計用材料數(shù)據(jù)手冊編委會.航空發(fā)動機(jī)設(shè)計用材料數(shù)據(jù)手冊(第3冊)[M].北京:航空工業(yè)出版社,2008:166-170.

    The Editorial Board of Material Data Handbook for Design of Aero-engines.Material data handbook for design of aero-engines(Volume 3)[M].Beijing:Aviation Industry Press,2008:166-170.

    [17]趙志平,李有堂,黨延祖.金屬材料疲勞設(shè)計及裂紋擴(kuò)展規(guī)律研究與應(yīng)用[J].石油礦場機(jī)械,2010,39(9):38-41.

    ZHAO Zhi-ping,LI You-tang,DANG Yan-zu.Study and application on fatigue crack propagation law and design of metal material[J].Oil Field Equipment,2010,39(9):38-41.

    [18]倪向貴,李新亮,王秀喜.疲勞裂紋擴(kuò)展規(guī)律Paris公式的一般修正及應(yīng)用[J].壓力容器,2006,23(12):8-15.

    NI Xiang-gui,LI Xin-liang,WANG Xiu-xi.General modification and application of the Paris law for fatigue crack propagation[J].Pressure Vessel Technology,2006,23(12):8-15.

    Finite element simulation for the growth of the sealing labyrinth crack in turbine disc

    REN Yuan1,ZHANG Cheng-cheng1,GAO Jing-yun1,LI Meng-guang2
    (1.Shanghai Center of Research for Commercial Aircraft Engine Engineering Techniques,AVIC(Aviation Industry Corporation of China)Commercial Aircraft Engine Co.,Ltd.,Shanghai 200240,China;2.Multiangle Virtual Technology Incorporation,Shanghai 200233,China)

    In order to investigate the growth mechanism of sealing labyrinth crack,the propagation process,which started from the tooth tip and stops when the crack was about to penetrate through the labyrinth seal ring,was numerically simulated.The finite element model of the turbine disc with a sealing labyrinth crack was established using submodeling method,and the stress intensity factors along the crack front were calculated with M integral.After determining the local direction and extension distance of each node along the sealing labyrinth crack front,the new crack front was fitted with spline curve,and then self-adaption mesh generation was used for updating the finite elements in the vicinity of the crack.The results of numerical simulation indicated that the top of the root fillet could be regarded as the demarcation point between the slow crack-growth stage and the fast crack-growth stage.The sealing labyrinth crack grew with a low velocity as a penetrating crack before the demarcation point,and it began to transfer to a surface crack with a much higher average velocity after the demarcation point.The crack growth life before the demarcation point was several times longer than the life after the demarcation point. From the results of numerical simulation,it can also be concluded that both increasing the radius of root fillet and reducing the space between adjacent sealing labyrinths are helpful for resisting the growth of the sealing labyrinth crack.

    sealing labyrinth crack;turbine disc;crack growth;stress intensity factor;fatigue life

    V 231.95;O 346.2

    A

    1006-754X(2016)02-0152-08

    10.3785/j.issn.1006-754X.2016.02.008

    2015-10-19.本刊網(wǎng)址·在線期刊:http://www.journals.zju.edu.cn/gcsjxb

    上海市“引進(jìn)技術(shù)消化與吸收”專項資金資助(12XI-04);上??莆蒲杏媱濏椖浚?3DJ1400200).

    任遠(yuǎn)(1982—),男,四川南充人,工程師,博士,從事航空發(fā)動機(jī)結(jié)構(gòu)壽命與可靠性技術(shù)研究,E-mail:renyuan@acae. com.cn.http://orcid.org//0000-0001-8644-1234

    猜你喜歡
    裂紋有限元模型
    一半模型
    裂紋長度對焊接接頭裂紋擴(kuò)展驅(qū)動力的影響
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
    Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
    3D打印中的模型分割與打包
    磨削淬硬殘余應(yīng)力的有限元分析
    預(yù)裂紋混凝土拉壓疲勞荷載下裂紋擴(kuò)展速率
    基于SolidWorks的吸嘴支撐臂有限元分析
    箱形孔軋制的有限元模擬
    上海金屬(2013年4期)2013-12-20 07:57:18
    精品久久久久久电影网| 女的被弄到高潮叫床怎么办| 久久久国产一区二区| √禁漫天堂资源中文www| 国产成人91sexporn| 日本猛色少妇xxxxx猛交久久| 自拍欧美九色日韩亚洲蝌蚪91| 男女国产视频网站| 一区二区三区四区激情视频| 久久精品久久精品一区二区三区| 国语对白做爰xxxⅹ性视频网站| av国产久精品久网站免费入址| 国产精品嫩草影院av在线观看| 麻豆乱淫一区二区| 免费黄网站久久成人精品| 久久精品国产亚洲av天美| 国产又色又爽无遮挡免| 亚洲国产日韩一区二区| 久久久国产精品麻豆| 亚洲精品美女久久av网站| 伊人久久精品亚洲午夜| 日韩电影二区| 国产又色又爽无遮挡免| 国产精品一区二区三区四区免费观看| 97超碰精品成人国产| 久久久欧美国产精品| 天堂俺去俺来也www色官网| 99久久人妻综合| 亚洲三级黄色毛片| 欧美变态另类bdsm刘玥| 亚洲欧美成人精品一区二区| 国产一区有黄有色的免费视频| 妹子高潮喷水视频| 水蜜桃什么品种好| 最后的刺客免费高清国语| 亚洲国产最新在线播放| 国产片特级美女逼逼视频| 2021少妇久久久久久久久久久| 中文字幕精品免费在线观看视频 | 九九在线视频观看精品| 九九在线视频观看精品| 黄片无遮挡物在线观看| 亚州av有码| 激情五月婷婷亚洲| 国产精品三级大全| 国产高清有码在线观看视频| 一级毛片我不卡| 亚洲成人手机| 搡老乐熟女国产| 自拍欧美九色日韩亚洲蝌蚪91| 日韩熟女老妇一区二区性免费视频| 美女脱内裤让男人舔精品视频| 两个人的视频大全免费| 美女大奶头黄色视频| 国产极品天堂在线| 一区在线观看完整版| 少妇人妻 视频| 一本色道久久久久久精品综合| 2018国产大陆天天弄谢| 18禁在线播放成人免费| 日日爽夜夜爽网站| 在现免费观看毛片| 亚洲精品日本国产第一区| 日本欧美国产在线视频| 插阴视频在线观看视频| 黄片播放在线免费| 国产精品三级大全| 综合色丁香网| 国产在线免费精品| 国产精品不卡视频一区二区| 久久毛片免费看一区二区三区| 十八禁高潮呻吟视频| 18禁动态无遮挡网站| 免费看av在线观看网站| 999精品在线视频| 一个人免费看片子| 亚洲欧美色中文字幕在线| 久久久久久久亚洲中文字幕| 晚上一个人看的免费电影| 婷婷色av中文字幕| 久久99热6这里只有精品| 久久精品夜色国产| 2021少妇久久久久久久久久久| 国产乱来视频区| 亚洲国产精品国产精品| 日日摸夜夜添夜夜添av毛片| 人体艺术视频欧美日本| .国产精品久久| 久久久久久久久久成人| 男女无遮挡免费网站观看| 午夜精品国产一区二区电影| 亚洲欧美日韩另类电影网站| 丰满少妇做爰视频| 国精品久久久久久国模美| 一二三四中文在线观看免费高清| 久久ye,这里只有精品| 精品久久久久久久久亚洲| 国产熟女欧美一区二区| 国产高清不卡午夜福利| 亚洲第一区二区三区不卡| 日韩成人伦理影院| 在线播放无遮挡| 亚洲精品美女久久av网站| 大陆偷拍与自拍| 国产有黄有色有爽视频| 高清黄色对白视频在线免费看| av播播在线观看一区| 91久久精品国产一区二区三区| 99热6这里只有精品| 久久久欧美国产精品| 亚洲av中文av极速乱| 青春草亚洲视频在线观看| 99热这里只有是精品在线观看| 日日摸夜夜添夜夜添av毛片| 婷婷色综合大香蕉| 亚洲精品av麻豆狂野| 亚洲人成网站在线播| 夫妻性生交免费视频一级片| 99精国产麻豆久久婷婷| 少妇丰满av| 午夜免费鲁丝| 水蜜桃什么品种好| 国产精品国产三级国产专区5o| 丝瓜视频免费看黄片| 欧美一级a爱片免费观看看| 蜜桃在线观看..| 精品熟女少妇av免费看| 多毛熟女@视频| 91在线精品国自产拍蜜月| 久久久久久人妻| 色网站视频免费| 久久久久视频综合| 天堂俺去俺来也www色官网| 不卡视频在线观看欧美| 久久久久久久久大av| 亚洲综合精品二区| 国产精品不卡视频一区二区| 99热国产这里只有精品6| 欧美日韩国产mv在线观看视频| 久久久a久久爽久久v久久| 日本黄大片高清| 欧美另类一区| 婷婷成人精品国产| 精品国产露脸久久av麻豆| 亚洲av福利一区| 国产 精品1| 欧美老熟妇乱子伦牲交| 亚洲色图 男人天堂 中文字幕 | a级毛片免费高清观看在线播放| 丰满饥渴人妻一区二区三| 久久久久久久久久成人| 一本色道久久久久久精品综合| 简卡轻食公司| 国产欧美亚洲国产| 成人国语在线视频| 国产69精品久久久久777片| www.av在线官网国产| 精品一区二区免费观看| 蜜桃久久精品国产亚洲av| 人体艺术视频欧美日本| 久久久久人妻精品一区果冻| 九九爱精品视频在线观看| 亚洲成色77777| 日韩中文字幕视频在线看片| 亚洲精品中文字幕在线视频| 免费观看在线日韩| 国国产精品蜜臀av免费| 九九爱精品视频在线观看| 天堂中文最新版在线下载| 2018国产大陆天天弄谢| 狠狠精品人妻久久久久久综合| 色网站视频免费| 九色成人免费人妻av| 在线免费观看不下载黄p国产| 亚洲av成人精品一区久久| 亚洲丝袜综合中文字幕| 80岁老熟妇乱子伦牲交| 看十八女毛片水多多多| 免费观看av网站的网址| 国产精品国产三级专区第一集| 亚洲精品日韩av片在线观看| 黄片无遮挡物在线观看| 五月玫瑰六月丁香| 一边亲一边摸免费视频| 最近中文字幕高清免费大全6| 青春草亚洲视频在线观看| xxxhd国产人妻xxx| 老司机亚洲免费影院| 婷婷色av中文字幕| 久久久精品区二区三区| 欧美 日韩 精品 国产| 少妇被粗大的猛进出69影院 | 欧美97在线视频| 国产在线视频一区二区| 性高湖久久久久久久久免费观看| 国产有黄有色有爽视频| 18在线观看网站| 五月伊人婷婷丁香| 国产精品国产三级国产av玫瑰| 久久久久精品性色| .国产精品久久| 国产亚洲欧美精品永久| 婷婷成人精品国产| 男女免费视频国产| 少妇猛男粗大的猛烈进出视频| 久久久久久久亚洲中文字幕| 久久精品国产自在天天线| 一本一本综合久久| 大香蕉久久网| 国产男人的电影天堂91| 亚洲人成网站在线播| 国产成人精品无人区| 天美传媒精品一区二区| 成年女人在线观看亚洲视频| 99热国产这里只有精品6| 日韩亚洲欧美综合| 大香蕉久久网| a级毛片免费高清观看在线播放| 成年av动漫网址| 亚洲精品aⅴ在线观看| 天美传媒精品一区二区| 日本欧美视频一区| 久久精品国产自在天天线| 午夜福利影视在线免费观看| h视频一区二区三区| 26uuu在线亚洲综合色| 日韩亚洲欧美综合| 亚洲人成77777在线视频| 99热网站在线观看| 丝袜在线中文字幕| 亚洲激情五月婷婷啪啪| 2021少妇久久久久久久久久久| 国产精品99久久久久久久久| 国产精品免费大片| 免费黄频网站在线观看国产| 亚洲精品,欧美精品| 日韩免费高清中文字幕av| 在线播放无遮挡| 久久久久久伊人网av| 久久人人爽人人片av| 久久久久久久久久人人人人人人| 黑人欧美特级aaaaaa片| 国产淫语在线视频| 久久精品久久久久久噜噜老黄| 亚洲国产色片| 成人亚洲欧美一区二区av| 亚洲精品国产色婷婷电影| 亚洲国产精品国产精品| 亚洲美女黄色视频免费看| 春色校园在线视频观看| 性高湖久久久久久久久免费观看| 伦理电影免费视频| 亚洲三级黄色毛片| 亚洲精品日韩在线中文字幕| 久久久久久久精品精品| 国产在线一区二区三区精| a级片在线免费高清观看视频| 国产探花极品一区二区| 你懂的网址亚洲精品在线观看| 久久国产精品大桥未久av| 满18在线观看网站| 亚洲精品,欧美精品| 国产精品 国内视频| 99视频精品全部免费 在线| 美女大奶头黄色视频| 国产精品成人在线| 一级,二级,三级黄色视频| 一级,二级,三级黄色视频| 制服丝袜香蕉在线| 欧美日本中文国产一区发布| 色5月婷婷丁香| 久久av网站| 国产一区二区三区综合在线观看 | 亚洲人成网站在线播| 午夜免费男女啪啪视频观看| 日韩av不卡免费在线播放| 精品久久久久久久久av| 免费日韩欧美在线观看| 欧美丝袜亚洲另类| 国产色婷婷99| 亚洲av中文av极速乱| 人人澡人人妻人| 精品亚洲乱码少妇综合久久| 满18在线观看网站| 制服人妻中文乱码| 欧美xxxx性猛交bbbb| 亚洲三级黄色毛片| 国产免费福利视频在线观看| 日韩中文字幕视频在线看片| 欧美最新免费一区二区三区| 欧美人与性动交α欧美精品济南到 | 免费人妻精品一区二区三区视频| 欧美日韩视频精品一区| av又黄又爽大尺度在线免费看| 曰老女人黄片| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产精品偷伦视频观看了| 久久人人爽av亚洲精品天堂| 在线亚洲精品国产二区图片欧美 | 18禁在线无遮挡免费观看视频| 国产精品麻豆人妻色哟哟久久| 成人午夜精彩视频在线观看| 青春草国产在线视频| 一级毛片aaaaaa免费看小| 国产黄色免费在线视频| 18禁在线播放成人免费| 久久人人爽人人片av| 午夜福利,免费看| 搡女人真爽免费视频火全软件| 国产在视频线精品| 国产淫语在线视频| 日本色播在线视频| 夜夜骑夜夜射夜夜干| 黄色一级大片看看| 91在线精品国自产拍蜜月| 日韩成人伦理影院| 搡老乐熟女国产| 亚洲精品久久成人aⅴ小说 | 七月丁香在线播放| 久久人人爽人人片av| 久久鲁丝午夜福利片| 日本wwww免费看| 免费日韩欧美在线观看| 久久精品国产亚洲av涩爱| 久久久久精品久久久久真实原创| 亚洲国产精品成人久久小说| 91精品三级在线观看| 人人妻人人澡人人看| 99久久精品国产国产毛片| 精品人妻在线不人妻| 免费人妻精品一区二区三区视频| 亚洲婷婷狠狠爱综合网| 国产精品国产av在线观看| 哪个播放器可以免费观看大片| 91精品伊人久久大香线蕉| 另类精品久久| 国产精品嫩草影院av在线观看| 满18在线观看网站| 欧美日韩视频精品一区| 一级片'在线观看视频| 中文精品一卡2卡3卡4更新| av一本久久久久| 一级毛片电影观看| 伦理电影免费视频| 视频中文字幕在线观看| 亚洲精品日本国产第一区| 成人综合一区亚洲| 国产黄频视频在线观看| 亚洲婷婷狠狠爱综合网| 精品一品国产午夜福利视频| 国产午夜精品久久久久久一区二区三区| 久久国产亚洲av麻豆专区| 亚洲欧美色中文字幕在线| 美女国产视频在线观看| 涩涩av久久男人的天堂| 极品人妻少妇av视频| 九九爱精品视频在线观看| 精品国产乱码久久久久久小说| 国产无遮挡羞羞视频在线观看| av在线观看视频网站免费| 少妇人妻精品综合一区二区| 18在线观看网站| 少妇人妻精品综合一区二区| 国产极品天堂在线| 欧美3d第一页| 日本色播在线视频| 狠狠精品人妻久久久久久综合| 观看av在线不卡| 国产精品久久久久久久久免| 国产极品粉嫩免费观看在线 | 人妻制服诱惑在线中文字幕| 日日摸夜夜添夜夜爱| 麻豆成人av视频| 国产日韩一区二区三区精品不卡 | 免费不卡的大黄色大毛片视频在线观看| 亚洲精品一区蜜桃| 色吧在线观看| 精品亚洲乱码少妇综合久久| 色吧在线观看| 夜夜爽夜夜爽视频| 午夜视频国产福利| 久久久国产一区二区| 亚洲第一av免费看| 国产精品三级大全| 日韩在线高清观看一区二区三区| 久久午夜综合久久蜜桃| 在线精品无人区一区二区三| 免费观看性生交大片5| 国产在线视频一区二区| 毛片一级片免费看久久久久| 免费高清在线观看日韩| 婷婷色综合www| av.在线天堂| 我要看黄色一级片免费的| 你懂的网址亚洲精品在线观看| 亚洲欧美清纯卡通| 高清黄色对白视频在线免费看| 久久国内精品自在自线图片| 免费观看av网站的网址| 简卡轻食公司| 久久久午夜欧美精品| 亚洲欧美日韩另类电影网站| 免费人妻精品一区二区三区视频| 91精品国产九色| 国产一区二区在线观看日韩| 亚洲av在线观看美女高潮| 少妇人妻久久综合中文| 亚洲第一av免费看| 男女高潮啪啪啪动态图| 亚洲,欧美,日韩| xxxhd国产人妻xxx| 不卡视频在线观看欧美| 爱豆传媒免费全集在线观看| 韩国av在线不卡| 亚洲欧美成人精品一区二区| 亚洲av中文av极速乱| av天堂久久9| 久久久久久久久久久丰满| 国产欧美亚洲国产| 一本色道久久久久久精品综合| 简卡轻食公司| 午夜免费观看性视频| 超色免费av| 精品久久久精品久久久| 夜夜爽夜夜爽视频| 国产黄频视频在线观看| 一二三四中文在线观看免费高清| 中文字幕人妻熟人妻熟丝袜美| 成人毛片a级毛片在线播放| 欧美另类一区| 婷婷成人精品国产| 一级毛片黄色毛片免费观看视频| 亚洲国产av新网站| 国产精品国产三级国产专区5o| 我的老师免费观看完整版| 大片免费播放器 马上看| 亚洲精品美女久久av网站| 精品国产一区二区久久| 久久ye,这里只有精品| av女优亚洲男人天堂| 男人添女人高潮全过程视频| 欧美+日韩+精品| 亚洲精品自拍成人| 精品久久蜜臀av无| 国产在视频线精品| 精品一区二区三卡| 天天躁夜夜躁狠狠久久av| 中国国产av一级| 一区二区三区乱码不卡18| 亚洲国产精品一区三区| 国产午夜精品一二区理论片| 精品久久蜜臀av无| 久久精品久久精品一区二区三区| 亚洲精品中文字幕在线视频| 日本免费在线观看一区| 人人澡人人妻人| 午夜视频国产福利| 国产亚洲午夜精品一区二区久久| 国产色婷婷99| 一本色道久久久久久精品综合| 999精品在线视频| 母亲3免费完整高清在线观看 | 国产乱人偷精品视频| 日本色播在线视频| 麻豆精品久久久久久蜜桃| av网站免费在线观看视频| 美女cb高潮喷水在线观看| 波野结衣二区三区在线| 亚洲精品日本国产第一区| 欧美亚洲 丝袜 人妻 在线| 亚洲av不卡在线观看| 国产成人免费观看mmmm| 久久这里有精品视频免费| 青青草视频在线视频观看| 亚洲经典国产精华液单| 亚洲精品日本国产第一区| av在线观看视频网站免费| 精品99又大又爽又粗少妇毛片| 97超视频在线观看视频| 国产又色又爽无遮挡免| 午夜福利视频精品| 国产成人精品无人区| 天堂俺去俺来也www色官网| 亚洲欧洲国产日韩| 日本欧美视频一区| 欧美日韩视频高清一区二区三区二| 国产白丝娇喘喷水9色精品| 日本欧美国产在线视频| 国产精品成人在线| 国产亚洲精品第一综合不卡 | 精品少妇久久久久久888优播| 国产成人午夜福利电影在线观看| 国产一区亚洲一区在线观看| av在线app专区| 十分钟在线观看高清视频www| 国产深夜福利视频在线观看| 免费高清在线观看视频在线观看| 色婷婷久久久亚洲欧美| 人人妻人人澡人人爽人人夜夜| 亚洲情色 制服丝袜| 黄色欧美视频在线观看| 久久久久精品久久久久真实原创| 蜜桃国产av成人99| 国产黄色免费在线视频| 国产亚洲精品久久久com| 日本-黄色视频高清免费观看| 国产在视频线精品| 国产av码专区亚洲av| 一区二区三区乱码不卡18| 国产男女超爽视频在线观看| 热re99久久国产66热| 人人妻人人澡人人爽人人夜夜| 女性被躁到高潮视频| 制服诱惑二区| 午夜影院在线不卡| 人人妻人人添人人爽欧美一区卜| 欧美三级亚洲精品| 少妇熟女欧美另类| 久久久久久伊人网av| 国产精品国产av在线观看| 午夜福利视频精品| 一级,二级,三级黄色视频| 亚洲人成网站在线播| 2018国产大陆天天弄谢| 岛国毛片在线播放| 亚洲国产色片| 日韩伦理黄色片| 亚洲丝袜综合中文字幕| 久久国产亚洲av麻豆专区| 边亲边吃奶的免费视频| 国产欧美日韩一区二区三区在线 | 尾随美女入室| 夜夜骑夜夜射夜夜干| 在线观看免费日韩欧美大片 | 伦理电影免费视频| 国产男女超爽视频在线观看| 亚洲精品第二区| 亚洲欧洲国产日韩| 少妇的逼水好多| 亚洲高清免费不卡视频| 国产黄色免费在线视频| 日日啪夜夜爽| 久久影院123| 成人亚洲欧美一区二区av| 一边摸一边做爽爽视频免费| 亚洲av.av天堂| 国产精品不卡视频一区二区| 一本久久精品| av播播在线观看一区| 亚洲国产毛片av蜜桃av| 精品少妇久久久久久888优播| 男女免费视频国产| 亚洲精品国产av蜜桃| 国产亚洲午夜精品一区二区久久| 亚洲经典国产精华液单| 久久人人爽人人片av| 国产成人免费观看mmmm| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 欧美人与性动交α欧美精品济南到 | 免费久久久久久久精品成人欧美视频 | av专区在线播放| 午夜福利视频精品| 亚洲av国产av综合av卡| av在线app专区| 免费大片18禁| 三级国产精品欧美在线观看| 纯流量卡能插随身wifi吗| 男女边吃奶边做爰视频| 久久久久久久久久人人人人人人| 天堂俺去俺来也www色官网| 日韩三级伦理在线观看| 免费大片黄手机在线观看| 精品人妻熟女av久视频| 亚洲国产av新网站| 中国美白少妇内射xxxbb| 亚洲四区av| 3wmmmm亚洲av在线观看| 人妻 亚洲 视频| 最后的刺客免费高清国语| 日本av免费视频播放| 男女无遮挡免费网站观看| 毛片一级片免费看久久久久| 日产精品乱码卡一卡2卡三| 亚洲国产毛片av蜜桃av| 亚洲欧美日韩卡通动漫| 男女边摸边吃奶| 涩涩av久久男人的天堂| 精品少妇久久久久久888优播| 成人毛片a级毛片在线播放| 老熟女久久久| 亚洲情色 制服丝袜| 永久免费av网站大全| 亚洲精品乱码久久久v下载方式| 在线观看国产h片| 大话2 男鬼变身卡| 亚洲不卡免费看| 9色porny在线观看| 免费观看在线日韩| 18禁观看日本| 亚洲欧洲精品一区二区精品久久久 | 亚洲婷婷狠狠爱综合网| 黑丝袜美女国产一区| 老司机影院成人| 国产免费一级a男人的天堂| 51国产日韩欧美| videos熟女内射| 中文精品一卡2卡3卡4更新| 少妇的逼水好多| 中文字幕最新亚洲高清| 亚洲性久久影院| 久久人人爽av亚洲精品天堂| 国产精品不卡视频一区二区| 人人妻人人添人人爽欧美一区卜| 国产精品.久久久| 丰满少妇做爰视频|