• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    有限元-點插值耦合法大地電磁二維正演模擬

    2015-06-27 05:54:49李俊杰嚴家斌
    石油物探 2015年4期
    關(guān)鍵詞:網(wǎng)格法電阻率電磁

    李俊杰,嚴家斌

    (1.浙江省水利水電勘測設(shè)計院,浙江杭州310002;2.中南大學地球科學與信息物理學院,有色資源與地質(zhì)災害探查湖南省重點實驗室,湖南長沙410083)

    有限元-點插值耦合法大地電磁二維正演模擬

    李俊杰1,嚴家斌2

    (1.浙江省水利水電勘測設(shè)計院,浙江杭州310002;2.中南大學地球科學與信息物理學院,有色資源與地質(zhì)災害探查湖南省重點實驗室,湖南長沙410083)

    點插值法(PIM)作為一種典型的全域弱式無網(wǎng)格法,該方法在地質(zhì)建模時將物性加載到只與坐標有關(guān)的高斯積分點上,因此處理復雜模型時較常規(guī)網(wǎng)格方法便利,但缺點是計算效率低。將有限元法(FEM)與PIM耦合,形成FE-PIM,用于大地電磁二維正演模擬。利用Galerkin法代入插值法構(gòu)造的形函數(shù)并結(jié)合高斯積分公式推導了大地電磁二維無網(wǎng)格化總體矩陣表達式,簡述了背景網(wǎng)格積分與邊界條件的加載技術(shù),理論模型的數(shù)值計算驗證了FE-PIM算法的正確性、高效性及其在處理復雜模型上的便利性。

    點插值法;全域弱式無網(wǎng)格法;大地電磁;有限元-點插值法

    大地電磁二維正演多采用有限差分法(finite deference method,FDM)[1-2]、有限元法(finite element method,FEM)[3-6]等網(wǎng)格方法。對于地質(zhì)體相對規(guī)則、地形簡單的模型,FDM利用較小的網(wǎng)格就能獲得較高的計算精度,效率高;而對于復雜的地質(zhì)體,就需要精細網(wǎng)格才能完成高精度的模擬。FEM適用于復雜邊界與復雜電性結(jié)構(gòu)模型的計算,但高維問題中非結(jié)構(gòu)化三角單元生成程序?qū)崿F(xiàn)困難。無單元Galerkin法(element-free Galerkin method,EFGM)[7]是一種基于節(jié)點的全域弱式無網(wǎng)格方法,其形函數(shù)構(gòu)造不依賴預定義的網(wǎng)格單元,選擇離散的節(jié)點便可求Helmholtz方程的解。相對于低階(low-order)有限元分析,EFGM只需少量的節(jié)點就可獲得光滑的解,且場分量的導數(shù)同樣光滑連續(xù),對于大地電磁輔助分量的求解非常有利,可以獲得與主分量(水平電場或磁場)相同的計算精度,對于節(jié)點的參數(shù)化處理也較單元分析容易。EFGM能處理網(wǎng)格方法難以處理的問題[8-10],有廣闊的應(yīng)用前景。李俊杰等[11]綜述了無網(wǎng)格法進展并推導了大地電磁二維問題的無網(wǎng)格化系統(tǒng)矩陣表達式;賈曉峰等[12-13]將EFGM用于疊前地震模擬,并討論了幾種吸收邊界方法在地震模擬中與EFGM的結(jié)合,研究結(jié)果表明,EFGM精度較FDM精度高,且具有較好的穩(wěn)定性;王月英[14]分析了基函數(shù)階次對EFGM用于地震波正演計算精度的影響;馮德山等[15]、戴前偉等[16]將EFGM用于雷達波場的二維正演模擬,詳細推導了探地雷達無單元法波動方程,并提出了適用于EFGM雷達正演模擬的透射吸收邊界、Sarma吸收邊界的改進方法。嚴家斌等[17]將EFGM成功應(yīng)用于大地電磁二維正演模擬,并分析了支持域無量綱尺寸與高斯點數(shù)量對無網(wǎng)格法正演計算精度與效率的影響。無網(wǎng)格點插值法(point interpolation method,PIM)[18]是另一種簡單高效的全域弱式無網(wǎng)格法,與EFGM的區(qū)別主要在于形函數(shù)采用過點插值的方法進行構(gòu)造,邊界條件加載便利。李俊杰等[19]將PIM用于復雜模型的MT二維正演模擬,凸顯了PIM較常規(guī)網(wǎng)格方法在地質(zhì)建模上的優(yōu)越性。

    我們將PIM與高效的FEM相耦合,形成有限元-點插值法(finite element-point interpolation method,FE-PIM),從大地電磁場的二維正演理論出發(fā),采用Galerkin法推導了FE-PIM系統(tǒng)矩陣表達式,數(shù)值計算結(jié)果驗證了FE-PIM的有效性及其在計算復雜模型時高效、便利的特性。

    1 大地電磁二維變分問題

    研究地下二維電性結(jié)構(gòu),取y軸垂直向上,z軸為走向,x軸向右且與y軸、z軸垂直,求解域Ω為矩形,4個頂點依次以A,B,C,D順時針編號,Γ1為地質(zhì)體邊界。大地電磁二維正演滿足(1)式所示的變分問題[20]:

    (1a)

    uAB=1

    (1b)

    δF(u)=0

    (1c)

    TE模式為:

    u=Ez

    (2)

    TM模式為:

    u=Hz

    (3)

    式中:ω為角頻率;μ為磁導率;σ為電導率;ε為介電常數(shù);Ez,Hz分別為電場及磁場平行于異常體走向的水平分量。

    2 大地電磁二維變分問題的無網(wǎng)格解法

    PIM以積分點支持域內(nèi)的場節(jié)點利用多項式基插值構(gòu)造形函數(shù),背景網(wǎng)格用于計算總體矩陣表達式中包含的求解域及求解域邊界的積分項。圖1 為PIM支持域、高斯點、背景網(wǎng)格與場節(jié)點示意圖。由于網(wǎng)格積分多采用高斯積分法,故積分點也稱高斯點。

    圖1 全域弱式無網(wǎng)格法支持域、高斯點、背景網(wǎng)格與場節(jié)點圖示

    2.1 支持域

    如圖1所示,支持域形狀常為矩形或圓形,對于任一高斯點XQ,支持域尺寸d可表示為:

    (4)

    式中:α為支持域無量綱尺寸,其大小影響無網(wǎng)格法的計算精度[17];dc為支持域內(nèi)高斯點XQ的平均結(jié)點間距,有:

    (5)

    式中:A為支持域面積;n為包含在A中的節(jié)點總數(shù)。本文采用矩形支持域,x,y方向的支持域尺寸為:

    (6)

    式中:dcx與dcy分別為x,y方向的節(jié)點間距;αx與αy為對應(yīng)的無量綱尺寸。為便于程序設(shè)計,一般令αx=αy=α,本文取α=1。

    2.2 FE-PIM離散系統(tǒng)方程的構(gòu)造

    PIM計算效率低但處理復雜模型便利[19],FEM求解高效,由于PIM與FEM系統(tǒng)矩陣采用相同的方法構(gòu)成,故PIM與FEM具有很好的耦合條件。

    如圖2所示,FE-PIM的基本原理是將求解域Ω劃分為實線表示的FEM區(qū)域Ω1及虛線表示的PIM背景單元區(qū)域Ω2,異常體置于Ω2,無窮遠邊界單元設(shè)為Ω1。將變分符號δ代入(1)式并將FEM區(qū)域單元離散化,有:

    (7)

    式中:e表示FEM的子單元;CD為求解域底邊界;Γe為邊界單元。

    圖2 FE-PIM計算模型二有限單元與背景網(wǎng)格分布圖示

    將FEM區(qū)域的積分項即δF1(u)離散,有:

    (8)

    (9)

    式中:U為全部節(jié)點的場值向量;KFEM為FEM區(qū)域計算的總體系統(tǒng)矩陣。求得δF1(u)的離散表達式后,將δF2(u)離散,將u表示為PIM形函數(shù)向量與場節(jié)點向量之積的形式,有

    (10)

    式中:Φ為用PIM構(gòu)造的形函數(shù)矩陣[18-19];n為支持域內(nèi)的節(jié)點數(shù);u為支持域內(nèi)n個節(jié)點的場向量值。對(10)式進行變分運算,有:

    (11)

    將(10)式與(11)式代入(1)式的第1項并引入總體編號體系,總體矩陣最終變?yōu)?

    (12)

    式中:M1與Mn表示PIM區(qū)域第一個節(jié)點與最后一個節(jié)點的編號;KPIM為背景單元域PIM計算的總體系統(tǒng)矩陣。KPIM表達式中包含對求解域Ω2的積分,積分利用高斯積分法求解,有:

    (13)

    (14)

    由于δUT的任意性,(14)式成立的條件變?yōu)?

    (15)

    (15)式即為用FE-PIM構(gòu)造的系統(tǒng)矩陣,由于FEM與PIM構(gòu)造的總體矩陣KFEM及KPIM均具有帶狀、稀疏、對稱的性質(zhì),故K也具備該特性。

    線性方程組KU=0的求解需加載(1)式所示的邊界條件,PIM與FEM邊界條件均可直接加載,計算時先找出邊界節(jié)點在總體剛度矩陣中對應(yīng)的KII,KII對應(yīng)的行與列設(shè)置為0,其值設(shè)置為1,最后將方程右端項對應(yīng)的值PI設(shè)置為邊界上的預設(shè)場值即可。PIM邊界條件可精確加載,該特性是PIM較EFGM最大的優(yōu)勢。

    3 正演計算

    3.1 算法驗證

    PIM與FEM在一維介質(zhì)大地電磁(MT)正演中具有較高的精度[19],FE-PIM為兩者的耦合,先通過圖3所示的COMMEMI 2D-0模型[21]來驗證文中算法的正確性。此模型異常體頂部與地表重合,縱向規(guī)模50km,橫向規(guī)模20km,電阻率為1Ω·m,其左、右兩側(cè)電阻率分別為10Ω·m和2Ω·m,深度大于50km區(qū)域為超導層,電阻率為0,求解域TE模式取120km×110 km,TM模式取120km×70km,場節(jié)點橫縱向等間距分布,節(jié)點間距為1km。

    圖3 COMMEMI 2D-0模型

    表1列出了PIM與FEM計算COMMEMI 2D-0模型的視電阻率結(jié)果。從表1可以看出,兩種方法計算結(jié)果與COMMEMI解析解基本吻合,驗證了數(shù)值算法的有效性。此外,表1中FEM計算結(jié)果也與文獻[22]一致。

    表1 COMMEMI 2D-0模型數(shù)值方法視電阻率計算結(jié)果

    3.2 有效性及優(yōu)勢

    為進一步驗證FE-PIM有效性及其優(yōu)勢,設(shè)計了如圖4所示的二維地質(zhì)模型:模型二異常體為截面方形二度體,異常體邊長400m,埋深800m;模型三為水平橢圓,長半軸300m,短半軸200m,埋深800m;模型四為等腰直角三角形地壘,斜邊長800m;模型五為出露于地表的條帶狀低阻體,異常體截面呈平行四邊形,下底長400m,上、下底間距800m,異常體右邊界與地面的夾角呈45°。

    模型背景電阻率1000Ω·m,異常體電阻率100Ω·m。場節(jié)點等間距分布于求解域,TE模式3321(41×81)個場節(jié)點,TM模式1681(41×41)個場節(jié)點,橫、縱向節(jié)點間距均為200m,FEM與FE-PIM采用相同的節(jié)點分布。

    圖5顯示了頻率為100Hz時PIM與FEM模型二視電阻率計算結(jié)果。從圖5可以看出,求解域兩側(cè)邊界的視電阻率值約為1000Ω·m,低阻異常在里程0附近最顯著,TE模式視電阻率值約780Ω·m,TM模式約920Ω·m。兩種模式下FE-PIM,PIM與FEM的計算結(jié)果基本一致,驗證了FE-PIM二維算法的正確性。

    圖4 二維地質(zhì)模型

    圖5 頻率為100Hz時模型二視電阻率計算結(jié)果

    圖6為不同模型的FE-PIM視電阻率計算結(jié)果。從圖6a和圖6b可以看出,TE模式低阻異常中心與模型中心一致,能明顯地反映出橢圓低阻體的存在,但異常規(guī)模比實際低阻體大,視電阻率為740~1080Ω·m。TM模式低阻異常呈直立狀并無限向下延伸,異常的水平寬度反映了實際低阻體水平方向的尺寸,視電阻率為800~1060Ω·m。從圖6c可以看出,模型四山脊位置在視電阻率斷面圖上顯示為高阻極大值區(qū)域,兩側(cè)的山谷區(qū)則呈現(xiàn)低阻極小值,且頻率越高此特性越明顯,視電阻率變化范圍為900~1950Ω·m,此現(xiàn)象表明地形對大地電磁測深法影響規(guī)律與直流電阻率法存在差異。斷面圖高阻區(qū)的橫向規(guī)模約800m,與山脊橫向規(guī)模相當,驗證了FE-PIM的有效性。從圖6d 可以看出,模型五的視電阻率斷面圖左右非對稱,低阻區(qū)域呈“上窄下寬”條帶狀傾斜分布,傾向與地質(zhì)體相同,視電阻率為100~1100Ω·m,較好地反映出了異常體的物性及空間賦存狀態(tài)。

    圖6 不同模型的視電阻率FE-PIM計算結(jié)果

    FE-PIM延續(xù)PIM地質(zhì)建模便利特點的同時也兼顧了計算效率,表2列出了17個頻點PIM,FE-PIM與FEM的計算耗時。從表2可以看出,PIM耗時約為FEM耗時的8~9倍,FE-PIM耗時相對PIM耗時呈著降低,無網(wǎng)格區(qū)域(10×4)的FE-PIM與PIM相比,TE模式下效率約提高87%,TM模式下效率約提高83%,隨著FEM區(qū)域的增大,FE-PIM耗時逐漸接近FEM耗時。

    表2 模型二不同數(shù)值方法計算耗時

    4 結(jié)束語

    1) 將FE-PIM應(yīng)用于大地電磁二維正演模擬,通過把求解區(qū)域劃分為有限元區(qū)域及背景單元區(qū)域,利用Galerkin法和高斯積分公式導出了大地電磁二維FE-PIM耦合的系統(tǒng)方程離散表達式,介紹了邊界條件直接加載的技巧,COMMEMI 2D-0 模型的數(shù)值計算驗證了文中PIM與FEM算法的正確性。

    2) 采用節(jié)點規(guī)則分布的FE-PIM計算了截面為方形、橢圓、平行四邊形及三角地壘二維模型的大地電磁場響應(yīng),視電阻率斷面圖較好地反映了地質(zhì)體產(chǎn)生的電磁異常響應(yīng)及空間賦存狀態(tài),驗證了耦合算法的正確性,凸顯了FE-PIM處理復雜地質(zhì)模型的優(yōu)越性。

    3) 比較分析了PIM,FE-PIM及FEM的計算效率,PIM計算耗時明顯高于FEM計算耗時;FE-PIM計算耗時受異常體區(qū)域的影響較大,當異常體區(qū)域較小時,FE-PIM計算耗時遠小于PIM計算耗時而略高于FEM計算耗時,但FE-PIM綜合了FEM計算高效及PIM處理復雜模型便利的優(yōu)點,是一種優(yōu)越的數(shù)值方法。

    [1] 胡祥云,霍光譜,高銳,等.大地電磁各向異性二維模擬及實例分析[J].地球物理學報,2013,56(12):4268-4277 Hu X Y,Huo G P,Gao R,et al.The magnetotelluric anisotropic two-dimensional simulation and case analysis[J].Chinese Journal of Geophysics,2013,56(12):4268-4277

    [2] 董浩,魏文博,葉高峰,等.基于有限差分正演的帶地形三維大地電磁反演方法[J].地球物理學報,2014,57(3):939-952 Dong H,Wei W B,Ye G F,et al.Study of three-dimensional magnetotelluric inversion including surface topography based on finite-difference method[J].Chinese Journal of Geophysics,2014,57(3):939-952

    [3] Key K,Weiss C.Adaptive finite-element modeling using unstructured grids:the 2D magnetotelluric example[J].Geophysics,2006,71(6):291-299

    [4] 劉長生,湯井田,任政勇,等.基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬[J].中南大學學報(自然科學版),2010,41(5):1855-1859 Liu C S,Tang J T,Ren Z Y,et al.Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes[J].Journal of Central South University(Science and Technology),2010,41(5):1855-1859

    [5] Ren Z Y,Tang J T.3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method[J].Geophysics,2010,75(1):7-17

    [6] 柳建新,孫麗影,童孝忠,等.基于自適應(yīng)有限元的起伏地形MT二維正演模擬[J].地球物理學進展,2012,27(5):2016-2023 Liu J X,Sun L Y,Tong X Z,et al.Undulating terrain 2D MT forward modelling with adaptive finite element[J].Progress in Geophysics,2012,27(5):2016-2023

    [7] Belytschko T,Lu Y Y,Gu L.Element-free Galerkin methods[J].International Journal for Numerical Methods in Engineering,1994,37(2):229-256

    [8] Belytschko T,Lu Y Y,Gu L.Fracture and crack growth by element-free Galerkin methods[J].Modelling and Simulation in Materials Science and Engineering,1994,2(3):519-534

    [9] Li D M,Bai F N,Cheng Y M,et al.A novel complex variable element-free Galerkin method for two-dimensional large deformation problems[J].Computer Methods in Applied Mechanics and Engineering,2012,233-236(1):1-10

    [10] Liu L C,Dong X H,Li C X.Adaptive finite element-element-free Galerkin coupling method for bulk metal forming processes[J].Journal of Zhejiang University-Science A,2009,10(3):353-360

    [11] 李俊杰,嚴家斌.無網(wǎng)格法進展及其在地球物理學中的應(yīng)用[J].地球物理學進展,2014,29(1):452-461 Li J J,Yan J B.Developments of meshless method and applications in geophysics[J].Progress in Geophysics,2014,29(1):452-461

    [12] 賈曉峰,胡天躍,王潤秋.復雜介質(zhì)中地震波模擬的無單元法[J].石油地球物理勘探,2006,41(1):43-48 Jia X F,Hu T Y,Wang R Q.A mesh-free algorithm of seismic wave simulation in complex medium[J].Oil Geophysical Prospecting,2006,41(1):43-48

    [13] 賈曉峰,胡天躍,王潤秋.無單元法用于地震波波動方程模擬與成像[J].地球物理學進展,2006,21(1):11-17 Jia X F,Hu T Y,Wang R Q.Wave equation modeling and imaging by using element-free method[J].Progress in Geophysics,2006,21(1):11-17

    [14] 王月英.地震波正演模擬中無單元Galerkin法初探[J].地球物理學進展,2007,22(5):1539-1544 Wang Y Y.Study of element-free Galerkin method in the seismic forward modeling[J].Progress in Geophysics,2007,22(5):1539-1544

    [15] 馮德山,王洪華,戴前偉.基于無單元Galerkin法探地雷達正演模擬[J].地球物理學報,2013,56(1):298-308 Feng D S,Wang H H,Dai Q W.Forward simulation of ground penetrating radar based on the element-free Galerkin method[J].Chinese Journal of Geophysics,2013,56(1):298-308

    [16] 戴前偉,王洪華.基于隨機介質(zhì)模型的GPR無單元法正演模擬[J].中國有色金屬學報,2013,23(9):1-9 Dai Q W,Wang H H.Element free method forward modeling of GPR based on random medium model[J].The Chinese Journal of Nonferrous Metals,2013,23(9):1-9

    [17] 嚴家斌,李俊杰.無網(wǎng)格法在大地電磁正演計算中的應(yīng)用[J].中南大學學報(自然科學版),2014,45(10):3513-3520 Yan J B,Li J J.Magnetotelluric forward calculation by meshless method[J].Journal of Central South University(Science and Technology),2014,45(10):3513-3520

    [18] Liu G R,Gu Y T.A point interpolation method for two-dimensional solids[J].International Journal for Numerical Methods in Engineering,2001,50(4):937-951

    [19] 李俊杰,嚴家斌.無網(wǎng)格點插值法大地電磁二維正演數(shù)值模擬[J].石油物探,2014,53(5):617-626 Li J J,Yan J B.Magnetotelluric two-dimensional for-

    ward numerical modeling by meshfree point interpolation method[J].Geophysical Prospecting for Petroleum,2014,53(5):617-626

    [20] 戴前偉,張富強,楊坤平,等.電導率分塊線性變化的二維高頻大地電磁法有限元正演模擬[J].物探化探計算技術(shù),2012,34(5):552-558 Dai Q W,Zhang F Q,Yang K P,et al.The finite element method for modeling 2-D high-frequency electromagnetic with continuous variation of conductivity within each block[J].Computing Techniques for Geophysical and Geochemical Exploration,2012,34(5):552-558

    [21] Zhdanov M S,Varentsov I M,Weaver J T,et al.Methods for modelling electromagnetic fields results from COMMEMI-the international project on the comparison of modelling methods for electromagnetic induction[J].Journal of Applied Geophysics,1997,37(3):133-271

    [22] 許建榮.起伏地形條件下大地電磁測深二維正反演研究及應(yīng)用[D].長沙:中南大學地球科學與信息物理學院,2010 Xu J R.Research and applications of 2-D MT forward modeling and inversion with topography[D].Changsha:Central South University of Geosciences and Info-Physics,2010

    (編輯:顧石慶)

    Magnetotelluric two-dimensional forward modelling by finite element-point interpolation coupling method

    Li Junjie1,Yan Jiabin2

    (1.ZhejiangDesignInstituteofWaterConservancyandHydroelectricPower,Hangzhou310002,China; 2.KeyLaboratoryofNonferrousResourcesandGeologicalHazardDetection,SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha410083,China)

    Point interpolation method (PIM) is a typical global weak-form meshfree numerical calculation method.The physical properties of PIM are loaded on Gauss integral points which are only related with coordinates during geological modeling.Therefore,PIM is more convenient while dealing complex model than conventional grid method,but the former computation efficiency is low.We couple FEM and PIM into finite element-point interpolation method (FE-PIM) for magnetotelluric two-dimensional forward.Firstly,magnetotelluric two-dimensional discrete system matrix is deduced through substituting the shape function constructed by interpolation method and combining Galerkin method with Gauss integral formula.Then,the principle of background grid integral and the loading technique of boundary conditions are briefly characterized.Finally,the effectiveness and high efficiency of the FE-PIM algorithm and the convenience for complex models are proved by several numerical models.

    point interpolation method,global weak-form meshfree method,magnetotelluric,finite element-point interpolation method

    2014-12-02;改回日期:2015-02-15。

    李俊杰(1989—),碩士,助理工程師,現(xiàn)主要從事大地電磁無網(wǎng)格化正演模擬研究工作。

    國家自然科學基金項目(40874055)和湖南省自然科學基金項目(14JJ2012)共同資助。

    P631

    A

    1000-1441(2015)04-0477-08

    10.3969/j.issn.1000-1441.2015.04.015

    猜你喜歡
    網(wǎng)格法電阻率電磁
    雷擊條件下接地系統(tǒng)的分布參數(shù)
    科技風(2020年13期)2020-05-03 13:44:08
    三維多孔電磁復合支架構(gòu)建與理化表征
    角接觸球軸承的優(yōu)化設(shè)計算法
    科學與財富(2019年3期)2019-02-28 07:33:42
    基于遺傳算法的機器人路徑規(guī)劃研究
    基于GIS的植物葉片信息測量研究
    掌握基礎(chǔ)知識 不懼電磁偏轉(zhuǎn)
    三維電阻率成像與高聚物注漿在水閘加固中的應(yīng)用
    隨鉆電阻率測井的固定探測深度合成方法
    海洋可控源電磁場視電阻率計算方法
    粉煤灰摻量對水泥漿體電阻率與自收縮的影響
    91精品三级在线观看| 久久青草综合色| 国产欧美日韩一区二区精品| 国产免费福利视频在线观看| 大陆偷拍与自拍| 国产精品九九99| 久久精品国产亚洲av香蕉五月 | 美女扒开内裤让男人捅视频| bbb黄色大片| 在线观看免费日韩欧美大片| 久久99一区二区三区| 国产免费现黄频在线看| 久久久水蜜桃国产精品网| 国产一卡二卡三卡精品| 成年动漫av网址| 免费少妇av软件| 亚洲精品久久成人aⅴ小说| 波多野结衣av一区二区av| 成年av动漫网址| 久久久久久免费高清国产稀缺| 久久国产精品影院| 99热全是精品| 国产精品国产三级国产专区5o| 黄网站色视频无遮挡免费观看| 黄网站色视频无遮挡免费观看| 水蜜桃什么品种好| 久热爱精品视频在线9| 久久久久网色| 水蜜桃什么品种好| 精品少妇一区二区三区视频日本电影| 国产精品99久久99久久久不卡| 又黄又粗又硬又大视频| 一本久久精品| 精品卡一卡二卡四卡免费| 91精品国产国语对白视频| 久久人人爽av亚洲精品天堂| 亚洲av电影在线观看一区二区三区| 国产高清videossex| 青春草亚洲视频在线观看| 成人影院久久| 50天的宝宝边吃奶边哭怎么回事| 黑人欧美特级aaaaaa片| 9色porny在线观看| 久久久久久久久免费视频了| av天堂久久9| 亚洲欧美日韩高清在线视频 | 国产欧美日韩一区二区精品| 国产黄频视频在线观看| 十八禁人妻一区二区| 汤姆久久久久久久影院中文字幕| 最新的欧美精品一区二区| 18在线观看网站| 中文字幕人妻熟女乱码| 国产成人av教育| 高清av免费在线| 高清在线国产一区| 亚洲精品久久午夜乱码| 王馨瑶露胸无遮挡在线观看| 午夜91福利影院| 日韩免费高清中文字幕av| 欧美日韩av久久| 国产精品免费视频内射| 动漫黄色视频在线观看| 国产欧美日韩一区二区三区在线| 亚洲精品一区蜜桃| 日本欧美视频一区| 欧美变态另类bdsm刘玥| 99国产精品免费福利视频| 大香蕉久久网| videosex国产| 俄罗斯特黄特色一大片| 亚洲一码二码三码区别大吗| 一级,二级,三级黄色视频| 男人操女人黄网站| 久久久国产成人免费| 性少妇av在线| 久久久久视频综合| 高清在线国产一区| 精品国产乱子伦一区二区三区 | 亚洲视频免费观看视频| 在线永久观看黄色视频| 80岁老熟妇乱子伦牲交| 亚洲精品美女久久久久99蜜臀| 女人久久www免费人成看片| 国产在线视频一区二区| 视频区欧美日本亚洲| 亚洲av日韩在线播放| 亚洲精品自拍成人| 99精品久久久久人妻精品| 久久久久久亚洲精品国产蜜桃av| 一本色道久久久久久精品综合| 午夜激情久久久久久久| 国产高清videossex| 丰满迷人的少妇在线观看| 啦啦啦在线免费观看视频4| 在线观看舔阴道视频| 欧美日韩精品网址| 天天躁夜夜躁狠狠躁躁| 亚洲成人免费av在线播放| 久久亚洲国产成人精品v| 欧美日本中文国产一区发布| 国产免费现黄频在线看| cao死你这个sao货| 大香蕉久久成人网| 夜夜夜夜夜久久久久| 正在播放国产对白刺激| av视频免费观看在线观看| 考比视频在线观看| 老汉色av国产亚洲站长工具| 美女脱内裤让男人舔精品视频| 视频在线观看一区二区三区| 黄色片一级片一级黄色片| 日本猛色少妇xxxxx猛交久久| 久久ye,这里只有精品| netflix在线观看网站| 黄色视频在线播放观看不卡| 999久久久精品免费观看国产| 亚洲国产成人一精品久久久| 国产一区有黄有色的免费视频| 成年动漫av网址| 19禁男女啪啪无遮挡网站| 免费女性裸体啪啪无遮挡网站| a在线观看视频网站| 国产成人免费观看mmmm| 亚洲综合色网址| 精品国内亚洲2022精品成人 | 欧美97在线视频| 亚洲久久久国产精品| 少妇裸体淫交视频免费看高清 | 国产亚洲一区二区精品| 午夜免费成人在线视频| 日韩视频在线欧美| 考比视频在线观看| 国产一区二区三区综合在线观看| 国产伦人伦偷精品视频| 午夜精品久久久久久毛片777| 婷婷丁香在线五月| 黑人欧美特级aaaaaa片| 丝袜美足系列| 精品少妇一区二区三区视频日本电影| 日韩制服丝袜自拍偷拍| 久久精品国产亚洲av香蕉五月 | 美国免费a级毛片| 亚洲专区国产一区二区| videos熟女内射| 亚洲久久久国产精品| 午夜免费观看性视频| 王馨瑶露胸无遮挡在线观看| 欧美午夜高清在线| 51午夜福利影视在线观看| 中文字幕人妻熟女乱码| 在线精品无人区一区二区三| 蜜桃在线观看..| 亚洲人成77777在线视频| 国产色视频综合| 亚洲精品乱久久久久久| 久久 成人 亚洲| 热re99久久精品国产66热6| 亚洲成人免费av在线播放| 考比视频在线观看| 99精国产麻豆久久婷婷| av超薄肉色丝袜交足视频| 亚洲人成电影免费在线| 免费久久久久久久精品成人欧美视频| 亚洲 欧美一区二区三区| 免费在线观看视频国产中文字幕亚洲 | av福利片在线| 国产成人免费观看mmmm| 国产精品 欧美亚洲| 久久热在线av| 淫妇啪啪啪对白视频 | 国产亚洲精品一区二区www | 另类亚洲欧美激情| 日韩三级视频一区二区三区| 国产精品自产拍在线观看55亚洲 | 99精品久久久久人妻精品| 性少妇av在线| 欧美中文综合在线视频| 免费在线观看日本一区| 咕卡用的链子| 一二三四在线观看免费中文在| 每晚都被弄得嗷嗷叫到高潮| 久久精品国产a三级三级三级| 天天躁狠狠躁夜夜躁狠狠躁| 蜜桃国产av成人99| 免费在线观看影片大全网站| 久久人人爽人人片av| 丝瓜视频免费看黄片| 国产免费av片在线观看野外av| 欧美日本中文国产一区发布| 久久精品国产亚洲av高清一级| 亚洲熟女精品中文字幕| 亚洲成国产人片在线观看| 亚洲av电影在线观看一区二区三区| 亚洲精品国产av成人精品| 国产精品久久久久久精品古装| 欧美av亚洲av综合av国产av| 亚洲精品国产精品久久久不卡| 天天操日日干夜夜撸| tocl精华| 黄网站色视频无遮挡免费观看| 午夜精品国产一区二区电影| 男女午夜视频在线观看| 亚洲国产精品999| 国产成人精品无人区| 国产精品秋霞免费鲁丝片| 男女国产视频网站| 超碰成人久久| 精品亚洲成a人片在线观看| 女性生殖器流出的白浆| 黄片大片在线免费观看| 狂野欧美激情性xxxx| 国产免费一区二区三区四区乱码| 亚洲精品国产av蜜桃| 十八禁网站网址无遮挡| 午夜福利在线观看吧| 亚洲欧洲日产国产| 国产亚洲精品一区二区www | 国产成人免费观看mmmm| 一级毛片精品| 老汉色∧v一级毛片| 精品国产一区二区三区四区第35| 一个人免费在线观看的高清视频 | 啦啦啦中文免费视频观看日本| 老司机靠b影院| 宅男免费午夜| 母亲3免费完整高清在线观看| 午夜久久久在线观看| 成年人午夜在线观看视频| 国产精品久久久av美女十八| 天天添夜夜摸| 18禁黄网站禁片午夜丰满| 精品国产乱子伦一区二区三区 | 1024视频免费在线观看| av网站在线播放免费| 99国产精品一区二区三区| av在线老鸭窝| 色94色欧美一区二区| 菩萨蛮人人尽说江南好唐韦庄| 亚洲中文日韩欧美视频| 成人av一区二区三区在线看 | 亚洲av电影在线观看一区二区三区| 国产精品99久久99久久久不卡| 91成人精品电影| 午夜福利在线免费观看网站| 美女主播在线视频| 久久 成人 亚洲| 婷婷色av中文字幕| 人人妻人人爽人人添夜夜欢视频| 在线观看免费日韩欧美大片| 精品国产乱子伦一区二区三区 | 建设人人有责人人尽责人人享有的| 波多野结衣av一区二区av| 久久中文字幕一级| 欧美黄色淫秽网站| 午夜91福利影院| 精品亚洲乱码少妇综合久久| 日韩视频一区二区在线观看| 99热网站在线观看| 无限看片的www在线观看| 无遮挡黄片免费观看| 黄频高清免费视频| 免费久久久久久久精品成人欧美视频| 在线 av 中文字幕| 久久国产精品男人的天堂亚洲| 高清黄色对白视频在线免费看| 精品一区二区三区四区五区乱码| 91精品国产国语对白视频| 丝袜脚勾引网站| 午夜福利乱码中文字幕| 永久免费av网站大全| 熟女少妇亚洲综合色aaa.| 亚洲欧美精品综合一区二区三区| 在线观看免费视频网站a站| 欧美精品一区二区大全| 黑人猛操日本美女一级片| 欧美精品亚洲一区二区| 欧美午夜高清在线| 91国产中文字幕| 脱女人内裤的视频| 老熟妇仑乱视频hdxx| 国产精品.久久久| 久久久精品国产亚洲av高清涩受| 精品高清国产在线一区| 亚洲少妇的诱惑av| 成年美女黄网站色视频大全免费| 中文精品一卡2卡3卡4更新| 手机成人av网站| 久久精品国产亚洲av香蕉五月 | 国产成人欧美| 国产亚洲欧美在线一区二区| 色婷婷av一区二区三区视频| 黑人操中国人逼视频| 国产野战对白在线观看| 久久人妻福利社区极品人妻图片| 成年女人毛片免费观看观看9 | 国产精品一区二区在线观看99| 一本色道久久久久久精品综合| 女性被躁到高潮视频| 久热这里只有精品99| 嫁个100分男人电影在线观看| 高潮久久久久久久久久久不卡| 午夜视频精品福利| 99精国产麻豆久久婷婷| 中文字幕色久视频| 老司机影院成人| 国产欧美日韩一区二区三 | 精品国内亚洲2022精品成人 | 欧美日韩视频精品一区| 啦啦啦视频在线资源免费观看| av在线app专区| 国产免费福利视频在线观看| 伦理电影免费视频| 国产成人欧美| 天天添夜夜摸| 一级片'在线观看视频| av免费在线观看网站| 宅男免费午夜| 国产熟女午夜一区二区三区| 日韩 亚洲 欧美在线| av视频免费观看在线观看| 99香蕉大伊视频| netflix在线观看网站| 美女午夜性视频免费| 亚洲精品日韩在线中文字幕| 91老司机精品| www.自偷自拍.com| 在线看a的网站| 2018国产大陆天天弄谢| 水蜜桃什么品种好| 十八禁网站免费在线| 精品久久久久久电影网| 欧美性长视频在线观看| 91麻豆精品激情在线观看国产 | 午夜激情久久久久久久| 国产亚洲av高清不卡| 一本一本久久a久久精品综合妖精| 69精品国产乱码久久久| 美女福利国产在线| 天天影视国产精品| 日韩人妻精品一区2区三区| 黄片播放在线免费| 91字幕亚洲| 亚洲欧美精品自产自拍| √禁漫天堂资源中文www| 亚洲精品国产区一区二| 久久女婷五月综合色啪小说| 高清av免费在线| 久久久精品区二区三区| 亚洲精品美女久久av网站| 一本色道久久久久久精品综合| 侵犯人妻中文字幕一二三四区| 亚洲,欧美精品.| 久久这里只有精品19| 色婷婷久久久亚洲欧美| av在线老鸭窝| 汤姆久久久久久久影院中文字幕| 精品国产乱子伦一区二区三区 | 亚洲伊人久久精品综合| 老司机深夜福利视频在线观看 | 午夜久久久在线观看| 在线观看www视频免费| 五月开心婷婷网| 日韩大片免费观看网站| 欧美日韩av久久| 老汉色av国产亚洲站长工具| 日本vs欧美在线观看视频| 超碰成人久久| 久久av网站| 日韩大片免费观看网站| 亚洲熟女毛片儿| 亚洲av电影在线进入| www.熟女人妻精品国产| 纯流量卡能插随身wifi吗| 性色av乱码一区二区三区2| 午夜久久久在线观看| 国产无遮挡羞羞视频在线观看| 99精国产麻豆久久婷婷| 免费av中文字幕在线| 日韩中文字幕欧美一区二区| 大香蕉久久网| 久久久国产一区二区| 丰满饥渴人妻一区二区三| 叶爱在线成人免费视频播放| 精品一区在线观看国产| 99国产极品粉嫩在线观看| 丝袜人妻中文字幕| 色播在线永久视频| 日韩欧美一区视频在线观看| 99精品久久久久人妻精品| 我要看黄色一级片免费的| 国产精品国产av在线观看| 精品国产乱码久久久久久小说| 亚洲精品一卡2卡三卡4卡5卡 | 制服诱惑二区| 中文字幕精品免费在线观看视频| 欧美老熟妇乱子伦牲交| 国产成人免费观看mmmm| 国产视频一区二区在线看| 女性生殖器流出的白浆| 国产成人啪精品午夜网站| 欧美+亚洲+日韩+国产| 男女下面插进去视频免费观看| 免费在线观看日本一区| 日韩电影二区| 亚洲一卡2卡3卡4卡5卡精品中文| 又黄又粗又硬又大视频| 精品乱码久久久久久99久播| 久久久久久亚洲精品国产蜜桃av| 99热全是精品| 亚洲人成77777在线视频| 中文字幕色久视频| 国产精品久久久久久精品古装| 成年美女黄网站色视频大全免费| 久久99热这里只频精品6学生| 久久狼人影院| 久久精品久久久久久噜噜老黄| 黄色视频在线播放观看不卡| 女人久久www免费人成看片| 精品一区在线观看国产| 日本五十路高清| 两性夫妻黄色片| 一本大道久久a久久精品| 高清视频免费观看一区二区| 国产亚洲欧美精品永久| 欧美日韩黄片免| 永久免费av网站大全| 99久久99久久久精品蜜桃| 蜜桃在线观看..| 午夜激情久久久久久久| 男人舔女人的私密视频| 亚洲av欧美aⅴ国产| 色播在线永久视频| 免费在线观看完整版高清| 两人在一起打扑克的视频| 国产亚洲一区二区精品| 91av网站免费观看| 9色porny在线观看| 一区福利在线观看| 99久久精品国产亚洲精品| 老汉色∧v一级毛片| 欧美精品一区二区大全| 男人舔女人的私密视频| 91av网站免费观看| 久久青草综合色| 欧美变态另类bdsm刘玥| 一本色道久久久久久精品综合| 久久久久久久国产电影| 又大又爽又粗| 9191精品国产免费久久| 高清视频免费观看一区二区| 熟女少妇亚洲综合色aaa.| 精品亚洲乱码少妇综合久久| 国产欧美日韩一区二区精品| 9热在线视频观看99| 精品国产一区二区三区久久久樱花| 国产一区二区 视频在线| 欧美在线黄色| 午夜免费成人在线视频| 制服诱惑二区| 老司机影院成人| 国产精品自产拍在线观看55亚洲 | 人人妻人人爽人人添夜夜欢视频| 在线av久久热| 婷婷成人精品国产| 久久久久久久久久久久大奶| 99国产精品99久久久久| 亚洲,欧美精品.| 久久久精品94久久精品| 国产成人系列免费观看| 日本一区二区免费在线视频| 建设人人有责人人尽责人人享有的| 国产精品av久久久久免费| 国产av国产精品国产| 午夜福利,免费看| 久久久久精品国产欧美久久久 | 成年美女黄网站色视频大全免费| 一级毛片精品| 国产一卡二卡三卡精品| 免费久久久久久久精品成人欧美视频| 精品亚洲乱码少妇综合久久| 两性夫妻黄色片| 久久人人爽人人片av| 最近最新中文字幕大全免费视频| 在线观看舔阴道视频| 老汉色av国产亚洲站长工具| 成人18禁高潮啪啪吃奶动态图| 久久狼人影院| 色综合欧美亚洲国产小说| 老熟妇仑乱视频hdxx| 国产精品偷伦视频观看了| a级毛片黄视频| 99精国产麻豆久久婷婷| 日日摸夜夜添夜夜添小说| 99精品欧美一区二区三区四区| 久久久久精品国产欧美久久久 | 最近中文字幕2019免费版| 美女福利国产在线| 日本黄色日本黄色录像| 国产精品一区二区免费欧美 | 国产激情久久老熟女| 国产xxxxx性猛交| 国产精品亚洲av一区麻豆| 亚洲伊人色综图| 18禁黄网站禁片午夜丰满| 女人被躁到高潮嗷嗷叫费观| 国产精品偷伦视频观看了| 老熟妇乱子伦视频在线观看 | 女人精品久久久久毛片| 日韩制服丝袜自拍偷拍| 精品熟女少妇八av免费久了| 亚洲熟女毛片儿| 最新的欧美精品一区二区| 久热爱精品视频在线9| 纵有疾风起免费观看全集完整版| 超碰成人久久| 免费观看av网站的网址| 欧美人与性动交α欧美软件| 亚洲中文av在线| 一区二区三区四区激情视频| 国产一卡二卡三卡精品| 亚洲欧洲精品一区二区精品久久久| 午夜福利一区二区在线看| 午夜成年电影在线免费观看| 美女中出高潮动态图| 亚洲国产欧美在线一区| av在线老鸭窝| 日本vs欧美在线观看视频| 久久精品国产a三级三级三级| 人人妻人人澡人人看| 男人舔女人的私密视频| 久久精品久久久久久噜噜老黄| 亚洲精品国产av蜜桃| 桃红色精品国产亚洲av| 秋霞在线观看毛片| 最新的欧美精品一区二区| 久久人人爽人人片av| 又大又爽又粗| 一边摸一边做爽爽视频免费| 国产成人影院久久av| 国产日韩欧美视频二区| 国产黄色免费在线视频| 亚洲精品一卡2卡三卡4卡5卡 | 精品一区在线观看国产| 欧美变态另类bdsm刘玥| 大片电影免费在线观看免费| 日韩视频在线欧美| 欧美另类亚洲清纯唯美| 丝袜喷水一区| 欧美一级毛片孕妇| 日韩视频在线欧美| 亚洲avbb在线观看| 久久久精品免费免费高清| 正在播放国产对白刺激| 久久天堂一区二区三区四区| 一级毛片电影观看| 亚洲成人国产一区在线观看| 欧美大码av| 免费在线观看完整版高清| 91老司机精品| 亚洲国产日韩一区二区| 亚洲 欧美一区二区三区| 五月开心婷婷网| 天天躁日日躁夜夜躁夜夜| 丝袜美腿诱惑在线| 在线十欧美十亚洲十日本专区| 欧美在线一区亚洲| 99热全是精品| 欧美日韩视频精品一区| 久久人人97超碰香蕉20202| 亚洲国产av影院在线观看| 亚洲成人免费av在线播放| 久久久久视频综合| 精品人妻在线不人妻| 国产人伦9x9x在线观看| 在线精品无人区一区二区三| 国产成人免费无遮挡视频| 亚洲性夜色夜夜综合| av不卡在线播放| 国产福利在线免费观看视频| 手机成人av网站| 亚洲一区二区三区欧美精品| 999久久久精品免费观看国产| 99精品久久久久人妻精品| 午夜免费成人在线视频| 中文字幕色久视频| 80岁老熟妇乱子伦牲交| 亚洲国产精品一区三区| 岛国毛片在线播放| 一本久久精品| 国产精品 国内视频| 亚洲精品一卡2卡三卡4卡5卡 | 青青草视频在线视频观看| 在线天堂中文资源库| 国产成人a∨麻豆精品| 51午夜福利影视在线观看| 国产精品一区二区精品视频观看| 汤姆久久久久久久影院中文字幕| 激情视频va一区二区三区| 免费在线观看黄色视频的| 女人高潮潮喷娇喘18禁视频| 精品欧美一区二区三区在线| 久久久欧美国产精品| 看免费av毛片| 免费高清在线观看日韩| 男人爽女人下面视频在线观看| 最新的欧美精品一区二区| 每晚都被弄得嗷嗷叫到高潮| 国产精品秋霞免费鲁丝片| 久久国产精品男人的天堂亚洲| 欧美国产精品一级二级三级| 中文字幕人妻丝袜一区二区| 色播在线永久视频| 大型av网站在线播放| 日韩中文字幕欧美一区二区|