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

        ?

        軸向分布式藥包激發(fā)地震波場模型*

        2021-07-30 02:54:00王仲琦
        爆炸與沖擊 2021年7期
        關(guān)鍵詞:振動模型

        徐 謙,王仲琦

        (北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室,北京 100081)

        炸藥震源是油氣勘探中人工激發(fā)地震波的主要震源[1]。炸藥起爆時產(chǎn)生大量高溫、高壓氣體直接作用于巖土介質(zhì)上,在緊鄰炸藥區(qū)域造成劇烈的破壞并形成爆腔,隨著爆炸沖擊波發(fā)展,超壓值迅速降低直至低于巖土介質(zhì)破壞強度,在巖土介質(zhì)中形成地震波[2]。其中,炸藥震源性質(zhì)會影響氣體爆壓大小及作用時間,不同類型的炸藥會形成不同的爆腔,此外,藥包外形還會影響爆腔的幾何形狀。由于爆腔性質(zhì)決定了地震波的特性,因此可以認(rèn)為炸藥震源對地震波特征具有重要的影響作用。

        為了研究炸藥震源形成地震波的過程,Jeffreys[3]最早建立一維空間內(nèi)空腔振動模型,給出球面脈沖下全空間空腔問題的解答。Sharpe[4]對球形封閉空間內(nèi)爆炸產(chǎn)生的子波進(jìn)行研究,獲得了爆炸空腔壁面壓力的彈性波解析解。Blake[5]針對點源爆炸理論模型在Sharpe的研究基礎(chǔ)上獲得了非泊松體的彈性波的解析解。肖建華等[6-7]基于球面波在各向同質(zhì)介質(zhì)的傳播特征,獲得炸藥震源形成球面波振幅和波形;丁樺等[8]通過等效載荷模型開展爆破振動研究,發(fā)現(xiàn)其理論結(jié)果與實際情況相符;林大超等[9]針對炸藥在巖土中的爆炸過程求解出爆腔中的壓力解。這些球形空腔震源模型可以建立爆炸壓力和彈性空腔半徑與地震波場之間的關(guān)系,但是這些方法簡化了炸藥震源激發(fā)地震波的過程,無法建立炸藥震源特征與地震波場幅頻特征之間的直接關(guān)系。Yu 等[10]針對這一問題發(fā)展了空腔震源模型,建立了炸藥震源激發(fā)地震波場理論模型,該方法描述了炸藥震源的形成地震波的作用過程,建立了炸藥震源初始參數(shù)與地震波場之間的關(guān)系,能夠?qū)η蛐嗡幇珗龅卣鸩ㄌ卣鬟M(jìn)行描述。

        雖然球形空腔震源模型已經(jīng)可以描述炸藥爆炸地震波場的形成及變化,但是在地震勘探中會采用多段柱形藥包作為震源,由于球形藥包和柱形藥包形成的地震波場特征并不相同,所以傳統(tǒng)球形空腔震源模型在柱形藥包形成的地震波場描述中并不適用。Heelan[11]針對有限長度的柱形藥包激發(fā)地震波子波進(jìn)行分析,得到柱形藥包的地震波遠(yuǎn)場解。這種方法將炸藥爆炸壓力沿著藥柱長度方向同時施加到激發(fā)介質(zhì)上,是對無限長柱形藥包形成地震波場進(jìn)行直接求解,并不能精確描述有限長柱形藥包的地震波場特征。Starfield 等[12]將長柱形藥包視為多個短柱藥包的疊加,進(jìn)而對柱形藥包端部起爆時的爆炸應(yīng)力場進(jìn)行了求解,計算結(jié)果與實驗測試結(jié)果基本一致。龍源等[13]通過實驗和數(shù)值模擬方法獲得了柱形藥包形成空腔機器發(fā)展過程。李鵬毅等[14]在于成龍等[15]的結(jié)果上通過球形藥包疊加方法,提出了有限長柱形藥包的地震波場模型。胡立新等[16]基于彈性波動理論,通過現(xiàn)場實驗獲得了延遲震源地震波場與藥包類型、藥量以及藥柱間隔之間的關(guān)系。黃文堯等[17]和牟杰[18]對柱形藥包的實驗研究發(fā)現(xiàn)了間隔式柱形藥包起爆的延遲時間對地震波場的影響。這些方法往往基于實驗或者數(shù)值模擬方法對分布式藥包形成的地震波特征進(jìn)行探索,并不能明確描述軸向分布式藥包與地震波場之間的理論關(guān)系。為了能夠?qū)崿F(xiàn)精細(xì)勘探的目的,達(dá)到通過改變炸藥震源激發(fā)方案實現(xiàn)對地震波特征控制的目標(biāo),必須建立軸向分布式藥包激發(fā)地震波場模型。

        本文首先分析炸藥震源激發(fā)地震波的整個過程,概略描述炸藥震源激發(fā)地震波場的不同階段;考慮到藥包的空間結(jié)構(gòu),提出軸向分布式藥包激發(fā)地震波場的計算方法;并通過現(xiàn)場實驗數(shù)據(jù)驗證該計算模型的適用性,并進(jìn)一步分析軸向分布式藥包在巖土中激發(fā)形成地震波場的演變過程。

        1 炸藥震源激發(fā)地震波過程分析

        炸藥從爆炸到最終在遠(yuǎn)處形成彈性波的過程伴隨著一系列的化學(xué)與物理變化過程。由于各種耗散機制的作用,能量在傳遞過程中不斷衰減,由最初具有很強能量的沖擊波演變?yōu)樽罱K的彈性波。該過程分為4個不同階段,分別是:流體動力學(xué)階段、巖土介質(zhì)粉碎階段、動態(tài)膨脹階段和彈性波傳播階段。同時,衰減過程中介質(zhì)在強爆炸作用下會發(fā)生一些不可逆變形,在爆源附近巖土介質(zhì)在巨大的能量作用下表現(xiàn)出流體性質(zhì)。在爆炸瞬間,爆室內(nèi)部推出一個與炸藥形狀相同的波陣面,隨著沖擊波的發(fā)展,沖擊波在向外傳播過程中應(yīng)力峰值急速衰減,直至應(yīng)力值小于巖土極限破壞強度,此時巖土介質(zhì)從流體受力狀態(tài)轉(zhuǎn)化為彈塑性受力形態(tài),應(yīng)力波的這種衰減一直持續(xù)到其峰值小于某個值后,巖土介質(zhì)從塑性狀態(tài)轉(zhuǎn)化為彈性狀態(tài),由此在能量傳遞方向上炸藥震源在巖土中爆炸時依次形成爆炸空腔區(qū)、塑性區(qū)和彈性區(qū)。但是由于柱形藥包與球形藥包在幾何形狀上的區(qū)別,導(dǎo)致各個區(qū)域的發(fā)展過程有所變化,柱形藥包近場各區(qū)域與藥柱形狀相似,隨著爆心距增加,柱形藥包與球形藥包形成的地震波場相似。軸向分布式藥包形成的地震波場逐漸由橢球形向球形發(fā)展,直至最后形成與球形藥包相似的波陣面,各區(qū)域的形狀對比如圖1所示。

        圖1 軸向分布式藥包在巖土介質(zhì)中爆炸分區(qū)Fig.1 Explosive zoning of an axially distributed charge in geotechnical media

        2 軸向分布式炸藥震源激發(fā)地震波場模型

        2.1 球形炸藥激發(fā)地震波場模型

        球形空腔震源模型能夠描述炸藥震源與地震波場之間的關(guān)系。在該模型中,炸藥起爆瞬間產(chǎn)生大量高溫、高壓氣體直接作用于巖土介質(zhì)上形成爆腔。由于空腔中的氣體壓力遠(yuǎn)大于巖土的強度,炸藥附近的巖土介質(zhì)在巨大的壓力下呈現(xiàn)流體性質(zhì),爆腔可以視為在不可壓縮流體介質(zhì)進(jìn)行擴張。隨著爆炸沖擊波在介質(zhì)中的傳播和衰減,在距離震源一定距離處,沖擊波壓力衰減至介質(zhì)破壞強度以下時,此時激發(fā)介質(zhì)呈現(xiàn)彈性性質(zhì)。于成龍等[15]假設(shè)爆轟過程和爆腔瞬間形成,在土介質(zhì)不可壓縮并忽略性質(zhì)改變條件下,建立了球形藥包爆腔預(yù)測準(zhǔn)靜態(tài)模型。通過該模型可以獲得球形藥包在巖土中激發(fā)時的爆腔半徑和塑性區(qū)半徑。爆腔半徑為:

        其中:

        塑性區(qū)半徑為

        式中:a0為球形藥包半徑;p0為炸藥初始爆壓;γ為炸藥膨脹指數(shù);φ為土介質(zhì)黏聚力;α為土介質(zhì)內(nèi)摩擦角;σc為土介質(zhì)抗壓強度;σt為土介質(zhì)抗拉強度;G為拉梅系數(shù),在彈性介質(zhì)中為剪切模量。

        通過球?qū)ΨQ和線性理論,引入線性徑向應(yīng)變和環(huán)向應(yīng)變,質(zhì)點位移u滿足運動方程:

        該方程可以描述質(zhì)點在黏滯阻尼下的受迫振動,根據(jù)方程可以解出振動速度:

        其中:

        式中:ρs為土介質(zhì)初始密度,σ 為應(yīng)力,p為壓力。在上述方程中,振動的頻率為

        振動幅值為:

        2.2 軸向分布式藥包激發(fā)地震波場計算方法

        軸向分布式藥包由等長度的柱形藥包沿著藥柱軸向等間距布置組成。將每段柱形藥包等效為多個球形藥包疊加,則軸向分布式藥包可以視為多級連續(xù)球形藥包間隔疊加。由于各段柱形藥包空間位置不同,不同球形等效藥包與空間某點相對位置有差異,每個藥包引起的指定位置處的振動方向也不同,所以需要對各個球形藥包引起的振動進(jìn)行分解,在x和y方向進(jìn)行分別疊加后,最終獲得合速度。各段藥柱需要滿足以下條件:(1)藥包總體積和所有球形藥包總體積相等;(2)藥包總長度和所有球形藥包直徑之和相等;(3)爆轟過程的連續(xù)性,各段柱形藥包的等效球形藥包間距為零。

        如圖2所示,第m段的第n個球形等效藥包爆炸地震波在A點的振動速度為u˙mn(r,t),x方向速度分量為u˙mnx(r,t),y方向速度分量為u˙mny(r,t),則

        圖2 軸向分布式激發(fā)方案Fig.2 Axial distributed excitation

        式中:Lmnx、Lmny為第m段的第n個球形等效藥包與A點的橫向距離和縱向距離,rmn為第m段的第n個球形等效藥包與A點的實際距離。

        那么第m段藥柱在A處的振動速度在x方向分量u˙mx(r,t)和y方向分量u˙my(r,t)分別為

        整個軸向分布式藥包在A點處的振動速度在x方向分量u˙x(r,t)和y方向分量u˙y(r,t)分別為

        則合速度為

        軸向分布式藥包激發(fā)地震波幅頻特征與各個藥包之間間隔和延時激發(fā)時間有密切關(guān)系,選擇合適的激發(fā)間隔可以獲得更大的能量。通過調(diào)整各段柱形藥包激發(fā)的延遲時間,使第一段藥柱激發(fā)的沖擊波的壓力峰值可以與第2段藥柱激發(fā)沖擊波峰值重疊,就可以達(dá)到最佳疊加效果。由于各段藥柱從激發(fā)到?jīng)_擊波達(dá)到峰值時間相同,所以認(rèn)為最佳激發(fā)間隔等于各段藥柱傳爆時間與兩段藥柱之間的沖擊波傳播時間之和,即:

        式中:l為每段藥柱長度,h為藥柱間距,D為炸藥爆速,c為波速。

        2.3 軸向分布式藥包激發(fā)地震波場模型驗證

        由于特殊的幾何結(jié)構(gòu),軸向分布式藥包在土介質(zhì)中激發(fā)的地震波在軸向和徑向存在差異。在距離藥包的近處,應(yīng)力波波陣面會以各段藥包的基本形狀向外擴散,隨著幾何發(fā)散以及地震波在土中的傳播,地震波特征在軸向和徑向上表現(xiàn)并不相同。為了驗證軸向分布式藥包激發(fā)地震波場計算方法的有效性,本文采用了3種不同的激發(fā)方案對軸向分布式藥包激發(fā)地震波場的理論計算結(jié)果與數(shù)值模擬計算結(jié)果在軸向和徑向分別進(jìn)行對比,詳細(xì)的激發(fā)參數(shù)如表1所示。

        表1 軸向分布式炸藥激發(fā)方案Table 1 Excitation scheme of axially distributed explosives

        將理論模型的計算結(jié)果與有限元軟件Autodyn 計算結(jié)果進(jìn)行對比。

        在Autodyn 模型中,土介質(zhì)采用理想彈塑性模型,其狀態(tài)方程為

        式中:μ=(ρ/ρs)?1,ρ 為土介質(zhì)密度;k為體積模量。

        土介質(zhì)的強度模型采用von Mises強度模型:

        式中:σy為屈服強度,G為剪切模量。土介質(zhì)參數(shù)見表2。

        表2 土介質(zhì)參數(shù)Table 2 Parameters of soil

        TNT炸藥采用理想爆轟產(chǎn)物標(biāo)準(zhǔn)JWL方程描述

        在藥包近區(qū)采用Euler 方法進(jìn)行單元網(wǎng)格劃分,起爆點選取在藥包中心位置,Euler 單元網(wǎng)格介質(zhì)尺寸為2 000 mm×1000 mm,劃分為500×250 個網(wǎng)格。Lagrange單元網(wǎng)格介質(zhì)尺寸為30000 mm×30000 mm,網(wǎng)格劃分為250×250。在藥柱軸向和徑向30 m 范圍內(nèi)分別設(shè)立觀察點,并在兩個方向上分別對比對應(yīng)觀察點處數(shù)值計算與模型計算結(jié)果。

        表3 TNT炸藥特性參數(shù)Table 3 Parameters of TNT

        三種不同激發(fā)方案激發(fā)的地震波最大振動速度的理論計算結(jié)果和數(shù)值模擬結(jié)果對比如圖3所示。在距離震源中心點較近位置處軸向的最大振動速度明顯大于徑向最大振動速度,這是因為軸向上各個藥包激發(fā)的彈性波在疊加過程中不存在x軸方向上的抵消作用,并且軸向上與藥包的距離比徑向與炸藥的距離更近。對比理論計算結(jié)果與數(shù)值模擬結(jié)果在相同爆心距位置處的最大振動速度,在軸向上,爆心距12 m 以內(nèi)差值較大,誤差在3%~7%之間;兩者差值隨著爆心距增加而減小,當(dāng)爆心距大于22 m 時,差值在3%以內(nèi);在徑向上,兩者差值在5%以內(nèi)??梢园l(fā)現(xiàn),在軸向上兩者差距較大,這是因為端部效應(yīng)導(dǎo)致理論計算結(jié)果偏大。隨著爆心距增加,幾何結(jié)構(gòu)的影響逐漸降低,距離較遠(yuǎn)處差距逐漸減小。由于土介質(zhì)吸收衰減作用,軸向上的振動速度會迅速衰減,隨著爆心距增加,相同爆心距處,軸向與徑向上的振動速度會逐漸接近直至相等。

        圖3 軸向分布式藥包最大振動速度Fig.3 Maximum vibration velocity at different locations of axial distributed excitation:

        3 軸向分布式藥包激發(fā)地震波場特征實驗

        為了對本文模型進(jìn)行驗證,開展軸向分布式藥包激發(fā)地震波場實驗,對比實驗與等效模型計算的質(zhì)點振動速度之間的差距。實驗在地面空曠、平整的區(qū)域開展,以減少復(fù)雜地形及構(gòu)筑物對地震波數(shù)據(jù)采集工作的影響,增加實驗數(shù)據(jù)的準(zhǔn)確性。軸向分布式藥包激發(fā)方案如表1所示。根據(jù)現(xiàn)場打孔取樣的勘察結(jié)果,距離地面1~5 m 為粉土,5~10 m 為粉質(zhì)黏土,粉質(zhì)黏土的詳細(xì)參數(shù)如表2所示。本次測試振動速度的傳感器指標(biāo)如表4所示。

        表4 檢波器部分參數(shù)Table 4 Parameters of sensors

        傳感器布設(shè)在地面,距離炮孔最近傳感器的水平距離為10 m,其后間隔4 m 呈直線依次布設(shè)6個傳感器,共布置12個傳感器,埋深為0.3 m,檢波器布置如圖4所示。

        圖4 傳感器布置Fig.4 Sensor layout

        對不同激發(fā)方案的質(zhì)點振動峰值進(jìn)行統(tǒng)計,并與軸向分布式藥包激發(fā)地震波場模型的理論計算結(jié)果進(jìn)行對比,得到不同激發(fā)方案激發(fā)時質(zhì)點振動速度峰值衰減曲線,如圖5所示。其中實線為軸向分布式藥包模型計算的理論振動速度最大值,單點為對應(yīng)位置處振動速度檢波器實測質(zhì)點振動速度的最大值。表5列出不同位置處最大振動速度誤差對比結(jié)果。

        圖5 不同激發(fā)方案地震波的峰值粒子振動速度衰減規(guī)律Fig.5 Peak particle vibration velocities of seismic waveswith different excitation schemes

        表5 不同位置(x)處峰值粒子速度實驗結(jié)果相與計算結(jié)果的相對誤差Table 5 Table6 Relativeerror of experimental resultsof peak paritcle velocity to the calculational onesat different position x

        通過對比發(fā)現(xiàn)實驗結(jié)果與等效模型計算結(jié)果在3 kg×2和2 kg×3兩種激發(fā)方案中誤差較小,1.5 kg×4的誤差較大。這主要是地層分層結(jié)構(gòu)造成的,實驗位置處5~10 m 為物理參數(shù)相同的粉質(zhì)黏土,當(dāng)采用1.5 kg×4的激發(fā)方案時,第一段和最后一段炸藥起爆的爆腔及塑性區(qū)有一部分處于其他地層中,不同地層的巖土參數(shù)并不相同,導(dǎo)致炸藥形成的爆腔及地震波特性也不同,所以計算結(jié)果與實際結(jié)果會有差異。但是完全處于粉質(zhì)黏土層的3 kg×2和2 kg×3的兩種激發(fā)方案中的振動速度與理論值相差在5%以內(nèi),這說明軸向分布式藥包激發(fā)地震波模型可以描述工程應(yīng)用中炸藥地震波的特性,符合工程應(yīng)用需要。

        4 討 論

        軸向分布式藥包的幾何結(jié)構(gòu)決定了震源近場區(qū)域的地震波場特征在軸向與徑向并不相同,但是隨著爆心距的增加,其軸向與徑向上的最大振動速度呈趨近狀態(tài)。為了研究軸向分布式藥包激發(fā)地震波場的發(fā)展?fàn)顩r,以相同藥量的球形藥包激發(fā)地震波場作為參考,求解與球形藥包最大質(zhì)點振動速度相等的分布式藥包爆心距的軸向距離及徑向距離,求解示意圖見圖6。以3 kg×2 TNT在土中的爆炸過程為例,對在50 m×50 m 范圍內(nèi)分布式炸藥地震波軸向和徑向發(fā)展過程進(jìn)行分析,結(jié)果如圖7所示,其中每段柱形藥包長徑比為15。

        圖6 球形藥包和軸向分布式藥包地震波場Fig.6 Seismic wave field of spherical explosive and axially-distributed explosive

        從圖7可以看出,在10 m 范圍內(nèi)地震波振幅值迅速衰減,其中軸向分布式藥包在軸向分方向上的振動速度大于徑向的振動速度,主要原因是徑向上不同位置藥包形成的地震波會有一部分能量被抵消,而軸向上的各個藥包形成的地震波可以無損疊加。當(dāng)?shù)刃Ь嚯x為18 m 時,軸向分布式藥包激發(fā)的地震波振幅等值線在軸向距離與徑向距離上的差異為5%,這時的等效距離為藥柱長度的9.8倍。另外通過對比發(fā)現(xiàn),當(dāng)藥量相同時隨著距離增加,球形藥包與分布式柱形藥包在相同距離處質(zhì)點最大振動相同,可以認(rèn)為在距離爆源中心一定距離處分布式柱形藥包與球形藥包形成的地震波最大振動速度是相同的。

        圖7 相同質(zhì)點振動速度時球形藥包距離與柱形藥包軸向、徑向距離對比Fig.7 Comparison of the distance between spherical explosives and different directions of axially distributed cylindrical explosives at the same particle vibration velocity

        同樣的,3 kg×2軸向分布式藥包與6 kg 球形藥包在相同巖土介質(zhì)條件下,頻率衰減過程如圖8所示。

        圖8 軸向分布式柱形藥包和球形藥包形成地震波主頻對比Fig.8 Comparison of frequency of spherical charge and axially distributed charge

        可以發(fā)現(xiàn)軸向分布式柱形藥包產(chǎn)生地震波頻率大于球形裝藥產(chǎn)生的頻率,因為軸向分布式藥包主頻與各個等效球形藥包直徑相關(guān),這些等效藥包產(chǎn)生地震波主頻大于球形藥包主頻。軸向分布式炸藥藥包激發(fā)地震波在軸向上頻率大于在徑向上頻率,這是由于炸藥在軸向傳爆需要一定時間,地震波傳播到相同爆心距時,徑向已經(jīng)產(chǎn)生一部分衰減,而軸向則尚未衰減,所以兩者會產(chǎn)生差距,這種頻率差距會隨著距離增加而減小,當(dāng)爆心距到達(dá)15 m 時兩者差距為3%,達(dá)到一定距離后兩個方向上的地震波主頻相等。

        5 結(jié) 論

        為了明確柱軸向分布式藥包與地震波場之間的關(guān)系,本文以球形空腔震源模型為基礎(chǔ),通過多個球形藥包疊加等效替代柱形藥包的方法,提出軸向分布式藥包激發(fā)地震波場的計算方法,建立了軸向分布式藥包激發(fā)地震波場模型。該模型克服了傳統(tǒng)球形藥包震源模型和柱形藥包震源模型無法精確描述對地震勘探中分布式藥包地震波場特性的不足,為指導(dǎo)精細(xì)勘探生產(chǎn)需要提供理論基礎(chǔ)。

        采用數(shù)值模擬和現(xiàn)場實驗兩種方法驗證了理論模型的適用性。其中,數(shù)值模擬計算結(jié)果表明,徑向上,理論模型與數(shù)值模型計算結(jié)果誤差在5%以內(nèi);軸向上,理論模型與數(shù)值模型計算結(jié)果誤差3.4%以內(nèi)?,F(xiàn)場實驗結(jié)果表明,爆心距大于14 m 時理論模型地震波振動速度誤差都在10%以內(nèi),計算誤差隨著距離增加而減小,距離大于24 m 時誤差小于6%。對比實驗結(jié)果發(fā)現(xiàn),軸向分布式藥包激發(fā)地震波場模型可以對軸向分布式藥包在巖土中激發(fā)的地震波場進(jìn)行描述。

        軸向分布式藥包在巖土中激發(fā)時,在炸藥附近,其形成的地震波場存在方向性,即軸向與徑向上的波場存在差異。在爆心距相同條件下,軸向上振動速度大于徑向上振動速度,兩者差異隨著爆心距的增加而減小,當(dāng)爆心距為藥柱長度的16~21倍時,軸向與徑向振動速度差異在5%以內(nèi),此時軸向分布式藥包形成的地震波場與球形藥包形成地震波場在能量上基本一致,但是軸向分布式藥包形成的地震波場主頻高于球形藥包地震波場主頻。

        猜你喜歡
        振動模型
        一半模型
        振動的思考
        噴水推進(jìn)高速艇尾部振動響應(yīng)分析
        重要模型『一線三等角』
        This “Singing Highway”plays music
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        振動攪拌 震動創(chuàng)新
        中國公路(2017年18期)2018-01-23 03:00:38
        中立型Emden-Fowler微分方程的振動性
        3D打印中的模型分割與打包
        FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
        亚洲av毛片在线免费观看 | 91精品91久久久久久| 亚洲女同性恋在线播放专区| 在线a亚洲视频播放在线播放| 三年中文在线观看免费大全| 精品一区二区三区四区国产| 天天做天天爱天天综合网2021| 国产激情对白一区二区三区四| 亚洲国产一区二区三区,| 亚洲天堂av福利在线| 免费a级毛片无码av| 亚洲乱码av中文一区二区| 国产国拍亚洲精品午夜不卡17| 亚洲av综合色区久久精品| 国产精品第一区亚洲精品| 欧美高清视频手机在在线| 天堂…在线最新版资源| 久久精品性无码一区二区爱爱| 久久开心婷婷综合中文| 风韵丰满熟妇啪啪区老熟熟女 | 精品人无码一区二区三区| 91爱爱视频| 亚洲国产高清一区av| 亚洲色精品三区二区一区| 日本无遮挡吸乳呻吟视频| 国产亚洲女人久久久久久| 亚洲精品国产一区二区免费视频 | 亚洲一区二区三区香蕉| 亚洲综合性色一区| 久久天堂av综合合色| 亚洲av中文无码字幕色本草| 成人亚洲性情网站www在线观看| 香蕉久久夜色精品国产| 国产性感丝袜在线观看| 国产激情久久久久影院老熟女| 国产精品久久无码不卡黑寡妇| 中文字幕专区一区二区| 国产欧美性成人精品午夜| 精品视频一区二区三三区四区| 国产麻豆放荡av激情演绎| 狂猛欧美激情性xxxx大豆行情|