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

基于超梁降階模型的復雜梁式結構動力學模型修正及確認

2017-11-06 02:23:51王文勝魏豪杰侯中華李一帆
固體火箭技術 2017年5期
關鍵詞:有限元優化結構

王文勝,魏豪杰,侯中華,梅 群,李一帆

(河南科技大學 土木工程學院工程力學系,洛陽 471023)

2016-11-01;

2016-12-08。

國家自然科學基金項目(11402077)。

王文勝(1983—),男,博士,主要從事復雜結構動力模型降階及結構優化設計。E-mailwswang@live.cn

基于超梁降階模型的復雜梁式結構動力學模型修正及確認

王文勝,魏豪杰,侯中華,梅 群,李一帆

(河南科技大學 土木工程學院工程力學系,洛陽 471023)

提出了一種基于超梁降階模型的復雜梁式結構動力學模型修正技術。依據超梁降階模型理論,并考慮可能的修正參數,將復雜梁式結構降階為具有較高計算精度的超梁降階模型;而后,利用基于靈敏度分析及優化算法的模型修正技術,并根據模態測試數據對模型參數進行修正,獲得能夠準確反映實際結構動力學特征的修正模型。最后,通過多種動力問題分析對修正模型進行評估,并與實際結構的響應結果進行對比。以典型蒙皮加筋圓柱殼為例實現這一過程。結果表明,所提出的模型修正方法具有可行性和較高的計算效率。

超梁降階模型;梁式結構;動力學;模型修正;優化

0 引言

隨著計算機技術和科學的快速發展,過去的幾十年里有限元理論和方法取得了巨大的成功,被廣泛應用在航空航天、交通和建筑結構等工程領域,并發揮著越來越重要的作用。有限元模擬中最重要的問題是其求解精度問題,在模擬實際工程結構時,有限元方法對實際結構進行離散化,有時不能準確反映實際結構的真實力學特性,分析結果與試驗測量值之間的誤差較大。如賀爾銘等[1]研究發現某機翼有限元模型固有頻率分析結果和試驗測量數值之間的差異多數超過10%,甚至某些彎曲頻率的誤差能夠達到70%。馬雙超等[2]發現某航空發動機機匣的頻率分析結果與試驗值誤差均在60%以上。

Mottershead等[3]認為有限元模型計算誤差的原因主要有:(1)在有限元建模過程中,因為線性化假設以及邊界條件的近似所引起的有限元模型的結構誤差;(2)有限元方法將實際連續無限維結構劃分為具有有限自由度的模型所造成的階次誤差,不能模擬具有無限自由度的真實結構;(3)由于對初始模型不精確的簡化、近似以及加工制作等原因造成結構材料和幾何參數的不確定所引起的模型參數誤差。建立一個能夠準確反映實際結構力學特性的有限元模型,是進行結構靜動力響應計算、優化設計、損傷識別和健康檢測等的重要前提條件。因此,模型修正技術越來越受到研究人員的重視[4-6]。

模型修正技術興起于20世紀70年代,其目的是建立一個能夠重現試驗測得的所有模態參數或者頻率響應函數的有限元模型。常用到的模型修正方法主要有矩陣型修正法和參數型修正法[7-8]。前者直接對結構的剛度和質量矩陣進行修改,修正速度快,但修正丟失了原結構總剛度陣和質量陣的帶狀稀疏分布特征,無明顯物理意義。后者以結構的物理參數(如密度、彈性模量和截面積等)作為修正對象,保留了矩陣的帶狀分布特征,具有明確的物理意義,逐漸成為模型修正的主要方法[9]。

參數型模型修正中,一種較常見的方法是利用優化技術將模型修正問題轉化為優化問題來處理。將結構的動力響應(試驗測得)等作為優化目標或者約束條件,通過優化方法獲取設計變量(修正參數)的最優取值,經過多次迭代,可使有限元模型與實際結構的動力特性差異最小化。該方法近幾年逐漸成為研究熱點[10-12]。在優化修正過程中,每次迭代都會涉及到結構特征值求解和靈敏度分析[13]。對于實際工程結構,由于其結構復雜性,所建立的有限元模型通常規模比較大,具有幾十萬甚至上千萬的自由度,并包含很多不確定參數需要修正,結構特征值求解和靈敏度分析是非常耗時的。因此,基于此類模型進行結構模型的優化修正是幾乎不能的。Weng[14]、Zhu[15]、Papadimitriou[16]等研究了基于子結構法建立復雜結構的降階模型,并在降階模型基礎上實現了結構動力模型修正。分析認為子結構方法在模型修正中主要有以下優點:(1)將原結構用許多子結構代替,降低了模型規模,有利于加快結構特征值求解和靈敏度分析;(2)子結構具有獨立性,當不確定修正參數只在結構局部范圍內,只需要一個或多個子結構被修正,且僅修正后的子結構需要重新分析,而其它保持不變。但子結構劃分數量的多少對于模型修正計算精度有較大影響,合理的子結構劃分對于大多數學者還是一項挑戰,且隨著結構約束條件的改變子結構方法不可避免的需要重建。

