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

基于Morris靈敏度分析的彈性-黏彈性層合板阻尼性能代理模型構(gòu)建

2021-08-21 03:05:14吳云鵬
噪聲與振動控制 2021年4期
關(guān)鍵詞:模態(tài)分析模型

吳云鵬,鄭 輝

(上海交通大學(xué)振動、沖擊、噪聲研究所,上海200240)

黏彈性材料由于同時表現(xiàn)出類似固體和液體的性能,具有優(yōu)良的阻尼性能,被廣泛應(yīng)用于汽車、航空、航天等行業(yè)裝備的減振降噪設(shè)計[1]。但是黏彈性材料的力學(xué)性能尤其是剛度性能較差,因此實際工程應(yīng)用中通常會使用夾層結(jié)構(gòu),即在剛度較大的各向同性金屬層或各向異性復(fù)合材料層間嵌入黏彈性阻尼層,以強化結(jié)構(gòu)的力學(xué)性能。Kerwin 等[2]針對各向同性黏彈性材料夾層結(jié)構(gòu)建立了力學(xué)和聲學(xué)解析模型。為了進一步改善彈性-黏彈性層合結(jié)構(gòu)的力學(xué)或聲學(xué)性能,還需要進一步對阻尼層的內(nèi)部結(jié)構(gòu)進行設(shè)計,而本文中則通過對黏彈性層開孔并填充高強度材料的方式提高其力學(xué)性能,這種阻尼層內(nèi)部結(jié)構(gòu)的復(fù)雜性以及各項異性材料的存在,一方面使得彈性-黏彈性層合板的分析和優(yōu)化設(shè)計仍然以有限元分析為主[3–9],另一方面也使得這類結(jié)構(gòu)的優(yōu)化設(shè)計所涉及的參數(shù)數(shù)量進一步增加。

對于彈性-黏彈性夾層結(jié)構(gòu)的力學(xué)和聲學(xué)性能優(yōu)化,學(xué)者們已開展諸多研究。Teng 和Hu[3]通過研究黏彈性夾層梁阻尼性能發(fā)現(xiàn),阻尼層厚度對結(jié)構(gòu)整體的阻尼性能影響較大。而Arvin 等[4]針對各向異性黏彈性梁分析了其阻尼性能隨阻尼層厚度的變化規(guī)律。也有不少學(xué)者從參數(shù)靈敏度的角度出發(fā)研究黏彈性層合板的關(guān)鍵影響參數(shù),并對計算和優(yōu)化方法進行改進。Akoussan 等[5]結(jié)合靈敏度分析和基于梯度的優(yōu)化方法進行了黏彈性層合結(jié)構(gòu)的局部優(yōu)化設(shè)計,并通過算法優(yōu)化降低了運算量。但是這些工作大多圍繞局部優(yōu)化展開的。而通過全局優(yōu)化對黏彈性阻尼材料進行優(yōu)化設(shè)計的難點在于,全局優(yōu)化算法運算量遠高于局部優(yōu)化算法,并且運算量顯著受到待優(yōu)化問題維數(shù)的影響,隨著輸入變量數(shù)目的增加會迅速增長。為了克服全局優(yōu)化算法運算量過大的問題,一方面有學(xué)者針對全局優(yōu)化算法進行改進,如楊洋等[6]通過對遺傳算法進行優(yōu)化有效降低了收斂所需的代數(shù)。另一方面,關(guān)于使用代理模型對黏彈性材料的相關(guān)性能進行擬合與預(yù)測的研究也在不斷深入。臧獻國[7]使用響應(yīng)面代理模型結(jié)合模態(tài)響應(yīng)法對轎車車身結(jié)構(gòu)進行了優(yōu)化設(shè)計。而魏星等[8]則發(fā)現(xiàn)通過使用代理性較好的代理模型,可以使得有限元模型的調(diào)用次數(shù)降低至直接優(yōu)化的1/4 至1/6。但這些使用代理模型進行全局優(yōu)化設(shè)計的研究并非通過參數(shù)靈敏度分析篩選高靈敏度參數(shù)構(gòu)建代理模型,對于復(fù)雜結(jié)構(gòu)構(gòu)建的代理模型維數(shù)依然很大。

