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

艦載機壁板剪切后屈曲承載能力預測與試驗驗證

2019-04-22 10:43:38劉存張磊楊衛平
航空學報 2019年4期
關鍵詞:模態有限元結構

劉存,張磊,楊衛平

航空工業第一飛機設計研究院 強度設計研究所,西安 710089

加筋壁板是艦載飛機的重要承載結構,在飛機著艦撞擊構型下,蒙皮承受的剪應力會出現大于臨界剪應力的情況,此時蒙皮開始失穩,由于壁板加強筋的存在,蒙皮失穩后還能夠繼續承受增加的外載荷,但其內力重新分配,失穩波紋數隨之增多并趨于規則,最終可發生剪切后屈曲破壞。目前工程上對加筋壁板的剪切屈曲特性的計算仍在初始屈曲階段,而需要考慮大擾度和塑性效應的工程后屈曲計算方法還不成熟。因此,采用有限元MSC.NASTRAN軟件進行加筋壁板剪切后屈曲承載能力預測十分必要。

國內外學者對加筋壁板剪切后屈曲承載能力進行了研究。張國凡等[1]應用ABAQUS UMAT子程序,建立考慮漸進損傷的剪切破壞分析模型預估了復合材料加筋壁板在剪切載荷下的失效模式和失效載荷。Mallela和Upadhyay[2]采用有限元對承受面內剪切載荷下的復合材料加筋板進行了線性屈曲分析,但并沒有進一步的后屈曲研究。Ambur等[3]采用ABAQUS建立了帶初始幾何缺陷的有限元模型,數值模擬了受面內剪切載荷下的復合材料加筋板后屈曲行為,并與試驗進行了對比分析。田偉偉等[4]通過全場光學形貌掃描的方法得到了剪切載荷下鋁合金加筋壁板屈曲模態的演變過程、光學離面位移數字結果,但僅限于試驗測量,未通過有限元預測其破壞載荷。劉存和趙謀周[5]論述了加筋壁板面內剪切試驗方法,通過試驗得到了加筋壁板剪切后屈曲破壞載荷,僅用有限元分析了應力分布,未預測其后屈曲承載能力。王平安等[6]利用試驗得到了剪切載荷下的復合材料加筋壁板屈曲特性,但僅通過有限元模擬與試驗初始屈曲載荷進行了對比分析。馮宇等[7]用工程算法和有限元法對剪切載荷下的復合材料加筋壁板蒙皮初始屈曲載荷進行了計算,與試驗結果吻合較好,但未涉及有限元預測其后屈曲承載能力。李愛環和支晗[8]通過PATRAN & NASTRAN分析剪切載荷下不同加筋尺寸對壁板蒙皮屈曲的影響,得到了屈曲載荷與腹板高度變化曲線,但沒有分析其后屈曲承載能力。

在安全設計的前提下盡可能地減輕結構重量是飛機設計工程師一直努力的方向,而準確預測結構的承載能力直接關系飛機的安全性和減重設計。文中基于MSC.PATRAN軟件對加筋壁板剪切試驗件進行仿真建模,運用NASTRAN SOL105進行線性屈曲求解,采用一致缺陷模態法,運用FILED場將屈曲結果作為擾動引入后屈曲模型,考慮了加筋壁板后屈曲的大擾度和塑性效應,調用MARC非線性求解器對結構進行后屈曲計算,得到載荷施加點全過程的載荷-位移曲線、后屈曲模態、極限載荷及破壞模式,通過試驗對比驗證了該方法的準確性,在工程上為同類型加筋壁板剪切后屈曲承載能力的預測提供了方法和技術支持。

1 計算模型

1.1 加筋壁板剪切試驗件

加筋壁板剪切試驗件的考核區長度為700 mm,試驗件的考核區由4個相同的長桁單元組成,長桁單元剖面如圖1所示??己藚^四周伸出過渡區用于與試驗夾具相連,過渡區厚度為10.5 mm,寬度為150 mm,拉力為P1。試驗件材料選用7150-T7751,其基本力學性能參數[9]如下:彈性模量Ec=73 700 MPa,屈服應力σ0.2t=524 MPa,極限強度σb=565 MPa,泊松比μ=0.33。