火箭、火車車廂等大型復雜結構的共同特征是縱向尺度顯著大于橫向。設計部門對這些結構進行設計分析時都將其處理作梁式結構,采用各種均勻近似方法將其降階為梁模型,同時根據在已有分析或實驗中獲得的結構前若干階頻率,對模型的參數進行調整,并反復這一試錯過程,建立一個較精確的梁模型。利用這一高度降階的模型,可在初步設計階段估計結構的靜動力性能,確定和分配作用在結構各組成部分上的載荷,考慮結構方案修改的影響。降階梁模型是目前最為廣泛應用的運載火箭降階動力學模型。基于梁平截面假設和位移插值函數的超梁降階模型理論[17-20]主要適用于復雜梁式結構的動力學模型降階,其模型降階減縮基具有清晰的物理意義,有利于降階模型的進一步修改,模型建立靈活,在一定程度上克服了子結構方法的不足。

本文基于超梁降階模型理論研究復雜梁式結構的動力學模型修正問題。主要包括以下內容:基于超梁降階模型理論及考慮可能的模型修正參數,建立由若干梁超單元組成的較高精度大規模降階梁模型;在此基礎上建立以結構彎曲振動特性為修正目標的優化模型,利用優化技術實現其動力學模型的修正;基于修正模型確認準則[21-22]對修正后的模型進行多種動力分析和確認。以典型蒙皮加筋圓柱殼為例,對本文所提出的模型修正方法的可行性和計算效率進行驗證。

1 超梁降階模型理論

針對復雜梁式結構,文獻[17-20]中提出了一種基于梁平截面假設和位移插值函數的超梁降階模型理論。方法是首先將梁式結構用垂直于軸線的平面劃分為若干個梁段,垂直于軸線的平面與軸線的交點定義為主節點。梁彎曲理論假定變形前垂直于軸線的截面變形后仍然保持為平面。如圖1所示,根據梁式結構的特點,在發生非局部變形時,截面的變形可近似為剛體位移,可采用平截面假定來描寫。根據這一假定,結構模型上某一梁段i內任一點j和其在結構軸線上的投影點之間可建立如下位移映射關系:

Uj=Rjqji

(1)

其展開形式為

(2)

式中uj為第j點的位移向量;qji為投影點的準剛體模態;Rj為j點的位移轉換矩陣;ycj、zcj為第j點坐標與投影點坐標在y、z方向的差值。

如圖2所示,取某一梁段i為研究對象,假設其左邊主節點的位移為uL=(u1,v1,w1,θx1,θy1,θz1)T,右邊主節點的位移為uR=(u2,v2,w2,θx2,θy2,θz2)T。依據梁位移插值函數,投影點的位移可用這兩個主節點的位移插值得到,即

(3)

三維梁結構的位移插值函數N可表示為

(4)

式中ξ=Δx/l,l為梁段的長度。

由式(1)及式(3)可建立梁段內任意節點位移與其兩個主節點位移之間的轉換關系矩陣,即

(5)

假設梁段內含有m個節點,則梁段內所有節點的位移向量Ui=(u1,…,um)T可表示為

(6)

假設梁段i的剛度陣和質量陣分別為ki、mi,可計算得到梁段降階后的剛度陣和質量陣分別為

(7)

(8)

位移轉換矩陣Ti大小為6m×12,因此式(7)和式(8)中凝聚得到kRi、mRi均為12×12的矩陣,與梁單元矩陣大小相當,分別定義為梁超單元的剛度陣和質量陣。kRi、mRi可按照序列組裝成一個超梁降階模型的總剛度陣和總質量陣,如果模型均勻,只需要建立一次梁超單元;反之,則需要建立若干個梁超單元。文獻[18]中給出了考慮剪切變形的梁超單元剛度陣的修正方法,提高了超梁降階模型的計算精度,本文不再做詳細介紹。