參數(shù)靈敏度分析對代理模型的優(yōu)化主要體現(xiàn)在兩方面,首先是通過降低問題維數(shù),可有效減少后續(xù)全局優(yōu)化的運算量并提高收斂性;其次在穩(wěn)健性分析中能夠評估結(jié)果的魯棒性。本文主要關(guān)注靈敏度分析的第一類作用,即降低問題維數(shù)的作用。

目前動力學(xué)和聲學(xué)靈敏度分析常用的方法為伴隨變量法,其本質(zhì)是基于偏導(dǎo)數(shù)計算的攝動分析法,該類方法的優(yōu)勢是運算量較小,缺點是僅能反映某個參數(shù)在局部范圍內(nèi)的靈敏度。而Morris方法能夠反映某個參數(shù)在全局范圍內(nèi)的靈敏度,缺點為運算量高于攝動法。

本文在嵌入開孔阻尼層的彈性-黏彈性層合板阻尼性能代理模型構(gòu)建中,通過靈敏度分析對模型進行降維簡化。首先分別使用攝動法和Morris方法對層合板的1 階模態(tài)阻尼比進行參數(shù)靈敏度分析,根據(jù)分析結(jié)果篩除相對不敏感參數(shù),之后使用多項式響應(yīng)面法對模型1 階模態(tài)阻尼比進行擬合,并用構(gòu)建的代理模型預(yù)測新采樣點的1階模態(tài)阻尼比并計算其均方根誤差,評價代理模型的代理性能,通過對比明確兩種靈敏度分析方法對構(gòu)建代理模型的簡化效果。

1 參數(shù)靈敏度分析

1.1 攝動法

攝動法的基礎(chǔ)是在感興趣的點(通常情況下是參數(shù)變化范圍的均值點,或是初始模型的參數(shù)值)周圍使用泰勒展開代替原函數(shù)關(guān)系,只需要做出1 階泰勒展開,就可以分析各個模型的均值與方差。模型的響應(yīng)值f(x) 在x0處的1階泰勒展開形式如下:

式中:x0i為展開點x0在第i個變量處的分量。利用泰勒展開式,可以對原函數(shù)的均值和方差進行估計,當展開點x0恰是輸入分布的均值時,有:

如果設(shè)計參數(shù)之間彼此獨立,則C轉(zhuǎn)化為對角矩陣,對角線上的元素為對應(yīng)設(shè)計參數(shù)輸入量的方差,此時響應(yīng)量的方差可以表示為各分量方差σi的線性組合:

因此可以通過響應(yīng)量在展開點附近的偏導(dǎo)數(shù)和設(shè)計參數(shù)的方差描述某個設(shè)計參數(shù)對響應(yīng)值的重要性。

在基于梯度運算的局部優(yōu)化算法中,這一類靈敏度分析方法被廣泛使用[4],原因是基于攝動法的靈敏度分析本身能夠提供局部優(yōu)化算法所需的梯度信息,同時在基于梯度的優(yōu)化算法中,隨著優(yōu)化過程進行,攝動點會發(fā)生改變,克服了攝動法只能反映局部靈敏度的劣勢。但在使用遺傳算法等優(yōu)化算法的全局優(yōu)化設(shè)計中,梯度信息變得不再必要,并且攝動法給出的靈敏度信息較為粗略。

1.2 Morris方法

由于不同輸入變量的取值范圍有較大差異,為便于描述,在本文中將采樣空間所有維度取值范圍放縮至[0,1]之間。基于軌跡的采樣策略的Morris方法所需的采樣次數(shù)為(M+1)×r次[10],其中M為靈敏度分析問題的維數(shù),r描述采樣過程中對每個維度的細分程度,采樣軌跡可以用維數(shù)為(M+1)×r的矩陣B′表示,該矩陣在每列中都包含一個從隨機點x?開始的點。下一個軌跡點通過將x?某個維度的坐標增加Δ獲得,其中:

式中:n為取值范圍從1到r-1的整數(shù)。

一個滿足上述條件的矩陣B′可通過式(5)計算:

式中:J為只包含元素1的矩陣,B為包含一個嚴格下三角矩陣,而其余部分為零的(M+1)×r階矩陣。起始點x?的必要條件為從該點沿任意維度增加Δ后仍在可行域G內(nèi)。

通過上述方法獲得的采樣點除起點外軌跡都是固定的,為了克服這一問題,Morris 提出了式(6)所示的更有隨機性的采樣方式[10]:

