999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于機(jī)構(gòu)隨機(jī)振動(dòng)分析的載荷譜非參數(shù)上限統(tǒng)計(jì)歸納方法對(duì)比研究

2022-07-30 08:26:00馮志杰肖慧婷楊永鋒
航空科學(xué)技術(shù) 2022年7期
關(guān)鍵詞:方法

馮志杰,肖慧婷,楊永鋒

1.航空工業(yè)航宇救生裝備有限公司航空防護(hù)救生技術(shù)航空科技重點(diǎn)實(shí)驗(yàn)室,湖北 襄陽(yáng) 441000

2.西北工業(yè)大學(xué),陜西 西安 710072

機(jī)構(gòu)在使用過(guò)程中會(huì)受到各種因素的影響,其中振動(dòng)環(huán)境會(huì)影響機(jī)構(gòu)關(guān)鍵部件的使用壽命,甚至影響整個(gè)機(jī)構(gòu)的安全性及可靠性。而載荷譜作為機(jī)構(gòu)關(guān)鍵部件疲勞壽命及環(huán)境可靠性研究的重要組成部分,能否對(duì)其進(jìn)行正確的數(shù)據(jù)處理與編制,對(duì)于保障整個(gè)機(jī)構(gòu)的安全性及可靠性具有十分重要的意義。

李雷[1]等建立了某型彈射座椅的有限元計(jì)算模型,提出了飛機(jī)彈射座椅振動(dòng)響應(yīng)的有限元計(jì)算方法。基于該方法,馮志杰[2]等對(duì)某型飛機(jī)彈射座椅椅載設(shè)備進(jìn)行了隨機(jī)振動(dòng)仿真研究,獲取了椅載設(shè)備兩個(gè)不同測(cè)點(diǎn)處的振動(dòng)載荷譜,并對(duì)其進(jìn)行數(shù)據(jù)歸納。當(dāng)選取少量測(cè)點(diǎn)對(duì)振動(dòng)數(shù)據(jù)進(jìn)行歸納時(shí),出現(xiàn)偶然誤差的概率較大,可能導(dǎo)致得到的載荷譜過(guò)考核或欠考核,因此需對(duì)部件多個(gè)測(cè)點(diǎn)響應(yīng)進(jìn)行歸納分析。豐志強(qiáng)[3]等針對(duì)飛機(jī)機(jī)載設(shè)備在振動(dòng)環(huán)境可靠性試驗(yàn)剖面的編制需要,依據(jù)標(biāo)準(zhǔn)對(duì)機(jī)載設(shè)備振動(dòng)試驗(yàn)數(shù)據(jù)進(jìn)行了歸納。田永衛(wèi)[4]等對(duì)某型飛機(jī)多個(gè)艙位和主要結(jié)構(gòu)部位進(jìn)行了振動(dòng)量值的飛行實(shí)測(cè),并對(duì)多個(gè)不同測(cè)點(diǎn)的試驗(yàn)數(shù)據(jù)采用參數(shù)歸納法進(jìn)行了數(shù)據(jù)處理。依據(jù)GJB 150A的準(zhǔn)則,參數(shù)歸納法僅適用于總體近似服從正態(tài)分布的數(shù)據(jù)集[5],但由于測(cè)點(diǎn)位置的選取可能會(huì)影響響應(yīng)數(shù)據(jù)集概率分布情況,因此當(dāng)測(cè)得的數(shù)據(jù)集總體不完全服從正態(tài)分布時(shí),參數(shù)歸納法便不再適用,應(yīng)采用非參數(shù)歸納法處理不同測(cè)點(diǎn)數(shù)據(jù)。