假設將具有s個節點的梁式結構(6s個自由度)劃分為p(p·s)個梁段,可建立一個具有p+1個主節點的超梁降階模型(6(p+1)個自由度,遠小于原結構的自由度數),其結構剛度陣和質量陣可表示為

(9)

(10)

在模型降階過程中,采用梁平截面假設限制了梁的變形模式,基于超梁降階模型只能求解整體振動模態,而這些模態往往在結構設計中是較重要的。

2 基于超梁降階模型的特征值和靈敏度分析

動力學模型修正問題是以最小化數值分析與試驗測量所得結構的模態特征之間的差異為主要目標,包括頻率和振型的修正。通過求解經典特征值方程可計算獲得結構的頻率和振型:

Kφl=λlMφl

(11)

式中K、M分別為原結構的剛度陣和質量陣;φl、λl分別為其第l階振型向量及相應的特征值。

原結構模型的規模通常較大,其結構剛度陣和質量陣規模也較大。在模型修正迭代過程中,需要多次特征值求解以及相應的靈敏度分析,基于原模型的求解通常需要很大的工作量。因此,文中的特征值和靈敏度分析均是基于超梁降階模型,旨在降低剛度和質量矩陣的規模以及提高優化迭代過程中特征值問題的求解效率。

在模型降階過程中,同時考慮到可能的模型修正參數,劃分后各個梁段內的修正參數是相互獨立的,且假設這些參數與梁段的剛度陣和質量陣呈線性關系,這在實際的模型修正中經常遇到。具體地,假設每一梁段具有Ni個修正參數Xi=(xi1,xi2,…,xiNi)T,梁段的剛度和質量矩陣與修正參數之間的關系分別為

(12)

(13)

式中kk、mk為常數矩陣,獨立于參數xik。

應當指出的是,位移轉換矩陣也獨立于修正參數,因此梁段降階后的梁超單元同樣可表示為修正參數的函數,如式(14)和式(15)所示。

(14)

(15)

則超梁降階模型的剛度和質量矩陣可分別表示為

(16)

(17)

進而可在超梁降階模型基礎上進行結構的特征問題分析:

KRφRl=λRlMRφRl

(18)

φRl、λRl分別表示利用超梁降階模型求出的第l階振型向量及相應的特征值。考慮剪切變形影響可使λRl與λl之間的誤差變得很小,而原結構的振型向量則可由通過位移變換獲得:

(19)

靈敏度分析通常指計算一個特定的響應量相對于一個修正參數的變化率。在超梁降階模型基礎上,由于梁段劃分時保證了修正參數的獨立性,原結構響應對于某一參數的靈敏度分析可通過計算包含該參數的梁超單元矩陣對參數的靈敏度獲得,而其他梁超單元對于該參數的靈敏度為零。特征值和振型向量對于修正參數的靈敏度可表示為[23]

(20)

(21)

式(21)中

(22)

由于超梁降階模型的特征值方程(式(18))求解規模比原始結構的小得多(式(11)),靈敏度矩陣Sλ、Sφ的計算可比基于原結構的分析更快,而特征值和靈敏度計算在整個模型優化修正過程中占主導地位。因此,基于超梁降階模型的動力學模型修正,將顯著提高計算效率,下面通過具體算例進行驗證。

3 動力學模型修正算例

算例為圖3所示自由-自由的蒙皮加筋柱殼。模型長12 000 mm,圓柱段外徑2000 mm,錐段最小外徑為1000 mm,蒙皮厚度為10 mm。沿結構軸向每200 mm布置1條環向加筋,筋條高30 mm,筋條厚30 mm。沿結構環向每12°布置一條軸向加筋,筋條高30 mm,筋條厚15 mm。模型采用鋁合金材料,彈性模量E=70 GPa,泊松比ν=0.3,材料密度ρ=2.7×10-3g/mm3。采用ANSYS殼單元(SHELL181)建立有限元模型,3660個節點,21 960個自由度。

沿圖3所示結構軸線方向,每200 mm劃分為一個梁段,基于超梁降階模型理論,將原結構模型降階為具有60個梁超單元(61個節點,366個自由度)的降階模型,如圖4所示。圓柱段模型均勻,只需建立一個梁超單元,錐段由于模型不均勻,建立了10個梁超單元。利用所建立的超梁降階模型進行模態分析,并采用模態置信因子MAC[24]對降階模型與原結構模型的模態相似度進行比較。表1給出前幾階整體模態分析的結果。