式中:矩陣D為M×M階對角矩陣,其對角線元素以相同的概率隨機取1或-1,矩陣P為M×M階單位矩陣的隨機置換矩陣。 給定一個輸入x={x1,x2,…,xN}∈G,并要求滿足xi+Δ<1,則xi對于函數(shù)值變化的貢獻可以通過式(7)描述。

在G內(nèi)隨機獲取多組采樣點,并對每個點都計算di(x),通過計算該組數(shù)據(jù)的平均值和標準差,可以分別從不同的角度評估某個輸入變量的重要性:平均值主要描述輸入變量xi自身的重要性;而標準差從另一方面描述變量的重要性,如果變量xi的評估指標di(x) 有著較高的標準差,則變量xi的重要性主要來自該變量的非線性行為或與其他輸入變量的相互作用。

由于存在非線性關(guān)系,一些變量的效果可能會在不同取值范圍內(nèi)相互抵消,為了避免這種情況出現(xiàn),Campolongo 等[11]在2007年引入了新評估方法,即使用函數(shù)di(x)絕對值|di(x)|的算術(shù)平均*作為新評估標準。新評估方法避免了某個參數(shù)在不同取值范圍內(nèi)由于效果符號改變造成互相抵消,但也忽略了有關(guān)輸入變量重要性的符號信息,需與先前的和組合使用,在不進行額外采樣的情況下可以分析每個輸入變量的實際重要性和這種貢獻性質(zhì)(線性、非線性或相互作用)的更完整信息。

相比于攝動法,Morris 方法的顯著優(yōu)勢在于能夠識別響應(yīng)值與某個輸入變量之間的非線性關(guān)系,并提供某個參數(shù)在一定取值范圍內(nèi)的整體貢獻。因此Morris方法對于存在非線性行為函數(shù)的靈敏度分析結(jié)果比攝動法更精準,如果在攝動點附近某個參數(shù)對目標函數(shù)的影響要大于或小于在其他區(qū)域內(nèi)的影響,則攝動法可能會高估或低估該參數(shù)的重要性,Morris方法則能綜合評價某個指標在取值范圍內(nèi)的重要性。此外如果模型存在明顯的非線性行為,或是某個參數(shù)和目標函數(shù)的相關(guān)性以線性關(guān)系為主,都可以通過Morris 分析進行區(qū)分,因此除了針對代理模型的元素數(shù)量進行簡化,通過Morris 靈敏度分析還可以針對代理模型中包含某個元素的線性項或非線性項進行簡化,這對于代理模型的建模和簡化過程更有指導(dǎo)作用。

1.3 參數(shù)靈敏度分析結(jié)果及討論

用有限元商業(yè)軟件COMSOL 計算參數(shù)靈敏度分析所需模型1階阻尼比。用于結(jié)果比較的模型為圖1 所示的阻尼層開孔的彈性-黏彈性層合板,層合板外形尺寸為200 mm×200 mm×1.2 mm,層合板的邊界條件為四邊固支。

圖1 阻尼層開孔的層合板示意圖

阻尼層料使用黏彈性材料SA-3C[9],其楊氏模量和阻尼損耗因子具有頻變特性,函數(shù)關(guān)系為

式中:E(f)為阻尼材料的楊氏模量,單位為MPa,f為頻率,單位Hz,β(f) 為材料的黏彈性損耗因子。由于彈性-黏彈性阻尼層合板的黏彈性阻尼材料層的楊氏模量和阻尼損耗因子具有頻變特性,需通過迭代方法[9]獲得夾層板1 階模態(tài)的特征頻率與阻尼損耗因子,迭代方式如圖2 所示。計算時設(shè)定特征頻率的最大容差為0.02 Hz。

圖2 阻尼比迭代計算流程

得到模型的1階模態(tài)特征頻率后,按式(9)計算模型1 階模態(tài)阻尼比,并將其作為彈性-黏彈性阻尼夾層板阻尼性能的評價指標。

式中:ζ1為模型的1 階模態(tài)阻尼比,Im(f1)為1 階特征頻率的虛部,|f1|為1階特征頻率的模。

將彈性-黏彈性層合板可在一定范圍內(nèi)變化的參數(shù)作為輸入變量,分別用攝動法和Morris 方法對這些變量進行參數(shù)靈敏度分析,采用攝動法計算時的采樣數(shù)量為8 組,采用Morris 方法時的采樣數(shù)量為24組。

