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

考慮加速度頻響函數不確定的有限元模型修正

2023-08-04 11:44:32彭珍瑞張亞峰張雪萍
哈爾濱工業大學學報 2023年8期
關鍵詞:有限元結構模型

彭珍瑞,張亞峰,張雪萍

(蘭州交通大學 機電工程學院,蘭州 730070)

用實驗數據作為精確的參考數據來改進有限元預測,稱為有限元模型修正[1]。當前多數模型修正方法屬于確定性方法,沒有考慮結構參數和響應的不確定性,降低了其工程實際應用價值。然而,大多數實際工程結構的材料參數、幾何尺寸和試驗數據的不確定性廣泛存在。因此,研究涉及不確定性分析的有限元模型修正非常必要[2-3]。

在模型修正過程中,使用代理模型代替有限元模型可減小計算量,是解決復雜工程結構優化問題的有效途徑之一[4]。常見的代理模型中,徑向基(Radial basis function,RBF)模型結構簡單,預測精度高,廣泛應用于土木、航天和機械等領域[5]。因此,本文選用RBF模型為模型修正的代理模型。

近年來,眾多學者研究了考慮不確定性的模型修正方法,并取得了一定的成果。不確定性有限元模型修正主要可分為概率和非概率方法,基于概率的模型修正一般為隨機模型修正方法,基于非概率的模型修正有模糊方法和區間分析方法[6]。概率方法需要豐富的試驗數據以得到較準確的統計信息,計算成本高[7];基于模糊數學的模糊方法在確定隸屬函數時容易引入額外的不確定性[8];區間分析方法對數據要求低,適用于缺乏數據信息的情況[9]。Fang等[10]將區間響應面模型用于區間模型修正,實現參數和響應的區間估計。駱勇鵬等[11]提出一種基于響應面的靈敏度模態區間分析方法,通過比較相對模態區間靈敏度判斷結構響應對參數的敏感程度。姜東等[12]采用區間擴張理論與靈敏度方法,用于不確定性和初始有限元模型誤差均較小情況的有限元模型修正。Deng等[13]提出了一種利用一階攝動法和RBF神經網絡進行結構參數區間識別的兩步確定性修正方法。Wang等[14]提出了一種改進的參數攝動法來預測不確定結構的特征值區間,該方法考慮了計算效率,但它依賴于區間變量的初始值和不確定程度。Deng等[15]提出了一種結構不確定性度量區間識別策略,可以得到系統響應的精確區間估計。Zheng等[8]提出了一種利用泛灰數學和高斯過程回歸模型的兩步區間模型修正方法。上述區間模型修正研究均以結構振動的模態參數作為響應量,相比于模態參數,結構振動的頻響函數(Frequency response function,FRF)則包含了更充分的結構信息,能夠更加準確地反映結構的動力學特性。同時,采用FRF進行模型修正可以有效避免由模態分析所引入的誤差[16-17]。此外,小波變換是進行信號時頻研究的重要工具,可以聚焦信號的細節特征,已被廣泛應用于水文、地震勘探、生態保護和經濟管理等領域的預測研究[18]。本文對加速度FRF小波變換提取低頻小波系數作為響應特征量進行模型修正。

綜上所述,本文提出了基于加速度頻響函數的區間有限元模型修正方法。首先構造RBF代理模型并對RBF模型方差值進行了優選;然后利用區間重疊率和巴氏距離分別構造兩步和同步修正所需要的目標函數,用來度量兩個樣本區間分布的相似性和相異性;最后結合灰色數學方法,運用花朵授粉算法分別進行兩步和同步待修正參數區間的求解。通過支撐框架結構,空間桁架結構和簡支梁結構驗證所提方法的有效性。

1 小波變換

小波變換通過伸縮平移運算對信號多尺度細化,最終高頻處時間細分,低頻處頻率細分,凸顯問題某些方面的特征。對一維信號y(t)的小波變換,給定一個基函數,令

(1)

式中:a為伸縮因子且大于0,b為平移因子。為了計算方便,在實際分析和應用中需要將連續小波及小波系數進行離散化處理,則y(t)的小波變換為

(2)