由表1可知,超梁降階模型頻率分析結果與原結構相比,誤差大部分都在5%以內;MAC值也在0.8以上,表明超梁降階模型與原模型的模態相似性好。因此,基于超梁降階模型的結構動力分析能夠滿足模型修正的要求。觀察發現相對于彎曲振動,基于超梁降階模型的扭轉和軸向振動不僅頻率誤差較小,且模態相似度也較高,這是因為建立降階模型時,采用了平截面假定,能較好地模擬這些振動模式,在后面分析中,重點將關注彎曲振動的分析與修正。

階次原模型/Hz降階模型/Hz誤差/%MAC1階彎曲70.0867.194.120.851階扭轉118.96122.152.680.992階彎曲151.94154.871.920.911階軸向192.47198.823.290.992階扭轉236.75243.112.680.993階彎曲237.78250.955.540.83

在模型修正中,所模擬的“試驗”數據通常是通過

在有限元模型上故意引入一些缺陷,然后修正分析模型來識別這些“損害”[25]。文中修正所需的“試驗”數據是通過故意減少圖3所示結構某些區域的剛度后分析獲得。

表2給出了假定的剛度損失區域及損失值。在此基礎上,以圖4所示超梁降階模型為分析模型,并進行優化修正,使其模態分析結果與基于表2中數據獲得的“試驗”結果相吻合。

表2 假定的剛度損失區域及損失值

以模型的彎曲振動頻率和振型為優化目標,建立如下優化問題:

(23)

以圖4所示超梁降階模型中每個梁超單元的剛度作為優化參數,因此一共有60個優化設計變量。采用MMA[26]算法進行優化求解,優化終止準則設為兩次迭代的優化目標之差小于0.001。

模型修正前后的頻率和模態MAC值的對比如表3所示。

表3 修正前后模型的頻率和振型比較

在情況1中,Part3和Part4的剛度被故意減少了20%,Part6和Part7區域的剛度降低了30%。從表3可知,修正前的頻率和MAC值與試驗數據相比誤差均較大,修正后的數據得到了較大改善。不失一般性,在情況2中選擇不同區域減少剛度,具體為Part4和Part5減少30%,Part8和Part9減少20%。從表3可見,經過優化修正,分析獲得的結果與試驗數據吻合較好。

圖5給出了兩種情況下優化修正后每個梁超單元的剛度改變量,與表2中假定的剛度損失區域和損失值較好的相對應,修正實現了對假定“損害”的識別。由于基于超梁降階模型的分析結果與原結構模型分析數據之間一直存在一定誤差,結構其余部分的剛度也有稍微的改變,但對修正結果的影響可忽略。

需要注意的是,優化修正過程中的特征值和靈敏度分析,均是在圖4所示超梁降階模型基礎進行求解的。從計算效率上看,每次優化迭代求解結構前三階整體彎曲振動,在相同配置電腦上,均采用MATLAB軟件特征值求解程序,基于超梁降階模型所需要的計算時間為0.08 s,而原結構模型則需要10.8 s。因此,基于超梁降階模型的模型修正能顯著提高計算效率。

4 修正模型的確認

修正后模型不僅要滿足修正頻段范圍內的計算精度要求,還要能預測修正頻段范圍外的結構動力特性,或同時對修正后模型和試驗模型做一定修改后,修正后模型還能準確預測試驗模型的各項動力學特征。為了驗證修正后模型的可靠性,針對上述情況2的修正后降階模型和試驗模型,研究了一端固支約束條件下結構的頻率分析和瞬態沖擊分析,如圖6所示。

表4給出了一端固支約束條件下,基于修正后降階模型計算得到的前幾階整體振動分析的結果。由表4可知,修正后降階模型施加一端固支約束后,頻率分析結果與試驗結果一致,誤差在5%左右,滿足工程計算精度的要求。而一端固支降階模型的建立只需要在圖4所示自由-自由降階模型基礎上,在1號節點施加固支約束即可,不需要重新建立降階模型,有效提高了分析效率。

表4 修正后降階模型的一端固支模態分析結果

在圖6所示一端固支約束模型的自由端截面所有節點上施加y軸負方向的沖擊載荷(如圖7所示),分別基于修正后降階模型和試驗模型進行瞬態沖擊響應分析。

圖8~圖10分別給出了模型自由端節點(編號3605)在沖擊載荷作用下,沿y方向的位移、速度和加速度響應曲線。由圖可知,修正后降階模型的響應曲線與試驗模型的響應曲線幾乎一致,進一步證明了修正后模型的有效性。

