• <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麻豆精品| 国产精品国产三级国产专区5o| 干丝袜人妻中文字幕| 欧美日韩亚洲高清精品| 啦啦啦啦在线视频资源| av黄色大香蕉| 日本欧美视频一区| 热re99久久精品国产66热6| 久久精品夜色国产| 免费看av在线观看网站| 亚洲人成网站在线播| 丝袜喷水一区| 夜夜爽夜夜爽视频| 日韩视频在线欧美| 精品一区二区三卡| 少妇的逼水好多| 久热这里只有精品99| 久久亚洲国产成人精品v| 熟女av电影| 在线精品无人区一区二区三 | 国产精品成人在线| 美女主播在线视频| 高清日韩中文字幕在线| 激情五月婷婷亚洲| 色视频在线一区二区三区| 噜噜噜噜噜久久久久久91| 久久国产亚洲av麻豆专区| 最近手机中文字幕大全| 久久国内精品自在自线图片| 少妇 在线观看| 欧美丝袜亚洲另类| 香蕉精品网在线| 日日啪夜夜爽| 国产伦精品一区二区三区四那| 亚洲无线观看免费| 性高湖久久久久久久久免费观看| 欧美日韩一区二区视频在线观看视频在线| av在线蜜桃| 黄片无遮挡物在线观看| 搡女人真爽免费视频火全软件| 日韩大片免费观看网站| 精品久久久久久久久av| 亚洲国产精品一区三区| 国产免费一区二区三区四区乱码| 久久人人爽人人片av| 高清在线视频一区二区三区| 精品视频人人做人人爽| 久久久精品免费免费高清| 国产精品99久久久久久久久| 大陆偷拍与自拍| 亚洲av欧美aⅴ国产| 亚洲精品,欧美精品| 七月丁香在线播放| 国产精品福利在线免费观看| 在线观看一区二区三区| 日本欧美视频一区| 99久久精品一区二区三区| 一个人看视频在线观看www免费| 天堂中文最新版在线下载| 少妇人妻一区二区三区视频| 国产亚洲91精品色在线| 午夜免费鲁丝| 国产人妻一区二区三区在| 午夜精品国产一区二区电影| 国产欧美日韩精品一区二区| 日韩国内少妇激情av| 最近的中文字幕免费完整| 亚洲精品亚洲一区二区| 久久精品国产亚洲av天美| 深爱激情五月婷婷| 免费观看在线日韩| 国产欧美日韩一区二区三区在线 | kizo精华| 91狼人影院| 一个人看的www免费观看视频| 身体一侧抽搐| 搡女人真爽免费视频火全软件| 久久人人爽人人爽人人片va| 国产黄片视频在线免费观看| 人妻少妇偷人精品九色| 亚洲美女搞黄在线观看| 赤兔流量卡办理| 网址你懂的国产日韩在线| 国产黄色免费在线视频| 最近最新中文字幕免费大全7| 偷拍熟女少妇极品色| 久久国内精品自在自线图片| 晚上一个人看的免费电影| 国产欧美日韩精品一区二区| 久久国产精品大桥未久av | 一区二区三区精品91| 亚州av有码| 亚洲,一卡二卡三卡| 国产人妻一区二区三区在| 老熟女久久久| 日韩在线高清观看一区二区三区| 欧美日韩综合久久久久久| av网站免费在线观看视频| 观看av在线不卡| 久久久国产一区二区| 国产综合精华液| 97热精品久久久久久| 婷婷色av中文字幕| 国产乱人偷精品视频| 欧美人与善性xxx| 午夜福利网站1000一区二区三区| 国产精品偷伦视频观看了| 十分钟在线观看高清视频www | 国产精品无大码| 91精品伊人久久大香线蕉| 最近的中文字幕免费完整| 新久久久久国产一级毛片| 国产又色又爽无遮挡免| 国产av精品麻豆| av播播在线观看一区| 麻豆国产97在线/欧美| 啦啦啦视频在线资源免费观看| 色婷婷久久久亚洲欧美| 国产精品伦人一区二区| 精品亚洲乱码少妇综合久久| 蜜桃在线观看..| 黄色日韩在线| 在线观看免费日韩欧美大片 | 美女中出高潮动态图| 久久女婷五月综合色啪小说| 街头女战士在线观看网站| 国产欧美另类精品又又久久亚洲欧美| 91在线精品国自产拍蜜月| 午夜激情福利司机影院| 日本欧美视频一区| av福利片在线观看| 中文欧美无线码| 麻豆国产97在线/欧美| 国语对白做爰xxxⅹ性视频网站| 99热全是精品| 亚洲av成人精品一区久久| h日本视频在线播放| 日本与韩国留学比较| 大片电影免费在线观看免费| 久久6这里有精品| 国产精品99久久99久久久不卡 | 日本免费在线观看一区| 少妇人妻一区二区三区视频| 精品午夜福利在线看| 国产成人免费无遮挡视频| 五月天丁香电影| 国产一级毛片在线| 联通29元200g的流量卡| 网址你懂的国产日韩在线| 街头女战士在线观看网站| 午夜激情久久久久久久| 91精品伊人久久大香线蕉| 国产乱来视频区| av黄色大香蕉| 亚洲,欧美,日韩| 永久免费av网站大全| 亚洲综合色惰| 国国产精品蜜臀av免费| 免费观看无遮挡的男女| 久久精品夜色国产| 秋霞伦理黄片| 我的老师免费观看完整版| 黑人高潮一二区| 婷婷色综合www| 极品少妇高潮喷水抽搐| 精品少妇黑人巨大在线播放| 午夜福利网站1000一区二区三区| 亚洲伊人久久精品综合| 欧美性感艳星| 国产精品一区www在线观看| 久久精品国产鲁丝片午夜精品| 国产欧美亚洲国产| 日本色播在线视频| 91久久精品国产一区二区三区| 国产精品99久久99久久久不卡 | 男的添女的下面高潮视频| 久久久国产一区二区| 久久久色成人| 国产日韩欧美在线精品| 国产熟女欧美一区二区| av免费观看日本| 丰满迷人的少妇在线观看| 狂野欧美白嫩少妇大欣赏| 国产成人a区在线观看| 欧美日韩综合久久久久久| 国产精品av视频在线免费观看| 亚洲成人手机| 99热这里只有是精品在线观看| 日韩一区二区视频免费看| 寂寞人妻少妇视频99o| 亚洲一级一片aⅴ在线观看| av卡一久久| 女的被弄到高潮叫床怎么办| 亚洲精华国产精华液的使用体验| 国产女主播在线喷水免费视频网站| 国产欧美另类精品又又久久亚洲欧美| 嫩草影院入口| 午夜精品国产一区二区电影| 人妻一区二区av| 日本与韩国留学比较| 成人18禁高潮啪啪吃奶动态图 | 黑丝袜美女国产一区| 男女下面进入的视频免费午夜| 久久国产亚洲av麻豆专区| 女的被弄到高潮叫床怎么办| 国产伦理片在线播放av一区| 亚洲av国产av综合av卡| 免费久久久久久久精品成人欧美视频 | 人妻制服诱惑在线中文字幕| 亚洲美女黄色视频免费看| 亚洲精品久久午夜乱码| 国产一级毛片在线| 国产一区二区三区综合在线观看 | 多毛熟女@视频| 熟女人妻精品中文字幕| av视频免费观看在线观看| 插逼视频在线观看| 免费人妻精品一区二区三区视频| 国产精品一区二区在线不卡| 国产精品99久久99久久久不卡 | 日本色播在线视频| 国产高清有码在线观看视频| 亚洲成人一二三区av| 欧美人与善性xxx| 日韩中字成人| 久久6这里有精品| 国产91av在线免费观看| 看十八女毛片水多多多| 日韩不卡一区二区三区视频在线| 国产精品久久久久久av不卡| 97在线人人人人妻| av.在线天堂| 久久久久久伊人网av| 婷婷色综合www| 青春草视频在线免费观看| 三级国产精品欧美在线观看| 国产亚洲91精品色在线| 国产人妻一区二区三区在| 午夜免费观看性视频| freevideosex欧美| 一级二级三级毛片免费看| 日本黄色片子视频| 99re6热这里在线精品视频| 亚洲成人中文字幕在线播放| 综合色丁香网| 丝瓜视频免费看黄片| 看十八女毛片水多多多| 日韩不卡一区二区三区视频在线| 欧美zozozo另类| 高清在线视频一区二区三区| 最近中文字幕2019免费版| 99久久中文字幕三级久久日本| av国产久精品久网站免费入址| 涩涩av久久男人的天堂| 欧美三级亚洲精品| 国产欧美另类精品又又久久亚洲欧美| 老司机影院成人| 午夜免费观看性视频| 免费高清在线观看视频在线观看| 青青草视频在线视频观看| 国产v大片淫在线免费观看| 国产白丝娇喘喷水9色精品| 婷婷色综合大香蕉| 最后的刺客免费高清国语| 成人美女网站在线观看视频| 国产女主播在线喷水免费视频网站| 97在线视频观看| 三级国产精品片| 国产成人freesex在线| 精品国产乱码久久久久久小说| 少妇猛男粗大的猛烈进出视频| 国产久久久一区二区三区| 亚洲成人av在线免费| 国产欧美另类精品又又久久亚洲欧美| 国产片特级美女逼逼视频| a级毛色黄片| 在线精品无人区一区二区三 | 成年av动漫网址| 免费av不卡在线播放| av免费在线看不卡| 麻豆成人午夜福利视频| 人妻 亚洲 视频| 国产精品成人在线| 免费观看无遮挡的男女| 尤物成人国产欧美一区二区三区| 少妇高潮的动态图| av免费观看日本| 91精品伊人久久大香线蕉| av在线老鸭窝| 成人毛片60女人毛片免费| 激情五月婷婷亚洲| 久久久久久久大尺度免费视频| av福利片在线观看| 国产欧美另类精品又又久久亚洲欧美| 老司机影院成人| 蜜桃在线观看..| 观看av在线不卡| 亚洲国产最新在线播放| 精品视频人人做人人爽| av在线观看视频网站免费| av在线播放精品| av在线app专区| 91午夜精品亚洲一区二区三区| 国产综合精华液| 高清在线视频一区二区三区| 男女啪啪激烈高潮av片| 一区二区三区四区激情视频| 精品亚洲乱码少妇综合久久| 肉色欧美久久久久久久蜜桃| 2018国产大陆天天弄谢| 亚洲欧美日韩另类电影网站 | 少妇丰满av| 97热精品久久久久久| 亚洲三级黄色毛片| 欧美亚洲 丝袜 人妻 在线| 高清不卡的av网站| 亚洲精品一区蜜桃| 寂寞人妻少妇视频99o| 狂野欧美激情性xxxx在线观看| 国产精品伦人一区二区| 精品国产一区二区三区久久久樱花 | 国产精品秋霞免费鲁丝片| 精品一区二区三区视频在线| 亚洲精品一二三| 一区二区三区精品91| 日韩欧美精品免费久久| 大又大粗又爽又黄少妇毛片口| 中文欧美无线码| 一级毛片我不卡| av专区在线播放| 久久99精品国语久久久| 久久 成人 亚洲| 99久久人妻综合| 国产免费一级a男人的天堂| 超碰av人人做人人爽久久| 欧美高清性xxxxhd video| 亚洲av欧美aⅴ国产| 观看美女的网站| 久久精品国产鲁丝片午夜精品| 国产精品一区二区性色av| 观看美女的网站| 国产 一区 欧美 日韩| 午夜福利高清视频| av在线播放精品| 久久婷婷青草| 精品国产乱码久久久久久小说| 国产老妇伦熟女老妇高清| 日本黄色日本黄色录像| 国产一区亚洲一区在线观看| 日韩一区二区三区影片| 网址你懂的国产日韩在线| 亚洲国产欧美人成| 五月开心婷婷网| 97超视频在线观看视频| 一级黄片播放器| 人体艺术视频欧美日本| 日韩强制内射视频| 久久精品国产亚洲av天美| 日韩一本色道免费dvd| 18禁在线无遮挡免费观看视频| 美女中出高潮动态图| 亚洲真实伦在线观看| 我要看黄色一级片免费的| 在线精品无人区一区二区三 | 丝袜喷水一区| av在线观看视频网站免费| 性色avwww在线观看| 国产成人一区二区在线| 亚洲av男天堂| 97精品久久久久久久久久精品| 久久久久久久久久成人| 亚洲国产色片| 久久鲁丝午夜福利片| 国产一区亚洲一区在线观看| 久久精品久久久久久久性| 高清黄色对白视频在线免费看 | 偷拍熟女少妇极品色| av在线老鸭窝| 91aial.com中文字幕在线观看| 一级毛片我不卡| 国产精品国产av在线观看| 亚洲伊人久久精品综合| 久久久久国产精品人妻一区二区| 看免费成人av毛片| 国产亚洲午夜精品一区二区久久| 亚洲av不卡在线观看| 男的添女的下面高潮视频| 好男人视频免费观看在线| 午夜免费男女啪啪视频观看| 视频区图区小说| 日本黄色片子视频| 国产精品久久久久久久电影| 亚洲精品久久午夜乱码| 在线观看免费视频网站a站| 高清日韩中文字幕在线| 干丝袜人妻中文字幕| av在线蜜桃| 美女内射精品一级片tv| 午夜视频国产福利| 亚洲av男天堂| 美女高潮的动态| 男女免费视频国产| 久久99热这里只频精品6学生| 人妻 亚洲 视频| 最近最新中文字幕免费大全7| 亚洲国产成人一精品久久久| 中文字幕av成人在线电影| 精品久久久久久久末码| 欧美日本视频| 成人国产av品久久久| 美女xxoo啪啪120秒动态图| 国产色爽女视频免费观看| 国产真实伦视频高清在线观看| 国产精品熟女久久久久浪| 久久久久久久久久久免费av| 精品一区在线观看国产| 国产亚洲av片在线观看秒播厂| 少妇人妻一区二区三区视频| 亚洲欧美中文字幕日韩二区| 老司机影院毛片| 天美传媒精品一区二区| 日韩在线高清观看一区二区三区| 18禁裸乳无遮挡动漫免费视频| 亚洲国产成人一精品久久久| 国产黄色免费在线视频| 亚洲av免费高清在线观看| 亚洲成色77777| 人人妻人人添人人爽欧美一区卜 | 午夜福利在线在线| 午夜福利影视在线免费观看| 久久鲁丝午夜福利片| a级毛片免费高清观看在线播放| 国产大屁股一区二区在线视频| 免费观看a级毛片全部| 国产人妻一区二区三区在| 日本免费在线观看一区| 一二三四中文在线观看免费高清| 色5月婷婷丁香| 国产精品爽爽va在线观看网站| 欧美一级a爱片免费观看看| 99久久中文字幕三级久久日本| 久久久久久久国产电影| 国产精品福利在线免费观看| 丝袜脚勾引网站| 你懂的网址亚洲精品在线观看| 一级毛片我不卡| 高清在线视频一区二区三区| 国产 一区精品| 少妇 在线观看| 亚洲av欧美aⅴ国产| 免费观看的影片在线观看| 亚洲av二区三区四区| 日本av手机在线免费观看| kizo精华| 久久久精品94久久精品| 亚洲人成网站高清观看| 草草在线视频免费看| 建设人人有责人人尽责人人享有的 | 卡戴珊不雅视频在线播放| 日日摸夜夜添夜夜添av毛片| 亚洲国产欧美在线一区| 精品亚洲成a人片在线观看 | 欧美 日韩 精品 国产| 精品人妻偷拍中文字幕| 久久国内精品自在自线图片| 91狼人影院| 激情 狠狠 欧美| 欧美97在线视频| av免费在线看不卡| 熟妇人妻不卡中文字幕| 在线观看免费日韩欧美大片 | 亚洲av成人精品一二三区| 超碰97精品在线观看| 国产极品天堂在线| 亚洲人与动物交配视频| 三级国产精品欧美在线观看| 亚洲精品国产av成人精品| 欧美区成人在线视频| 少妇的逼好多水| 国产一区有黄有色的免费视频| 色婷婷av一区二区三区视频| 久久精品国产亚洲av天美| 精品久久久久久久久亚洲| 免费在线观看成人毛片| 欧美精品人与动牲交sv欧美| 久久久色成人| 制服丝袜香蕉在线| 一级毛片我不卡| 亚洲成色77777| 久久 成人 亚洲| 岛国毛片在线播放| 天堂8中文在线网| 内射极品少妇av片p| 日韩强制内射视频| 久热这里只有精品99| 2018国产大陆天天弄谢| 在线观看一区二区三区激情| 我的女老师完整版在线观看| 亚洲不卡免费看| 国内揄拍国产精品人妻在线| 超碰97精品在线观看| 久久99精品国语久久久| 狂野欧美激情性bbbbbb| 午夜福利高清视频| 久久久久精品性色| 成人毛片a级毛片在线播放| 在线免费十八禁| av在线蜜桃| 男人和女人高潮做爰伦理| 你懂的网址亚洲精品在线观看| 亚洲国产最新在线播放| 亚洲色图av天堂| 亚洲aⅴ乱码一区二区在线播放| 亚洲欧美日韩东京热| 国产av国产精品国产| 大香蕉久久网| 欧美zozozo另类| 岛国毛片在线播放| 亚洲aⅴ乱码一区二区在线播放| 成人高潮视频无遮挡免费网站| 在线观看美女被高潮喷水网站| 高清视频免费观看一区二区| 少妇裸体淫交视频免费看高清| 亚洲精品一二三| 青春草亚洲视频在线观看| 18+在线观看网站| 国产精品欧美亚洲77777| 国产av一区二区精品久久 | 日韩 亚洲 欧美在线| 久久热精品热| 在线免费观看不下载黄p国产| 欧美亚洲 丝袜 人妻 在线| 国产亚洲一区二区精品| 久久国产乱子免费精品| 美女视频免费永久观看网站| 国产一区二区在线观看日韩| 熟妇人妻不卡中文字幕| freevideosex欧美| 久久99热6这里只有精品| 久久99精品国语久久久| 在线免费观看不下载黄p国产| 亚洲精品久久久久久婷婷小说| 国产精品国产三级国产专区5o| 日本av免费视频播放| 亚洲美女搞黄在线观看| 99久久中文字幕三级久久日本| videos熟女内射| 日韩中文字幕视频在线看片 | 夫妻午夜视频| 蜜桃在线观看..| 一级二级三级毛片免费看| 久久久久久久久久成人| 高清不卡的av网站| 久久久久久久久久人人人人人人| 亚洲国产最新在线播放| 午夜福利在线观看免费完整高清在| 成年人午夜在线观看视频| 亚洲av国产av综合av卡| 超碰97精品在线观看| 欧美bdsm另类| 国产乱人偷精品视频| 一本—道久久a久久精品蜜桃钙片| 精品久久久噜噜| 精品久久国产蜜桃| 久久av网站| 黄色怎么调成土黄色| 看免费成人av毛片| 午夜老司机福利剧场| 寂寞人妻少妇视频99o| 简卡轻食公司| 汤姆久久久久久久影院中文字幕| 美女视频免费永久观看网站| 最新中文字幕久久久久| 国产 一区 欧美 日韩| 少妇的逼好多水| 亚洲精品一二三| 亚洲色图综合在线观看| av天堂中文字幕网| 99久国产av精品国产电影| 91久久精品电影网| 久久99热6这里只有精品| 国产一区二区在线观看日韩| 欧美日韩亚洲高清精品| www.色视频.com| 亚洲第一av免费看| 精品一区在线观看国产| 欧美激情国产日韩精品一区| 男男h啪啪无遮挡| 国产成人精品久久久久久| 97精品久久久久久久久久精品| 久久99蜜桃精品久久| 久久久久久久国产电影| 日韩中文字幕视频在线看片 | 国产男人的电影天堂91| 午夜视频国产福利| 91在线精品国自产拍蜜月| 九九在线视频观看精品| 精品一区二区三区视频在线| av在线观看视频网站免费| 久久久久久久大尺度免费视频| 一级毛片 在线播放| 99久久精品一区二区三区| 天美传媒精品一区二区| 大陆偷拍与自拍| 国产片特级美女逼逼视频| 午夜福利在线在线| 蜜桃亚洲精品一区二区三区| 亚洲综合色惰| 国产伦精品一区二区三区四那| 综合色丁香网| 久久国产精品大桥未久av | 国产无遮挡羞羞视频在线观看| 99热网站在线观看| 小蜜桃在线观看免费完整版高清|