具體剖面參數如表1所示。其中:t為長桁腹板厚度;t1和b分別為長桁自由凸緣的厚度和寬度;t2和b1分別為長桁底邊凸緣的厚度和寬度;h為長桁高度;R和R1為倒圓半徑。

為了監測結構在剪切載荷下的后屈曲變形,在加筋壁板蒙皮和長桁的關鍵部位布置了應變片,試驗件構型及應變測量點布置如圖2所示,共142個測量點,長桁上布置單片,蒙皮表面布置花片,且正反面鏡像分布。有限元模擬與試驗對比分析所取的關鍵部位應變均來源于圖2。

圖1 試驗件剖面圖Fig.1 Cross-section of specimen

表1 試驗件的剖面參數Table 1 Parameters of cross-section of specimen

參數蒙皮厚度/mm長桁參數/mmtt1t2bb1hRR1數值3.53.53.51.0282061.765

圖2 試驗件應變片布置圖Fig.2 Layout diagram of strain gauge of specimen

1.2 有限元模型

有限元建模時,蒙皮和長桁的模擬采用四節點殼單元,考慮整體加筋壁板,長桁和蒙皮共節點。為了更好地模擬蒙皮的局部屈曲,采用細化的有限元模型。與壓縮屈曲波相比,剪切屈曲波的波長較小,故在剪切載荷下需要足夠細密的有限元網格來模擬剪切失穩波,每半波至少采用5個節點。加筋壁板剪切后屈曲破壞過程包含材料的彈性和塑性行為,為了準確地預測和模擬結構這一行為,采用材料真實的本構關系。鋁合金7150-T7751板材[9]的應力-應變曲線如圖3所示。

在MSC.PATRAN中定義材料彈塑性本構關系時,先輸入彈性變形參數,定義材料的本構為彈性模型;再將該材料的本構模型指定為彈塑性,非線性應力應變數據通過FIELD調用圖3中逐點輸入的坐標。

采用類型為RBE2的多點約束單元模擬夾具與試驗件連接。模型中處于對角拉伸方向上的加載點和約束點為主動節點,加載點約束垂向(Z向)位移,約束點約束蒙皮面內的X和Y向位移以及繞Z軸的轉動位移。夾具與試驗件用2排螺栓連接,其中心線對應的有限元節點為從動節點,分別約束其X或Y向自由度,同時約束其繞Z軸的轉動位移。以便使模型達到試驗狀態下的約束。模型施加的載荷為強迫位移。有限元模型的加載和約束條件如圖4所示。

圖3 7150-T7751板材應力-應變曲線[9]Fig.3 Stress-strain curve of 7150-T7751 plate[9]

圖4 有限元模型加載及約束圖Fig.4 Loading and constrain conditions of finite element model

2 有限元模型的求解策略

建立能夠真實反映結構抗屈曲性能的計算模型是進行加筋壁板剪切后屈曲計算的前提。然而求解要能夠跟蹤整個結構的平衡路徑,能夠跨越屈曲分叉點或極限強度點,又能夠追蹤整個失穩過程中實際的載荷-位移關系而獲得結構失穩前后的全部信息,進而獲得較為準確的結構后屈曲承載能力,需要制定求解策略。

首先要基于計算模型采用SOL105進行線性屈曲分析,獲取結構的屈曲模態。其次將1階屈曲模態乘以擾動系數,作為結構的初始缺陷,即一致模態缺陷。然后運用FILED將初始缺陷加入模型。在此基礎上采用基于Newton-Raphson迭代的弧長法(Arc-Length Method)進行非線性有限元計算。

2.1 線性屈曲分析

線性屈曲分析通過提取使線性系統剛度矩陣奇異的特征值來獲得結構的臨界失穩載荷及失穩模態。MSC.NASTRAN軟件SOL105為線性屈曲分析求解序列,適用于結構為小變形、單元應力必須是彈性(應力應變關系為線性)等范圍內的求解,故對于非線性材料和大變形結構,SOL105計算得到的線性屈曲載荷和實際的臨界載荷相差較大。根據最小勢能原理,在僅考慮材料線彈性情況時節點外載荷矩陣F與節點位移矩陣u的關系為