從計算時間上看,利用修正后降階模型完成沖擊載荷分析所需時間為9.30 s,而試驗模型則需要292.03 s,遠大于前者。同時,也比較了結構上其余節點的響應,位移、速度和加速度曲線,均吻合較好。

5 結論

(1)推導了基于超梁降階模型的特征值及靈敏度分析求解公式,建立了包含結構頻率和振型的優化修正目標,利用優化方法實現了復雜蒙皮加筋柱殼結構的動力學模型修正。計算結果,表明基于超梁降階模型的優化修正具有較高的計算精度和求解效率。

(2)通過一端固支約束條件下的頻率和瞬態沖擊響應分析,對修正后模型進行評估和確認,結果表明了本文提出的模型修正方法的有效性。

[1] 賀爾銘,陳熠,李玉龍,等.機翼雙梁模型的動力學修正及應用[J].應用力學學報,2013,30(3):367-372.

[2] 馬雙超,臧朝平,蘭海波.某航空發動機機匣的動力學模型修正[J].航空動力學報,2013,28(4):878-884.

[3] Mottershead J E,Friswell M I.Model updating in structural dynamics:a survey[J].Journal of Sound and Vibration,1993,167(2):347-375.

[4] 韓建平,駱勇鵬,鄭沛娟,等.基于響應面的剛構-連續組合梁橋有限元模型修正[J].工程力學,2013,30(12):85-90,106.

[5] 張根輩,臧朝平,王曉偉,等.螺栓連接框架結構的有限元模型修正[J].工程力學,2014,31(4):26-33.

[6] 楊春峰,路林華,張盛,等.基于SiPESC.OPT的結構動力模型修正研究[J].大連理工大學學報,2013,53(5):630-636.

[7] 朱安文,曲廣吉,高耀南,等.結構動力模型修正技術的發展[J].力學進展,2002,32(3):337-348.

[8] 李輝,丁樺.結構動力模型修正方法研究進展[J].力學進展,2005,35(2):170-180.

[9] Brownjohn J M W,Xia P Q,Hao H,et al..Civil structure condition assessment by FE model updating:methodology and case studies[J].Finite Element in Analysis and Design,2001,37(10):761-75.

[10] 鄭惠強,陳鵬程,宓為建,等.大型橋吊結構動力有限元模型修正[J].同濟大學學報,2001,29(12):1412-1415.

[11] 劉繼承,周傳榮.一個基于優化的有限元模型修正方法[J].振動與沖擊,2003,22(2):33-35.

[12] 王元清,姚南,張天申,等.基于最優化理論的多階段模型修正及其在橋梁安全評估中的應用[J].工程力學,2010,27(1):91-97.

[13] Bakira P G,Reynders E,Roeck G D.Sensitivity-based finite element model updating using constrained optimization with a trust region algorithm[J].Journal of Sound and Vibration,2007,305(1-2):211-25.

[14] Weng S,Xia Y,Xu Y,et al..Substructure based approach to finite element model updating[J].Computers and Structures,2011,89(9-10):772-782.

[15] Zhu D P,Dong X J,Wang Y.Substructure model updating through modal dynamic residual approach[C]//Proceedings of the 9th International Workshop on Structural Health Monitoring,Stanford,CA.2013.

[16] Papadimitriou C,Papadioti D C.Component mode synthesis techniques for finite element model updating[J].Computers and Structures,2013,126(15):15-28.

[17] Wang Wensheng,Cheng Gengdong,Li Quhao.Fast dynamic performance optimization of complicated beam type structure based on two new reduced physical models[J].Engineering Optimization,2013,45(7):835-850.

[18] Cheng Gengdong,Wang Wensheng.Fast dynamic analysis of complicated beam-type structure based on reduced super beam model[J].AIAA Journal,2014,52(5):952-963.

[19] 李取浩,王文勝,程耿東.基于梁理論的復雜梁式結構低階頻率快速求解方法[J].工程力學,2014,31(11):1-8.

[20] 王文勝,程耿東,郝鵬.基于超梁降階模型的蒙皮加筋圓柱殼頻率分析[J].固體火箭技術,2015,38(3):401-406.

[21] Link M,Friswell M I.Generation of validated structural dynamic models:results of a bechmark study utilising the GARTEUR SM-AG19 test-bed[J].Mechanical Systems and Signal Processing,2003,17(1):9-21.

[22] 畢司峰,鄧忠民.基于隨機抽樣與距離判別的GARTEUR模型修正與確認研究[J].航空學報,2013,34(12):2757-2767.