式中,WTy(a,b)為y(t)的小波系數。在一定的頻率點數下,小波分解層數過少會使每層小波系數過多,計算量過高,反之,分解層數過多會使每層小波系數過少而影響模型修正效果。通過對比多次試驗效果,本文設定基函數為db3,分解層數為4、5層。對加速度FRF進行小波變換,提取低頻小波系數作為構造RBF模型的輸出,也作為加速度FRF的響應特征量進行模型修正。

2 RBF模型的構造

2.1 RBF神經網絡的基本原理

在d維設計空間中生成一組數量為N的初始樣本點x=[x1,x2,…,xN]∈Rd,目標的真實響應值為y(x1),y(x2),…,y(xN),RBF插值函數的數學形式為

(3)

式中:N為樣本點個數,wk為權值系數,‖·‖為歐氏距離,φ(·)為基函數。

由s(xk)=y(xk)可得:Φ·W=Y,即W=Φ-1·Y,其中

式中:W=(w1,w2,…,wk)T為權系數向量,Y=[y(x1),y(x2),…,y(xk)]T為樣本點的真實響應向量。選用應用較廣的高斯插值基函數,表達式為

(4)

式中:r=‖x-xk‖為當前預測點與任意樣本點之間的歐氏距離,σ2為方差。

2.2 RBF模型的構造及驗證

建立一個嚴格RBF模型,需要先設定高斯插值基函數的方差值,方差直接影響模型預測精度。通過鯨魚優化算法(Whale optimization algorithm,WOA)優選最佳方差值。WOA以泡泡網搜索策略為靈感,主要步驟為包圍獵物、氣泡網攻擊和搜索獵物[19]。

采用拉丁超立方抽樣方法抽取初始待修正參數±20%區間內的樣本,將其分為訓練集和測試集。建立目標函數:

(5)

(6)

3 灰色數學方法

實際工程結構中所能獲得的樣本信息非常有限,且具有不確定性。灰色數學方法需要的樣本量小,主觀依賴性小,具有效率高和精度高的優點[9]。在工程問題中,不確定信息可以使用灰色數學方法來獲得實際值的區間估計,其基本原理如下。

對一組初始不規則數據:

Y={y(c),c=1,2,…,n}

(7)

將數據從小到大排列,得到:

Y(0)={y(0)(c),c=1,2,…,n}

(8)

把Y(0)的元素按順序相加,得到一個新的Y(1)序列:

Y(1)={y(1)(c),c=1,2,…,n}={y(0)(1),y(0)(1)+

y(0)(2),…,y(0)(1)+y(0)(2)+…+y(0)(n)}

(9)

定義:

(10)

4 區間分析理論

不確定度量用以表征理論分析響應和測量響應之間的區間一致性,反映兩個區間分布的相似性和相異性,它直接影響模型修正的精度。

4.1 區間重疊率

(11)

定義區間P和區間Q的區間重疊率(Interval overlap ratio,IOR)[15]為

(12)

式中:運算符∪、∩分別為并集、交集;len(P∩Q)、len((P∪Q)-(P∩Q))分別為表征區間P和區間Q之間的相似程度和不相似程度。由式(12)可知,在P=Q的情況下,IOR(P|Q)=1。對于區間模型修正問題,利用IOR(P|Q)可以度量理論分析數據和測量數據之間區間分布的相似性。

4.2 區間巴氏距離

巴氏距離用來度量兩個概率分布之間的相似性,巴氏距離越小,表示兩個概率分布越相似[20]。在高斯分布下U=N(mU,θU)和V=N(mV,θV)的巴氏距離為

(13)

式中:θ=(θU+θV)/2。本文將式(13)中概率分布的均值和方差定義為區間分布的區間中點和半徑。將高斯分布下的巴氏距離推廣應用于區間有限元模型修正,度量兩個區間分布的相似性,稱為區間巴氏距離,表示為

(14)

5 模型修正過程

花朵授粉算法(Flower pollination algorithm,FPA)具有異花授粉(全局尋優)和自花授粉(局部尋優)兩大特性。異花授粉過程可以逃離任何局部景觀,進而探索更大的空間;自花授粉過程使相似的解被更頻繁地選擇,本質上是指數收斂,收斂速度更高,性能優于遺傳算法和粒子群算法[21]。本文運用FPA分別進行兩步和同步的區間有限元模型修正。