KTu=F

(1)

式中:KT為切線剛度矩陣。

2.2 幾何缺陷假定

由于制造工藝等因素,加筋壁板不可避免地存在一定的幾何缺陷,如初始彎曲、初始偏心等。幾何缺陷的實質是以附加應力的形式促使部分結構截面提前進入屈服,影響其承載能力。缺陷的存在使得結構的失穩一般都呈現為彈塑性狀態,對其承載能力的分析應該是計及幾何缺陷的彈塑性疊加幾何非線性的雙重非線性計算問題。

一致缺陷模態法[10]認為結構的最低階臨界點所對應的屈曲模態為結構的最低階屈曲模態,結構按該模態變形將處于勢能最小狀態,所以對于實際結構來說,在載荷施加的最初階段即有沿著該模態變形的趨勢,如果結構的缺陷分布形式恰好與最低階屈曲模態相吻合,這將對其受力性能產生最不利影響。一致缺陷模態法就是用最低階屈曲模態來模擬結構的最不利幾何缺陷的分布,并認為最低階屈曲模態所對應的臨界荷載就是該結構的最小臨界荷載。采用這種缺陷形式主要面臨2個問題:① 如何確定初始缺陷的模式;② 如何確定缺陷的最大值。Xu和Soares[11]利用有限元方法計算了5種不同構型的試驗件在傳感器測量的真實缺陷值(實測缺陷法)和認為假定缺陷(一致缺陷模態法)下的加筋壁板結構的承載能力,并與試驗結果進行對比,結果表明引入實測缺陷計算得到的破壞載荷和破壞模式與試驗結果吻合更好。萬春華等[12]采用考慮幾何和材料雙重非線性的弧長法分析了初始幾何缺陷對加筋結構后屈曲的影響,結果表明:當初始缺陷較小時,對結構的承載能力影響很小,但當初始缺陷增大至蒙皮厚度0.1倍以上時,載荷-位移曲線向下偏移。另外,相同模式、不同尺度的初始幾何缺陷會得到不同的破壞模式。推薦在無實測結構初始幾何缺陷情況下,以線性屈曲特征值法計算得到的1階屈曲模態為初始缺陷模式,初始缺陷的最大位移值取蒙皮厚度的0.005~0.1倍。

文中通過SOL105線性屈曲特征值分析,得到加筋壁板剪切特征值屈曲模態。進行非線性屈曲分析前,將加筋壁板1階屈曲模態矢量歸一化,乘以缺陷基矢量10-3得到缺陷偏移矢量,通過SPCD模型數據卡來施加強迫節點位移,以更新單元節點坐標的形式將考核區的缺陷引入完善加筋壁板結構中。

2.3 非線性屈曲分析

在增量加載過程中,用包含加載過程中所有非線性影響的剛度矩陣來評定屈曲特征值,求得的失穩載荷更接近結構的真實臨界載荷值。通過增量定義的非線性屈曲分析,可以同時考慮材料和幾何非線性對結構穩定性的影響,選中材料本構關系后,在增量加載分析中激活LGDISP等選項,程序會自動地在幾何非線性的基礎上計算材料非線性對剛度矩陣的貢獻[13]。

在考慮材料和幾何雙重非線性的后屈曲分析中[14],式(1)中的切線剛度矩陣KT應改為

KT=Kep+Kgs+Kld-KR

(2)

式中:Kep為結構彈塑性矩陣;Kgs為幾何剛度矩陣;Kld為大位移剛度矩陣;KR為載荷矯正矩陣。

在每一子步的切線剛度求解完成后,采用基于Newton-Raphson迭代的弧長法進行子步位移增量計算,Newton-Raphson迭代公式為

φ(un)=KTun-F

(3)

(4)

un+1=un+Δunn=1,2,…

(5)

式中:φ(un)為第n步的剩余載荷;KTn為第n步的切線剛度矩陣;un和un+1為第n和n+1步的位移向量;Δun為位移增量。

由于Newton-Raphson方法無法計及剛度系數等于零的情況(通常對應于極限載荷點),而對于非線性后屈曲求解來說,主要目的在于探究結構的極限載荷點和追蹤平衡路徑,而弧長法在增量步中加入弧長約束來實現剛度系數過零的力學狀態。具體弧長迭代公式為