[23] Mottershead J E,Link M,Friswell M.The sensitivity method in finite element model updating:A tutorial[J].Mechanical Systems and Signal Processing,2011,25(7):2275-2296.

[24] 邱吉寶,向樹紅,張正平.計算結構動力學[M].合肥:中國科學技術大學出版社,2009.

[25] Modak S V,Kundra T K,Nakra B C.Comparative study of model updating methods using simulated experimental data[J].Computers and Structures,2002,80(5-6):437-47.

[26] Svanber K.The method of moving asymptotes-a new method for structural optimization[J].International Journal for Numerical Methods in Engineering,1987,24(2):359-373.

Dynamicmodelupdatingandvalidationofcomplicatedbeam-typestructuresbasedonreducedsuperbeam

WANG Wen-sheng,WEI Hao-jie,HOU Zhong-hua,MEI Qun, LI Yi-fan

(1.Department of Engineering Mechanics,School of Civil Engineering,Henan University of Science and Technology,Luoyang,Henan 471023,China)

A method based on reduced super beam for dynamic model updating of complicated beam-type structures was proposed in this paper.Firstly,a reduced super beam model with high accuracy was constructed for the complicated beam-type structure,considering the possible updating parameters.Then, the model parameters of reduced super beam model were updated using the tested modal data based on sensitivity analysis and optimization algorithm, and then an updated model that can represent the dynamic behavior of the practical structure accurately was obtained.Finally,application in a variety of dynamic problem analysis was performed on the updated model and the updating results were assessed by comparing with the practical structure.A typical stiffened cylindrical shell model,serving as the illustrative example,was employed for model updating and validation.The results demonstrate that the proposed method is valid and efficient.

reduced super beam;beam-type structure;dynamics;model updating;optimization

V414.3

A

1006-2793(2017)05-0620-07

10.7673/j.issn.1006-2793.2017.05.016

(編輯:薛永利)

猜你喜歡
有限元優化結構
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 亚洲AV电影不卡在线观看| 精品国产99久久| 国产在线视频二区| 丁香六月激情综合| 亚洲高清中文字幕| 久草青青在线视频| 亚洲一区二区日韩欧美gif| 国产精品伦视频观看免费| 色婷婷在线播放| 玖玖精品在线| 国产精品密蕾丝视频| 日韩欧美综合在线制服| 99re在线视频观看| 亚洲最新地址| 四虎在线观看视频高清无码| 亚洲欧洲日产无码AV| 最新加勒比隔壁人妻| 日本不卡视频在线| 真实国产乱子伦视频| 18禁影院亚洲专区| 亚洲成a人片| 国产一区亚洲一区| 欧美色综合网站| 波多野结衣亚洲一区| 欧美精品一二三区| 亚洲乱强伦| 在线观看亚洲人成网站| 国产在线自揄拍揄视频网站| 国产成人永久免费视频| 久久国产免费观看| 精品国产黑色丝袜高跟鞋 | 欧美全免费aaaaaa特黄在线| 色综合中文字幕| 免费人成视网站在线不卡| 久久国产精品国产自线拍| 欧美a在线| 毛片网站观看| 久久综合亚洲色一区二区三区| 国产在线观看91精品亚瑟| 欧美综合区自拍亚洲综合绿色 | 波多野结衣AV无码久久一区| 综合久久五月天| 国产一在线观看| 国产激情在线视频| 亚洲v日韩v欧美在线观看| 亚洲成人网在线观看| 中文字幕在线播放不卡| 亚洲欧美国产五月天综合| 亚洲伊人久久精品影院| 欧美伊人色综合久久天天| 97se亚洲综合在线| v天堂中文在线| 无码网站免费观看| 欧美视频免费一区二区三区| 国产精品99一区不卡| 日韩少妇激情一区二区| 精品国产免费人成在线观看| 3344在线观看无码| 亚洲成A人V欧美综合| 黑人巨大精品欧美一区二区区| 大香网伊人久久综合网2020| 亚洲国产综合精品一区| 少妇精品在线| 国产精品护士| 亚洲第一成年网| 国产精品午夜福利麻豆| 精品91视频| 欧美亚洲网| 天堂成人av| 国产精品久线在线观看| 国产午夜在线观看视频| 日本亚洲欧美在线| 四虎永久免费地址| 久久永久精品免费视频| 亚洲天堂网2014| 国产精品综合久久久| 国产精品一区在线麻豆| 手机成人午夜在线视频| 国产精品免费久久久久影院无码| 69av免费视频| 永久天堂网Av| 亚洲一区二区三区在线视频|