實際層合板的參數(shù)變化來自兩方面,首先是通過尺寸設(shè)計人為控制的尺寸參數(shù),其次為加工精度導(dǎo)致的材料性能差異。可以人為控制的尺寸參數(shù)包括阻尼層厚度、阻尼層中面至層合板中面的偏移量、阻尼層開孔邊長和相鄰孔間距,其中,阻尼層厚度x1變化范圍為0.2 mm~0.24 mm,阻尼層中面到層合板中面的偏移量x2變化范圍為-0.2 mm~0.2 mm。阻尼層開孔邊長d變化范圍為20 mm~90 mm,相鄰兩孔的間距l(xiāng)變化范圍為88 mm~90 mm。

金屬層和開孔處的填充物為鋁合金,其楊氏模量E1在67.9 GPa~69.9 GPa 的范圍內(nèi)浮動,密度ρ1為2 740 kg/m3,泊松比v1為0.3。并且認為阻尼層材料的密度ρ2變化范圍為1 425 kg/m3~1 475 kg/m3,泊松比v2變化范圍為0.46~0.52。

首先設(shè)定所有可變參數(shù)的取值均為取值范圍的中值,代入有限元模型驗證迭代程序的可行性,經(jīng)計算,模型的1階特征頻率為210.67 Hz+16.830 Hz,依據(jù)式(8)獲得黏彈性材料的楊氏模量和阻尼損耗因子,并將新的材料參數(shù)代入有限元模型計算其1 階模態(tài),計算結(jié)果如圖3所示。

圖3 模型的1階模態(tài)有限元計算結(jié)果

通過對比可以發(fā)現(xiàn),最后一次迭代前后的1 階特征頻率已經(jīng)非常接近,能夠滿足最大容差要求,可以將該結(jié)果應(yīng)用于后續(xù)的靈敏度分析。

攝動法的靈敏度分析結(jié)果如圖4 所示。由圖4結(jié)果可以發(fā)現(xiàn),高靈敏度參數(shù)分別為阻尼層厚度x1、阻尼材料泊松比v2、開孔邊長d和開孔間距l(xiāng)。參數(shù)x1靈敏度較高的結(jié)論符合文獻[3]中關(guān)于阻尼層厚度能夠顯著影響阻尼性能的描述,而開孔邊長對阻尼性能影響較大的原因在于開孔的大小會直接影響阻尼層的質(zhì)量占比,從而對阻尼性能產(chǎn)生較大影響,因此也符合對計算結(jié)果的預(yù)期。

圖4 攝動法靈敏度分析結(jié)果

Morris 方法靈敏度分析結(jié)果如圖5 所示。與攝動法分析結(jié)果的主要不同在于:基于Morris 方法分析結(jié)果表明,阻尼層中面偏移量x2也有較高的靈敏度,而基于攝動法分析結(jié)果則表明該參數(shù)的靈敏度很低。通過進一步觀察基于Morris靈敏度分析結(jié)果的縱軸坐標可以發(fā)現(xiàn),參數(shù)x2的靈敏度主要來自于自身的非線性行為或與其他參數(shù)相互作用,而當使用函數(shù)絕對值|di(xi)|作為Morris 靈敏度分析的評估指標時,其參數(shù)靈敏度進一步提高,這表明x2可能存在對目標函數(shù)值的貢獻度相互抵消的情況。

圖5 Morris方法靈敏度分析結(jié)果

2 代理模型的代理性能比較

2.1 代理模型構(gòu)建

代理模型的主要工程意義在于,對于需要經(jīng)過高精度有限元分析或?qū)嶒灢拍塬@得的數(shù)據(jù),通過使用結(jié)果與有限元分析結(jié)果相近而所需運算量卻遠小于有限元計算的近似模型進行預(yù)測,可以在滿足一定精度要求的同時顯著減少運算量。

本文將要構(gòu)建的代理模型為多項式響應(yīng)面代理模型,是工程上較為常用的一種代理模型方法,其一般形式為

式中:y為采樣點獲取的數(shù)據(jù),β0、βi、βii、βij為待定系數(shù),xi、xj為輸入變量,ε為近似誤差。

通過最小二乘法可獲得待定系數(shù)表達式:

式中:β為待定系數(shù)矢量,X為設(shè)計變量矩陣。