本文對(duì)于機(jī)構(gòu)部件載荷譜獲取方法進(jìn)行了研究,首先對(duì)機(jī)構(gòu)進(jìn)行隨機(jī)振動(dòng)分析,然后提取機(jī)構(gòu)關(guān)鍵部件上不同測(cè)點(diǎn)的響應(yīng)數(shù)據(jù),再通過(guò)非參數(shù)上限統(tǒng)計(jì)歸納方法處理響應(yīng)數(shù)據(jù)。基于有效性和無(wú)偏性分析,評(píng)價(jià)了不同歸納方法的載荷譜歸納結(jié)果,最終得出較優(yōu)的機(jī)構(gòu)關(guān)鍵部件載荷譜歸納結(jié)果,為后續(xù)關(guān)于部件的可靠性及壽命的研究提供幫助。

1 機(jī)構(gòu)部件動(dòng)力學(xué)響應(yīng)獲取方法研究

若直接將機(jī)構(gòu)安裝位置載荷譜作為關(guān)鍵部件安裝位置處的載荷譜對(duì)部件進(jìn)行動(dòng)力學(xué)研究[6-7],則不能考慮從激勵(lì)輸入點(diǎn)到關(guān)鍵部件之間的傳力路徑,因此提出更符合實(shí)際的機(jī)構(gòu)關(guān)鍵部件載荷譜獲取方法[8-12]。

將機(jī)構(gòu)安裝位置處的載荷譜作為機(jī)構(gòu)載荷譜,對(duì)機(jī)構(gòu)進(jìn)行隨機(jī)振動(dòng)仿真分析,提取關(guān)鍵部件安裝位置附近不同測(cè)點(diǎn)處的響應(yīng),對(duì)不同測(cè)點(diǎn)處的響應(yīng)進(jìn)行歸納包絡(luò),將歸納后的曲線作為關(guān)鍵部件的載荷譜,以探究關(guān)鍵部件在振動(dòng)環(huán)境下的動(dòng)力學(xué)響應(yīng)及壽命等。機(jī)構(gòu)關(guān)鍵部件載荷譜獲取過(guò)程如圖1所示。

圖1 機(jī)構(gòu)關(guān)鍵部件載荷譜獲取流程Fig.1 Load spectrum acquisition process of key components of mechanism

隨機(jī)振動(dòng)的計(jì)算過(guò)程包括兩個(gè)階段[13]:(1)計(jì)算系統(tǒng)響應(yīng)特性,在時(shí)域內(nèi)用脈沖響應(yīng)函數(shù)h(t)表示,在頻域內(nèi)用復(fù)頻響函數(shù)H(ω)描述,其物理意義為系統(tǒng)響應(yīng)與激勵(lì)之比;(2)通過(guò)系統(tǒng)響應(yīng)特性和隨機(jī)激勵(lì)計(jì)算功率譜密度矩陣Sxx(ω)。

通過(guò)隨機(jī)振動(dòng)分析得到關(guān)鍵部件不同測(cè)點(diǎn)的響應(yīng)后,對(duì)結(jié)果進(jìn)行歸納。以往在對(duì)機(jī)構(gòu)關(guān)鍵部件響應(yīng)數(shù)據(jù)進(jìn)行處理時(shí),直接用仿真或?qū)崪y(cè)的某一個(gè)點(diǎn)或某兩個(gè)點(diǎn)的包絡(luò)作為關(guān)鍵部件的載荷譜[1-2]。取點(diǎn)數(shù)目較少時(shí),得到的部件載荷譜不夠準(zhǔn)確,因此需對(duì)關(guān)鍵部件多個(gè)測(cè)點(diǎn)響應(yīng)進(jìn)行歸納研究。

2 振動(dòng)數(shù)據(jù)歸納方法研究