結構待修正參數具有不確定性,將其描述為實數域R內的一個閉區間:

(15)

1)兩步修正參數區間中點和半徑。第1步求解結構待修正參數的區間中點。建立目標函數:

(16)

(17)

2)同步修正參數區間中點和半徑。首先抽取200個參數區間樣本,通過已構造的RBF模型預測樣本響應特征量;然后使用灰色數學方法估計響應特征量區間中點和半徑;最后運用FPA直接以式(14)為目標函數同時求解待修正參數區間中點和半徑。其中,在每次同時求解參數區間中點和半徑迭代過程中,同樣抽取200個區間樣本參與有限元模型修正。值得說明,式(14)完全可以替換式(17),實現兩步修正參數區間中點和半徑,由于篇幅所限,不再詳細說明。

6 算 例

6.1 數值算例1

如圖2所示支撐框架結構,采用空間剛架單元模擬該結構的有限元模型。該結構共有54個節點,144個自由度,70個梁單元和6個邊界支撐。其中,1~6為邊界上的固定節點。箭頭表示激勵位置,?表示測點位置。選擇結構彈性模量E和材料密度d的區間中點及半徑為模型待修正參數,初始有限元模型E和d的區間中點分別為231 GPa和7 020 kg/m3,試驗模型E和d的區間中點分別為210 GPa和7 800 kg/m3,區間半徑分別為3 GPa和100 kg/m3。初始有限元模型和試驗模型E和d的區間中點的相對誤差分別為10%和-10%。

圖2 支撐框架結構

本算例中選擇500個頻率點,設定小波變換分解層數為5。根據所述方法,通過WOA優化確定最佳方差值來建立RBF模型,得到RBF模型的方差值為0.089 3。然后計算得到RMSE為9.29×10-5,說明構造的RBF模型預測精度較高。為進一步評估RBF模型的預測精度,圖3給出了測試集樣本第2和第10響應特征量的有限元模型計算值和RBF模型預測值,從圖3可以看到,有限元模型計算值和RBF模型預測值幾乎全部重合,可以代替有限元模型。

圖3 測試集樣本的預測值和真實值

6.1.1 兩步修正

以試驗模型參數區間中點和半徑抽取50個樣本進行計算有限元模型加速度FRF,通過小波變換并提取第5層低頻小波系數作為仿真試驗響應特征量,使用灰色數學方法評估響應特征量區間。然后按照所述兩步修正過程,通過FPA第1步迭代求解結構參數區間中點,第2步迭代求解待修正參數的區間半徑,修正后的結構參數結果見表1。

表1 結構修正前、后的參數區間

從表1可以看出,修正后的相對誤差較小,表明所提兩步進行區間模型修正取得了很好的效果,對結構參數區間的修正是精確的,也表明了所提方法可以有效地解決試驗小樣本的區間模型修正問題。為進一步驗證修正效果,使用表1中修正后的結構參數區間抽取50個樣本分別計算修正后的有限元模型加速度FRF(稱為RBF模型FRF)和響應特征量。其中,FRF區間收斂如圖4所示,計算修正后有限元模型第4和第12響應特征量樣本與試驗模型響應特征量數據比較的散點如圖5所示。

圖4 頻響函數區間收斂

圖5 響應特征量區間散點圖

從圖4、5可以看出,修正后的有限元模型和試驗模型的FRF區間吻合良好,修正后的有限元模型響應特征量計算樣本和試驗模型響應特征量數據也吻合良好。綜上所述,驗證了所提兩步進行區間有限元模型修正的有效性。

6.1.2 同步修正

按照所述同步修正過程,使用與兩步修正相同的仿真試驗數據,通過FPA同時迭代求解結構參數區間中點和半徑,修正后的結果見表2。

表2 結構修正前、后的參數區間

從表2可以看出,同步修正與兩步修正具有相同的修正效果。同樣,為進一步驗證修正效果,使用表2中修正后的結構參數區間抽取50個樣本分別計算修正后的有限元模型FRF和響應特征量。其中,FRF區間收斂如圖6所示,計算修正后有限元模型第4和第12響應特征量樣本與試驗模型響應特征量數據比較的散點如圖7所示。從圖6、7可以看出,修正后的有限元模型和試驗模型的FRF區間及響應特征量區間樣本同樣吻合良好。綜上所述,驗證了所提同步進行區間有限元模型修正的有效性。