首先根據(jù)本項目中的輸入變量的取值范圍進行采樣,獲得采樣點中輸入變量的取值后生成設(shè)計變量矩陣,并將采樣點處獲得的有限元分析結(jié)果記錄為矢量y,之后根據(jù)式(11)獲得待定系數(shù)矢量β,將結(jié)果代入式(10)即可獲得代理模型表達式。

2.2 代理模型的性能比較

為了評估基于Morris靈敏度分析和基于攝動法的靈敏度分析建立的代理模型的代理性能,分別使用根據(jù)Morris分析方法獲得的高靈敏度參數(shù)(x1、v2、d和l)和根據(jù)攝動法獲得的高靈敏度參數(shù)(x1、x2、v2、d和l)對層合板的一階模態(tài)阻尼比進行二項式響應(yīng)面建模,建模點數(shù)量均為40 組。獲得代理模型后,分別使用兩組代理模型預(yù)測20組新采樣點的1階模態(tài)阻尼比,將預(yù)測結(jié)果與有限元計算結(jié)果進行對比,并計算預(yù)測結(jié)果的均方根誤差,即:

其中:nt為預(yù)測采樣點的數(shù)目,此處取20。yi為待預(yù)測點的真實函數(shù)值,通過有限元計算獲得。為待預(yù)測點的預(yù)測值,通過代理模型計算獲得。

將分別用攝動法和Morris方法進行靈敏度分析構(gòu)建的代理模型的代理性能與不進行靈敏度分析、直接用全部參數(shù)進行代理模型建模的情況進行對比,結(jié)果如表1所示。

表1 代理模型的代理性能對比

結(jié)果表明,雖然通過攝動法篩選獲得的高靈敏度參數(shù)數(shù)量更少,模型的形式也更為簡單,但是模型的均方根誤差顯著高于另外兩種情況,這表示模型的代理性嚴重不足。而通過Morris方法構(gòu)建的代理模型卻能保持良好的代理性能,其均方根誤差低于無靈敏度分析直接建立代理模型的情況。考慮到基于攝動法和Morris方法靈敏度分析的主要區(qū)別在于兩者對于參數(shù)x2的重要性判斷不同,為了確認參數(shù)x2的重要性,現(xiàn)將所有輸入變量的數(shù)值設(shè)定為各自取值范圍的中值,隨后使x2在取值范圍內(nèi)由低到高逐步變化,并計算模型1 階模態(tài)阻尼比。之后針對攝動法分析結(jié)果中有較高靈敏度的參數(shù)d進行同樣的操作,將兩者對阻尼比變化的影響程度進行對比,計算結(jié)果如圖6所示。

圖6 不同參數(shù)對模型阻尼比的影響

不難發(fā)現(xiàn),通過調(diào)節(jié)阻尼層中面的偏移量,層合板的1 階模態(tài)阻尼比變化范圍為0.061 1 至0.079 6,變化幅度超過30%,而且阻尼比與x2的關(guān)系曲線為單峰曲線,符合x2存在對目標函數(shù)值的貢獻度相互抵消情形的預(yù)判。而通過調(diào)節(jié)阻尼層開孔的邊長實現(xiàn)的阻尼比變化范圍為0.077 4至0.081 3,變化幅度約為5%。這表明阻尼層中面的偏移量x2對層合板1階模態(tài)阻尼比的影響程度要遠高于阻尼層開孔邊長d。攝動法對于參數(shù)x2的靈敏度給出了錯誤的判斷,主要原因為:在攝動點附近,隨著參數(shù)x2的改變,模型1 階模態(tài)阻尼比的變化較為平緩,攝動法將攝動點處阻尼比的變化情況作為對整體范圍內(nèi)變化情況的預(yù)估,低估了x2的重要性。而Morris 方法則能夠?qū)γ總€維度取多個點進行全面考察,避免對參數(shù)靈敏度做出錯誤判斷。

3 結(jié)語