GJB 150A 的18A[5]中規(guī)定了5 種數(shù)據(jù)集上限的統(tǒng)計(jì)分析方法,包含兩種參數(shù)上限統(tǒng)計(jì)估計(jì)方法和三種非參數(shù)上限統(tǒng)計(jì)估計(jì)方法。兩種參數(shù)上限統(tǒng)計(jì)估計(jì)方法分別是正態(tài)單邊容差上限法(NTL)和正態(tài)預(yù)測(cè)上限法(NPL);三種非參數(shù)上限統(tǒng)計(jì)估計(jì)方法分別是包絡(luò)上限法(ENV)、無(wú)驗(yàn)前分布容差上限法(DFL)和經(jīng)驗(yàn)容差上限法(ETL)。參數(shù)上限統(tǒng)計(jì)法適用于滿足正態(tài)分布或者經(jīng)轉(zhuǎn)化后滿足正態(tài)分布的數(shù)據(jù)集,而非參數(shù)上限統(tǒng)計(jì)法適用于不滿足正態(tài)分布或者轉(zhuǎn)化后不滿足正態(tài)分布的數(shù)據(jù)集。在對(duì)部件不同測(cè)點(diǎn)響應(yīng)進(jìn)行數(shù)據(jù)歸納時(shí),測(cè)點(diǎn)位置的選取會(huì)影響數(shù)據(jù)集分布,因此采用非參數(shù)上限統(tǒng)計(jì)方法歸納部件不同測(cè)點(diǎn)數(shù)據(jù)集。

2.1 ENV法

ENV 法的計(jì)算方法為選取數(shù)據(jù)集中最大估計(jì)值作為最大上限,以譜曲線一個(gè)頻段上的譜線作為研究對(duì)象,所有樣本在該段譜線上的容差上限為

式中:N表示數(shù)據(jù)集的樣本容量,xi表示第i個(gè)樣本的譜線值。

該方法計(jì)算簡(jiǎn)單,忽略了樣本的分布特性,無(wú)法給出超過(guò)該值的概率。當(dāng)估計(jì)集有異常的情況下,包絡(luò)上限法可能過(guò)于保守,對(duì)于譜線帶寬也很敏感。

2.2 DFL法

DFL法計(jì)算方法同ENV法一致,取數(shù)據(jù)集的最大值作為上限估計(jì)值,但不同于ENV 法的是,DFL 法引入了置信度和覆蓋率。

式中:β表示覆蓋率,γ表示置信度。用xβ代表數(shù)據(jù)集樣本的真實(shí)上限,則數(shù)據(jù)集樣本小于xβ的概率為β,即

若要使樣本中的最大值xmax小于xβ,則數(shù)據(jù)集中每一個(gè)樣本都應(yīng)小于xβ。每個(gè)樣本之間相互獨(dú)立,故有

xmax≥xβ的概率為

工程中通常使用的置信度至少為50%,覆蓋率為95%[11]。若要使估計(jì)上限置信度γ≥50%且β= 95%,由式(5)計(jì)算得N≥14;若要保持置信度大小不變,取覆蓋率較大時(shí)(β>0.95),計(jì)算所需的樣本數(shù)量較大。

2.3 ETL法

經(jīng)驗(yàn)容差上限一般表示為ETL(β)。假設(shè)共有N個(gè)測(cè)量點(diǎn),每個(gè)測(cè)量點(diǎn)的輸出響應(yīng)譜曲線有M段譜線,則總數(shù)據(jù)集中共包含NM個(gè)樣本

式中:j= 1,2,…,M,i= 1,2,…,N,mj表示第j段譜線的樣本均值。在M段譜線上構(gòu)造的正則化估計(jì)集為

將正則化估計(jì)集{u}中所有數(shù)據(jù)按從升序排列,數(shù)據(jù)集中第k個(gè)元素以u(píng)(k)來(lái)表示。構(gòu)造{u}的經(jīng)驗(yàn)分布函數(shù)FN

取覆蓋率為β時(shí),有

采用經(jīng)驗(yàn)分布函數(shù)FN代替實(shí)際分布函數(shù)F,取覆蓋率為β,由式(9)計(jì)算得k值,用u(k)作為容差上限uβ的估計(jì)值。每段譜線上的經(jīng)驗(yàn)容差上限為

取置信度為50%,此時(shí)xβj表示在第j段譜線上,EBL(β)以50%的置信度超過(guò)所有值的100β%。若選擇大于xβj的值,置信度會(huì)增加。當(dāng)樣本容量N達(dá)到一定要求時(shí),不同頻率分辨率帶寬上的正則化估計(jì)集應(yīng)服從同一分布。

