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

    高階數(shù)值流形方法的速度公式

    2016-08-09 08:31:13蘇海東
    長江科學(xué)院院報 2016年7期
    關(guān)鍵詞:構(gòu)形流形結(jié)點

    蘇海東

    ?

    高階數(shù)值流形方法的速度公式

    蘇海東

    (長江科學(xué)院材料與結(jié)構(gòu)研究所,武漢430010)

    高階數(shù)值流形方法可以顯著提高結(jié)構(gòu)計算精度,但目前在涉及大位移的動力分析中往往得到精度很差、甚至不正確的速度結(jié)果?;谄矫嫒切螖?shù)學(xué)網(wǎng)格和一階多項式覆蓋函數(shù),通過一個剛體桿件旋轉(zhuǎn)算例探討其中的原因,得出必須考慮構(gòu)形坐標(biāo)變化對速度的影響,并提出高階流形法的3種速度處理方法及相應(yīng)的高階速度公式。該方法對一些在結(jié)點處增加廣義自由度的類似方法(如廣義有限元)的幾何非線性問題分析也具有一定的參考價值。

    數(shù)值流形方法;高階多項式覆蓋函數(shù);大位移;速度公式;廣義自由度

    doi:10.11988/ckyyb.20150343

    1 研究背景

    作為計算力學(xué)的經(jīng)典和普遍使用的方法,有限元法廣泛應(yīng)用于實際工程計算和科學(xué)研究中已長達(dá)數(shù)十年,其特點是以待求的場變量作為網(wǎng)格結(jié)點上的未知數(shù),在結(jié)構(gòu)應(yīng)力分析中表現(xiàn)為位移自由度。近年來在有限元法的基礎(chǔ)上,數(shù)值流形方法[1]、擴展有限元[2]、廣義有限元[3]等新方法,通過在結(jié)點處增加廣義自由度來提高計算精度,在結(jié)構(gòu)的小變形應(yīng)力分析中取得了滿意的效果[4-6]。

    在上述新方法中,數(shù)值流形方法率先進(jìn)入了大位移或大變形的幾何非線性問題研究領(lǐng)域。研究表明,在結(jié)點上引入除常數(shù)位移場變量以外的多項式函數(shù)后(被稱為高階數(shù)值流形方法),這種廣義自由度對計算結(jié)果的影響逐漸顯現(xiàn)出來。作者在高階數(shù)值流形法的大變形靜力分析中曾指出[7],有必要在應(yīng)力累積中考慮構(gòu)形坐標(biāo)的變化,并相應(yīng)地提出了高階初應(yīng)力公式,通過懸臂梁大變形靜力分析的計算結(jié)果得到了驗證。

    進(jìn)一步討論,如果分析大位移條件下的高階數(shù)值流形法的動力問題,在速度計算中是否也需要進(jìn)行特殊處理呢?對于此問題,目前的研究均未涉及到,本文通過一個簡單例子展開討論。

    為使讀者對數(shù)值流形方法中的廣義自由度有所了解,以下首先簡要介紹高階數(shù)值流形方法。

    2 高階數(shù)值流形方法簡介

    數(shù)值流形方法[1](以下簡稱流形法)是我國留美學(xué)者石根華博士利用現(xiàn)代數(shù)學(xué)中流形分析的有限覆蓋技術(shù)建立起來的一種新的數(shù)值分析方法,具有廣闊的發(fā)展前景。流形法引入2套覆蓋網(wǎng)格:一是用于構(gòu)造物理場近似解的數(shù)學(xué)網(wǎng)格;另一個是表示材料邊界、用于定義積分區(qū)域的物理網(wǎng)格。2套網(wǎng)格相互獨立,只要求前者在空間上完全包容后者。

    流形法對各個覆蓋獨立定義覆蓋函數(shù),通過聯(lián)絡(luò)函數(shù)(又稱權(quán)函數(shù))將各個獨立覆蓋函數(shù)連起來形成總體的覆蓋函數(shù)。規(guī)定覆蓋Ui上的局部覆蓋函數(shù)為ui(x)(x∈Ui),它可以是多項式級數(shù)或其他級數(shù),隨著級數(shù)項數(shù)的增加及階數(shù)的提高,計算精度得以提高。這些覆蓋函數(shù)用權(quán)函數(shù)聯(lián)系起來,其含義是加權(quán)平均,權(quán)函數(shù)滿足

    因此總體覆蓋函數(shù)為

    式中n為覆蓋總數(shù)。

    現(xiàn)有的流形法多采用有限元網(wǎng)格來定義數(shù)學(xué)網(wǎng)格,與有限元任一結(jié)點相連的所有網(wǎng)格形成一個數(shù)學(xué)覆蓋,局部覆蓋函數(shù)定義于結(jié)點上形成結(jié)點覆蓋。這樣,任一原始的有限元網(wǎng)格就是各結(jié)點覆蓋的重疊部分,在此重疊區(qū)域(即有限元網(wǎng)格)內(nèi),覆蓋的權(quán)函數(shù)取為有限單元的形函數(shù)。本文以三角形數(shù)學(xué)網(wǎng)格為例,其形函數(shù)為面積坐標(biāo)表示的一階多項式Li(下文統(tǒng)一用上標(biāo)表示結(jié)點序號,注意與指數(shù)相區(qū)別,用下標(biāo)n-1,n,n+1表示構(gòu)形序號)。

    考慮在各個結(jié)點覆蓋上采用一階多項式的局部覆蓋函數(shù),即

    式中的a至f為多項式系數(shù),即待求未知數(shù),可以看出,這些未知數(shù)已沒有通常意義下的位移自由度的物理含義,因此被稱為廣義自由度。

    三角形網(wǎng)格內(nèi)的位移函數(shù)如式(3)所示,結(jié)點覆蓋函數(shù)用三角形網(wǎng)格的形函數(shù)連接成為二階多項式函數(shù),即

    式中Te為網(wǎng)格內(nèi)的插值函數(shù)(三角形網(wǎng)格的形函數(shù)與各單項式函數(shù)1,x,y之乘積)為以ai至fi表示的網(wǎng)格自由度向量。

    與有限元法一樣,流形法通過最小勢能原理建立平衡方程來進(jìn)行材料的變形和應(yīng)力分析,不同之處僅在于位移近似函數(shù)的表示方式。文獻(xiàn)[8]指出,流形法的動力分析采用了經(jīng)典的NewMark公式,只是將參數(shù)定為β=1/2和γ=1。其動力影響反映在慣性力引起的剛度矩陣(即通常所指的質(zhì)量矩陣)和荷載向量上,后者包含速度[1]。第n步計算完畢后,n+1步的初始速度為[1]

    式中:Vn為第n步的速度;Δ un為第n步的增量位移;Δt為時間步長。

    對于采用一階多項式覆蓋函數(shù)的高階流形法而言,其結(jié)點的速度覆蓋按式(4)計算后也具有式(2)的形式,同時,網(wǎng)格內(nèi)的速度分布具有式(3)的形式,只是系數(shù)不同。

    高階流形法由于在結(jié)點處引入了帶有廣義自由度的多項式覆蓋函數(shù),因而構(gòu)形坐標(biāo)的變化會對材料體的特征量(如應(yīng)力、速度)產(chǎn)生影響。作者已在文獻(xiàn)[7]中對高階流形法的初應(yīng)力公式進(jìn)行了修正,取得了很好的效果。以下針對速度問題,以剛性桿件的勻速旋轉(zhuǎn)為例進(jìn)行討論。

    3 高階流形法速度公式存在的問題

    如圖1所示,采用二維流形元網(wǎng)格計算剛性桿件,桿件長10 m,矩形截面的高度和寬度均為1 m。按平面應(yīng)力計算,彈性模量取大值以模擬剛性。桿件的左端中點為固定約束,從圖示的位置開始按1 rad/s的角速度順時針旋轉(zhuǎn),計算步長為0.01 s。

    圖1 剛性桿件及其流形元網(wǎng)格Fig.1 A rigid bar and its NMM meshes

    在文獻(xiàn)[9]編制的靜力分析程序的基礎(chǔ)上增加動力計算部分,分別取覆蓋函數(shù)為0階(即常規(guī)的位移自由度)和式(2)的一階多項式,計算圖1的桿件轉(zhuǎn)動,得到的桿件自由端中點的坐標(biāo)如表1所示。

    表1 原始方法的桿件自由端中點坐標(biāo)Table 1 Coordinates of the middle point at the free end of bar(conventional method)

    表1中可見,0階計算結(jié)果與理論值(x= 10cos(t),y=-10sin(t))很接近,說明對于剛體旋轉(zhuǎn)運動,用0階計算就已足夠準(zhǔn)確。而在同樣條件下,一階計算結(jié)果與理論值沒有可比性,旋轉(zhuǎn)速度明顯快于理論值。這表明,當(dāng)引入一階廣義自由度時,高階流形法的速度存在問題。

    為便于分析,下面以一維問題為例進(jìn)行討論。見圖2所示的一個一維網(wǎng)格。在第n-1步時結(jié)點1的初始位置為x1n-1,結(jié)點2的初始位置為x2n-1,網(wǎng)格內(nèi)任一點的坐標(biāo)為xn-1。位移計算完畢后,按式(4)計算得到的結(jié)點1和結(jié)點2的速度覆蓋函數(shù)即第n步的初始速度)分別為:

    式中a-1n,a-2n,b-1n,b-2n為構(gòu)形變化前的多項式系數(shù),則一維網(wǎng)格內(nèi)任意點的速度為

    式中N1,N2分別為對應(yīng)于結(jié)點1和結(jié)點2的一維網(wǎng)格形函數(shù)。

    圖2 從第n-1步移至第n步的一維網(wǎng)格Fig.2 Schematic diagram of one-dimensional mesh moving from the(n-1)thstep to the nthstep

    如圖2所示,在進(jìn)行第n步計算前,根據(jù)計算得到的位移量Δu將材料體從第n-1構(gòu)形更新至第n構(gòu)形,附著在材料體上的數(shù)學(xué)網(wǎng)格結(jié)點隨材料體移動至新的結(jié)點位置x1n和x2n,網(wǎng)格內(nèi)的任一點由xn-1移至xn。若按原始的高階流形法速度公式,結(jié)點1和結(jié)點2的第n步速度覆蓋函數(shù)分別為其中,在構(gòu)形更新過程中的多項式系數(shù)保持不變,即??紤]到剛體運動中的形函數(shù)N1和

    作為第n步計算的初始速度。然而表1表明,這種速度處理是不正確的。

    實際上,材料的構(gòu)形坐標(biāo)發(fā)生了變化,即

    從以上分析可知,由坐標(biāo)移動產(chǎn)生的速度誤差正是造成表1中的一階速度計算結(jié)果不正確的原因,因此需要對高階流形法的速度進(jìn)行修正。

    4 高階流形法速度公式的修正

    高階流形法速度公式的修正考慮以下3種方法。

    4.1方法1——結(jié)點覆蓋上的速度修正

    比較式(8)與式(6)可見,欲使Vn=V-n,則

    考慮到式(11)中的Δu隨網(wǎng)格內(nèi)的材料點位置而變化,作為一種近似,在上述公式中可取結(jié)點處的Δu1和Δu2,即結(jié)點上的覆蓋多項式系數(shù)修正為

    推廣到二維問題,設(shè)結(jié)點i(i=1,2,3)的第n步速度覆蓋函數(shù)為

    式中ain至fi

    n等為系數(shù),三角形網(wǎng)格內(nèi)的速度分布為

    其中Lin為在第n步構(gòu)形下的三角形網(wǎng)格形函數(shù)。

    與一維情況的式(10)類似,相對于V-n的速度差值為

    設(shè)三角形網(wǎng)格第i結(jié)點的位移增量為Δui和Δvi,相應(yīng)于一維公式(12),二維第n步的速度公式為

    式中a-i

    n至-fin為構(gòu)形變化前的多項式系數(shù)。

    重新計算如圖1所示的剛性桿件旋轉(zhuǎn)問題,計算結(jié)果見表2,可見相對表1有很大改善,表明了覆蓋修正公式的作用。但與理論值相比還有一定的誤差,如t=1 s時,x坐標(biāo)相差約7 cm,y坐標(biāo)相差約5 cm。顯而易見,這是由速度覆蓋公式中近似取結(jié)點處的Δu和Δv所造成的,而準(zhǔn)確的Δu和Δv在網(wǎng)格內(nèi)應(yīng)隨坐標(biāo)變化。

    4.2方法2——速度荷載修正

    仍考慮如圖2所示的一維網(wǎng)格。網(wǎng)格內(nèi)每一點的位移變化為

    表2 結(jié)點覆蓋速度修正的桿件自由端中點坐標(biāo)Table 2 Coordinates of the middle point at the free end of bar(velocity modification of nodal covers)

    其中,系數(shù) Δain和 Δbi

    n是將它們均與Δ u1和Δ u2有關(guān),這表明速度修正量+已不局限于結(jié)點覆蓋本身,而與一維網(wǎng)格內(nèi)的另一結(jié)點覆蓋有關(guān)。對于處在相鄰網(wǎng)格的交界點,在不同網(wǎng)格內(nèi)的速度修正量是不同的,這將造成修正后的速度不連續(xù)。

    同理,對于二維三角形網(wǎng)格,將

    及形函數(shù)L1

    n,L2n,L3

    n的表達(dá)式代入式(15)可得

    同樣,上式中的Δain,Δbin,Δci

    n與網(wǎng)格各結(jié)點的位移增量有關(guān),對于處在相鄰網(wǎng)格的交界處,在不同網(wǎng)格內(nèi)的速度修正量也是不同的,因此,無法像方法1那樣對某一結(jié)點的速度覆蓋函數(shù)進(jìn)行簡單修正,需要尋找其他方法。

    考察流形法的基本控制方程,由最小勢能原理推導(dǎo)得出[1],其中考慮慣性力的勢能表達(dá)式(仍以一維問題為例)為

    式中:a=2u-2V,為加速度;ρ為密度;A為流形Δt2Δt元面積??疾炱渲械牡?項∫ρu2VdA

    AΔt,一般的做法是[1],將網(wǎng)格內(nèi)的位移u=uTeTTe及速度V=TeVe代入,取極值后得到關(guān)于單元速度的荷載表達(dá)式為式中:Ve為網(wǎng)格結(jié)點速度的廣義自由度向量;M= ∫Aρ TTeTedA為單元質(zhì)量矩陣。

    實際上,式(22)為單元的積分表達(dá)式,只要求計算Fv的積分,并未要求將V寫成包含結(jié)點速度的廣義自由度向量Ve的形式??梢詢H將u=uTeTTe代入并取極值得這就給網(wǎng)格內(nèi)的速度修正(而不是結(jié)點速度修正)提供了依據(jù)。

    修正后的速度荷載為

    由于在計算速度荷載Fv時進(jìn)行了速度修正,因而得到的增量位移Δ是比較準(zhǔn)確的。但式(4)中的Vn仍采用了構(gòu)形變化前的系數(shù),并沒有進(jìn)行修正,因此計算出來的Vn+1仍然不準(zhǔn)確,進(jìn)而在計算下一步(第n+1步)的Fv時產(chǎn)生誤差而生成不平衡荷載。

    下面考慮一種方法,在本步(第n步)式(24)的Fv中對速度增加一個修正量αΔV(α為系數(shù)),成為Vn-ΔV+αΔV,使得計算出的增量位移Δ un包含對Vn的修正,從而預(yù)先消除上述不平衡荷載。由于流形法計算要求每步的增量位移很小,這樣,相鄰計算步的速度變化量也很小,因此,基于Vn+1≈Vn的假設(shè),式(4)可看成即2 V對應(yīng)于,可見,由修正量αΔV引起的

    nFv荷載,其計算出的位移值Δ un,是對Vn的修正量的2倍。因此,系數(shù)α應(yīng)取為0.5。

    再計算剛性桿件旋轉(zhuǎn),計算結(jié)果見表3,可見相對表2又有很大改善,與理論值更接近。

    表3 速度荷載修正的桿件自由端中點坐標(biāo)Table 3 Coordinates of the middle point at the free end ofbar(modification of velocity loads)

    4.3方法3——網(wǎng)格速度修正

    對上述思路進(jìn)一步擴展,認(rèn)為同一結(jié)點在不同的網(wǎng)格中的速度不同,即速度不再連續(xù)。計算程序不再以結(jié)點的形式記錄速度,而以網(wǎng)格內(nèi)的速度信息來記錄,即使是連接相鄰網(wǎng)格的同一結(jié)點,它在各網(wǎng)格內(nèi)記錄的速度信息也是不同的。

    對于二維三角形網(wǎng)格情況,按式(20)計算網(wǎng)格內(nèi)的結(jié)點速度修正量,對網(wǎng)格內(nèi)的各結(jié)點Vn直接進(jìn)行修正(對于同一結(jié)點,在不同網(wǎng)格內(nèi)的修正值不一樣),然后代入式(23)計算速度荷載。該步計算完畢后,基于網(wǎng)格內(nèi)的速度值按式(4)計算新的速度。

    剛性桿件旋轉(zhuǎn)的計算結(jié)果見表4,在保留小數(shù)點后三位有效數(shù)字時,與0階的計算結(jié)果完全相同,說明已經(jīng)完全排除了線性項誤差的影響。

    表4 網(wǎng)格速度修正的桿件自由端中點坐標(biāo)Table 4 Coordinates of the middle point at the free end of bar(modification of mesh velocities)

    5結(jié)語

    對于采用一階多項式覆蓋函數(shù)的高階流形法,本文提出了3種速度修正方法:

    方法1(結(jié)點覆蓋修正)最簡單,但計算精度較差;方法2(速度荷載修正)的計算精度基本令人滿意,并能保持速度場的連續(xù)性,但要求相鄰計算步的速度差別小;方法3(網(wǎng)格速度修正)計算精度最高,但速度在相鄰網(wǎng)格間有間斷,不再保持連續(xù)性。3種方法各有優(yōu)缺點,在實際計算中可視情況選用。二階或更高階多項式覆蓋函數(shù)情況,可采用類似方法。另外,本文的研究是針對采用整體坐標(biāo)來表示多項式覆蓋函數(shù)時所產(chǎn)生的高階速度問題,若采用局部化的覆蓋函數(shù)形式[10],也許會減弱甚至避免上述問題,這還需進(jìn)一步研究。

    在結(jié)點處增加廣義自由度的各種方法中,數(shù)值流形方法已在大位移和大變形的計算中領(lǐng)先一步,本文方法對于其他具有廣義自由度的方法(如廣義有限元)也具有一定的參考價值。

    致謝:感謝數(shù)值流形方法的發(fā)明者石根華博士的指導(dǎo)。感謝國家自然科學(xué)基金(10772034)的資助。

    [1] 石根華.?dāng)?shù)值流形方法與非連續(xù)變形分析[M].裴覺民,譯.北京:清華大學(xué)出版社,1997.

    [2]DUARTE C A,BABUSKA I,ODEN J T.Generalized Finite Element Methods for Three-dimensional Structural Mechanics Problems[J].Computer&Structures,2000,77:215-232.

    [3]DAUX C,MOES N,DOLBOW J,et al.Arbitrary Branched and Intersecting Cracks with the Extended Finite Element Method[J].International Journal for Numerical Methods in Engineering,2000,48:1741-1760.

    [4] 張湘?zhèn)ィ聽帢s,呂文閣.?dāng)?shù)值流形方法研究及應(yīng)用進(jìn)展[J].力學(xué)進(jìn)展,2010,40(1):1-12.

    [5]李錄賢,王鐵軍.?dāng)U展有限元法(XFEM)及其應(yīng)用[J].力學(xué)進(jìn)展,2005,35(1):5-20.

    [6]李錄賢,劉書靜,張慧華.廣義有限元方法研究進(jìn)展[J].應(yīng)用力學(xué)學(xué)報,2009,26(1):96-108.

    [7]蘇海東,崔建華,謝小玲.高階數(shù)值流形方法的初應(yīng)力公式[J].計算力學(xué)學(xué)報,2010,27(2):270-274.

    [8] WANG C Y,CHUANG C C,SHENG J.Time Integration Theories for the DDA Method with Finite Element Meshes [C]∥SALAMI M R,BANKS D.Proceedings of the First International Conference on Analysis of Discontinuous Deformation.Analysis(DDA)and Simulations of Discontinuous Media.Berkeley,CA:TSI Press,Albuquerque,NM,1996:263-287.

    [9]蘇海東,謝小玲,陳琴.高階數(shù)值流形方法在結(jié)構(gòu)靜力分析中的應(yīng)用研究[J].長江科學(xué)院院報,2005,22(5):74-77.

    [10]林紹忠,祁勇峰,蘇海東.?dāng)?shù)值流形法中覆蓋函數(shù)的改進(jìn)形式及其應(yīng)用[J].長江科學(xué)院院報,2006,23(6):55-58.

    (編輯:王慰)

    Velocity Equations for High-order Numerical Manifold Method

    SU Hai-dong
    (Department of Material and Structure,Yangtze River Scientific Research Institute,Wuhan430010,China)

    Computational accuracy of structure deformation can be improved greatly by the high-order Numerical Manifold Method(NMM).However,poor accuracy or even incorrect velocity results were obtained in the dynamic analysis involved in large displacement.Based on 2-D triangular mathematical meshes and 1-order polynomial cover functions,the reason of the above cases is discussed through an example of rotation of a rigid bar in this paper.Three treatments and the corresponding equations for high-order velocities are presented for the first time,reflecting the change of configuration coordinates under large displacement.The high-order numerical manifold method is useful to other methods such as Generalized Finite Element Method(GFEM)which introduces generalized freedoms at nodes when solving geometric nonlinear problems.

    Numerical Manifold Method(NMM);high-order polynomial cover function;large displacement;velocity equation;generalized degree of freedom

    O33;O34

    A

    1001-5485(2016)07-0121-05

    2015-04-22;

    2015-06-04

    國家自然科學(xué)基金項目(10772034)

    蘇海東(1968-),男,湖北武漢人,教授級高級工程師,博士,主要從事水工結(jié)構(gòu)數(shù)值分析工作和計算方法研究,(電話)027-82927167(電子信箱)suhd@mail.crsri.cn。

    猜你喜歡
    構(gòu)形流形結(jié)點
    雙星跟飛立體成像的構(gòu)形保持控制
    緊流形上的Schr?dinger算子的譜間隙估計
    通有構(gòu)形的特征多項式
    迷向表示分為6個不可約直和的旗流形上不變愛因斯坦度量
    Nearly Kaehler流形S3×S3上的切觸拉格朗日子流形
    Ladyzhenskaya流體力學(xué)方程組的確定模與確定結(jié)點個數(shù)估計
    對一個幾何構(gòu)形的探究
    基于多故障流形的旋轉(zhuǎn)機械故障診斷
    基于Raspberry PI為結(jié)點的天氣云測量網(wǎng)絡(luò)實現(xiàn)
    甲骨文構(gòu)形法研究
    ponron亚洲| 99九九线精品视频在线观看视频| 亚洲最大成人av| 亚洲色图av天堂| 国产精品一二三区在线看| 露出奶头的视频| av免费在线看不卡| 国产av在哪里看| 小说图片视频综合网站| 午夜精品在线福利| 国产精品久久视频播放| 国产精品女同一区二区软件| 深夜精品福利| 亚洲av.av天堂| 少妇丰满av| 婷婷精品国产亚洲av在线| 国产视频一区二区在线看| 午夜福利在线观看免费完整高清在 | 99热这里只有是精品在线观看| 免费在线观看成人毛片| 青春草视频在线免费观看| 国产毛片a区久久久久| 国产爱豆传媒在线观看| 国产精品一区二区三区四区久久| 国产亚洲91精品色在线| 亚洲av.av天堂| 亚洲av一区综合| 给我免费播放毛片高清在线观看| 精品人妻一区二区三区麻豆 | 久久久久久久久久黄片| 亚洲精品久久国产高清桃花| 内地一区二区视频在线| 日韩精品中文字幕看吧| 91av网一区二区| 亚洲国产欧美人成| 亚洲国产精品成人久久小说 | 变态另类成人亚洲欧美熟女| 国内揄拍国产精品人妻在线| 成人午夜高清在线视频| 国产v大片淫在线免费观看| 国产av麻豆久久久久久久| 变态另类成人亚洲欧美熟女| 国产精品久久久久久精品电影| 美女xxoo啪啪120秒动态图| 午夜精品在线福利| 国产探花在线观看一区二区| 国产一级毛片七仙女欲春2| 在线观看一区二区三区| 两性午夜刺激爽爽歪歪视频在线观看| 18禁在线无遮挡免费观看视频 | 国产精品,欧美在线| 亚洲精品一卡2卡三卡4卡5卡| 久久人妻av系列| 给我免费播放毛片高清在线观看| 久久久久久国产a免费观看| 久久久午夜欧美精品| 国产白丝娇喘喷水9色精品| 观看免费一级毛片| 美女高潮的动态| 男人和女人高潮做爰伦理| 久久九九热精品免费| 久久精品国产鲁丝片午夜精品| 免费看日本二区| 国产成人影院久久av| 国产不卡一卡二| 欧美激情久久久久久爽电影| 99热这里只有精品一区| 99久国产av精品国产电影| 丰满乱子伦码专区| 久久精品国产自在天天线| 国产成人精品久久久久久| 欧美+日韩+精品| 国产精品一区二区免费欧美| 人妻久久中文字幕网| 能在线免费观看的黄片| 校园春色视频在线观看| 大香蕉久久网| 国产av一区在线观看免费| 久久久久国内视频| 日日啪夜夜撸| 久久99热6这里只有精品| 日本爱情动作片www.在线观看 | 亚洲精品456在线播放app| 欧美成人a在线观看| 91久久精品国产一区二区三区| 久久亚洲精品不卡| 欧美一区二区精品小视频在线| 欧美绝顶高潮抽搐喷水| 日韩欧美在线乱码| 日韩av不卡免费在线播放| h日本视频在线播放| 一区二区三区免费毛片| 97人妻精品一区二区三区麻豆| 中国美白少妇内射xxxbb| 成人鲁丝片一二三区免费| 成人一区二区视频在线观看| 91在线观看av| 亚洲成人精品中文字幕电影| 在线观看美女被高潮喷水网站| 国产中年淑女户外野战色| 99久久精品一区二区三区| 亚洲中文字幕一区二区三区有码在线看| 久久这里只有精品中国| 亚洲精品久久国产高清桃花| 成人美女网站在线观看视频| 一进一出抽搐gif免费好疼| 波多野结衣高清作品| 日韩欧美精品免费久久| 精品一区二区三区av网在线观看| 亚洲精品日韩在线中文字幕 | 日本欧美国产在线视频| 99精品在免费线老司机午夜| 亚洲成人精品中文字幕电影| 少妇猛男粗大的猛烈进出视频 | 日本成人三级电影网站| 国产高清视频在线观看网站| 国产人妻一区二区三区在| 免费人成在线观看视频色| 99久久久亚洲精品蜜臀av| 亚洲av成人av| 女人被狂操c到高潮| 伦理电影大哥的女人| 日韩欧美一区二区三区在线观看| 18禁裸乳无遮挡免费网站照片| 国产久久久一区二区三区| 成年版毛片免费区| 国产伦精品一区二区三区视频9| 亚洲婷婷狠狠爱综合网| 乱码一卡2卡4卡精品| 免费av不卡在线播放| 久久这里只有精品中国| 男插女下体视频免费在线播放| 国产大屁股一区二区在线视频| 国产成人freesex在线 | АⅤ资源中文在线天堂| 国产黄a三级三级三级人| 亚洲无线在线观看| 嫩草影院新地址| 一本久久中文字幕| 中文字幕精品亚洲无线码一区| 女的被弄到高潮叫床怎么办| 久久久精品大字幕| 免费av观看视频| .国产精品久久| 亚洲成人久久性| 亚洲真实伦在线观看| 91av网一区二区| 97在线视频观看| 最后的刺客免费高清国语| 免费看日本二区| 色5月婷婷丁香| 国产一区二区在线观看日韩| 最近在线观看免费完整版| 国产av麻豆久久久久久久| 亚洲真实伦在线观看| 成人午夜高清在线视频| 免费人成在线观看视频色| 人人妻,人人澡人人爽秒播| 一级毛片久久久久久久久女| 久久久久精品国产欧美久久久| 三级国产精品欧美在线观看| 我要搜黄色片| 亚洲第一电影网av| 国产乱人视频| 91久久精品国产一区二区成人| 欧美国产日韩亚洲一区| 精品一区二区三区视频在线观看免费| 男女边吃奶边做爰视频| 久久草成人影院| h日本视频在线播放| 亚洲一区二区三区色噜噜| 黑人高潮一二区| 最近的中文字幕免费完整| 国产成人一区二区在线| av免费在线看不卡| 欧美色欧美亚洲另类二区| 尾随美女入室| 久久精品国产清高在天天线| 最近2019中文字幕mv第一页| 一个人看视频在线观看www免费| 国产精品无大码| 精品免费久久久久久久清纯| 99热这里只有是精品50| 日韩欧美精品v在线| 国产高清有码在线观看视频| 国产av一区在线观看免费| 日本-黄色视频高清免费观看| 国产av不卡久久| 少妇人妻精品综合一区二区 | 亚洲av二区三区四区| 国产男人的电影天堂91| 最新在线观看一区二区三区| 九九在线视频观看精品| 亚洲内射少妇av| 欧美色视频一区免费| 天堂影院成人在线观看| 午夜精品国产一区二区电影 | 黄色一级大片看看| 国产 一区精品| 国产高清视频在线观看网站| 精品欧美国产一区二区三| 中文字幕人妻熟人妻熟丝袜美| 搡老熟女国产l中国老女人| 伦理电影大哥的女人| av国产免费在线观看| 国产精品日韩av在线免费观看| 美女xxoo啪啪120秒动态图| 免费搜索国产男女视频| videossex国产| 一进一出抽搐动态| 九九爱精品视频在线观看| av黄色大香蕉| 97超碰精品成人国产| 亚洲最大成人中文| 三级男女做爰猛烈吃奶摸视频| 悠悠久久av| 97碰自拍视频| 99在线人妻在线中文字幕| 亚洲精品久久国产高清桃花| 欧美极品一区二区三区四区| 亚洲成av人片在线播放无| 不卡一级毛片| 一进一出好大好爽视频| а√天堂www在线а√下载| 韩国av在线不卡| 蜜臀久久99精品久久宅男| 天美传媒精品一区二区| 久久久色成人| 欧美极品一区二区三区四区| 免费看光身美女| 亚洲欧美日韩无卡精品| 久久久久久久午夜电影| 乱系列少妇在线播放| 国产色婷婷99| 国产老妇女一区| 久久国产乱子免费精品| 国产精品人妻久久久影院| 久久午夜亚洲精品久久| 天堂av国产一区二区熟女人妻| 1000部很黄的大片| 日韩欧美一区二区三区在线观看| 露出奶头的视频| 在线看三级毛片| 亚洲人成网站在线观看播放| 亚洲国产高清在线一区二区三| 婷婷六月久久综合丁香| 美女内射精品一级片tv| 亚洲国产精品国产精品| 国语自产精品视频在线第100页| 91av网一区二区| 女人十人毛片免费观看3o分钟| 欧美色欧美亚洲另类二区| 亚洲中文日韩欧美视频| 伦理电影大哥的女人| 免费看a级黄色片| av免费在线看不卡| 亚洲国产精品sss在线观看| 亚洲内射少妇av| 久久精品影院6| 春色校园在线视频观看| 久久精品国产亚洲av涩爱 | 午夜日韩欧美国产| 欧美bdsm另类| 韩国av在线不卡| 日韩国内少妇激情av| 国产不卡一卡二| 国产中年淑女户外野战色| 免费看光身美女| 亚洲国产欧美人成| 国产精品伦人一区二区| 嫩草影视91久久| 欧美xxxx黑人xx丫x性爽| 亚洲精品日韩在线中文字幕 | 赤兔流量卡办理| 亚洲精品在线观看二区| 国产伦精品一区二区三区视频9| 日韩欧美免费精品| videossex国产| 免费在线观看成人毛片| 中文字幕免费在线视频6| 国产成人freesex在线 | 午夜视频国产福利| 91在线观看av| 亚洲美女黄片视频| 国产精品一区二区三区四区久久| 午夜福利视频1000在线观看| 国产单亲对白刺激| 国产激情偷乱视频一区二区| 久久精品国产亚洲av涩爱 | 亚洲国产欧美人成| 久久99热6这里只有精品| 免费在线观看影片大全网站| 欧美性感艳星| 欧美最黄视频在线播放免费| 国产精品无大码| 国产片特级美女逼逼视频| 国产亚洲精品综合一区在线观看| 搡女人真爽免费视频火全软件 | 色尼玛亚洲综合影院| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲综合色惰| 韩国av在线不卡| 中国美白少妇内射xxxbb| 国产成人福利小说| 蜜桃亚洲精品一区二区三区| 在线天堂最新版资源| av卡一久久| 精品久久久久久久久久免费视频| 成人特级黄色片久久久久久久| 99久久精品热视频| 欧美另类亚洲清纯唯美| 国产极品精品免费视频能看的| 久久久久久国产a免费观看| 精品一区二区三区av网在线观看| 国产高潮美女av| a级毛片a级免费在线| 国产av一区在线观看免费| 国产精品1区2区在线观看.| 亚洲自偷自拍三级| 如何舔出高潮| 日韩av不卡免费在线播放| 欧美一区二区国产精品久久精品| 特大巨黑吊av在线直播| 寂寞人妻少妇视频99o| 69人妻影院| av黄色大香蕉| 国产午夜精品久久久久久一区二区三区 | 免费电影在线观看免费观看| 欧美极品一区二区三区四区| 老熟妇仑乱视频hdxx| 香蕉av资源在线| 久久6这里有精品| 国产精品美女特级片免费视频播放器| 亚洲内射少妇av| 草草在线视频免费看| 欧美一区二区亚洲| 亚洲av免费高清在线观看| 干丝袜人妻中文字幕| 国产高潮美女av| 久久婷婷人人爽人人干人人爱| 欧美一区二区国产精品久久精品| 高清毛片免费观看视频网站| 中文字幕免费在线视频6| 亚洲精华国产精华液的使用体验 | 青春草视频在线免费观看| 波多野结衣高清无吗| 日日干狠狠操夜夜爽| 色尼玛亚洲综合影院| 亚洲欧美成人综合另类久久久 | 久久人人爽人人爽人人片va| 男人狂女人下面高潮的视频| 少妇人妻一区二区三区视频| 少妇猛男粗大的猛烈进出视频 | 亚洲人成网站在线播| 蜜臀久久99精品久久宅男| 卡戴珊不雅视频在线播放| 免费看光身美女| 日韩欧美三级三区| 一级黄片播放器| 在线天堂最新版资源| 色噜噜av男人的天堂激情| 成熟少妇高潮喷水视频| 天天躁夜夜躁狠狠久久av| 麻豆成人午夜福利视频| 国产精品1区2区在线观看.| 欧美zozozo另类| 国产精品野战在线观看| 在线观看一区二区三区| 51国产日韩欧美| 老熟妇乱子伦视频在线观看| 亚洲精品色激情综合| 久久精品国产鲁丝片午夜精品| 亚洲精品亚洲一区二区| 日韩欧美免费精品| 日韩精品中文字幕看吧| 国产男人的电影天堂91| 可以在线观看毛片的网站| 国产乱人视频| 午夜激情欧美在线| 小说图片视频综合网站| 久久久色成人| 久久精品国产亚洲av涩爱 | 中文资源天堂在线| 我的女老师完整版在线观看| 国产 一区 欧美 日韩| 亚洲精品在线观看二区| 性欧美人与动物交配| 亚洲,欧美,日韩| 国产高清激情床上av| 免费搜索国产男女视频| 蜜桃久久精品国产亚洲av| 两性午夜刺激爽爽歪歪视频在线观看| 免费不卡的大黄色大毛片视频在线观看 | 亚洲精华国产精华液的使用体验 | 插逼视频在线观看| 一进一出抽搐gif免费好疼| 卡戴珊不雅视频在线播放| 日韩欧美国产在线观看| 国产成人精品久久久久久| 国产高清有码在线观看视频| 日本三级黄在线观看| 一级黄色大片毛片| 欧美日韩乱码在线| 欧洲精品卡2卡3卡4卡5卡区| .国产精品久久| 深爱激情五月婷婷| 少妇熟女欧美另类| 国产激情偷乱视频一区二区| 国产成人91sexporn| 性色avwww在线观看| 亚洲,欧美,日韩| 亚洲av熟女| 97超碰精品成人国产| 我要搜黄色片| 午夜免费男女啪啪视频观看 | 丝袜喷水一区| 日本精品一区二区三区蜜桃| 亚洲久久久久久中文字幕| 男女做爰动态图高潮gif福利片| 国产精品一二三区在线看| 人人妻人人澡欧美一区二区| av女优亚洲男人天堂| 国产午夜精品论理片| 精品免费久久久久久久清纯| 欧美高清成人免费视频www| 99九九线精品视频在线观看视频| 中文字幕av在线有码专区| 深爱激情五月婷婷| 国产综合懂色| 精品乱码久久久久久99久播| 国产久久久一区二区三区| 婷婷精品国产亚洲av在线| 人妻少妇偷人精品九色| 九九在线视频观看精品| 岛国在线免费视频观看| 搡老熟女国产l中国老女人| 日韩制服骚丝袜av| 十八禁国产超污无遮挡网站| 日韩欧美 国产精品| 美女xxoo啪啪120秒动态图| 国内精品一区二区在线观看| 精品久久国产蜜桃| 免费看av在线观看网站| 在线看三级毛片| 国产精品亚洲美女久久久| 少妇的逼好多水| 99热这里只有精品一区| 国产亚洲精品久久久久久毛片| 亚州av有码| 最新在线观看一区二区三区| 毛片女人毛片| 国产成人aa在线观看| 五月玫瑰六月丁香| 性插视频无遮挡在线免费观看| 在线免费观看不下载黄p国产| 看免费成人av毛片| 色噜噜av男人的天堂激情| 99久久精品国产国产毛片| 大又大粗又爽又黄少妇毛片口| 两个人的视频大全免费| 综合色av麻豆| 国产麻豆成人av免费视频| 少妇熟女aⅴ在线视频| 欧美又色又爽又黄视频| 免费看美女性在线毛片视频| 国产v大片淫在线免费观看| 国产精品亚洲一级av第二区| 免费不卡的大黄色大毛片视频在线观看 | av国产免费在线观看| 在线观看午夜福利视频| 久久久国产成人精品二区| 日本黄大片高清| 国产精品日韩av在线免费观看| 亚洲精品456在线播放app| 在线观看66精品国产| 国产精品日韩av在线免费观看| 女人被狂操c到高潮| 国产高清不卡午夜福利| 最近视频中文字幕2019在线8| 欧美日本亚洲视频在线播放| 国模一区二区三区四区视频| 免费av观看视频| 变态另类成人亚洲欧美熟女| 超碰av人人做人人爽久久| 国产亚洲精品av在线| 成年女人永久免费观看视频| 亚洲中文日韩欧美视频| 一进一出好大好爽视频| 人妻久久中文字幕网| 中文资源天堂在线| 亚洲美女黄片视频| 亚洲成a人片在线一区二区| 亚洲欧美日韩东京热| 最新中文字幕久久久久| 亚洲熟妇中文字幕五十中出| 夜夜夜夜夜久久久久| 啦啦啦观看免费观看视频高清| 亚洲真实伦在线观看| 69av精品久久久久久| 日韩欧美国产在线观看| 亚洲电影在线观看av| 欧美在线一区亚洲| 九九爱精品视频在线观看| 12—13女人毛片做爰片一| 高清日韩中文字幕在线| 国产亚洲91精品色在线| 中出人妻视频一区二区| 久久精品国产亚洲av涩爱 | 成人永久免费在线观看视频| 99久久精品热视频| 亚洲不卡免费看| 一进一出抽搐gif免费好疼| 校园人妻丝袜中文字幕| 亚洲性久久影院| 菩萨蛮人人尽说江南好唐韦庄 | 午夜免费激情av| 国产熟女欧美一区二区| 天美传媒精品一区二区| 91在线观看av| АⅤ资源中文在线天堂| 性色avwww在线观看| 老熟妇乱子伦视频在线观看| 狂野欧美激情性xxxx在线观看| 久久久久精品国产欧美久久久| 搡女人真爽免费视频火全软件 | 中文字幕免费在线视频6| .国产精品久久| 99riav亚洲国产免费| 小蜜桃在线观看免费完整版高清| 国产精品日韩av在线免费观看| 免费黄网站久久成人精品| 天天躁夜夜躁狠狠久久av| 国产欧美日韩一区二区精品| 国产伦精品一区二区三区四那| 国产精品福利在线免费观看| 久久久久九九精品影院| 久久久久国产精品人妻aⅴ院| 91久久精品国产一区二区三区| 亚洲丝袜综合中文字幕| 午夜福利在线在线| 欧美xxxx黑人xx丫x性爽| 婷婷精品国产亚洲av| 97在线视频观看| 久久精品国产亚洲av天美| 精品人妻视频免费看| 禁无遮挡网站| 日韩,欧美,国产一区二区三区 | a级毛片免费高清观看在线播放| 看片在线看免费视频| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲成a人片在线一区二区| 可以在线观看毛片的网站| 女生性感内裤真人,穿戴方法视频| 日本黄色片子视频| 国产午夜福利久久久久久| 春色校园在线视频观看| 免费人成在线观看视频色| 日韩精品有码人妻一区| 中国美女看黄片| 国产免费男女视频| 精品免费久久久久久久清纯| 白带黄色成豆腐渣| 在线观看美女被高潮喷水网站| 亚洲最大成人中文| 色哟哟·www| 又爽又黄a免费视频| 午夜影院日韩av| 欧美成人一区二区免费高清观看| 麻豆一二三区av精品| 国产精品av视频在线免费观看| 色视频www国产| 国语自产精品视频在线第100页| 国产精品一区二区三区四区久久| 成人鲁丝片一二三区免费| 国产免费一级a男人的天堂| 欧美激情国产日韩精品一区| 亚洲精品成人久久久久久| 午夜精品在线福利| 搡老熟女国产l中国老女人| 日本五十路高清| 精品一区二区三区人妻视频| 久久午夜亚洲精品久久| 亚洲成人av在线免费| 麻豆国产av国片精品| 久久热精品热| 干丝袜人妻中文字幕| 在线观看一区二区三区| 一本精品99久久精品77| 亚洲精品色激情综合| 99久久精品一区二区三区| 最近视频中文字幕2019在线8| 亚洲最大成人中文| 亚洲va在线va天堂va国产| 老女人水多毛片| 亚洲精品日韩在线中文字幕 | 国产一区亚洲一区在线观看| 最近2019中文字幕mv第一页| 久久亚洲国产成人精品v| 男人舔女人下体高潮全视频| 免费观看的影片在线观看| 看免费成人av毛片| 听说在线观看完整版免费高清| 老师上课跳d突然被开到最大视频| 亚洲精品色激情综合| 男人狂女人下面高潮的视频| 中文资源天堂在线| 免费不卡的大黄色大毛片视频在线观看 | 观看免费一级毛片| 国产不卡一卡二| 亚洲精品456在线播放app| 日本-黄色视频高清免费观看| 天天躁日日操中文字幕| 啦啦啦韩国在线观看视频| 亚洲av.av天堂| 一级av片app|