本文將基于Morris 方法與基于傳統(tǒng)振-聲學(xué)優(yōu)化中常用的攝動法靈敏度分析應(yīng)用于開孔阻尼層的彈性-黏彈性層合板1 階模態(tài)阻尼比代理模型的構(gòu)建,以降低模型維數(shù)和加速后續(xù)將進行的黏彈性阻尼層參數(shù)優(yōu)化進程。以一彈性-黏彈性-彈性阻尼夾層板為例,對分別用這兩種靈敏度分析方法辨識高靈敏度參數(shù)后構(gòu)建的阻尼代理模型與不進行靈敏度分析、直接用全部參數(shù)構(gòu)建的代理模型進行了代理性能對比。結(jié)果表明:雖然通過攝動法獲得的高靈敏度參數(shù)數(shù)量更少,代理模型的形式也更為簡單,但是模型的新采樣點預(yù)測結(jié)果均方根誤差顯著高于另外兩種情況,表明模型的代理性能嚴重不足;而通過Morris方法構(gòu)建的代理模型卻能保持良好的代理性能,其新采樣點預(yù)測結(jié)果均方根誤差甚至低于無靈敏度分析直接建立代理模型的情況。研究結(jié)果可為嵌入式阻尼金屬或復(fù)合材料層合板的高阻尼結(jié)構(gòu)參數(shù)優(yōu)化代理模型構(gòu)建提供依據(jù)。

需要說明的是,雖然Morris 方法相較于攝動法不易對參數(shù)的重要性做出誤判,但運算量高于攝動法,本文采用Morris 方法進行參數(shù)敏感性分析時的采樣點數(shù)量達到攝動法的3倍。但是在阻尼優(yōu)化設(shè)計等問題中,需要直接使用代理模型代替有限元模型輸出結(jié)果,如果參數(shù)靈敏度分析結(jié)果出現(xiàn)錯誤,將會使得降維簡化后的代理模型失去代理性,進而使得整個優(yōu)化過程的結(jié)果失去意義,因此在代理模型構(gòu)建階段使用精度更高的靈敏度分析方法仍是有必要的。

猜你喜歡
模態(tài)分析模型
一半模型
隱蔽失效適航要求符合性驗證分析
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動化發(fā)展趨勢分析
3D打印中的模型分割與打包
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
由單個模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 四虎免费视频网站| 精品成人免费自拍视频| 色综合久久久久8天国| 亚洲精品波多野结衣| 国产理论最新国产精品视频| 久久精品国产国语对白| 一级毛片a女人刺激视频免费| 国产麻豆91网在线看| 国产丝袜无码精品| 久久香蕉国产线看观看精品蕉| 精品国产污污免费网站| 毛片三级在线观看| 国产精品久久久久久久久kt| 97se亚洲| 九九这里只有精品视频| 手机看片1024久久精品你懂的| 中文字幕va| 无码内射在线| 免费人成在线观看成人片| 久草热视频在线| 国产午夜精品一区二区三| 亚洲国产成人精品一二区| 四虎永久在线视频| 92午夜福利影院一区二区三区| 国产高清自拍视频| 亚洲天堂日韩av电影| 免费看黄片一区二区三区| 永久免费无码成人网站| 亚洲不卡影院| 久久人午夜亚洲精品无码区| 久久精品中文无码资源站| 色综合天天娱乐综合网| 亚洲aaa视频| 成AV人片一区二区三区久久| 欧美精品亚洲精品日韩专区| 毛片一级在线| 天堂在线亚洲| 丁香五月亚洲综合在线| 少妇精品在线| 狠狠亚洲婷婷综合色香| 国产又爽又黄无遮挡免费观看| 国产精品综合久久久| 高潮毛片免费观看| 精品国产Av电影无码久久久| 奇米影视狠狠精品7777| 超碰aⅴ人人做人人爽欧美| 亚洲色图另类| 老司机精品一区在线视频| 国产精品9| 91青草视频| 中文字幕在线一区二区在线| 久久a毛片| 国产成人精品午夜视频'| 日韩麻豆小视频| 精品国产美女福到在线不卡f| 精品国产一区91在线| 女人爽到高潮免费视频大全| 五月婷婷中文字幕| 国产av剧情无码精品色午夜| 九九久久精品国产av片囯产区| 97精品伊人久久大香线蕉| 91视频国产高清| 亚洲国产综合自在线另类| …亚洲 欧洲 另类 春色| 成人免费视频一区二区三区 | 国产第八页| 日韩av在线直播| 97亚洲色综久久精品| 欧美97欧美综合色伦图 | 欧美国产综合色视频| 久久久精品国产SM调教网站| 欧美在线导航| 国产精品视频3p| 在线免费不卡视频| 成人一级免费视频| 国产精品人人做人人爽人人添| 亚洲日韩AV无码一区二区三区人| 成人另类稀缺在线观看| 91精品小视频| 中文字幕免费播放| 欧美一区福利| 综合网天天|