ETL法是根據(jù)整個(gè)頻帶分布計(jì)算各頻率段譜線上的容差上限,因此可以消除一些偶然因素帶來(lái)的誤差,可以得到相對(duì)穩(wěn)定的結(jié)果,但其樣本量通常應(yīng)大于10[10]。

3 上限估計(jì)方法優(yōu)劣性評(píng)價(jià)

不同方法歸納得到的載荷譜必然存在優(yōu)良之分,故需要評(píng)價(jià)哪個(gè)估計(jì)方法相對(duì)較好。評(píng)價(jià)標(biāo)準(zhǔn)主要有無(wú)偏性評(píng)價(jià)和有效性評(píng)價(jià)。(1)無(wú)偏性評(píng)價(jià):對(duì)于真實(shí)容差上限θ,其估計(jì)量?滿足E=θ,則稱估計(jì)量?為真實(shí)容差上限θ的無(wú)偏估計(jì)量。(2)有效性評(píng)價(jià):設(shè)總體X~F(x;θ),對(duì)于真實(shí)容差上限θ的兩個(gè)無(wú)偏估計(jì)量?,,如果,則比有效。

本文通過(guò)對(duì)不同歸納方法得到的載荷譜結(jié)果的有效性和無(wú)偏性評(píng)價(jià),對(duì)比分析歸納結(jié)果的優(yōu)良性。取總樣本容量為N,將每段譜線樣本按樣本觀測(cè)值升序排列,取第95%N個(gè)樣本作為總體覆蓋率為95%的理論上限。基于無(wú)偏性和有效性,評(píng)價(jià)各方法得到的容差上限估計(jì)量的優(yōu)劣性。當(dāng)各容差上限估計(jì)量均值越接近理論上限,方差越小,則說(shuō)明該方法估計(jì)得到的載荷譜越接近真實(shí)工況。

4 實(shí)例分析

以某型彈射座椅為例,通過(guò)有限元分析獲得仿真結(jié)果,將椅載關(guān)鍵設(shè)備安裝位置不同測(cè)點(diǎn)處的響應(yīng)結(jié)果進(jìn)行載荷譜歸納。能準(zhǔn)確獲得人椅系統(tǒng)彈射座艙的各種參數(shù)是提高救生性能的關(guān)鍵[12],但在緊急狀況下彈射座椅幫助飛行員順利出艙的首要前提是保證彈射座椅關(guān)鍵設(shè)備在飛機(jī)正常飛行過(guò)程中的安全性和可靠性。采用程序控制器作為彈射座椅關(guān)鍵部件,討論不同數(shù)據(jù)歸納方法得到的載荷譜歸納結(jié)果的優(yōu)良性。

4.1 模型建立及響應(yīng)計(jì)算

以保證主結(jié)構(gòu)完整性、傳力路徑不變?yōu)樵瓌t,對(duì)某型飛機(jī)彈射座椅模型零部件進(jìn)行簡(jiǎn)化。該型飛機(jī)彈射座椅主要材料為鋁合金,其各項(xiàng)材料參數(shù)見表1。

表1 材料參數(shù)Table 1 Material parameters

座椅激勵(lì)譜是通過(guò)加速度傳感器實(shí)測(cè)得到,再根據(jù)GJB 150[5]中相關(guān)規(guī)定對(duì)試驗(yàn)數(shù)據(jù)歸納處理得到。以某型彈射座椅振動(dòng)耐久垂向載荷譜為例進(jìn)行計(jì)算分析,載荷曲線如圖2所示。

圖2 飛機(jī)座艙振動(dòng)耐久垂向加速度載荷譜Fig.2 Aircraft cockpit vibration endurance vertical acceleration load spectrum

