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

    采用粘彈性人工邊界單元時顯式算法穩(wěn)定性分析

    2020-11-14 06:49:58李述濤劉晶波陶西貴陳一村賈藝凡
    工程力學 2020年11期
    關(guān)鍵詞:粘彈性角點子系統(tǒng)

    李述濤,劉晶波,寶 鑫,陶西貴,陳一村,肖 蘭,賈藝凡

    (1. 清華大學土木工程系,北京 100084;2. 軍事科學院國防工程研究院,北京 100036)

    采用數(shù)值方法計算土-結(jié)構(gòu)地震反應(yīng)或近場波動問題時,需要從半無限介質(zhì)中切取出有限的近場區(qū)域進行計算,同時需在截斷邊界處設(shè)置人工邊界來模擬外行散射波向無窮遠的輻射效應(yīng)[1 ? 4]。較為常見的粘彈性人工邊界[5]計算精度較高但前處理操作較為復雜,需遍歷所有邊界節(jié)點施加彈簧-阻尼器系統(tǒng),在此基礎(chǔ)上發(fā)展了粘彈性人工邊界單元技術(shù)[6 ? 7],它是沿邊界法線向外延伸的一層單元,通過賦予該層單元等效物理參數(shù)即可模擬粘彈性人工邊界對于外行散射波的吸收作用。由于該方法具有良好的計算精度和魯棒性,且前處理過程簡單,在實際工程計算中得到較多的應(yīng)用[8 ? 12]。

    隨著人工邊界技術(shù)的不斷發(fā)展,針對人工邊界條件的穩(wěn)定性研究逐漸得到重視[13 ? 14],目前針對此類問題的主要研究方法包括:Trefethen 等[15]提出的GKS 定理—基于偏微分方程初邊值問題有限差分格式的穩(wěn)定性理論,給出初邊值問題多層線性差分格式穩(wěn)定的充要條件;Liao 等[16]研究了離散模型中穩(wěn)態(tài)波動的完備解,給出了多次透射邊界的穩(wěn)定性條件;Kamel 等[17]和關(guān)慧敏等[18]通過傳遞矩陣譜半徑分析積分格式的穩(wěn)定性。理論研究表明,如果不全面考慮逐步積分過程中人工邊界節(jié)點與內(nèi)部節(jié)點運動方程的耦合效應(yīng),穩(wěn)定性準則可能會失效。

    由于粘彈性人工邊界或粘彈性人工邊界單元本身均為物理上穩(wěn)定的系統(tǒng),將其引入計算系統(tǒng)中,并不影響整體的物理穩(wěn)定條件。但采用顯式時域逐步積分方法進行計算時,受粘彈性人工邊界單元粘性阻尼、剛度及幾何尺寸的影響,整體計算模型的數(shù)值積分穩(wěn)定性條件將發(fā)生改變,目前尚未給出明確而實用的含粘彈性人工邊界條件影響的顯式算法穩(wěn)定性準則,影響了數(shù)值分析時穩(wěn)定時間步長的判斷和選取,進一步限制了粘彈性人工邊界單元在顯式動力分析中的應(yīng)用。因此目前在引入粘彈性人工邊界單元進行波動問題分析時,多采用隱式的無條件穩(wěn)定的時域逐步積分算法,以避免顯式算法帶來的數(shù)值穩(wěn)定性問題。隱式算法需要求解聯(lián)立方程組,計算效率不高,并不適合大規(guī)模波動問題計算。隨著顯式時域逐步積分算法在工程結(jié)構(gòu)和大范圍場地分析問題中的廣泛應(yīng)用[19 ? 22],有必要對含有粘彈性人工邊界單元的系統(tǒng)進行顯式時域逐步積分算法穩(wěn)定性研究。

    本文考慮人工邊界和內(nèi)部單元節(jié)點的運動耦合效應(yīng),利用傳遞矩陣的譜半徑分析時域逐步積分格式(中心差分)的穩(wěn)定性,給出不同位置處局部節(jié)點系統(tǒng)顯式時域逐步積分算法穩(wěn)定條件的解析解,通過各子系統(tǒng)穩(wěn)定性條件及內(nèi)部介質(zhì)穩(wěn)定性條件的對比分析,給出考慮粘彈性人工邊界條件影響的整體耦合系統(tǒng)顯式時域逐步積分算法的穩(wěn)定性條件。

    1 粘彈性人工邊界單元等效物理參數(shù)及尺寸設(shè)置

    粘彈性人工邊界單元是在模型截斷邊界處沿法線向外延伸的一層單元(見圖1),人工邊界單元的厚度為h,質(zhì)量密度為0,單元最外層節(jié)點固定。通過賦予單元等效物理參數(shù)來模擬粘彈性人工邊界。二維粘彈性人工邊界單元等效彈性模量、等效泊松比和等效阻尼系數(shù)分別為[2]:

    圖1 人工邊界單元尺寸示意圖Fig. 1 Artificial boundary elements size diagram

    其中:設(shè)G為內(nèi)部介質(zhì)剪切模量; ρ為介質(zhì)質(zhì)量密度;CS和CP分別為介質(zhì)S 波和P 波波速;R為散射波源至邊界節(jié)點的距離; α=αN/αT, αT與 αN分別為切向與法向粘彈性人工邊界參數(shù),對于二維粘彈性人工邊界,其推薦值分別為0.5 和1,根據(jù)式(1),此時為0。

    盡管粘彈性人工邊界單元等效剛度矩陣的推導過程引入了邊界單元厚度h遠小于寬度L的假設(shè),如圖1 所示,劉晶波等[6]和谷音等[7]通過進一步的研究表明,邊界單元厚度的魯棒性良好,可靈活取值而對精度影響很小。由于顯式算法的穩(wěn)定條件受模型中的最小單元尺寸制約,在實際計算分析中,建議將邊界單元的寬度設(shè)定為與內(nèi)部單元尺寸一致,以排除邊界單元尺寸對整體穩(wěn)定性的影響,如圖2 所示。此時粘彈性人工邊界單元等效彈性模量、等效剪切模量可以寫為:

    圖2 適用于顯式算法的人工邊界單元Fig. 2 Artificial boundary elements adapted to explicit algorithms

    使用通用有限元軟件對粘彈性人工邊界單元賦予材料參數(shù)時,由于是各向同性介質(zhì),輸入等效彈性模量即可。另外還包括式(1)計算得到的等效阻尼系數(shù),以及泊松比(=0)和質(zhì)量密度(=0)。由于大多數(shù)軟件不支持將密度設(shè)為0,可采用近似0 的小數(shù)代替。

    2 顯式時域逐步積分算法穩(wěn)定性分析方法

    時域逐步積分算法穩(wěn)定性分析的目的是獲得時域逐步積分計算時滿足穩(wěn)定性要求的時間離散步長 ?t。最大的穩(wěn)定時間步長 ?t與單元的尺寸和系統(tǒng)的物理性質(zhì)有關(guān)。為此在進行穩(wěn)定性分析時,通常采用均勻介質(zhì)和均勻離散化網(wǎng)格模型,如圖3 所示。

    圖3 均勻離散化網(wǎng)格模型Fig. 3 Uniform discretization mesh model

    當實際力學模型介質(zhì)非均勻,單元尺寸不相同時,可以根據(jù)最小尺寸的單元或綜合考慮單元尺寸及其介質(zhì)性質(zhì)來確定整體計算模型的穩(wěn)定離散時間步長 ?t。算法的穩(wěn)定性通??刹捎靡韵聨追N分析方法:

    ① 不考慮自由邊界和人工邊界的影響,取無限計算區(qū)域模型,采用馮諾依曼方法進行穩(wěn)定性分析。這一方法可以獲得內(nèi)部系統(tǒng)的數(shù)值計算穩(wěn)定條件,但無法考慮含人工邊界條件時系統(tǒng)的穩(wěn)定性。

    ② 取實際的計算模型,考慮人工邊界的影響,通過對整體系統(tǒng)時域逐步積分方法的傳遞矩陣進行特征值分析,以獲得考慮耦合人工邊界條件影響的整體系統(tǒng)的數(shù)值計算穩(wěn)定性條件。但由于涉及到整體模型傳遞矩陣的建立和分析,這一方法在實際工作中通常是不可行的,特別是對大規(guī)模的計算問題。

    ③ 為對透射人工邊界的數(shù)值計算穩(wěn)定性進行有效分析,Kamel 等[17]和關(guān)慧敏等[18]提出了一種對含透射人工邊界條件子系統(tǒng)的傳遞矩陣進行特征值分析,以獲得波動問題中透射人工邊界穩(wěn)定性的分析方法。由于分析中采用的子系統(tǒng)為沿人工邊界切向(寬度方向)上若干排節(jié)點和沿法向(深度方向)上若干排節(jié)點組成的節(jié)點系,子系統(tǒng)自由度較多,僅能靠數(shù)值方法對人工邊界的穩(wěn)定性進行分析和判斷,未給出具有解析形式的更易于使用的人工邊界穩(wěn)定性準則。

    從以上有限元離散模型數(shù)值積分穩(wěn)定性準則的研究工作中可以發(fā)現(xiàn)以下兩個特點:1)算法的穩(wěn)定性與模型的截止頻率(系統(tǒng)最高階自振頻率)有關(guān)。2)截止頻率相應(yīng)的振型呈現(xiàn)局部節(jié)點系相鄰節(jié)點交錯振動的形態(tài)。由以上兩個特點可以判斷,整體系統(tǒng)的截止頻率可通過對局部子系統(tǒng)模型的分析得到,進一步可通過局部子系統(tǒng)模型的數(shù)值穩(wěn)定性分析,獲得整體系統(tǒng)的數(shù)值穩(wěn)定性條件。

    首先以一維剪切梁問題為例證明以上判斷。圖4 為一維剪切梁的離散模型及最高階振型示意圖,梁的剪切模量為G,密度為 ρ,橫截面面積為A,單元長度為L。取兩個子系統(tǒng)進行分析。子系統(tǒng)a:取振型兩節(jié)點之間的子系統(tǒng),兩端施加約束,由于約束是施加于振型節(jié)點(振幅為零)之上,因此并不改變剪切梁有限元模型的自振頻率,如圖5(a)所示;子系統(tǒng)b:取兩個單元,在兩端施加約束,此時子系統(tǒng)的自振頻率與原系統(tǒng)的截止頻率不相同,如圖5(b)所示。

    圖4 一維剪切梁振動示意圖Fig. 4 One-dimensional shear beam vibration diagram

    圖5 一維局部子系統(tǒng)Fig. 5 One-dimensional local subsystem

    根據(jù)子系統(tǒng)的剛度和質(zhì)量,可以建立子系統(tǒng)a的單自由度運動方程,當采用中心差分法進行時域逐步積分計算時,一維子系統(tǒng)a 的數(shù)值穩(wěn)定條件為:

    其中,Tn為子系統(tǒng)的自振周期,將式(3)中第三式代入式(4),得到子系統(tǒng)數(shù)值積分穩(wěn)定性條件為:

    式(5)即為一維連續(xù)介質(zhì)波動有限元分析時中心差分算法的穩(wěn)定性條件,可見采用子系統(tǒng)a 獲得的局部穩(wěn)定性條件即為連續(xù)介質(zhì)整體模型的穩(wěn)定性條件。

    一維子系統(tǒng)b 的剛度kb、質(zhì)量mb、自振頻率ωb和中心差分穩(wěn)定性條件分別為:

    對比式(5)和式(6),可以發(fā)現(xiàn)一維子系統(tǒng)b的穩(wěn)定性條件比實際整體模型的穩(wěn)定條件更為寬松,但仍然可以給出穩(wěn)定時間積分步長 ?t的一個上限估計。

    圖6 二維節(jié)點系統(tǒng)最高階振型圖Fig. 6 Highest-order modal pattern of two-dimensional nodes system

    同樣取兩個子系統(tǒng)進行分析,二維問題局部子系統(tǒng)由4 個單元構(gòu)成,如圖7(a)和圖7(b)所示。其中二維子系統(tǒng)a 的4 個角點為振型節(jié)點,位移為0,4 個邊點的位移條件可由振型和4 邊形等參元的性質(zhì)確定,子系統(tǒng)a 對應(yīng)的邊界條件為:

    圖7 二維局部子系統(tǒng)Fig. 7 Two-dimensional local subsystem

    二維子系統(tǒng)a 自振頻率和中心差分穩(wěn)定條件為:

    式(9)即為二維波動問題有限元分析時中心差分算法的穩(wěn)定性條件,可見同樣可以由子系統(tǒng)a的局部穩(wěn)定性分析獲得二維有限元模型數(shù)值算法的整體穩(wěn)定性條件。

    二維子系統(tǒng)b 由4 個單元構(gòu)成,在周邊節(jié)點施加固定約束,子系統(tǒng)b 運動方程為:

    同樣,二維子系統(tǒng)b 給出了整體有限元模型穩(wěn)定時間積分步長的上限估計。

    一維和二維子系統(tǒng)分析的結(jié)果表明:1)若能較為準確地判斷系統(tǒng)截止頻率對應(yīng)的振動模態(tài)(振型),則可以利用最高振型的特點和振型節(jié)點的分布規(guī)律,從整體系統(tǒng)中分離出由陣型節(jié)點所包圍的局部子系統(tǒng),對振型節(jié)點施加約束條件,然后對該子系統(tǒng)進行穩(wěn)定性分析,獲得的局部子系統(tǒng)穩(wěn)定性條件即為整體模型的數(shù)值穩(wěn)定條件;2)若不能準確判斷截止頻率對應(yīng)振動模態(tài)和相應(yīng)振型節(jié)點的位置,也可選取一個能體現(xiàn)整體有限元模型特征的最小的子系統(tǒng),對該子系統(tǒng)的邊界節(jié)點施加約束,通過分析獲得子系統(tǒng)的數(shù)值穩(wěn)定性條件,從而獲得整體系統(tǒng)穩(wěn)定性判據(jù)中各物理量的關(guān)系式,且該條件是整體系統(tǒng)穩(wěn)定性條件的一個逼近。再通過擴展數(shù)值算例進行修正,確定合適的系數(shù),即可得到整體系統(tǒng)的數(shù)值積分穩(wěn)定性條件。

    當采用粘彈性人工邊界單元時,內(nèi)部系統(tǒng)的數(shù)值穩(wěn)定條件已知,但難以確定人工邊界區(qū)截止頻率對應(yīng)的振型,因而僅能采用類似于b 型的子系統(tǒng)進行數(shù)值計算的穩(wěn)定性分析。

    對非均質(zhì)子系統(tǒng)進行穩(wěn)定性分析時,可根據(jù)顯式逐步積分算法格式,將系統(tǒng)運動方程寫成如下形式:

    積分格式的穩(wěn)定性問題與外力向量Qp無關(guān)。如果滿足下列兩條件,則積分格式(12)是穩(wěn)定[23]:條件1: ρ(A)≤1 , ρ(A)是傳遞矩陣A的譜半徑,即ρ(A)=max|λi|;條件2:如果A具有多重特征值,則該特征值的模小于1。因而可將子系統(tǒng)的穩(wěn)定性分析歸結(jié)為傳遞矩陣A的形成及譜半徑計算,按上述條件即可建立穩(wěn)定性準則。

    3 人工邊界子系統(tǒng)穩(wěn)定性分析

    采用粘彈性人工邊界單元時,人工邊界區(qū)的子系統(tǒng)有兩種形式,一種在側(cè)面或底面切取,如圖8 所示,另一種在角點處切取,切取邊界可設(shè)置固定邊界(見圖9)或自由邊界(見圖10)。

    圖8 側(cè)邊固定邊界子系統(tǒng)Fig. 8 Side fixed boundary subsystem

    圖9 角點固定邊界子系統(tǒng)Fig. 9 Corner fixed boundary subsystem

    圖10 角點自由邊界子系統(tǒng)Fig. 10 Corner fixed boundary subsystem

    3.1 側(cè)邊固定邊界子系統(tǒng)穩(wěn)定性分析

    側(cè)邊固定邊界子系統(tǒng)由兩個粘彈性人工邊界單元和兩個內(nèi)部介質(zhì)單元組成,如圖8 所示。設(shè)內(nèi)部介質(zhì)單元的彈性模量為E,剪切模量為G,密度為 ρ,泊松比為μ,且無阻尼;粘彈性人工邊界單元彈性模量為,阻尼系數(shù)為,泊松比為0,密度為0,單元邊長均為L。四節(jié)點正方形平面單元的集中質(zhì)量矩陣、剛度矩陣、阻尼矩陣分別見文獻[6]和文獻[24],按節(jié)點編號進行矩陣組裝后得到子系統(tǒng)中5 號節(jié)點運動方程如下:

    觀察式(13)可知,x和y方向運動方程是解耦的且形式相同,因此可只對其中一個方向的運動方程進行穩(wěn)定性分析。顯式時域逐步積分格式(中心差分[25]如下:

    式(19)即為側(cè)邊局部人工邊界子系統(tǒng)穩(wěn)定性條件的解析表達形式,觀察可知該子系統(tǒng)的穩(wěn)定性條件與內(nèi)部介質(zhì)壓縮波速、泊松比、單元尺寸和模型大小有關(guān)。

    3.2 角點固定邊界子系統(tǒng)穩(wěn)定性分析

    角點處固定邊界子系統(tǒng)中由3 個粘彈性人工邊界單元和1 個內(nèi)部介質(zhì)單元組成,如圖9 所示。將四節(jié)點正方形平面單元質(zhì)量、阻尼和剛度矩陣按節(jié)點編號組裝后得到5 號節(jié)點的運動方程如下:

    由于該運動方程的坐標耦聯(lián),需對整體方程進行穩(wěn)定性分析。將式(1)代入式(20),然后根據(jù)式(14)將式(20)展開,寫成如下傳遞矩陣的形式:

    其中:

    式(24)即為角點處局部系統(tǒng)穩(wěn)定性條件的解析表達形式,觀察可知影響角點局部人工邊界子系統(tǒng)穩(wěn)定性的參數(shù)與側(cè)邊相同,但各自貢獻不同。

    3.3 角點自由邊界子系統(tǒng)穩(wěn)定性分析

    為進一步研究角點處自由邊界子系統(tǒng)的穩(wěn)定性,如圖10 所示。研究對象為編號為1、2、4、5 的四個節(jié)點。子系統(tǒng)中既有單元內(nèi)部節(jié)點,也包含邊界處節(jié)點。將四節(jié)點正方形平面單元質(zhì)量、阻尼和剛度矩陣按節(jié)點編號組裝后,可得到子系統(tǒng)的剛度、質(zhì)量和阻尼矩陣,均為8×8 階對稱矩陣。將整體運動方程按照式(14)中心差分格式展開,得到16×16 階傳遞矩陣。該傳遞矩陣無法得出解析形式的穩(wěn)定性條件,可代入?yún)?shù)計算數(shù)值解。考慮穩(wěn)定性條件受系統(tǒng)截止頻率影響,自由邊界子系統(tǒng)截止頻率要低于固定邊界子系統(tǒng),由此判斷后者穩(wěn)定性條件應(yīng)比前者嚴格,下一節(jié)將進行驗證。

    4 穩(wěn)定性條件比較

    為比較3 種人工邊界子系統(tǒng)的數(shù)值計算穩(wěn)定條件,采用兩組參數(shù)進行分析,具體參數(shù)見表1。

    表1 兩組參數(shù)值(無量綱)Table 1 Two sets of parameter values (dimensionless)

    側(cè)邊固定和角點固定子系統(tǒng)的穩(wěn)定條件可由式(19)和式(24)獲得,角點自由邊界子系統(tǒng)的穩(wěn)定條件可采用數(shù)值方法獲得。3 種子系統(tǒng)譜半徑的計算結(jié)果見圖11。當譜半徑小于等于1 時,數(shù)值計算滿足穩(wěn)定性條件,由圖11 可以直觀地比較兩組參數(shù)條件下3 種子系統(tǒng)的穩(wěn)定條件。

    表2 中對3 種局部子系統(tǒng)和內(nèi)部系統(tǒng)滿足穩(wěn)定性條件的臨界時間步長進行了定量比較。內(nèi)部系統(tǒng)的穩(wěn)定性條件( ?t=L/CP)未考慮粘彈性人工邊界單元對數(shù)值算法的影響,最為寬松;側(cè)邊和角點子系統(tǒng)考慮了粘彈性人工邊界單元質(zhì)量、等效剛度和阻尼對穩(wěn)定性的影響,穩(wěn)定性條件要比無人工邊界單元時更為嚴格。四周固定的角點處邊界節(jié)點只共享1/4 內(nèi)部單元質(zhì)量,其子系統(tǒng)截止頻率最高,因此穩(wěn)定性條件最為嚴格。以上3 種子系統(tǒng)數(shù)值積分穩(wěn)定條件的對比表明,采用粘彈性人工邊界單元時,系統(tǒng)數(shù)值積分的穩(wěn)定條件由角點區(qū)控制。

    圖11 三種子系統(tǒng)的譜半徑對比Fig. 11 Comparison of spectral radius of three subsystems

    表2 穩(wěn)定性條件比較(無量綱)Table 2 Comparison of stability conditions (dimensionless)

    5 算例驗證

    5.1 均勻半空間模型

    為驗證以上穩(wěn)定性分析的準確性,按第一組參數(shù)建立有限元近場模型,模型尺寸為320×160,內(nèi)部介質(zhì)密度為2500,剪切波速為500,泊松比為0.3,網(wǎng)格尺寸為2×2,模型側(cè)邊和底邊最外層單元是粘彈性人工邊界單元,如圖12 所示。采用持時為0.2 的脈沖作為動力荷載,施加于模型中點處,時程曲線如圖13 所示。分別按照內(nèi)部介質(zhì)數(shù)值穩(wěn)定條件(?t=0.0021)、側(cè)邊子系統(tǒng)穩(wěn)定條件(?t=0.00163)、固定邊界角點子系統(tǒng)穩(wěn)定條件(?t=0.00623),采用固定時間步長的顯示時域逐步積分算法對整體模型進行計算。

    圖12 均勻半空間算例模型圖Fig. 12 Homogeneous half apace example model diagram

    圖13 動力荷載時程曲線Fig. 13 Dynamic load time history curve

    按照內(nèi)部介質(zhì)數(shù)值穩(wěn)定條件(?t=0.0021)計算時,由于不滿足側(cè)邊(底邊)子系統(tǒng)穩(wěn)定條件,因此P 波傳播到距離波源最近的底邊節(jié)點時發(fā)生失穩(wěn),如圖14 所示;按照側(cè)邊子系統(tǒng)穩(wěn)定條件(?t=0.00163)計算時,由于不滿足角點處子系統(tǒng)穩(wěn)定條件,波動傳播到模型角點處時發(fā)生失穩(wěn),如圖15所示,其中U為介質(zhì)中的位移,以上兩種穩(wěn)定條件均無法完成整體有限元模型的計算。

    圖14 按內(nèi)部穩(wěn)定條件計算時底邊失穩(wěn)狀態(tài)(0.11 時刻)Fig. 14 The unstable state of the bottom edge when calculated according to internal stability conditions (time 0.11)

    按照角點處子系統(tǒng)穩(wěn)定條件(?t=0.000623)計算時,可順利完成整體模型的動力顯式計算,粘彈性人工邊界單元很好模擬了外行波向無窮遠的輻射,結(jié)果如圖16 所示。此外,通過進一步的計算分析發(fā)現(xiàn),當整體穩(wěn)定條件取值略大于角點子系統(tǒng)穩(wěn)定條件時(例如?t=0.00085),模型角點處也發(fā)生失穩(wěn),失穩(wěn)狀態(tài)與圖15 相同。

    圖15 按側(cè)邊子系統(tǒng)穩(wěn)定條件計算時角點失穩(wěn)狀態(tài)(0.19 時刻)Fig. 15 The unstable state of the bottom edge when calculated according to the stability conditions of side subsystem (time 0.19)

    圖16 按角點子系統(tǒng)穩(wěn)定條件計算結(jié)果Fig. 16 Results calculated by the stability conditions of the corner subsystem

    以上計算分析驗證了采用粘彈性人工邊界單元時,整體模型顯式數(shù)值積分算法的穩(wěn)定性由角點區(qū)域控制,式(24)給出的角點處子系統(tǒng)穩(wěn)定條件是整體模型數(shù)值穩(wěn)定的充分條件。

    5.2 成層半空間模型

    為滿足實際場地計算需要,對成層半空間算例進行驗證。成層半空間有限元模型尺寸為320×160,上半部分內(nèi)部介質(zhì)密度2000,剪切波速300,泊松比0.27,下半部分內(nèi)部介質(zhì)密度為2500,剪切波速為500,泊松比為0.3,整體模型網(wǎng)格尺寸為2×2,模型側(cè)邊和底邊最外層單元是粘彈性人工邊界單元,如圖17 所示。脈沖荷載施加于模型中心,如圖13 所示。

    表3 給出了成層半空間局部子系統(tǒng)和內(nèi)部系統(tǒng)滿足穩(wěn)定性條件的臨界時間步長。模型中的上層介質(zhì)只有側(cè)邊子系統(tǒng),無角點子系統(tǒng),下層介質(zhì)兩者均存在。

    圖17 成層半空間算例模型圖Fig. 17 Layered half space example model diagram

    表3 成層半空間穩(wěn)定性條件比較(無量綱)Table 3 Comparison of stability conditions in layered half space (dimensionless)

    經(jīng)比較發(fā)現(xiàn),上層介質(zhì)中滿足內(nèi)部系統(tǒng)穩(wěn)定條件的時間步長(?t=0.0037)大于側(cè)邊子系統(tǒng)穩(wěn)定時間步長(?t=0.0028)。而二者均比下層介質(zhì)內(nèi)部系統(tǒng)穩(wěn)定的時間步長 (?t=0.0021)要寬松。因此,當采用上層介質(zhì)穩(wěn)定性條件對整體模型計算時,首先發(fā)生失穩(wěn)的是下層介質(zhì)的內(nèi)部系統(tǒng),如圖18所示。

    圖18 按上層介質(zhì)穩(wěn)定條件計算時下層介質(zhì)內(nèi)部系統(tǒng)失穩(wěn)狀態(tài)(0.1 時刻)Fig. 18 The unstable state of the internal system in lower media when calculated according to the stability conditions in upper media (time 0.1)

    采用下層介質(zhì)內(nèi)部系統(tǒng)穩(wěn)定條件(?t=0.0021)計算時,整體系統(tǒng)失穩(wěn)狀態(tài)與圖14 相同;采用下層介質(zhì)側(cè)邊子系統(tǒng)穩(wěn)定條件計算時(?t=0.00163),失穩(wěn)狀態(tài)和圖15 相同。

    按照下層角點子系統(tǒng)穩(wěn)定條件(?t=0.000623)計算時,可順利完成整體模型的動力顯式計算,分層介質(zhì)中,粘彈性人工邊界單元也可很好地模擬外行波向無窮遠的輻射,結(jié)果如圖19 所示。成層半空間算例進一步驗證了采用粘彈性人工邊界單元時,整體模型顯式數(shù)值積分算法的穩(wěn)定性仍然由角點區(qū)域控制。

    圖19 按下層介質(zhì)角點子系統(tǒng)穩(wěn)定條件計算結(jié)果Fig. 19 Results calculated by the stability conditions of the corner subsystem in lower media

    6 結(jié)論與展望

    本文將整體模型數(shù)值穩(wěn)定性問題合理轉(zhuǎn)移到若干局部子系統(tǒng)中,充分考慮粘彈性人工邊界單元和內(nèi)部單元節(jié)點間的運動耦合效應(yīng),通過傳遞矩陣譜半徑分析方法推導出各局部子系統(tǒng)顯式時域逐步積分(中心差分)數(shù)值穩(wěn)定條件的解析解和數(shù)值解。通過計算軟件驗證理論分析的可靠性。具體結(jié)論如下:

    (1) 對于大規(guī)模數(shù)值計算問題,可選取局部的子系統(tǒng)并對其進行顯式時域逐步積分算法的穩(wěn)定性分析,該穩(wěn)定性條件與整體系統(tǒng)穩(wěn)定性條件相同或近似。

    (2) 采用粘彈性人工邊界單元時,受人工邊界單元質(zhì)量、剛度和阻尼影響,整體系統(tǒng)的穩(wěn)定性條件與內(nèi)部介質(zhì)的穩(wěn)定性條件不同,前者的穩(wěn)定性條件更為嚴格,需使用更小的積分步長以滿足整體系統(tǒng)的數(shù)值穩(wěn)定。

    (3) 采用粘彈性人工邊界單元時,整體模型顯式數(shù)值積分算法的穩(wěn)定性由角點區(qū)域控制,本文給出了角點處子系統(tǒng)數(shù)值積分的穩(wěn)定性條件,該穩(wěn)定性條件是整體模型數(shù)值積分穩(wěn)定的充分條件。此外,本文給出的穩(wěn)定性條件是以正方形平面單元為對象推導的,同樣適用于矩形平面單元,由于系統(tǒng)穩(wěn)定條件受最小單元尺寸影響,可使用矩形單元的最小邊長作為參數(shù)計算穩(wěn)定性條件。

    下一步展望:

    (1) 對于三維粘彈性人工邊界單元,同樣可以利用本文提出的傳遞矩陣譜半徑分析方法對顯式時域逐步積分算法的穩(wěn)定性條件進行分析,三維問題涉及的局部子結(jié)構(gòu)更為特殊,傳遞矩陣的生成和特征值求解更加復雜,需進一步開展研究。

    (2) 相比隱式算法,顯式算法的解耦特性對于求解大范圍復雜工程場地問題更有優(yōu)勢。本文的研究成果為在顯式算法中合理使用粘彈性人工邊界單元提供了理論依據(jù)??稍诖嘶A(chǔ)上進一步開展分析和研究,以改善使用粘彈性人工邊界單元時顯式算法的穩(wěn)定性,提高大范圍工程場地問題的計算效率。

    猜你喜歡
    粘彈性角點子系統(tǒng)
    不對中轉(zhuǎn)子系統(tǒng)耦合動力學特性研究
    二維粘彈性棒和板問題ADI有限差分法
    GSM-R基站子系統(tǒng)同步方案研究
    時變時滯粘彈性板方程的整體吸引子
    駝峰測長設(shè)備在線監(jiān)測子系統(tǒng)的設(shè)計與應(yīng)用
    基于FAST角點檢測算法上對Y型與X型角點的檢測
    不可壓粘彈性流體的Leray-α-Oldroyd模型整體解的存在性
    基于邊緣的角點分類和描述算法
    電子科技(2016年12期)2016-12-26 02:25:49
    基于圓環(huán)模板的改進Harris角點檢測算法
    車載ATP子系統(tǒng)緊急制動限制速度計算
    国国产精品蜜臀av免费| 老司机影院成人| 国产无遮挡羞羞视频在线观看| 伊人久久国产一区二区| 成年人免费黄色播放视频| 日本-黄色视频高清免费观看| 国产一区二区三区av在线| av又黄又爽大尺度在线免费看| 91在线精品国自产拍蜜月| 女性被躁到高潮视频| 成年av动漫网址| 高清视频免费观看一区二区| 亚洲av国产av综合av卡| 九色成人免费人妻av| 精品国产国语对白av| 性色av一级| 日韩电影二区| 99精国产麻豆久久婷婷| 国产欧美日韩综合在线一区二区| 黄色怎么调成土黄色| 最近2019中文字幕mv第一页| 男人添女人高潮全过程视频| 99精国产麻豆久久婷婷| 少妇人妻久久综合中文| 高清不卡的av网站| 韩国高清视频一区二区三区| 人妻 亚洲 视频| 五月开心婷婷网| 只有这里有精品99| 久久ye,这里只有精品| 在现免费观看毛片| 一个人免费看片子| 久久鲁丝午夜福利片| 亚洲精品自拍成人| 国产亚洲午夜精品一区二区久久| av不卡在线播放| 欧美 日韩 精品 国产| 日韩精品免费视频一区二区三区 | 免费黄频网站在线观看国产| 观看美女的网站| 精品国产乱码久久久久久小说| 美女国产视频在线观看| 观看美女的网站| 国产成人aa在线观看| 亚洲国产毛片av蜜桃av| 国产熟女欧美一区二区| 18禁裸乳无遮挡动漫免费视频| 亚洲精品av麻豆狂野| a级毛色黄片| 国产一区有黄有色的免费视频| 免费大片黄手机在线观看| 蜜臀久久99精品久久宅男| 国产男女内射视频| 曰老女人黄片| 国产淫语在线视频| 精品福利永久在线观看| 久久人人爽人人爽人人片va| 一本大道久久a久久精品| 韩国精品一区二区三区 | 国产免费一区二区三区四区乱码| 熟妇人妻不卡中文字幕| 老司机影院毛片| 看非洲黑人一级黄片| 日本爱情动作片www.在线观看| 少妇 在线观看| 一个人免费看片子| 中文字幕另类日韩欧美亚洲嫩草| 一本色道久久久久久精品综合| 在线观看免费高清a一片| 啦啦啦啦在线视频资源| 色5月婷婷丁香| 一级,二级,三级黄色视频| 在线观看三级黄色| 国产 精品1| 国产成人精品一,二区| 亚洲成人av在线免费| 国产欧美日韩综合在线一区二区| 久久国产亚洲av麻豆专区| 婷婷成人精品国产| 黄色配什么色好看| 精品一区二区三区四区五区乱码 | 欧美变态另类bdsm刘玥| 国产免费福利视频在线观看| 免费不卡的大黄色大毛片视频在线观看| 热re99久久国产66热| 晚上一个人看的免费电影| 亚洲成av片中文字幕在线观看 | 久久这里只有精品19| 亚洲,一卡二卡三卡| 蜜桃国产av成人99| av有码第一页| 99热网站在线观看| 国产亚洲精品久久久com| 男人舔女人的私密视频| 国产男女超爽视频在线观看| 精品国产一区二区久久| 一级毛片我不卡| av黄色大香蕉| 亚洲欧美一区二区三区黑人 | 老司机影院成人| 亚洲国产看品久久| 永久网站在线| 久久ye,这里只有精品| 一边摸一边做爽爽视频免费| 国产伦理片在线播放av一区| 国产一区二区在线观看日韩| 午夜91福利影院| 亚洲图色成人| 欧美xxⅹ黑人| 国产男人的电影天堂91| 日韩 亚洲 欧美在线| av视频免费观看在线观看| 国产精品 国内视频| 国产激情久久老熟女| 曰老女人黄片| 久久久欧美国产精品| 午夜日本视频在线| 啦啦啦中文免费视频观看日本| 欧美亚洲日本最大视频资源| 国产欧美亚洲国产| a级片在线免费高清观看视频| 这个男人来自地球电影免费观看 | 色吧在线观看| 性色avwww在线观看| 乱人伦中国视频| 精品人妻熟女毛片av久久网站| 欧美精品一区二区免费开放| 国精品久久久久久国模美| 女人被躁到高潮嗷嗷叫费观| 免费大片黄手机在线观看| 视频区图区小说| 国产免费视频播放在线视频| 极品少妇高潮喷水抽搐| 最近最新中文字幕大全免费视频 | 欧美激情国产日韩精品一区| 欧美日韩视频精品一区| 综合色丁香网| 精品99又大又爽又粗少妇毛片| 看免费成人av毛片| 亚洲伊人久久精品综合| av国产久精品久网站免费入址| 精品人妻一区二区三区麻豆| 国产成人精品久久久久久| 国产精品人妻久久久影院| 亚洲四区av| 亚洲精品久久成人aⅴ小说| 欧美精品一区二区大全| 亚洲情色 制服丝袜| 婷婷成人精品国产| 国产在线免费精品| 七月丁香在线播放| 丝袜脚勾引网站| av视频免费观看在线观看| 这个男人来自地球电影免费观看 | 日韩不卡一区二区三区视频在线| 夜夜骑夜夜射夜夜干| 国产色婷婷99| 大香蕉久久成人网| av国产久精品久网站免费入址| 免费少妇av软件| 伊人亚洲综合成人网| 两个人看的免费小视频| 丝瓜视频免费看黄片| 99精国产麻豆久久婷婷| 亚洲欧洲精品一区二区精品久久久 | 在线看a的网站| 国国产精品蜜臀av免费| 亚洲国产欧美日韩在线播放| 肉色欧美久久久久久久蜜桃| 国产日韩一区二区三区精品不卡| 国产一级毛片在线| 日本免费在线观看一区| 狂野欧美激情性xxxx在线观看| 久久久久久久大尺度免费视频| 久久精品人人爽人人爽视色| 欧美精品国产亚洲| 久久精品久久久久久久性| 久久精品国产亚洲av涩爱| 国产麻豆69| 高清欧美精品videossex| 乱人伦中国视频| 精品人妻熟女毛片av久久网站| 国产永久视频网站| 国产精品一国产av| 日韩av免费高清视频| av电影中文网址| 久久久国产欧美日韩av| 这个男人来自地球电影免费观看 | 少妇的丰满在线观看| 丝袜脚勾引网站| 尾随美女入室| 亚洲欧美成人综合另类久久久| 人妻人人澡人人爽人人| kizo精华| 色婷婷av一区二区三区视频| 日韩 亚洲 欧美在线| 亚洲,欧美,日韩| 观看av在线不卡| 大片免费播放器 马上看| 丝袜喷水一区| 亚洲av国产av综合av卡| 啦啦啦视频在线资源免费观看| 极品少妇高潮喷水抽搐| 汤姆久久久久久久影院中文字幕| 9191精品国产免费久久| 午夜免费观看性视频| 国产精品无大码| 好男人视频免费观看在线| 又黄又爽又刺激的免费视频.| 亚洲精品aⅴ在线观看| 有码 亚洲区| 女性被躁到高潮视频| 久久久国产欧美日韩av| 90打野战视频偷拍视频| av福利片在线| 2018国产大陆天天弄谢| 国产成人免费无遮挡视频| 国产亚洲一区二区精品| 久久这里有精品视频免费| 国产精品久久久久久久电影| 91aial.com中文字幕在线观看| 国产精品免费大片| 春色校园在线视频观看| 精品99又大又爽又粗少妇毛片| 亚洲精华国产精华液的使用体验| 满18在线观看网站| 精品一区二区免费观看| 久久97久久精品| 丁香六月天网| 一区在线观看完整版| 好男人视频免费观看在线| 美女中出高潮动态图| 日本av免费视频播放| 18禁在线无遮挡免费观看视频| 日韩欧美一区视频在线观看| 99香蕉大伊视频| 高清欧美精品videossex| 免费看不卡的av| 国产极品粉嫩免费观看在线| 亚洲欧美一区二区三区国产| 精品国产一区二区三区久久久樱花| 2018国产大陆天天弄谢| 亚洲图色成人| 国产精品成人在线| 2021少妇久久久久久久久久久| 热99久久久久精品小说推荐| 成人国产麻豆网| 国产老妇伦熟女老妇高清| 人妻一区二区av| 国产麻豆69| 最近2019中文字幕mv第一页| 美女脱内裤让男人舔精品视频| 有码 亚洲区| 午夜精品国产一区二区电影| 永久网站在线| 大片免费播放器 马上看| 免费看不卡的av| 精品久久久久久电影网| 99热这里只有是精品在线观看| 国产免费现黄频在线看| 亚洲内射少妇av| 国产一级毛片在线| 日韩人妻精品一区2区三区| 亚洲 欧美一区二区三区| 国产精品麻豆人妻色哟哟久久| 一区二区三区精品91| 婷婷色综合www| 99re6热这里在线精品视频| 亚洲美女视频黄频| 欧美3d第一页| 中文天堂在线官网| 人妻人人澡人人爽人人| 亚洲,欧美精品.| 色婷婷av一区二区三区视频| 丰满饥渴人妻一区二区三| 插逼视频在线观看| 国产高清国产精品国产三级| 亚洲精品一区蜜桃| 久久久亚洲精品成人影院| av卡一久久| 少妇的逼水好多| 亚洲精品乱码久久久久久按摩| 人人澡人人妻人| av在线老鸭窝| 亚洲欧美成人综合另类久久久| 久久精品夜色国产| 寂寞人妻少妇视频99o| 青青草视频在线视频观看| 一个人免费看片子| 亚洲美女视频黄频| 中文字幕制服av| 久久免费观看电影| 亚洲综合精品二区| 日韩,欧美,国产一区二区三区| 国产国拍精品亚洲av在线观看| 欧美日韩亚洲高清精品| 母亲3免费完整高清在线观看 | 国产精品久久久久久精品电影小说| www.av在线官网国产| 伦理电影免费视频| 国产色爽女视频免费观看| 波多野结衣一区麻豆| 国产老妇伦熟女老妇高清| 亚洲中文av在线| 亚洲欧洲日产国产| 国产爽快片一区二区三区| 亚洲精品日韩在线中文字幕| 美女国产高潮福利片在线看| 90打野战视频偷拍视频| 日本爱情动作片www.在线观看| 人妻 亚洲 视频| 七月丁香在线播放| 自线自在国产av| 午夜日本视频在线| 人妻少妇偷人精品九色| 国产熟女欧美一区二区| av卡一久久| 色5月婷婷丁香| 亚洲欧洲日产国产| 乱码一卡2卡4卡精品| 91aial.com中文字幕在线观看| 男人爽女人下面视频在线观看| 亚洲经典国产精华液单| 久久青草综合色| 久久毛片免费看一区二区三区| 哪个播放器可以免费观看大片| 亚洲高清免费不卡视频| 内地一区二区视频在线| 中文字幕人妻丝袜制服| 亚洲国产最新在线播放| 国产免费现黄频在线看| 99久久综合免费| 捣出白浆h1v1| 久久久久久久精品精品| 中国美白少妇内射xxxbb| 国产亚洲午夜精品一区二区久久| 久久狼人影院| xxx大片免费视频| 国产日韩欧美亚洲二区| 极品人妻少妇av视频| 18禁动态无遮挡网站| av在线观看视频网站免费| 在线免费观看不下载黄p国产| av.在线天堂| 久久午夜福利片| 岛国毛片在线播放| 91精品伊人久久大香线蕉| 国产探花极品一区二区| 人人澡人人妻人| 免费看不卡的av| 亚洲欧洲日产国产| 久久久精品免费免费高清| 一边亲一边摸免费视频| 日本欧美视频一区| 爱豆传媒免费全集在线观看| 国产色爽女视频免费观看| 亚洲伊人久久精品综合| av在线app专区| 日韩制服骚丝袜av| 韩国精品一区二区三区 | 国产精品麻豆人妻色哟哟久久| 精品一区二区三区四区五区乱码 | 丝袜人妻中文字幕| 51国产日韩欧美| 国产女主播在线喷水免费视频网站| 欧美人与性动交α欧美软件 | 亚洲欧美一区二区三区国产| 久久久亚洲精品成人影院| 国产精品99久久99久久久不卡 | 亚洲,欧美,日韩| 亚洲一区二区三区欧美精品| 在线观看美女被高潮喷水网站| 欧美亚洲日本最大视频资源| 亚洲精品美女久久久久99蜜臀 | 赤兔流量卡办理| 欧美精品国产亚洲| 亚洲丝袜综合中文字幕| 欧美另类一区| 欧美激情 高清一区二区三区| 欧美精品av麻豆av| 欧美老熟妇乱子伦牲交| 伦理电影免费视频| 久久毛片免费看一区二区三区| 男女国产视频网站| 国产伦理片在线播放av一区| 黄色怎么调成土黄色| 国产一区二区三区综合在线观看 | 欧美日韩综合久久久久久| 好男人视频免费观看在线| 成人亚洲欧美一区二区av| 精品少妇内射三级| 肉色欧美久久久久久久蜜桃| 国内精品宾馆在线| 满18在线观看网站| 侵犯人妻中文字幕一二三四区| 男的添女的下面高潮视频| 在现免费观看毛片| 亚洲精品视频女| 侵犯人妻中文字幕一二三四区| 亚洲精品久久久久久婷婷小说| 亚洲人与动物交配视频| 只有这里有精品99| 中文字幕制服av| 国产一区二区在线观看日韩| 日日啪夜夜爽| 亚洲精品乱久久久久久| 天天影视国产精品| 一级a做视频免费观看| 午夜免费观看性视频| 伦理电影免费视频| 久久av网站| 亚洲国产精品一区三区| www.av在线官网国产| 卡戴珊不雅视频在线播放| 久久久久久久久久久免费av| 丰满少妇做爰视频| 一级毛片黄色毛片免费观看视频| 五月天丁香电影| 交换朋友夫妻互换小说| 99热全是精品| 最黄视频免费看| 18禁动态无遮挡网站| 精品国产一区二区久久| av在线老鸭窝| 亚洲欧洲日产国产| 91精品伊人久久大香线蕉| 最近中文字幕2019免费版| 99九九在线精品视频| 最近中文字幕2019免费版| 婷婷色综合www| 最近中文字幕2019免费版| 国产在视频线精品| 国产女主播在线喷水免费视频网站| 最近的中文字幕免费完整| 亚洲精品久久久久久婷婷小说| 亚洲av综合色区一区| 中文字幕人妻熟女乱码| 韩国精品一区二区三区 | 精品一区二区三卡| 伦理电影免费视频| 精品人妻熟女毛片av久久网站| 亚洲,一卡二卡三卡| 久久国产精品大桥未久av| 夜夜骑夜夜射夜夜干| 交换朋友夫妻互换小说| 建设人人有责人人尽责人人享有的| 如日韩欧美国产精品一区二区三区| 国产精品蜜桃在线观看| 99热网站在线观看| 午夜免费观看性视频| 中国国产av一级| 免费看光身美女| 在线观看免费高清a一片| videosex国产| 国产亚洲精品久久久com| 免费人成在线观看视频色| 有码 亚洲区| 国产精品不卡视频一区二区| 少妇猛男粗大的猛烈进出视频| videosex国产| 亚洲av免费高清在线观看| 母亲3免费完整高清在线观看 | 久久精品国产鲁丝片午夜精品| 亚洲一级一片aⅴ在线观看| 午夜福利乱码中文字幕| 欧美日韩视频高清一区二区三区二| 久久精品人人爽人人爽视色| 国内精品宾馆在线| 亚洲精品久久久久久婷婷小说| 18禁动态无遮挡网站| 日本vs欧美在线观看视频| 美女中出高潮动态图| 国产视频首页在线观看| 国产一区二区三区综合在线观看 | 亚洲精品视频女| 超碰97精品在线观看| 亚洲欧洲国产日韩| 夜夜爽夜夜爽视频| 视频中文字幕在线观看| 亚洲成人av在线免费| 下体分泌物呈黄色| 午夜av观看不卡| av线在线观看网站| 国产亚洲欧美精品永久| 日韩在线高清观看一区二区三区| 国产精品久久久久久av不卡| 亚洲av日韩在线播放| 一边摸一边做爽爽视频免费| 亚洲美女视频黄频| 男女无遮挡免费网站观看| 成人18禁高潮啪啪吃奶动态图| 国产探花极品一区二区| 一本久久精品| 超色免费av| 男人舔女人的私密视频| 少妇的逼好多水| 人人妻人人爽人人添夜夜欢视频| 最后的刺客免费高清国语| 男的添女的下面高潮视频| 超碰97精品在线观看| 亚洲av.av天堂| 精品国产一区二区三区四区第35| 观看美女的网站| 最近中文字幕2019免费版| 人妻少妇偷人精品九色| 国产一区二区在线观看日韩| 国产成人精品一,二区| 欧美性感艳星| 最近中文字幕2019免费版| 少妇人妻精品综合一区二区| 最新中文字幕久久久久| 亚洲国产看品久久| 免费观看性生交大片5| 国产有黄有色有爽视频| 99热这里只有是精品在线观看| 少妇人妻 视频| 午夜福利乱码中文字幕| 久久久久久久精品精品| 久久人妻熟女aⅴ| 激情视频va一区二区三区| 人人妻人人澡人人看| 亚洲丝袜综合中文字幕| 免费观看性生交大片5| 国产有黄有色有爽视频| 激情视频va一区二区三区| 亚洲国产精品专区欧美| 欧美变态另类bdsm刘玥| 国产成人a∨麻豆精品| 22中文网久久字幕| 黄色 视频免费看| 一本色道久久久久久精品综合| 亚洲国产精品999| videossex国产| 亚洲精品日韩在线中文字幕| 国产深夜福利视频在线观看| 男女边吃奶边做爰视频| 亚洲经典国产精华液单| 日韩 亚洲 欧美在线| 亚洲av成人精品一二三区| 久久婷婷青草| 黄色毛片三级朝国网站| 成人毛片60女人毛片免费| 婷婷成人精品国产| 久久人人爽人人爽人人片va| 久久久久精品人妻al黑| 男人爽女人下面视频在线观看| 成人黄色视频免费在线看| 国产高清三级在线| 青春草亚洲视频在线观看| 国产免费现黄频在线看| 亚洲欧美日韩另类电影网站| 日韩精品有码人妻一区| 久久久精品免费免费高清| 卡戴珊不雅视频在线播放| 91午夜精品亚洲一区二区三区| 内地一区二区视频在线| 国国产精品蜜臀av免费| 久久久久久久久久人人人人人人| 日韩 亚洲 欧美在线| 成人国产麻豆网| 两个人免费观看高清视频| 精品卡一卡二卡四卡免费| 美女大奶头黄色视频| 欧美日韩亚洲高清精品| 日韩av免费高清视频| videosex国产| 久久久久国产精品人妻一区二区| 国产麻豆69| 18禁国产床啪视频网站| 一本久久精品| av免费观看日本| 美女中出高潮动态图| 久久午夜综合久久蜜桃| 精品少妇黑人巨大在线播放| 日日摸夜夜添夜夜爱| 国产有黄有色有爽视频| 亚洲欧美中文字幕日韩二区| 国产黄色免费在线视频| 欧美3d第一页| 精品国产乱码久久久久久小说| 十八禁高潮呻吟视频| 黑人高潮一二区| 成年动漫av网址| 精品第一国产精品| 亚洲,欧美,日韩| 亚洲欧美日韩卡通动漫| 免费av不卡在线播放| www日本在线高清视频| 人成视频在线观看免费观看| 久久久久精品性色| 18在线观看网站| 国产亚洲一区二区精品| 26uuu在线亚洲综合色| 十八禁高潮呻吟视频| 中文精品一卡2卡3卡4更新| 亚洲丝袜综合中文字幕| 国产亚洲最大av| 亚洲成av片中文字幕在线观看 | 一本—道久久a久久精品蜜桃钙片| 青春草视频在线免费观看| 两个人看的免费小视频| 欧美xxⅹ黑人| 中文字幕亚洲精品专区| 午夜久久久在线观看| 免费av不卡在线播放| 亚洲成人一二三区av| 又粗又硬又长又爽又黄的视频| 大香蕉久久网| 亚洲精品成人av观看孕妇| 亚洲一级一片aⅴ在线观看| 岛国毛片在线播放| 亚洲欧美日韩另类电影网站| 国产精品人妻久久久影院| 精品酒店卫生间| 日韩一区二区三区影片|