(6)

式中:A0為常量;Δλi為載荷因子;ΔI為控制弧長。

弧長法分析屈曲問題不僅考慮剛度奇異失穩點附近的平衡,而且是通過追蹤整個失穩過程中實際的載荷-位移關系,獲得結構失穩前后的全部信息?;¢L參數選擇得當,可完整跟蹤包括屈曲非穩定段的載荷-位移曲線。后屈曲模型采用NLPARM卡定義非線性分析,載荷被分為52個等增量,采用ITER方法控制切線剛度修正,矩陣修正之前的迭代次數為40,每個載荷增量的總迭代限為1 000。采用NLPCI卡定義非線性靜態分析中弧長增量求解策略,弧長法類型為MRIKS,最小弧長比為0.25,最大弧長比為4.0,期望收斂的迭代次數為40,每步最大迭代次數為500。采用弧長法進行每個載荷步的位移及切線剛度求解,每個子步重復迭代直到滿足收斂條件,最后根據每個載荷步內的計算結果繪制載荷-位移曲線,依據曲線得到加筋壁板結構的雙重非線性后屈曲極限載荷。

3 后屈曲承載能力預測與試驗驗證

3.1 屈曲和破壞載荷

對幾何參數如表1所示的3件加筋壁板剪切試驗件進行了試驗,屈曲載荷試驗結果與有限元結果的比較如表2所示,破壞載荷試驗結果與有限元結果的比較如表3所示,同時給出了試驗載荷平均值。表中ε為有限元計算值相對于試驗值的誤差,即

ε=(PFEA-PTEST)/PTEST×100%

(7)

式中:PTEST為試驗值;PFEA為有限元計算值。

可見,采用有限元計算的加筋壁板剪切后屈曲的初始屈曲載荷和破壞載荷與試驗值誤差均在5%以內,有限元仿真計算得到的初始屈曲載荷與試驗載荷平均值誤差為1.25%,破壞載荷平均值誤差為2.4%。

表2屈曲載荷試驗結果與有限元結果比較

Table2Comparisonofbucklingloadsbetweentestresultsandcalculationresults

試驗件序號屈曲載荷/kN試驗值計算值誤差ε/%113461375 2.1214181375-3.0313111375 4.8平均值13581375 1.25

表3破壞載荷的試驗與有限元結果比較

Table3Comparisonoffailureloadsbetweentestresultsandcalculationresults

試驗件序號破壞載荷/kN試驗值計算值誤差ε/%1152515894.12154415892.93158515890.2平均值155115892.4

3.2 載荷-位移曲線

利用SOL105模塊進行加筋壁板剪切試驗件有限元模型的線性屈曲分析,計算得到的屈曲模態如圖5所示,1階失穩特征值的臨界失穩系數為0.98,計算得到屈曲載荷為1 375 kN,即在98%的設計載荷時加筋壁板發生初始屈曲。此外,加筋壁板的蒙皮發生局部屈曲,關于對角拉伸線呈現反對稱分布的剪切屈曲波,單側出現由4個屈曲半波構成的2個完整屈曲波形,4根長桁均未發生局部屈曲,處于波形的節點上,起到了隔波的作用。

提取加筋壁板剪切試驗件考核區結構各點彈性屈曲位移,按照一致缺陷模態法將所得到位移按比例因子10-3進行折算,利用MSC.PATRAN中位移的FIELD定義模式施加,引入擾動的加筋壁板等值線圖與初始屈曲模態形貌完全一致。在求解設置中開啟大變形(LGDISP)考慮幾何非線性,引入材料的彈塑性曲線,在弧長法中設置合適的載荷步和弧長增量,進行非線性迭代計算。提取載荷及位移的計算結果,繪制載荷-位移曲線如圖6所示。

圖5 加筋壁板的初始屈曲模態Fig.5 Initial buckling modes of stiffened panel