完成有限元模型處理后,先對(duì)模型進(jìn)行頻率響應(yīng)分析獲取頻響函數(shù),再對(duì)其進(jìn)行隨機(jī)響應(yīng)分析,獲取測(cè)點(diǎn)處振動(dòng)響應(yīng)。在程序控制器安裝位置附近各取900個(gè)節(jié)點(diǎn)作為響應(yīng)輸出測(cè)點(diǎn),座椅簡(jiǎn)化模型及程序控制器安裝位置如圖3所示。

圖3 彈射座椅簡(jiǎn)化模型及程控器安裝位置Fig.3 Ejection seat simplified model and installation position of program controller

4.2 正態(tài)性檢驗(yàn)

當(dāng)譜密度曲線樣本的每個(gè)頻率分辨率帶寬數(shù)據(jù)集服從正態(tài)分布時(shí),加速度均方根值(RMS)也近似服從正態(tài)分布[10]。各測(cè)試點(diǎn)的加速度RMS值可由式(12)計(jì)算

式中:Gj和Δfj分別為第j個(gè)譜線上的譜密度和頻率帶寬。

加速度RMS值的分布特性可以反映信號(hào)的采樣點(diǎn)數(shù)、平均值和方差,因此可以通過(guò)對(duì)不同測(cè)點(diǎn)RMS值的分布研究譜密度曲線樣本是否服從正態(tài)分布。分別繪制程控器安裝位置不同測(cè)點(diǎn)處振動(dòng)響應(yīng)的加速度均方根值Q-Q 圖和直方圖,如圖4和圖5所示。

圖4 程控器各測(cè)加速度點(diǎn)RMS值Q-Q圖Fig.4 Quantile-Quantile diagram of RMS value of each accel eration point measured by the program controller

圖5 程控器各測(cè)點(diǎn)加速度RMS值直方圖Fig.5 Histogram of acceleration RMS value of each measuring point of the program controlle

通過(guò)程控器安裝位置附近不同測(cè)點(diǎn)處的加速度RMS值Q-Q圖和直方圖,可判斷該譜密度曲線集不完全滿足正態(tài)分布,因此需對(duì)其進(jìn)行非參數(shù)上限統(tǒng)計(jì)歸納。

4.3 歸納結(jié)果

對(duì)于ENV、DFL 法,計(jì)算方法較為簡(jiǎn)單,按照第2 節(jié)所述連接各段譜線上最大值,最終得到載荷歸納譜;ETL方法需將各個(gè)測(cè)試點(diǎn)的所有頻段上的譜線數(shù)據(jù)提取出來(lái)作為估計(jì)集,依據(jù)式(6)將所有樣本值進(jìn)行正則化,得到正則化估計(jì)集,對(duì)該正則化估計(jì)集取覆蓋率為β= 95%,利用式(9)~式(11)計(jì)算得到每段譜線寬上的經(jīng)驗(yàn)容差上限值,終得到95%覆蓋率下的載荷歸納譜。采用ENV、DFL和ETL方法繪制得到的歸納曲線如圖6~圖8所示。

圖6 ENV法歸納譜曲線Fig.6 The spectrum curve induced by ENV method

圖7 DFL法歸納譜曲線Fig.7 The spectrum curve induced by DFL method

圖8 ETL法歸納譜曲線Fig.8 The spectrum curve induced by ETL method

圖6和圖7中的載荷歸納譜一樣,原因是ENV法和DFL法計(jì)算容差上限估計(jì)值的方法是一樣的。兩種方法的區(qū)別在于DFL法引入了參數(shù)置信度和覆蓋率,可以通過(guò)樣本容量和其中一個(gè)參數(shù)計(jì)算另一個(gè)參數(shù)。在本文后面討論中將ENV 法和DFL 法歸并討論。對(duì)比ENV/DFL 法和ETL 法的歸納結(jié)果(見圖9)可看出,由于ETL法考慮了整個(gè)頻率范圍內(nèi)的分布,消除了一部分由取點(diǎn)位置隨機(jī)性等因素產(chǎn)生的譜線值波動(dòng),導(dǎo)致其在峰值點(diǎn)的上限估計(jì)值均小于ENV/DFL法峰值點(diǎn)的上限估計(jì)值。為了驗(yàn)證計(jì)算的準(zhǔn)確性,對(duì)比參考文獻(xiàn)[2]中的試驗(yàn)結(jié)果與ETL 法歸納結(jié)果(見圖10)可看出,試驗(yàn)結(jié)果與歸納結(jié)果幅值及趨勢(shì)基本一致。計(jì)算得歸納譜加速度RMS=4.8844g,試驗(yàn)測(cè)得加速度RMS=5.0893g,計(jì)算得誤差為4.02%。

