亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        爆炸地沖擊載荷計算方法綜述

        2023-11-02 12:59:38周岸峰李道奎周仕明
        工程爆破 2023年5期
        關(guān)鍵詞:特征頻率峰值沖擊

        周岸峰,李道奎,周仕明,周 旋,崔 達

        (1.國防科技大學(xué)空天科學(xué)學(xué)院,長沙 410073;2.空天任務(wù)智能規(guī)劃與仿真湖南省重點實驗室,長沙 410073)

        爆炸地沖擊主要是指地下或地面發(fā)生爆炸時周圍土石介質(zhì)經(jīng)沖擊波的壓縮和推動產(chǎn)生以加速度、速度、位移為表征的強烈地運動,也稱爆炸地震或爆炸沖擊振動[1-2]。與天然地震相比,爆炸地震對土石介質(zhì)的作用以球?qū)ΨQ壓縮為主,主要成分為縱波,震源時間函數(shù)表現(xiàn)為單脈沖形式且能量急速釋放的物理過程[3],具有幅值大、衰減快、頻率高、持時短等特點[2]。在隧道、核電站和人防工程等建筑設(shè)施的設(shè)計中,通常需要考慮爆破施工、爆炸事故以及導(dǎo)彈攻擊等極端環(huán)境[4-5];為了對結(jié)構(gòu)在極端環(huán)境下的動態(tài)響應(yīng)與損傷進行預(yù)測和評估,必須實現(xiàn)爆炸地沖擊載荷的準(zhǔn)確計算。

        目前,國內(nèi)外研究爆炸地沖擊載荷的主要方法有理論分析、實驗研究和數(shù)值仿真[6]。理論分析通常采用比較理想化的近似模型,只能處理一些簡單問題。爆炸實驗具有瞬時性和破壞性,對測量精度要求非常高,一些超常規(guī)的爆炸實驗因國際公約無法進行。數(shù)值仿真則具有很強的靈活性,可以處理許多復(fù)雜的問題,甚至模擬實驗達不到的條件,且方便進行參數(shù)靈敏度分析,從而指導(dǎo)理論分析和實驗研究。因此,大部分研究將上述方法相互結(jié)合、相互驗證,提出了許多關(guān)于爆炸地沖擊載荷的計算方法,包括爆炸地沖擊峰值載荷的預(yù)測[2,7],反應(yīng)譜的設(shè)計[8-9]以及時域過程的模擬[10-15]。

        由于現(xiàn)有爆炸地沖擊載荷的計算方法種類繁多,各方法的適用條件和應(yīng)用范圍不盡相同,為此對現(xiàn)有研究成果進行梳理和匯總,對爆炸地沖擊峰值載荷的經(jīng)驗公式進行比較分析,對爆炸地沖擊反應(yīng)譜的設(shè)計以及時域過程的模擬進行歸納總結(jié),從而為進一步的相關(guān)研究提供參考與建議。

        1 爆炸地沖擊峰值載荷的預(yù)測

        爆炸地沖擊峰值載荷主要指峰值位移Dm、峰值速度vm和峰值加速度Am,其大小與爆炸當(dāng)量、爆心距、爆源埋深、爆炸傳播的地質(zhì)條件等參數(shù)密切相關(guān)。目前常用經(jīng)驗公式和人工智能算法來表示上述參數(shù)之間的復(fù)雜關(guān)系,從而實現(xiàn)對爆炸地沖擊峰值載荷的預(yù)測。

        1.1 經(jīng)驗公式

        爆炸地沖擊峰值載荷的經(jīng)驗公式一般為特定的基函數(shù)形式,以爆炸當(dāng)量Q和爆心距R為主要變量,將其它參數(shù)的影響作用歸納到比例系數(shù)k和衰減系數(shù)α中去,然后通過實驗數(shù)據(jù)識別出k和α。

        在礦產(chǎn)開采、隧道開挖等國民經(jīng)濟建設(shè)領(lǐng)域,峰值速度是工程爆破安全評估的重要準(zhǔn)則[16],相關(guān)經(jīng)驗公式最多,一些比較經(jīng)典、常用的峰值速度經(jīng)驗公式[2]如表1所示。

        表1 峰值速度經(jīng)驗公式

        美國陸軍編寫的《常規(guī)武器防護設(shè)計原理》[17]給出的地沖擊峰值載荷的經(jīng)驗計算公式為

        (1)

        式中:C為地沖擊耦合系數(shù),封閉爆炸時C=1,其余情況取值可參考文獻[1];cp為縱波波速;α為衰減系數(shù);g為重力加速度。cp、α的取值如表2所示。

        表2 巖土特性參數(shù)

        《美國空軍防護結(jié)構(gòu)設(shè)計與分析手冊》[10]給出了不同條件下核爆地沖擊峰值載荷的估算方法

        (2)

        式中:kd、kv和ka為與地質(zhì)條件和爆炸方式相關(guān)的常系數(shù)(見表3)。

        表3 公式(2)相關(guān)系數(shù)

        此外,Murphy和Lahoud[18]研究了99起地下核實驗的地沖擊振動記錄,分析了爆心距為0.25~600 km之間的500多個觀測點的數(shù)據(jù),采用以下基函數(shù)形式來擬合峰值載荷經(jīng)驗公式

        (3)

        式中:Y為待擬合的物理量(峰值位移、峰值速度或峰值加速度);k、m、n為待定常系數(shù),具體如表4所示。

        表4 Murphy-Lahoud經(jīng)驗公式相關(guān)系數(shù)

        綜合比較上述經(jīng)驗公式可知,雖然各個公式的構(gòu)造形式接近,但由于不同學(xué)者對地質(zhì)條件歸類、實驗數(shù)據(jù)處理和基函數(shù)定義的方法不盡相同,導(dǎo)致各經(jīng)驗公式的系數(shù)取值存在差異,實際計算結(jié)果的離散性較大。因此,在一些實際工程中,仍然需要根據(jù)爆炸現(xiàn)場的實測數(shù)據(jù),對傳統(tǒng)經(jīng)驗公式的常系數(shù)重新進行識別,或增加其它變量進行多元回歸分析,從而提出符合特定工程環(huán)境的爆炸地沖擊峰值載荷預(yù)測公式[19-20]。

        1.2 人工智能算法

        爆炸地沖擊峰值載荷相關(guān)參數(shù)之間具有很強的非線性,因而采用具體的函數(shù)形式來描述這種非線性關(guān)系是十分困難的。隨著計算機技術(shù)的發(fā)展,人工智能(artificial intelligence,AI)算法由于非常適合處理復(fù)雜的非線性關(guān)系,也逐漸應(yīng)用于爆破工程領(lǐng)域[7]。研究表明,AI算法在預(yù)測爆炸地沖擊峰值速度時比經(jīng)驗公式和多元回歸模型的性能都要好,但由于算法本身的局限性、訓(xùn)練樣本來源的單一性以及建模的差異性(見表5),故AI算法還只能針對特定的工程爆破問題進行峰值速度的預(yù)測。

        表5 爆炸地沖擊峰值載荷預(yù)測AI算法

        未來的研究還需要:①理清影響爆炸地沖擊峰值載荷的主要因素,如爆炸當(dāng)量(或炸藥質(zhì)量)、爆心距、地質(zhì)條件、埋深、炮孔設(shè)置等,建立適用于不同工程問題的AI算法預(yù)測模型;②搜集不同地質(zhì)條件的爆炸現(xiàn)場測量數(shù)據(jù),建立一個更廣泛、更具多樣性的數(shù)據(jù)庫,以對AI算法預(yù)測模型進行充分的訓(xùn)練,擴大其適用范圍。

        2 爆炸地沖擊反應(yīng)譜的設(shè)計

        與峰值載荷這類單一判據(jù)相比,反應(yīng)譜由于考慮了結(jié)構(gòu)自身固有動力特性對地沖擊響應(yīng)的放大作用,以此進行隔震設(shè)計一般來說更加安全。目前常用于爆炸地沖擊反應(yīng)譜設(shè)計的方法主要是分量譜法、三系數(shù)法和特征頻率法[8,9,15]。其中,分量譜法是美國《抗核武器結(jié)構(gòu)設(shè)計手冊》提出的方法,繪制圖譜的過程十分繁瑣,需要給定參數(shù)較多且難以確定,實際使用存在一定困難,而三系數(shù)法和特征頻率法相對簡單方便,下面主要介紹這兩種方法。

        2.1 三系數(shù)法

        三系數(shù)法是根據(jù)爆炸地沖擊峰值加速度、峰值速度、峰值位移,并基于大量實測加速度信號的歸一化反應(yīng)譜計算結(jié)果,確定加速度放大系數(shù)Ka、速度放大系數(shù)Kv、位移放大系數(shù)Kd,并按下式確定最大加速度響應(yīng)Sa、最大速度響應(yīng)Sv和最大位移響應(yīng)Sd:

        Sa=KaAm,Sv=KvVm,Sd=KdDm

        (4)

        然后根據(jù)Sa、Sv、Sd在三對數(shù)坐標(biāo)圖上分別作加速度、速度、位移坐標(biāo)的等值線(見圖1中實線),從而構(gòu)成爆炸震動現(xiàn)場環(huán)境的反應(yīng)譜包絡(luò)線。

        圖1 三系數(shù)法Fig.1 Three-coefficient method

        三系數(shù)法形式簡潔,構(gòu)建方便,但由于加速度、速度、位移峰值的經(jīng)驗公式結(jié)果在運動學(xué)上難以協(xié)調(diào),容易導(dǎo)致反應(yīng)譜形狀畸變,出現(xiàn)三角形圖譜和沿頻率軸偏移等現(xiàn)象(見圖1中虛線)。

        2.2 特征頻率法

        特征頻率法是根據(jù)爆炸地沖擊加速度的衰減規(guī)律、加速度放大系數(shù)及反應(yīng)譜的特征頻率,通過量綱分析及統(tǒng)計得出的反應(yīng)譜設(shè)計方法,無需速度、位移參數(shù)即可在三對數(shù)坐標(biāo)系中同時繪制加速度、速度和位移反應(yīng)譜,其中加速度反應(yīng)譜(見圖2中實線)可以使用分段函數(shù)表示為[15]

        圖2 特征頻率法Fig.2 Specific frequency method

        式中:f為頻率;F1和F2分別定義為第一和第二特征頻率。特征頻率計算公式和加速度放大系數(shù)可參考文獻[8]和文獻[9],此處不一一列舉。

        年鑫哲等[15]認為式(5)確定的設(shè)計反應(yīng)譜在高頻部分存在較大誤差,以此合成爆炸地沖擊時程曲線時高頻分量將被不合理地放大;當(dāng)頻率較高(剛度較大)時,系統(tǒng)的最大動力響應(yīng)是趨近激勵峰值的,即動力放大系數(shù)應(yīng)趨近于1(見圖2中虛線)。因此,對式(5)設(shè)計反應(yīng)譜曲線的水平階段(即圖2中f>F3時)進行如下修正:

        (6)

        式中:F3和F4為按照F1、F2的方式定義的第三和第四特征頻率。

        綜上可知,頻率特征值法由于控制了特征頻率F1和F2,且只需要計算峰值加速度,故克服了三系數(shù)法容易產(chǎn)生畸變的缺點[8]。但現(xiàn)有的特征頻率經(jīng)驗公式僅僅來源于少量實測數(shù)據(jù)的回歸分析,置信度不高;不同學(xué)者采用的經(jīng)驗公式基函數(shù)不同,導(dǎo)致計算結(jié)果差異巨大,通用性不強;而且文獻[15]也沒有給出計算特征頻率值的詳細說明。因此,為建立適用范圍更廣的爆炸地沖擊反應(yīng)譜模型,還需結(jié)合爆炸現(xiàn)場的實測數(shù)據(jù)對反應(yīng)譜的圖譜特征進行合理的統(tǒng)計分析,一方面可以對現(xiàn)有特征頻率的經(jīng)驗計算公式進行修正和改進,另一方面可以嘗試采用人工智能算法建立合適的特征頻率的預(yù)測模型。

        3 爆炸地沖擊時程的模擬

        爆炸地沖擊的峰值載荷和反應(yīng)譜,應(yīng)用于工程設(shè)計時十分方便快捷,但不能分析結(jié)構(gòu)在地沖擊作用下的動態(tài)響應(yīng)規(guī)律,而爆炸實驗研究的成本太高、危險過大;因此,準(zhǔn)確模擬地沖擊載荷的時程曲線對結(jié)構(gòu)的動力響應(yīng)分析和減振設(shè)計十分必要。目前,計算爆炸地沖擊時程過程的方法主要有經(jīng)驗曲線、數(shù)值仿真技術(shù)以及非平穩(wěn)隨機過程方法。

        3.1 經(jīng)驗曲線

        《常規(guī)武器防護設(shè)計原理》給出的爆炸地沖擊速度時程經(jīng)驗曲線公式為[17]

        (7)

        式中:t為時間;β為修正系數(shù),一般可取0.4,t0為沖擊波到達時間,t0=R/cp;式(7)對時間求導(dǎo)即為加速度時程曲線。

        嚴(yán)東晉等[11]對《美國空軍防護結(jié)構(gòu)設(shè)計與分析手冊》提出的兩類典型波形進行了曲線擬合,給出了對應(yīng)的加速度時程經(jīng)驗公式;其中,類型I加速度波形的表達式為

        (8)

        類型II加速度波形的表達式為

        (9)

        分析可知,文獻[17]的經(jīng)驗公式十分簡單,卻忽略了地沖擊載荷的爬升過程;而嚴(yán)東晉等[11]提出的經(jīng)驗公式雖然能快速確定不同工況條件下的地沖擊近似波形,并為隔震設(shè)計提供一定依據(jù),但由于實際地質(zhì)條件的復(fù)雜性,上述近似波形尚不能完全反應(yīng)實際爆炸地沖擊時域過程的復(fù)雜變化[26],容易丟失一些高頻信息。因此,將上述經(jīng)驗公式確定的近似地沖擊波作為結(jié)構(gòu)動力分析的激勵輸入,本身就存在較大誤差。

        3.2 數(shù)值仿真技術(shù)

        爆炸沖擊波在巖土介質(zhì)中的傳播過程是一類典型的非線性瞬態(tài)動力學(xué)問題。理論分析可以清晰地認識問題的本質(zhì),但只有在巖土動態(tài)本構(gòu)模型、載荷和邊界條件都比較簡單時,才可能得到一些解析解[26-27];隨著數(shù)學(xué)和計算機技術(shù)的飛速發(fā)展,有限元法、離散元法、物質(zhì)點法、光滑粒子流體動力學(xué)法等數(shù)值方法被廣泛應(yīng)用于爆炸沖擊問題的研究[6,28-32]。

        由于巖土中的爆炸問題比空氣和水中的爆炸復(fù)雜得多,為提高數(shù)值仿真精度,除了開發(fā)高精度的算法格式以外,還需要建立科學(xué)合理的巖土介質(zhì)動態(tài)本構(gòu)模型。目前,常用的巖土動態(tài)本構(gòu)模型有泡沫模型、帽蓋模型等彈塑性模型,以及HJC模型、TCK模型等損傷斷裂模型[12]。此外,北京理工大學(xué)的寧建國團隊[12]還給出了基于細觀力學(xué)理論的混凝土材料動態(tài)本構(gòu)模型以及微裂紋演化與匯合的本構(gòu)模型。然而,上述方法和本構(gòu)模型都建立在連續(xù)介質(zhì)力學(xué)的框架內(nèi),而實際地表的巖土介質(zhì)存在大量斷層、裂隙等特征,爆炸后巖土塊的旋轉(zhuǎn)以及相互間的滑移、擠壓和分離會導(dǎo)致巖土介質(zhì)位移場的間斷[5],因此,在未來的研究中,仍然需要努力解決以下問題:①多介質(zhì)界面相互作用的高精度計算[6,28,31];②爆炸沖擊波強間斷現(xiàn)象的高精度計算[29-30];③高速、高壓、高溫等極端條件下材料動態(tài)本構(gòu)關(guān)系以及損傷演化規(guī)律的數(shù)學(xué)表征[32]。

        3.3 非平穩(wěn)隨機過程方法

        由于爆轟產(chǎn)物、氣體、巖土等介質(zhì)相互作用的復(fù)雜性,即使進行多次相同條件的爆炸實驗,同一觀測點記錄的地沖擊波一般也不會完全重現(xiàn),具有一定的隨機性。因此,可采用隨機振動分析的方法,將爆炸地沖擊時程按非平穩(wěn)隨機過程進行模擬[13-15]。目前常使用包絡(luò)函數(shù)、相位差譜或時變功率譜表征隨機過程的非平穩(wěn)特性,分述如下。

        1) 基于包絡(luò)線函數(shù)的非平穩(wěn)隨機過程模擬。假定爆炸地沖擊加速度時程為調(diào)幅非平穩(wěn)隨機過程為

        A(t)=L(t)X(t)

        (10)

        式中:L(t)為幅值包絡(luò)線,是一確定性函數(shù);X(t)為一均值為0的平穩(wěn)高斯隨機過程。

        關(guān)于幅值包絡(luò)線函數(shù)L(t),其數(shù)學(xué)模型較多,目前常用的有單指數(shù)型[33]、改進的單指數(shù)型[14]、雙指數(shù)型[34]以及分段函數(shù)型[35]。其中,單指數(shù)型包絡(luò)線函數(shù)忽略了幅值的增長過程,實際應(yīng)用較少。改進的單指數(shù)型和雙指數(shù)型包絡(luò)線函數(shù)不僅彌補了單指數(shù)型的不足,而且常系數(shù)物理意義明確,計算簡單方便;但兩者的衰減系數(shù)與幅值峰值時刻相關(guān),當(dāng)保持幅值峰值時刻相等時,兩種包絡(luò)線在衰減階段會出現(xiàn)明顯差異,故只適用于特定工況。分段函數(shù)型的平穩(wěn)段持續(xù)時間和衰減系數(shù)則具有一定的可調(diào)性,適用范圍比較大,但這些參數(shù)還缺乏明確的計算方法,需要根據(jù)實驗數(shù)據(jù)進行確定。

        對于平穩(wěn)隨機過程X(t)的模擬,目前常用的是Shinozuka和Jan[36]提出的加權(quán)振幅諧波疊加法,其表達式為

        (11)

        Δω=(ωu-ωl)/N,ωk=ωl+(k-1/2)Δω

        (12)

        式中:Sxx(ω)為目標(biāo)功率譜密度函數(shù);N為三角級數(shù)項數(shù);ωk為各離散化的子頻帶的中心圓頻率;φk為[0~2π]之間服從均勻分布的隨機相位,且φk與φj相互獨立;ωu、ωl分別為Sxx(ω)在正頻率范圍內(nèi)的上、下限,即X(t)的頻帶寬度為(ωu-ωl)。

        顯然,基于諧波疊加法進行非平穩(wěn)隨機過程模擬的先決條件是給定合適的目標(biāo)功率譜。國勝兵等[14]和周鵬等[37]提出了能夠反應(yīng)裝藥量和爆距影響的爆炸地沖擊功率譜模型,但該模型部分系數(shù)的物理意義不明確,必須借助實測地沖擊信號進行確定且計算過程比較繁瑣。楊佑發(fā)等[13]和年鑫哲等[15]則根據(jù)反應(yīng)譜轉(zhuǎn)功率譜的方法,實現(xiàn)了爆炸地沖擊信號的模擬,但目標(biāo)反應(yīng)譜的設(shè)計仍然依賴于實測地沖擊信號,因而不能普遍應(yīng)用。

        2)基于相位差譜的非平穩(wěn)隨機過程模擬。一些強天然地震(持時較短,能量集中的強震)的研究[38-40]表明,采用包絡(luò)線函數(shù)來反映地震動的非平穩(wěn)性具有一定的隨意性,且不能表征頻率非平穩(wěn)特性;實際地震動的非平穩(wěn)性是由其相位差譜決定的,而相位差譜與相位譜具有一一對應(yīng)的關(guān)系,即地震動的非平穩(wěn)性包含在其相位差譜之中。因此,基于相位差譜進行非平穩(wěn)隨機過程模擬時,無需進行強制調(diào)幅處理,即

        (13)

        式中:各參數(shù)的物理意義與上文相同,但隨機相位之間的差值應(yīng)符合特定分布規(guī)律。隨機相位之間的差值定義為

        (14)

        目前,針對不同地區(qū)的天然強震,不同學(xué)者對大量實測數(shù)據(jù)進行統(tǒng)計分析,得到了相關(guān)相位差譜模型,如正態(tài)分布、對數(shù)正態(tài)分布、均勻分布、Beta分布等[39-40],并實現(xiàn)了天然強震的非平穩(wěn)隨機過程模擬。但對于爆炸沖擊振動,目前尚無相位差譜的統(tǒng)計分析模型,該方法能否準(zhǔn)確模擬爆炸沖擊振動有待進一步研究和驗證。

        3)基于時變功率譜的非平穩(wěn)隨機過程模擬。除相位差譜外,還可以采用時變功率譜來模擬地震動的時頻非平穩(wěn)特性,其數(shù)學(xué)表達式為

        (15)

        式中:Gxx(t,ωk)為t時刻ωk頻率點的時變功率譜值。

        由式(15)可知,建立合理的時變功率譜模型是模擬爆炸地沖擊時程的關(guān)鍵。為此,宋浩等[41-42]先后以Kanai-Tajimi模型和Kameda模型作為目標(biāo)時變功率譜對爆炸地沖擊加速度時程進行了模擬;結(jié)果表明該模擬信號的反應(yīng)譜與實際地沖擊時程的反應(yīng)譜整體趨勢比較一致,但局部存在明顯差異。因此,未來的研究應(yīng)該以大量的爆炸現(xiàn)場的實測信號為基礎(chǔ),深入分析爆炸沖擊振動的時頻特征,提出更符合其客觀物理過程的功率譜或時變功率譜模型。

        4 總結(jié)與展望

        對現(xiàn)有爆炸地沖擊載荷計算方法進行梳理發(fā)現(xiàn),當(dāng)前計算方法適用范圍有限、通用性較差,準(zhǔn)確性有待進一步提高。未來需開展以下研究工作:

        1)充分收集和整理已有的實驗數(shù)據(jù),建立一個廣泛的數(shù)據(jù)庫;并詳細研究爆炸當(dāng)量、爆心距、地質(zhì)條件、埋深、炮孔設(shè)置等因素的影響作用,提高人工智能預(yù)測模型的適用性。

        2)基于爆炸現(xiàn)場的實測數(shù)據(jù),對反應(yīng)譜的圖譜特征進行合理的統(tǒng)計分析,從而對現(xiàn)有的特征頻率法進行修正和改進;此外,還可以嘗試采用人工智能算法建立合適的特征頻率預(yù)測模型。

        3) 研究更高效、更精確的計算方法和更準(zhǔn)確的巖土介質(zhì)動態(tài)本構(gòu)模型,以實現(xiàn)對爆炸全過程的高精度仿真。同時,加強對實驗數(shù)據(jù)的統(tǒng)計和分析,提出符合爆炸地沖擊的客觀物理規(guī)律的時變功率譜模型或者相位差譜模型,進一步提高非平穩(wěn)隨機過程模擬的準(zhǔn)確性。

        猜你喜歡
        特征頻率峰值沖擊
        “四單”聯(lián)動打造適齡兒童隊前教育峰值體驗
        少先隊活動(2022年9期)2022-11-23 06:55:52
        瓷磚檢測機器人的聲音信號處理
        光學(xué)波前參數(shù)的分析評價方法研究
        基于振動信號特征頻率的數(shù)控車床故障辨識方法
        基于小波去噪和EMD算法在齒輪故障檢測中的應(yīng)用
        寬占空比峰值電流型準(zhǔn)PWM/PFM混合控制
        基于峰值反饋的電流型PFM控制方法
        奧迪Q5換擋沖擊
        奧迪A8L換擋沖擊
        一汽奔騰CA7165AT4尊貴型車換擋沖擊
        国产白丝网站精品污在线入口| 中文人妻AV高清一区二区| 成美女黄网站18禁免费| 亚洲国产av一区二区三| 亚洲一区二区视频免费看| 久久久精品国产亚洲av网深田| 99re6在线视频精品免费| 真多人做人爱视频高清免费| 国产精品日韩欧美一区二区区| 女优av福利在线观看| 看一区二区日本视频免费| 熟女肥臀白浆一区二区| 欧美精品videosex极品| 欧美性猛交xxxx乱大交3| 亚洲熟妇色xxxxx欧美老妇| 天堂av在线播放观看| 人妻精品视频一区二区三区| 人人爽久久涩噜噜噜av| 人妻无码aⅴ中文系列久久免费| 久久精品一区二区三区夜夜| 日本一区二区免费高清| 公和我做好爽添厨房中文字幕| 馬与人黃色毛片一部| 人妻少妇久久精品一区二区| 国产亚洲av成人噜噜噜他| 成人爽a毛片免费视频| 久久综合精品国产丝袜长腿| 中年人妻丰满AV无码久久不卡| 国内偷拍第一视频第一视频区| 中文字幕 亚洲精品 第1页| 99久久精品费精品国产一区二区| 99re在线视频播放| 亚州韩国日本区一区二区片| 一级老熟女免费黄色片| 人人妻人人添人人爽欧美一区| 久久亚洲精品成人av| 午夜免费福利一区二区无码AV | 日本精品一级二区三级| 波多野42部无码喷潮在线| 亚洲国产精品线路久久| 麻豆国产AV网站|