由加筋壁板加載端載荷-位移曲線可見,曲線在初始階段呈線性,加載到1 364 kN時,曲線斜率略微下降,此時結構出現初始屈曲,結構面內剪切剛度發生變化,進入后屈曲階段,采用非線性分析中基于曲線獲得的屈曲載荷與特征值線性屈曲分析結果1 375 kN基本一致,表明有限元分析模型對于出現屈曲前的線性加載階段的模擬準確。而實際試驗中剪切試驗件的初始屈曲載荷平均值為1 358 kN,相對誤差為1.25%。隨著載荷繼續增加,在出現拐點前曲線斜率小幅下降,最終曲線在1 589 kN時達到頂點,即結構的極限載荷。隨后,曲線斜率為負值,位移繼續增加,載荷下降,結構喪失承載能力。實際試驗中剪切試驗件的破壞載荷平均值為1 551 KN,相對誤差為2.4%。

非線性有限元仿真計算加筋壁板試驗件的初始屈曲載荷為破壞載荷的86.5%,而實際試驗的初始屈曲載荷平均值為破壞載荷平均值的87.5%??梢?.5 mm厚的加筋壁板剪切試驗件在發生剪切屈曲后仍能繼續承載,直到發生結構破壞。

圖6 加筋壁板加載端載荷-位移曲線Fig.6 Load-displacement curve of loading end

3.3 破壞模式

為進一步分析破壞過程,對加筋壁板剪切試驗件考核區蒙皮在加載過程中的位移變化過程進行分析。破壞模式是蒙皮先發生局部屈曲,隨著載荷增加,局部波形擴大并穿越長桁導致加筋板破壞。有限元模擬破壞時刻位移云圖如圖7所示,試驗第1#件的破壞形貌如圖8所示。文中僅對第1#件的試驗數據與有限元模擬結果進行對比分析,可見兩者的破壞部位和破壞模式一致。

結合試驗錄像過程及試驗件的最終破壞形貌對剪切屈曲波描述如下:在試驗件對角拉伸的兩個角區,一個向蒙皮側凸起,另一個則向長桁側凹陷,同時角區產生小的屈曲波;在對角拉伸線的兩側交替分布著凸起和凹陷的屈曲波,波長和波深隨載荷的增大而變化;隨著載荷的繼續增加,屈曲波產生的面外變形持續增加,在大撓度和塑性效應的作用下,長桁和蒙皮開裂,進而蒙皮破壞。從試驗件的破壞形貌和有限元仿真破壞變形圖可見,兩者最終破壞模式基本一致。

圖7 有限元模擬破壞時刻變形云圖Fig.7 Deformation morphology of FE simulation failure

圖8 剪切試驗件的試驗破壞形貌Fig.8 Failure morphology of shear specimen in test

3.4 破壞過程

剪切試驗件處于純剪切狀態,22#和122#應變片位于壁板考核區中心點上,且在加載對角線上,25#和125#應變片位于考核區中心點附近??己藚^中心點長桁側22#和蒙皮側122#應變片以及25#和125#應變片的應變隨加載步變化曲線的計算值與試驗值的對比如圖9所示,圖中,“試驗”為實測應變花合成最大剪應變(試驗-22代表22#、23#、24#應變片合成的最大剪應變,其他類同),“有限元”為有限元仿真最大剪應變。剪切試驗件角區長桁側43#和蒙皮側143#應變片的應變隨加載步變化曲線的計算值與試驗值的對比如圖10所示。