圖9 ENV/DFL法和ETL法歸納結(jié)果對(duì)比Fig.9 Comparison of inductive results between ENV/DFL and ETL

圖10 ETL法歸納結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.10 Comparison between ETL inductive results and experimental results

4.4 歸納結(jié)果優(yōu)劣性評(píng)價(jià)

對(duì)三種方法的歸納結(jié)果進(jìn)行無(wú)偏性和有效性評(píng)價(jià)。取總體樣本容量N為900,按升序排列各譜線值,取第855 個(gè)數(shù)據(jù)(95%覆蓋率)作為理論上限值。樣本容量分別設(shè)定為5、13、30、45、60、75 和90,每個(gè)容量下重復(fù)抽樣10 次,計(jì)算10次抽樣的載荷譜歸納結(jié)果的均值和方差。載荷譜線的加速度RMS值能反映信號(hào)的平均能量,峰值點(diǎn)的波動(dòng)大小能反映載荷對(duì)機(jī)構(gòu)部件的破壞性強(qiáng)弱,因此取載荷譜線的加速度RMS 值和峰值f= 62.78Hz 處的歸納結(jié)果作為討論對(duì)象進(jìn)行優(yōu)良性評(píng)價(jià)。首先分別計(jì)算各樣本容量10 次抽樣歸納結(jié)果的均值和方差,然后將計(jì)算所得均值除以理論上限進(jìn)行歸一化處理。計(jì)算得到的結(jié)果見表2~表5,對(duì)比不同歸納方法得到的峰值點(diǎn)和加速度RMS 值的歸一化均值及方差可以發(fā)現(xiàn):(1)無(wú)論是對(duì)于峰值點(diǎn)還是加速度RMS值,ENV/DFL 法得到的歸納結(jié)果隨著樣本容量的增加,歸一化平均值不斷增加,且方差波動(dòng)幅度較大,無(wú)法趨于穩(wěn)定;(2)ETL 法由于考慮了各測(cè)點(diǎn)每段譜線整體分布,其歸納結(jié)果消除了一部分波動(dòng)。所以隨樣本量增加,其歸納結(jié)果的歸一化均值變化幅度較小,方差逐漸趨于穩(wěn)定。

表2 不同抽樣次數(shù)下各歸納方法在f=62.78Hz 處的歸一化上限平均值Table 2 The normalized upper bound mean of each inductive method at f=62.78Hz under different sampling times

表3 不同抽樣次數(shù)下各歸納方法在f=62.78Hz處的上限值方差((g2/Hz)2?10?6)Table 3 The upper bound variance of each inductive method at f=62.78Hz under different sampling times

表4 不同抽樣次數(shù)下各歸納方法加速度RMS值的歸一化上限平均值Table 4 The normalized upper bound mean of each inductive method of acceleration RMS under different sampling times

表5 不同抽樣次數(shù)下各歸納方法加速度RMS值的上限值方差((g2/Hz)2?10?4)Table 5 The upper bound variance of each inductive method of acceleration RMS under different sampling times

綜上所述,ETL 方法得到的載荷譜歸納結(jié)果的無(wú)偏性和有效性均優(yōu)于ENV/DFL方法得到的歸納結(jié)果。

5 結(jié)束語(yǔ)