圖6 頻響函數區間收斂

圖7 響應特征量區間散點圖

6.2 數值算例2

如圖8所示空間桁架結構,該桁架模型由72個桿單元組成,共有20個節點和48個自由度。單元橫截面積為1 cm2,箭頭表示激勵位置,×表示測點位置。選擇結構彈性模量E1(直桿)、E2(斜桿)和材料密度d的區間中點及半徑為模型待修正參數。初始有限元模型E1、E2和d的區間中點分別為231、209 GPa和7 020 kg/m3,試驗模型E1、E2和d的區間中點分別設定為210、190 GPa和7 800 kg/m3,初始有限元模型和試驗模型E1、E2和d的區間中點的相對誤差分別為10 %、10 %和-10 %。設定兩種不同的試驗參數區間半徑工況來驗證所提模型修正方法的有效性。工況1為試驗模型參數E1、E2和d區間半徑分別設定為3、2 GPa和60 kg/m3;工況2為試驗模型參數E1、E2和d區間半徑分別設定為6、5 GPa和200 kg/m3。

圖8 空間桁架結構

本算例中選擇500個頻率點,設定小波變換分解層數為5。采用RBF模型的構造及驗證所述方法建立RBF模型,得到RBF模型的方差值為0.246 9,迭代情況如圖9所示。為進一步評估RBF模型的預測精度,圖10給出了測試集樣本第7和第15響應特征量的有限元模型計算值和RBF模型預測值,從圖10可以看到,有限元模型計算值和RBF模型預測值幾乎全部重合,可以代替有限元模型。

圖9 鯨魚優化算法迭代曲線

圖10 測試集樣本的預測值和真實值

6.2.1 兩步修正

以試驗模型參數區間中點和半徑抽取100個樣本進行計算有限元模型加速度FRF,通過小波變換并提取第5層低頻小波系數作為仿真試驗響應特征量,并使用灰色數學方法估計響應特征量區間。然后根據所述兩步修正方法進行區間模型修正。其中,工況1結構參數區間中點和半徑的迭代情況分別如圖11、12所示,修正結果見表3,工況2結構參數區間修正結果見表4。

表3 空間桁架結構修正前、后參數區間(工況1)

表4 空間桁架結構修正前、后參數區間(工況2)

圖11 花朵授粉算法求解參數區間中點

圖12 花朵授粉算法求解參數區間半徑

由表3、4可以看出,對兩種工況參數區間修正的相對誤差值都較小,表明所提兩步修正方法取得很好的效果,且在不同試驗響應區間下對參數區間的修正具有魯棒性。進一步驗證修正效果,使用表3、4中修正后的結構參數區間分別抽取100個樣本計算兩種工況修正后的有限元模型FRF和響應特征量。其中,工況2的FRF區間收斂如圖13所示,計算工況1修正后的有限元模型第8和第16響應特征量樣本與試驗模型響應特征量數據比較的散點如圖14所示。從圖13、14可以看出,修正后的有限元模型和試驗模型的FRF區間及響應特征量區間樣本吻合良好。

圖13 頻響函數區間收斂

圖14 響應特征量區間散點圖

6.2.2 同步修正

根據所提同步修正方法,使用與兩步修正相同的仿真試驗數據來同時求解參數區間中點和半徑。其中,工況1結構參數區間中點和半徑同步修正的迭代情況如圖15所示。兩種工況的結構參數修正結果分別見表5、6。

表5 空間桁架結構修正前、后參數區間(工況1)

表6 空間桁架結構修正前、后參數區間(工況2)

圖15 花朵授粉算法同時求解參數區間中點和半徑