從圖9可見,試驗件中心點蒙皮處(應變花22#、122#)的載荷-應變曲線均在起始階段呈現線性變化,在載荷增加到1 346 kN左右時,試驗與有限元分析得到的載荷-應變曲線均出現了明顯的拐折,據此推斷該處產生了較大的面外位移,即發生了局部屈曲。其中,屈曲點的判斷為背對背的應變平均值隨載荷變化曲線斜率發生明顯變化的拐折點。圖10所示的加筋壁板角區有限元與試驗載荷-應變曲線的屈曲載荷為1 485 kN,與中心點相比,明顯較晚,這與角區邊界支持剛度較強有關。圖9和圖10顯示有限元分析和試驗實測曲線在線性段吻合很好,結構發生屈曲后,進入材料非線性和幾何非線性交織狀態,有限元分析和試驗實測曲線變化趨勢基本一致,但存在一些偏差,這對采用有限元準確模擬加筋壁板剪切試驗件的后屈曲過程提出了更高要求。有限元仿真得到的破壞模式與試驗一致,關鍵部位的載荷變化曲線與試驗數據吻合較好,承載能力計算值與試驗值基本一致,達到了采用NASTRAN預測加筋壁板剪切后屈曲承載能力的目的。

圖9 加筋壁板中心有限元分析與試驗所得載荷-應變曲線對比Fig.9 Comparison of load-strain curves obtained by FEA and test at stiffened panel center point

圖10 加筋壁板角區有限元分析與試驗所得載荷-應變曲線對比Fig.10 Comparison of load-strain curves obtained by FEA and test at stiffened panel corner point

4 結 論

通過對加筋壁板受面內剪切載荷全過程的分析及試驗驗證,得到以下結論:

1) 采用MSC.NASTRAN模擬了加筋壁板受面內剪切的全過程,包括其載荷-位移曲線的初始線性、結構后屈曲承載能力和結構的破壞模式。

2) 基于文中有限元模擬的屈曲載荷、破壞載荷、破壞模式和典型位置的應變與試驗對比結果,印證了有限元預測模型的有效性。

3) 引入初始缺陷,考慮了雙重非線性的NASTRAN有限元預測模型,工程上可作為加筋壁板強度預估的一種有效手段。

猜你喜歡
模態有限元結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
國內多模態教學研究回顧與展望
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
磨削淬硬殘余應力的有限元分析
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 不卡无码h在线观看| 欧美一级一级做性视频| 欧美成人二区| 永久免费无码成人网站| 日韩黄色精品| 夜夜操天天摸| 国产区免费精品视频| 伊人色天堂| 亚洲av日韩av制服丝袜| 麻豆精品视频在线原创| 亚洲一区黄色| 亚洲人成网站日本片| 台湾AV国片精品女同性| 一区二区无码在线视频| 91在线无码精品秘九色APP| 婷婷午夜天| 九色综合视频网| 91口爆吞精国产对白第三集| 国产国拍精品视频免费看| 国产精品蜜芽在线观看| 一级香蕉视频在线观看| 无码福利视频| 99激情网| 久久精品国产一区二区小说| 国产91久久久久久| 性色一区| 欧美日韩精品一区二区视频| 国产精品丝袜在线| 午夜精品久久久久久久2023| 国产在线拍偷自揄拍精品| 国产午夜无码片在线观看网站| 99无码熟妇丰满人妻啪啪| 欧美一级片在线| 久久久四虎成人永久免费网站| 99偷拍视频精品一区二区| 成人毛片免费观看| 高清精品美女在线播放| 丰满人妻一区二区三区视频| 亚洲一区二区精品无码久久久| 91免费片| 国产天天射| 亚洲另类国产欧美一区二区| 国产亚洲精品yxsp| 91无码视频在线观看| 99久久精品久久久久久婷婷| 亚洲国产av无码综合原创国产| 日本午夜精品一本在线观看 | 国产精品人成在线播放| 亚洲一道AV无码午夜福利| 日韩中文无码av超清| 午夜福利视频一区| 日韩区欧美国产区在线观看| 久久婷婷色综合老司机| 日韩第九页| 国产99在线| 国产精品va免费视频| 国产浮力第一页永久地址| 久久久久亚洲Av片无码观看| 国产精品微拍| 亚洲精品欧美日韩在线| 成年A级毛片| 亚洲成a∧人片在线观看无码| 久久伊伊香蕉综合精品| 亚洲av片在线免费观看| 黄色网页在线观看| 夜夜操天天摸| 另类欧美日韩| 国产裸舞福利在线视频合集| 亚洲精品第一页不卡| 九九九精品成人免费视频7| 六月婷婷精品视频在线观看| 国产成人精品视频一区二区电影| 国产你懂得| 99re视频在线| 国产亚洲高清视频| 亚洲色成人www在线观看| 成人中文在线| 99精品免费在线| 亚洲综合婷婷激情| www中文字幕在线观看| 国产熟睡乱子伦视频网站| 毛片基地视频|