本文建立了一種機(jī)構(gòu)關(guān)鍵部件載荷譜獲取方法。首先對(duì)機(jī)構(gòu)進(jìn)行隨機(jī)振動(dòng)分析,獲取關(guān)鍵部件測(cè)點(diǎn)響應(yīng),然后對(duì)響應(yīng)數(shù)據(jù)集進(jìn)行歸納處理,最終得出部件的載荷譜。通過(guò)對(duì)比不同歸納方法的理論分析以及實(shí)例結(jié)果,得出以下結(jié)論:當(dāng)樣本容量N≤13時(shí)采用ENV/DFL方法得到的估計(jì)結(jié)果波動(dòng)較大,故采用ENV/DFL 法時(shí)樣本量不宜過(guò)小;ETL方法相對(duì)ENV/DFL方法得到的估計(jì)結(jié)果更穩(wěn)定,但當(dāng)樣本量較大時(shí)ETL法對(duì)于載荷譜波動(dòng)較小的頻段上限估計(jì)值較大,而波動(dòng)較大的峰值附近其上限估計(jì)值較小,因此采用ETL法時(shí)樣本量不宜過(guò)大。

猜你喜歡
方法
中醫(yī)特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數(shù)學(xué)教學(xué)改革的方法
化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學(xué)習(xí)方法
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡(jiǎn)單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲美女一区二区三区| 亚洲国产精品VA在线看黑人| 九九这里只有精品视频| 国产一区二区三区在线观看视频 | 麻豆精品视频在线原创| 国产精品密蕾丝视频| 色亚洲激情综合精品无码视频 | 国产拍揄自揄精品视频网站| 欧美日韩中文国产| 狠狠综合久久久久综| 亚洲国产午夜精华无码福利| 动漫精品中文字幕无码| 91在线视频福利| 国产亚洲精品无码专| 亚洲成在人线av品善网好看| 免费女人18毛片a级毛片视频| 亚洲无码高清视频在线观看| 成人另类稀缺在线观看| 国产成人高清亚洲一区久久| 狼友视频国产精品首页| 亚洲国产欧美中日韩成人综合视频| 中国黄色一级视频| 亚洲精品第一在线观看视频| 亚洲视频在线青青| 亚洲人成人无码www| 国产成人午夜福利免费无码r| 无码日韩视频| 亚洲精品欧美重口| 亚洲欧洲日韩久久狠狠爱| 欧美日本激情| 国产香蕉国产精品偷在线观看| 伊人久久福利中文字幕| 精品久久久久成人码免费动漫| 露脸真实国语乱在线观看| 国产嫖妓91东北老熟女久久一| 久久久久久高潮白浆| 成人韩免费网站| 国产精品视频免费网站| 一级毛片在线播放免费| A级毛片高清免费视频就| 久久综合伊人 六十路| 91久久青青草原精品国产| 国产午夜在线观看视频| 中文字幕日韩视频欧美一区| 午夜一区二区三区| 国产成人综合亚洲欧洲色就色| 国产成a人片在线播放| 国产福利一区视频| 免费啪啪网址| 亚洲日本精品一区二区| 亚洲欧洲国产成人综合不卡| 欧美一级高清免费a| 成年午夜精品久久精品| 成人综合网址| 高清乱码精品福利在线视频| 亚洲国产午夜精华无码福利| 亚洲午夜片| 亚洲美女一区| 亚洲激情区| 自慰高潮喷白浆在线观看| 欧美另类视频一区二区三区| 99视频在线看| 国产一级毛片yw| a免费毛片在线播放| 久久6免费视频| 欧美激情一区二区三区成人| 五月激情综合网| 天天色综合4| 成人韩免费网站| 国产亚洲美日韩AV中文字幕无码成人 | 色天堂无毒不卡| 中文字幕中文字字幕码一二区| 手机看片1024久久精品你懂的| 国产成人综合亚洲欧洲色就色| 日韩小视频在线观看| 亚洲日产2021三区在线| 日韩欧美国产精品| 毛片视频网| 视频二区欧美| 日韩欧美国产精品| 亚洲中文字幕无码mv| 亚洲国产成人综合精品2020 |