由表5、6可以看出,在兩種工況下同步修正參數區間的相對誤差值都較小,表明在不同試驗響應區間下對參數區間的修正具有魯棒性,說明所提同步修正方法與兩步修正方法具有相同的效果。再由圖11、12和圖15相比可以看出,同步修正參數區間迭代350次之前就已收斂,而兩步修正參數區間的第2步迭代收斂就需500次以上,因此在時間上同步修正相比兩步修正節省了許多。同樣進一步驗證修正效果,使用表5、6修正參數區間分別抽取100個樣本計算修正后兩種工況的FRF和響應特征量。工況2的FRF區間收斂如圖16所示,計算工況1修正后的有限元模型第8和第16響應特征量樣本與試驗模型響應特征量數據比較的散點如圖17所示。從圖16、17可以看出,修正后的有限元模型與試驗模型的FRF區間和響應特征量區間樣本同樣吻合良好。

圖16 頻響函數區間收斂

圖17 響應特征量區間散點圖

6.3 試驗算例

如圖18所示,對長2 000 mm,寬100 mm,高10 mm的簡支梁結構進行試驗研究。該梁材料為Q235鋼,將其劃分為20個單元和21個節點。選擇16節點為激勵位置,選擇5、8、11、14、和17節點為測量位置。選擇結構彈性模量E和材料密度d的區間中點及半徑為結構待修正參數。根據Q235鋼材料參數范圍,設置初始有限元模型E和d的區間中點分別為205 GPa和7 850 kg/m3。

圖18 簡支梁

測試試驗過程通過東方所的DASP軟件進行50次試驗,設置采樣通道為6,采樣頻率為4 000 Hz,采樣點數為4 096,直接采集加速度響應數據和加速度FRF數據。圖19給出了一次試驗的16節點激勵與11節點響應的加速度響應曲線,圖20給出了16節點激勵與11節點響應的初始有限元模型和實際測量的加速度FRF中值曲線。

圖20 加速度頻響函數

6.3.1 兩步修正

本試驗算例中選擇250個頻率點,設定分解層數為4。將16節點激勵與11節點響應的加速度FRF經過小波變換并提取第4層低頻小波系數作為試驗響應特征量,通過灰色數學方法評估試驗響應特征量區間。然后按照兩步修正過程,通過FPA第1步迭代求解結構參數區間中點,第2步迭代求解待修正參數的區間半徑,修正后結果見表7。

表7 簡支梁修正后的參數值

為驗證修正效果,使用表7中修正后的結構參數區間抽取50個樣本計算修正后的有限元模型加速度FRF。圖21給出了16節點激勵與11節點響應的修正前、后的有限元模型加速度FRF中值曲線和試驗加速度FRF中值曲線,從圖21可以看出修正后的有限元模型加速度FRF中值曲線和試驗加速度FRF中值曲線大致吻合。圖22給出了16節點激勵與11節點響應的加速度FRF區間收斂,從圖22可以看出修正后的有限元模型加速度FRF區間與試驗加速度FRF區間大致吻合。

圖21 加速度頻響函數

圖22 頻響函數區間收斂

為進一步驗證修正效果,測量不同位置的加速度FRF與修正后的有限元模型相對應位置的加速度FRF進行比較。圖23給出了16節點激勵與17節點響應的修正前、后有限元模型加速度FRF中值曲線和試驗加速度FRF中值曲線,從圖23可以看出修正后的有限元模型加速度FRF中值曲線與試驗加速度FRF中值曲線大致吻合。圖24給出了16節點激勵與17節點響應的修正后有限元模型加速度FRF區間收斂和試驗加速度FRF區間收斂。從圖24可以看出,修正后的有限元模型加速度FRF區間和試驗加速度FRF區間大致吻合。綜上所述,驗證了所提兩步進行區間有限元模型修正的有效性。

圖23 加速度頻響函數

圖24 頻響函數區間收斂

6.3.2 同步修正

按照所述同步修正過程,使用與兩步修正相同的試驗數據,通過FPA同時迭代求解結構參數區間中點和半徑,修正后的結果見表8。

表8 簡支梁結構修正后的參數值

同樣為驗證修正效果修正,使用表8中修正后的結構參數區間抽取50個樣本計算修正后的有限元模型加速度FRF。圖25、26分別給出了16節點激勵與11節點響應的有限元模型與試驗加速度FRF中值曲線和區間收斂,從圖25、26可以看出修正后的有限元模型與試驗加速度FRF中值曲線和區間大致吻合。

圖25 加速度頻響函數

圖26 頻響函數區間收斂

同樣為進一步驗證修正效果,測量不同位置的加速度FRF進行比較,圖27、28分別給出了16節點激勵與17節點響應的有限元模型與試驗加速度FRF中值曲線和區間收斂,可以看出修正后的有限元模型與試驗加速度FRF中值曲線和區間大致吻合。綜上所述,驗證了所提同步進行區間有限元模型修正的有效性。

圖27 加速度頻響函數

圖28 頻響函數區間收斂

7 結 論

1)考慮參數不確定性對于結構響應的影響,提出了以結構參數區間中點和半徑為修正目標的兩步和同步的區間模型修正方法,可以較好地修正結構參數的區間中點和半徑,且在不同試驗響應區間下對參數區間的修正都具有魯棒性。同時,可以有效地解決試驗小樣本的不確定性模型修正問題。

2)算例表明,兩步和同步修正都可以達到相同的效果,且在迭代求解過程中表征區間進行模型修正的參數樣本量較小,提高了模型修正效率。但在時間上同步修正相比兩步修正節省了許多。

3)在使用參數樣本參與不確定性模型修正時構造了RBF模型,并采用WOA優選了RBF模型的方差值,使所構造的RBF模型具有良好的擬合精度和預測能力,能代替有限元模型進行迭代計算,提高模型修正效率。

4)將加速度FRF經過小波變換,提取低頻小波系數作為區間模型修正的響應特征量,具有保留FRF特性的優點,且可以大量減少輸出響應數目,提高了模型修正計算效率。

猜你喜歡
有限元結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
3D打印中的模型分割與打包
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 97国产在线视频| 9啪在线视频| 欧美日韩高清在线| 夜夜爽免费视频| 欧美精品啪啪| 成年A级毛片| 国产浮力第一页永久地址| 日韩黄色精品| 国产精品不卡永久免费| 国产97视频在线观看| 熟女成人国产精品视频| 欧美精品在线免费| 成年女人a毛片免费视频| 国产午夜福利在线小视频| 欧美翘臀一区二区三区| 午夜福利免费视频| 国产精品入口麻豆| 无遮挡一级毛片呦女视频| 日本人妻一区二区三区不卡影院| 57pao国产成视频免费播放| 国产精品无码AV中文| 久久国产亚洲偷自| 国产乱人视频免费观看| 精品一区二区三区视频免费观看| 亚洲国产亚洲综合在线尤物| 久久黄色一级片| 成人午夜免费视频| 亚洲欧美一区二区三区麻豆| 国产精品流白浆在线观看| 99热这里只有精品免费| 亚洲色图欧美视频| 一本久道热中字伊人| 久久不卡国产精品无码| 国产成人亚洲无吗淙合青草| 啦啦啦网站在线观看a毛片| 国产成人久视频免费| 日本三区视频| 亚洲人成日本在线观看| 国产99免费视频| 欧美色综合网站| 久久成人国产精品免费软件| 亚洲国产日韩视频观看| 国产aⅴ无码专区亚洲av综合网| 国产www网站| 孕妇高潮太爽了在线观看免费| 亚洲制服中文字幕一区二区| av性天堂网| 中日韩欧亚无码视频| 久久国产热| 欧美成人二区| 精品一区二区无码av| 国产欧美视频综合二区| www.亚洲国产| 欧美日韩理论| 不卡的在线视频免费观看| 精品一区二区三区水蜜桃| 狠狠亚洲五月天| 91久久精品国产| 国产在线视频导航| 国产人在线成免费视频| 亚洲免费黄色网| 91九色国产porny| 亚洲国产中文综合专区在| 日韩少妇激情一区二区| 久久综合国产乱子免费| 亚洲成av人无码综合在线观看| 亚洲成人高清在线观看| 777国产精品永久免费观看| 久久久噜噜噜| 亚洲男人天堂2018| 97国产精品视频自在拍| 成人小视频网| 996免费视频国产在线播放| 国产美女无遮挡免费视频网站 | 亚洲精品大秀视频| 成人毛片免费在线观看| 免费人成在线观看成人片 | 一级成人a做片免费| 免费毛片a| 国产成人精品在线1区| 日韩精品成人网页视频在线| 中文字幕天无码久久精品视频免费|