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

    一種多自由度阻尼體系動(dòng)力問(wèn)題的高階分析方法

    2021-01-27 08:51:20李鴻晶楊筱朋梅雨辰
    工程力學(xué) 2021年1期
    關(guān)鍵詞:時(shí)段計(jì)算結(jié)果動(dòng)力

    李鴻晶,楊筱朋,梅雨辰

    (南京工業(yè)大學(xué)工程力學(xué)研究所,南京 211816)

    結(jié)構(gòu)動(dòng)力響應(yīng)分析是認(rèn)識(shí)動(dòng)態(tài)系統(tǒng)性態(tài)和行為的基礎(chǔ)性工作。傳統(tǒng)的動(dòng)力響應(yīng)分析方法可以大致分為兩種類(lèi)型:一類(lèi)是基于坐標(biāo)變換[1?2]的方法,如振型疊加法、頻域法等,一般適用于線(xiàn)彈性體系和經(jīng)典阻尼體系的動(dòng)力問(wèn)題;另一類(lèi)為直接積分法[3?6],其直接對(duì)體系運(yùn)動(dòng)微分方程進(jìn)行求解,以Newmark-β 法等逐步積分法為代表,可以解決非線(xiàn)性體系和非經(jīng)典阻尼體系的動(dòng)力分析問(wèn)題。這些傳統(tǒng)動(dòng)力分析方法中有些算法具備無(wú)條件穩(wěn)定及可控算法阻尼等特性[7],多采用低階格式,精度一般不超過(guò)二階。也有一些算法采用了高階格式,精度能達(dá)到三階以上,但所需的計(jì)算量及存儲(chǔ)量要增大不少[8]。有別于傳統(tǒng)動(dòng)力方法,精細(xì)積分法(precise time-step integration method,PTIM)[9?18]則利用體系動(dòng)力響應(yīng)積分解,通過(guò)計(jì)算矩陣指數(shù)獲得動(dòng)力響應(yīng)。由于結(jié)構(gòu)大型化、復(fù)雜化以及精細(xì)化分析等方面的需要,發(fā)展高效率的動(dòng)力分析方法始終是一種非?,F(xiàn)實(shí)的需求。

    微分求積 (differential quadrature, DQ)法[19?22]是近年發(fā)展起來(lái)的一種求解微分方程(組)的數(shù)值方法,F(xiàn)ung[23?25]首先將其應(yīng)用于初值問(wèn)題,其后一些學(xué)者也開(kāi)展了相關(guān)研究[26?32],逐漸形成了動(dòng)力DQ 分析方法。該方法數(shù)值穩(wěn)定性好,計(jì)算簡(jiǎn)便且精度高,然而由于其以位移響應(yīng)為基本未知量來(lái)離散動(dòng)力方程,所得數(shù)值格式為隱式格式,因而需要求解方程。另外,目前已提出的這類(lèi)動(dòng)力方法因穩(wěn)定性和精度方面的要求多采用不等步距的時(shí)間離散方案,求解均勻布置的采樣時(shí)刻處的響應(yīng)值時(shí)需要進(jìn)行插值,對(duì)其計(jì)算效率造成了一定的限制。與傳統(tǒng)動(dòng)力響應(yīng)DQ 分析方法不同,在文獻(xiàn)[33 ? 34]中,以單自由度體系為對(duì)象提出了一種計(jì)算Duhamel 積分的高精度算法,即積分求微方法(integral differentiation method, IDM)。該算法為高階方法,動(dòng)力響應(yīng)計(jì)算僅需要完成一次方陣與向量的乘法運(yùn)算,一次計(jì)算即可同時(shí)獲得多個(gè)時(shí)刻的響應(yīng)值。整個(gè)過(guò)程無(wú)需求解方程,且可采用等時(shí)間步距的計(jì)算方案,不必通過(guò)插值即可直接獲得各等距分布時(shí)點(diǎn)的響應(yīng)。積分求微法可視為非齊次動(dòng)力方程的一種數(shù)值解法,如果將其向多自由度體系推廣,則將引入矩陣指數(shù)。由于積分求微法無(wú)需進(jìn)行數(shù)值積分計(jì)算,與單自由度體系動(dòng)力響應(yīng)分析過(guò)程相比較,多自由度體系動(dòng)力分析僅需要處理矩陣指數(shù)的計(jì)算。

    本文將逆用DQ 原理,建立多自由度阻尼體系動(dòng)力時(shí)程分析的積分求微算法,以豐富動(dòng)力響應(yīng)分析的技術(shù)手段。

    1 多自由度體系狀態(tài)響應(yīng)

    動(dòng)態(tài)荷載向量p(t)激勵(lì)下的多自由度阻尼體系運(yùn)動(dòng)方程寫(xiě)為:

    狀態(tài)方程式(3)中的系數(shù)矩陣H 僅包含體系剛度、質(zhì)量和阻尼的特性,因而它是反映體系動(dòng)態(tài)特性的。對(duì)于線(xiàn)彈性體系動(dòng)力分析,體系特性矩陣H 將是時(shí)不變的常數(shù)矩陣。而方程的非齊次項(xiàng)中含有荷載向量p(t),它描述了環(huán)境對(duì)體系施加的動(dòng)態(tài)作用的特性。

    式(3)的解可采用 Duhamel 積分的形式表示為:

    式中, eHt為矩陣指數(shù)。

    形如式(5)的方程解中包含了矩陣指數(shù)及其積分式,對(duì)于實(shí)際動(dòng)力問(wèn)題,試圖利用式(5)獲得體系動(dòng)力響應(yīng)y(t)的解析解答幾乎是不可能的,但可以通過(guò)數(shù)值方法得到y(tǒng)(t)的近似結(jié)果。

    2 動(dòng)力響應(yīng)積分求微分析

    利用式(5)計(jì)算y(t)的難點(diǎn)在于含有矩陣指數(shù)項(xiàng)的卷積的計(jì)算。將該卷積做如下的分離:

    顯見(jiàn),求得動(dòng)力響應(yīng)y(t)的關(guān)鍵是處理好兩個(gè)問(wèn)題:一是矩陣指數(shù)T(t)的計(jì)算;二是計(jì)算出變上限定積分函數(shù)向量D(t)。兩者都要求建立起適合的數(shù)值算法。

    考慮動(dòng)力持時(shí)中的任意時(shí)段[ti, tj],該時(shí)段包含有ρ 個(gè)時(shí)步,其中ti、tj均選擇位于采樣點(diǎn)處。如圖1 所示,假定采取等距采樣,采樣周期為δ,則時(shí)段[ti, tj]的長(zhǎng)度應(yīng)為:

    圖 1 動(dòng)態(tài)響應(yīng)待求解時(shí)段Fig. 1 The solved time interval for dynamic response

    2.1 矩陣指數(shù)Tk 的計(jì)算

    根據(jù)矩陣指數(shù)的Taylor 級(jí)數(shù)定義,矩陣指數(shù)Tk可表示為:

    2.2 定積分 Dk 的計(jì)算

    由式(8)可見(jiàn),D(t)對(duì)時(shí)間t 的導(dǎo)數(shù)可表示為:

    觀(guān)察上述求解動(dòng)力響應(yīng)的過(guò)程,可以發(fā)現(xiàn)除了在確定系統(tǒng)特性矩陣H 時(shí)需要對(duì)質(zhì)量矩陣m 求逆,以及在確定轉(zhuǎn)換矩陣B 時(shí)需要對(duì)Vandermonde矩陣進(jìn)行求逆計(jì)算外,其余過(guò)程都不需要矩陣求逆或者求解線(xiàn)性方程組。而對(duì)于結(jié)構(gòu)體系可采用集中質(zhì)量模型,即質(zhì)量矩陣m 為對(duì)角陣,其逆矩陣可以方便地求出。至于Vandermonde 矩陣求逆計(jì)算,可以事先算出逆矩陣結(jié)果,考慮到實(shí)際動(dòng)力分析時(shí)待求解時(shí)段 [0, ?tij]內(nèi)包含的時(shí)步數(shù)ρ 不可能很大 (一般取 ρ=10~15),故可將轉(zhuǎn)換矩陣B 制成表格,動(dòng)力響應(yīng)分析時(shí)直接取用即可。鑒于此,按照式(32)計(jì)算積分矩陣Φ 時(shí)僅需做矩陣乘法運(yùn)算,這相當(dāng)于是一種顯式算法。

    3 算法步驟

    依據(jù)上述推導(dǎo),可將多自由度阻尼體系動(dòng)力時(shí)程積分求微法總結(jié)為下述算法:

    4 數(shù)值穩(wěn)定性分析

    對(duì)于線(xiàn)性體系動(dòng)力方法的數(shù)值穩(wěn)定性,只需研究如下的單自由度體系方程:

    阻尼一般對(duì)穩(wěn)定性起到有利作用,因此研究穩(wěn)定性時(shí)多不考慮阻尼的影響。以下分析都是在消除阻尼項(xiàng)后開(kāi)展的。觀(guān)察式(12)的動(dòng)力響應(yīng)數(shù)值格式,可得時(shí)段初始點(diǎn)與末點(diǎn)的遞推關(guān)系式:

    表 1 L 不同取值時(shí)的穩(wěn)定極限值Table 1 The limit of stability with different L

    由此可得出,本文方法一般為有條件穩(wěn)定方法,當(dāng)L=3、4、7、8 時(shí),穩(wěn)定極限較大,穩(wěn)定性相對(duì)較好,再綜合考慮計(jì)算量的因素,取L=4 比較合適。由于在精細(xì)積分法中λ 一般選20,即時(shí)步小段的長(zhǎng)度 η ≈ 9.54×10?7δ,因而只需要各節(jié)點(diǎn)間距滿(mǎn)足:

    即可穩(wěn)定地進(jìn)行計(jì)算。顯然,通常情況下該穩(wěn)定極限條件是極易得到滿(mǎn)足的,當(dāng)取L=4 時(shí)本文方法時(shí)步大小的取值幾乎不受穩(wěn)定性的限制。

    5 算例分析

    通過(guò)兩個(gè)算例闡釋本文方法的特性,第一個(gè)算例用來(lái)驗(yàn)證本文方法的可靠性,包括收斂性和精度;第二個(gè)算例重點(diǎn)闡釋本文方法的適用性。

    5.1 兩自由度體系算例

    選擇文獻(xiàn)[12]中的兩自由度體系算例,其特性矩陣分別為:

    采用本文的IDM 法計(jì)算該體系的位移響應(yīng)和速度響應(yīng),并以精確解為標(biāo)準(zhǔn)驗(yàn)證計(jì)算結(jié)果。為了清晰地展現(xiàn)本文方法的收斂性及精確性,將時(shí)段內(nèi)時(shí)步數(shù)ρ 取為10,時(shí)段長(zhǎng)度Δt 從1.0 s 開(kāi)始依次減為原來(lái)的1/2 時(shí),即分別取Δt=1 s、1/2 s、1/22s、1/23s 時(shí),計(jì)算體系在 t=0.2 s、0.4 s、0.6 s、0.8 s、1.0 s 處的響應(yīng)值。計(jì)算結(jié)果同精確解的最大絕對(duì)誤差Err 與Δt 的關(guān)系曲線(xiàn)如圖2 和圖3 所示,其中橫坐標(biāo)和縱坐標(biāo)皆采用以2 為底的對(duì)數(shù)坐標(biāo)。

    觀(guān)察圖2 和圖3 可以發(fā)現(xiàn),隨著時(shí)段長(zhǎng)度Δt 的減小,本文方法計(jì)算結(jié)果的誤差亦迅速減小。當(dāng)Δt 縮小到1/23s 時(shí),無(wú)論是位移響應(yīng)還是速度響應(yīng),其最大絕對(duì)誤差皆不超過(guò)10?14數(shù)量級(jí),計(jì)算結(jié)果已非常精確,具備較完美的收斂效果。從整個(gè)時(shí)步的取值過(guò)程來(lái)看,本文IDM 法的收斂速度非???,幾乎是時(shí)步長(zhǎng)度每縮小一倍計(jì)算結(jié)果的絕對(duì)誤差就會(huì)減小約2?10量級(jí)。進(jìn)一步觀(guān)察可以直觀(guān)地發(fā)現(xiàn)在對(duì)數(shù)坐標(biāo)系下最大絕對(duì)誤差與時(shí)步長(zhǎng)度近似呈線(xiàn)性關(guān)系,這符合數(shù)值計(jì)算原理中誤差與時(shí)步長(zhǎng)度的關(guān)系:

    式中:絕對(duì)誤差Err 對(duì)于函數(shù)可以用其范數(shù)來(lái)近似替代;q 為數(shù)值格式的收斂精度階數(shù)。

    圖 2 不同Δt 下體系位移響應(yīng)最大絕對(duì)誤差Fig. 2 The maximum absolute error of the displacement response of the system with different Δt

    圖 3 不同Δt 下體系速度響應(yīng)最大絕對(duì)誤差Fig. 3 The maximum absolute error of the velocity response of the system with different Δt

    為合理有效地估計(jì)本文IDM 方法的收斂階數(shù),分別選取整體誤差指標(biāo)——最大絕對(duì)誤差以及某個(gè)時(shí)刻(如t=0.4 s)時(shí)的局部絕對(duì)誤差兩個(gè)指標(biāo)作為式(41)中的誤差項(xiàng)Err,采用最小二乘法估算出體系位移響應(yīng)和速度響應(yīng)的收斂階數(shù)q,結(jié)果列于表2 和表3。

    表 2 響應(yīng)收斂階數(shù)q 估計(jì)(整體誤差指標(biāo))Table 2 Evaluation of convergence order q for responses based on whole error criterion

    表 3 響應(yīng)收斂階數(shù)q 估計(jì)(局部誤差指標(biāo))Table 3 Evaluation of convergence order q for responses based on local error criterion

    由表2 和表3 的統(tǒng)計(jì)結(jié)果,無(wú)論是采用整體誤差指標(biāo)還是局部誤差指標(biāo),各質(zhì)點(diǎn)的位移響應(yīng)和速度響應(yīng)的收斂階數(shù)皆在10~11,即該數(shù)值格式對(duì)于位移響應(yīng)和速度響應(yīng)的收斂精度至少為10 階。而加速度響應(yīng)是依據(jù)位移響應(yīng)和速度響應(yīng)由運(yùn)動(dòng)方程求得的,因而其具備與位移響應(yīng)、速度響應(yīng)相同的收斂精度。另一方面,從理論上來(lái)看,本文IDM 數(shù)值格式的收斂精度主要由矩陣指數(shù)的計(jì)算精度和DQ 法的逼近精度所共同決定。而矩陣指數(shù)的計(jì)算由精細(xì)積分實(shí)現(xiàn),其數(shù)值逼近的時(shí)間節(jié)點(diǎn)間距非常小,計(jì)算精度相較于非齊次項(xiàng)DQ 逼近要高許多。因此可以認(rèn)為,本文IDM法的收斂精度主要由DQ 法的計(jì)算精度決定,它直接取決于時(shí)段內(nèi)的時(shí)步數(shù)ρ 的取值。由DQ 原理可知,在各時(shí)段Δt 內(nèi)對(duì)其積分函數(shù)D(τ)逼近的絕對(duì)誤差函數(shù)為:

    從式(42)來(lái)看,其誤差函數(shù)是Δt 的ρ+1 階無(wú)窮小,即本文格式具有ρ+1 階代數(shù)精度。當(dāng)取ρ=10 時(shí),本文方法收斂精度階數(shù)的理論值為11。考慮到計(jì)算機(jī)舍入誤差及矩陣指數(shù)的影響,實(shí)際計(jì)算時(shí)收斂精度會(huì)略微小于理論值。即實(shí)際計(jì)算時(shí)若取ρ=10,本文方法至少能達(dá)到10 階收斂精度,這同表2 和表3 中估計(jì)的結(jié)果相符合。一般地,當(dāng)計(jì)算時(shí)段內(nèi)包含的時(shí)步數(shù)為ρ 時(shí),理論上本文方法有ρ+1 階代數(shù)精度,實(shí)際計(jì)算時(shí)至少能達(dá)到ρ 階精度。

    下面將本文IDM 法與傳統(tǒng)動(dòng)力方法進(jìn)行比較,選擇的傳統(tǒng)方法包括Newmark 常平均加速度法 (取 γ=0.5、β=0.25)、Wilson-θ 法 (取 θ=1.4)和精細(xì)積分法(PTIM)。將所有這些方法的計(jì)算結(jié)果列于表4,同時(shí)給出了精確解作為比較的標(biāo)準(zhǔn)。

    表 4 不同動(dòng)力方法計(jì)算的位移響應(yīng)結(jié)果Table 4 Computational results of displacement response by different dynamic methods

    續(xù)表 4

    觀(guān)察表4 可以發(fā)現(xiàn),當(dāng)時(shí)段長(zhǎng)度取為Δt=1.0 s時(shí),劃分10 個(gè)時(shí)步,本文積分求微法的計(jì)算結(jié)果與精確解已非常接近,最大相對(duì)誤差在2%以?xún)?nèi)。

    注意到本文IDM 法為高階方法,此時(shí)待求解時(shí)段[ti, tj]內(nèi)含有的10 個(gè)時(shí)步的步距均為δ=0.1 s,而傳統(tǒng)的 Newmark-β 法和 Wilson-θ 法在步距取為0.05 s 時(shí)位移響應(yīng)的最大相對(duì)誤差已超過(guò)了10%,在取為0.02 s 時(shí),也超過(guò)了2%。也就是說(shuō),在步距是這兩種傳統(tǒng)逐步積分方法5 倍的條件下,本文IDM 法的計(jì)算結(jié)果仍然比其要精確些,這直觀(guān)顯示了本文方法高精度的特點(diǎn)。

    5.2 彈簧-質(zhì)量體系算例

    以圖4 所示的串聯(lián)彈簧質(zhì)量體系為例,取自由度數(shù)N=100。為了使體系固有頻率覆蓋范圍廣,彈簧剛度一半設(shè)置為較剛,另一半設(shè)置相對(duì)較柔,其特性參數(shù)見(jiàn)表5。每個(gè)質(zhì)點(diǎn)都受到相同的水平簡(jiǎn)諧荷載激勵(lì):

    圖 4 100 自由度彈簧質(zhì)量體系Fig. 4 Schematic plot of 100-degree of freedom spring-mass system

    表 5 彈簧-質(zhì)量體系的基本特性Table 5 Properties of the spring-mass system

    該100 自由度體系的基本頻率為0.98 rad/s,最高階頻率為1094.92 rad/s,從低頻覆蓋到極高頻,便于全面地考察本文方法的數(shù)值穩(wěn)定性。取截?cái)囗?xiàng)數(shù)L=4、時(shí)段分段數(shù)ρ=10,采用本文IDM法計(jì)算體系每隔0.01 s 的動(dòng)力響應(yīng)。選用不同的時(shí)段長(zhǎng)度Δt,計(jì)算結(jié)果與顯式中心差分法(取Δt=0.001 s)進(jìn)行比較。圖5 和圖6 分別為最右端質(zhì)點(diǎn)(DOF 100)和第50 個(gè)質(zhì)點(diǎn)(DOF 50)的位移時(shí)程曲線(xiàn)。這里需要說(shuō)明的是,對(duì)于該體系中心差分法的穩(wěn)定極限為Δt≤0.0018 s,因而采用Δt=0.001 s可得到穩(wěn)定的計(jì)算結(jié)果,且它是直接獲得0.01 s 時(shí)間間隔響應(yīng)所允許的最大步距。

    圖 5 體系第100 自由度位移響應(yīng)對(duì)比Fig. 5 Comparison of the displacement responses of the DOF 100 in this system

    圖 6 體系第50 自由度位移響應(yīng)對(duì)比Fig. 6 Comparison of the displacement responses of the DOF 50 in this system

    觀(guān)察圖5 和圖6 的動(dòng)力響應(yīng)計(jì)算結(jié)果,IDM法在選擇大得多的步距情形下依然能夠獲得穩(wěn)定的計(jì)算結(jié)果。例如,將時(shí)步長(zhǎng)度放大至0.1 s 時(shí),利用IDM 法計(jì)算最右端質(zhì)點(diǎn)的位移響應(yīng),仍然沒(méi)有任何失穩(wěn)放大現(xiàn)象產(chǎn)生。此情形下其時(shí)段長(zhǎng)度已是中心差分法穩(wěn)定極限的50 倍以上,在較高自振頻率的參與下依然能得到穩(wěn)定可靠的結(jié)果,這充分顯示了本文方法良好的數(shù)值穩(wěn)定性。

    下面進(jìn)一步考察本文IDM 法的計(jì)算效率。為了更直觀(guān)地體現(xiàn)算法的執(zhí)行效率,表6 給出了一些情況下IDM 法與中心差分法(CDM)耗費(fèi)的計(jì)算時(shí)間及兩種方法計(jì)算時(shí)間的比值。其中所有的計(jì)算皆是在同一臺(tái)計(jì)算機(jī)上完成的,該計(jì)算機(jī)的架構(gòu)為英特爾Xeon(R)W3565 中央處理器,主頻3.20 GHz。各類(lèi)情況的計(jì)算時(shí)間均用符號(hào)CPU 來(lái)表示,且該時(shí)間是純計(jì)算時(shí)間,不包括數(shù)據(jù)輸出及單位換算的時(shí)間。

    表 6 不同方法CPU 時(shí)間對(duì)比Table 6 Comparison of CPU time with different methods

    觀(guān)察表6,對(duì)于IDM 法,當(dāng)時(shí)段長(zhǎng)度取為0.05 s 時(shí)其計(jì)算時(shí)間還不到中心差分法的40%,當(dāng)時(shí)步縮小到0.025 s 時(shí)也只有CDM 法的70%左右。通過(guò)對(duì)比各響應(yīng)時(shí)程計(jì)算結(jié)果,可以發(fā)現(xiàn)對(duì)于DOF 100 質(zhì)點(diǎn)的位移響應(yīng),IDM (Δt=0.05 s)和IDM (Δt=0.025 s)的計(jì)算結(jié)果與CDM 的結(jié)果基本一致,而對(duì)于DOF 50 質(zhì)點(diǎn),IDM (Δt=0.025 s)的計(jì)算結(jié)果與中心差分法基本一致,而IDM (Δt=0.05 s)的精度稍差一點(diǎn)。但即便如此,在IDM 法與CDM法精度大致接近的情況下,其計(jì)算時(shí)間依然比中心差分法要少30%左右,從而在時(shí)間復(fù)雜度上IDM 法要低于傳統(tǒng)的顯式中心差分法。另外,IDM(Δt=0.025 s)在整個(gè)持時(shí)內(nèi)的時(shí)間采樣點(diǎn)數(shù)量?jī)H為CDM(Δt=0.001 s)的 40%,雖然 IDM 法結(jié)構(gòu)特性矩陣的元素?cái)?shù)量是中心差分法的4 倍,但對(duì)于線(xiàn)性時(shí)不變體系而言?xún)H需存儲(chǔ)一次,且有近1/2 是零元素。相較于中心差分法,IDM 法憑借更少的時(shí)間采樣點(diǎn)數(shù),隨著荷載持時(shí)的延長(zhǎng)其空間復(fù)雜度的優(yōu)勢(shì)必將愈加明顯。綜上所述,本文IDM 法在計(jì)算復(fù)雜度(包括時(shí)間復(fù)雜度和空間復(fù)雜度)上要低于傳統(tǒng)的顯式CDM 法,由于CDM 已是目前效率較高的動(dòng)力分析方法,表明本文IDM 法在計(jì)算效率方面具有一定優(yōu)勢(shì)。

    6 討論與結(jié)論

    根據(jù)DQ 原理,函數(shù)在指定點(diǎn)處的導(dǎo)數(shù)值可以用域內(nèi)一系列離散點(diǎn)處函數(shù)值的加權(quán)和來(lái)近似表示。將DQ 原理用于偏微分方程求解時(shí),傳統(tǒng)上都是將域內(nèi)離散點(diǎn)處的函數(shù)值視為待求解的未知量,然后利用DQ 原理將各離散點(diǎn)處的導(dǎo)數(shù)值分別以函數(shù)值線(xiàn)性加權(quán)和替代并代入微分方程,從而將微分方程轉(zhuǎn)化為代數(shù)方程組,進(jìn)而獲得微分方程的近似解答。由于DQ 法高精度、低計(jì)算消耗的特點(diǎn),這種基于DQ 原理的微分方程數(shù)值解法僅需付出不大的計(jì)算工作量即可獲得較高精度的數(shù)值解。對(duì)于時(shí)間相關(guān)問(wèn)題,傳統(tǒng)的基于DQ原理的動(dòng)力響應(yīng)分析方法都是沿用了這種求解思路,即將各離散時(shí)刻的速度響應(yīng)和加速度響應(yīng)分別表示為時(shí)段(步)內(nèi)全部離散時(shí)刻位移響應(yīng)的加權(quán)和。但是像動(dòng)力響應(yīng)Duhamel 積分解這樣的在被積函數(shù)中含有指數(shù)函數(shù)乘積項(xiàng)的定積分式,在積分變量分離后其導(dǎo)數(shù)即為被積函數(shù)自身,所以這種類(lèi)型卷積的導(dǎo)數(shù)是非常容易求得的。同傳統(tǒng)的DQ 法應(yīng)用不同,本文的積分求微法中離散點(diǎn)處的導(dǎo)數(shù)值是在應(yīng)用DQ 原理前就已經(jīng)求得的,這樣只需通過(guò)權(quán)系數(shù)逆矩陣進(jìn)行變換即可求得位移響應(yīng)向量,即動(dòng)力響應(yīng)計(jì)算僅需執(zhí)行矩陣乘法運(yùn)算即可一次性地求出時(shí)段內(nèi)的全部位移響應(yīng)值。而權(quán)系數(shù)矩陣的求逆計(jì)算只涉及Vandermonde矩陣求逆,這可以事先求出并在程序中存儲(chǔ)。因此,本文的動(dòng)力響應(yīng)積分求微方法相當(dāng)于顯式方法,其關(guān)鍵在于對(duì)DQ 原理的逆用。

    經(jīng)過(guò)逆用DQ 原理后,動(dòng)力響應(yīng)分析只需解決矩陣指數(shù)計(jì)算即可完成。PTIM 是實(shí)現(xiàn)矩陣指數(shù)計(jì)算的一種優(yōu)秀算法,具體應(yīng)用到本文動(dòng)力響應(yīng)積分求微方法中,需要注意幾個(gè)問(wèn)題:

    首先,本文動(dòng)力響應(yīng)積分求微法是按照時(shí)段進(jìn)行求解的,每個(gè)待求時(shí)段內(nèi)都包含有若干個(gè)時(shí)步。根據(jù)矩陣指數(shù)加法定理,在等距時(shí)步采樣即待求時(shí)段內(nèi)各時(shí)步的步距相同的前提下(事實(shí)上動(dòng)力分析一般都會(huì)采取這種等時(shí)步方案),各個(gè)時(shí)步端點(diǎn)的響應(yīng)值才能表示成像式(19)那樣的矩陣指數(shù)遞推格式。這樣,當(dāng)采取等距時(shí)步方案時(shí)只需要計(jì)算出待求時(shí)段內(nèi)第一個(gè)時(shí)步對(duì)應(yīng)的矩陣指數(shù),然后通過(guò)遞推公式就可以計(jì)算出時(shí)段內(nèi)所有時(shí)步的矩陣指數(shù)。顯然,矩陣指數(shù)精細(xì)積分計(jì)算僅僅被限制在第一個(gè)時(shí)步,即整個(gè)時(shí)段動(dòng)力響應(yīng)分析只要求一次精細(xì)積分計(jì)算,而不必采用PTIM計(jì)算多個(gè)矩陣指數(shù)。這是實(shí)施本文動(dòng)力響應(yīng)積分求微法的一個(gè)要點(diǎn)。也就是說(shuō),在本文積分求微動(dòng)力分析方法中應(yīng)用PTIM 計(jì)算矩陣指數(shù)可以獲得較高的計(jì)算效率,關(guān)鍵在于動(dòng)力時(shí)程等距分布的采樣點(diǎn)同時(shí)也被選擇為DQ 離散節(jié)點(diǎn),這樣就只需要在第一個(gè)時(shí)步內(nèi)應(yīng)用PTIM,即在每個(gè)求解時(shí)段內(nèi)只需使用一次PTIM 完成矩陣指數(shù)計(jì)算即可。

    但是,使用DQ 法時(shí)一般會(huì)選擇不等距節(jié)點(diǎn)離散方案,采取等距節(jié)點(diǎn)方案會(huì)降低DQ 分析的精度。如果采取不等距時(shí)步方案,則勢(shì)必要求多次執(zhí)行矩陣指數(shù)的精細(xì)積分計(jì)算,而且還需要進(jìn)行插值才能獲得動(dòng)力分析結(jié)果,這將極大地增加計(jì)算工作量。所以,具體進(jìn)行動(dòng)力分析時(shí)就需要權(quán)衡計(jì)算的精度和效率。本文建議采取等時(shí)步方案,不僅因?yàn)檫@是動(dòng)力分析的通常選擇,而且可以將計(jì)算工作量控制在最低范圍。由于DQ 法的高精度屬性,即便采取等時(shí)步方案時(shí)也會(huì)取得相當(dāng)高的計(jì)算精度。這從本文給出的算例中可以看出。

    其次,由于PTIM 計(jì)算矩陣指數(shù)時(shí)分為2 個(gè)步驟,需要計(jì)算經(jīng)過(guò)細(xì)分的很小步長(zhǎng)內(nèi)的矩陣指數(shù),即首先需要計(jì)算出矩陣指數(shù) eHη,然后才能通過(guò)遞推得到需要求得的矩陣指數(shù) eHδ(δ 為時(shí)步的步距,本文方法中只要求對(duì) eHδ進(jìn)行精細(xì)積分計(jì)算)。在計(jì)算 eHη時(shí),本文引入了秦九韶算法,將冪函數(shù)求和(被截?cái)嗟腡aylor 級(jí)數(shù))計(jì)算轉(zhuǎn)化為遞推格式,每一步遞推計(jì)算都只做一次矩陣乘法運(yùn)算,有效地降低了計(jì)算工作量和數(shù)據(jù)存儲(chǔ)空間。截?cái)嚯A數(shù)L 越高,采用秦九韶算法的收益就會(huì)越大。

    最后,由于本文方法逆用DQ 原理,所以不要求進(jìn)行數(shù)值積分計(jì)算,但需要計(jì)算矩陣指數(shù)e?Hδ。在使用精細(xì)積分法時(shí),矩陣指數(shù) e?Hδ的計(jì)算可以連同 eHδ一起計(jì)算出來(lái)。只要將細(xì)分后的矩陣指數(shù) eHη的被截?cái)嗉?jí)數(shù)表達(dá)式中的奇數(shù)項(xiàng)和偶數(shù)項(xiàng)分列為2 組,然后分別使用秦九韶算法進(jìn)行計(jì)算,最后通過(guò)兩組結(jié)果相加得到 eHη,相減得到e?Hη。這樣算法可同時(shí)計(jì)算出矩陣指數(shù) eHη和e?Hη,而沒(méi)有增加額外的計(jì)算工作量。其后遞推計(jì)算 e?Hδ的過(guò)程與 eHδ完全一致,可以方便地計(jì)算出來(lái)。

    另外,在文獻(xiàn)[33]中針對(duì)單自由度體系實(shí)施積分求微分析時(shí),待求時(shí)段內(nèi)各時(shí)刻的響應(yīng)值構(gòu)成一個(gè)向量,因而可以只通過(guò)執(zhí)行一次矩陣與向量的乘法運(yùn)算就能計(jì)算出時(shí)段內(nèi)的全部響應(yīng)。而多自由度體系有所不同,因?yàn)榇髸r(shí)段內(nèi)每個(gè)時(shí)刻的響應(yīng)都是一個(gè)向量,故時(shí)段內(nèi)所有時(shí)刻的響應(yīng)實(shí)際上構(gòu)成了一個(gè)矩陣,所以時(shí)段內(nèi)全部響應(yīng)需要通過(guò)矩陣與矩陣的乘法運(yùn)算才能實(shí)現(xiàn)。

    關(guān)于本文積分求微法的穩(wěn)定性問(wèn)題,在文獻(xiàn)[33]中我們已經(jīng)證明其為無(wú)條件穩(wěn)定算法。但是對(duì)于多自由度體系,由于采用了PTIM 計(jì)算矩陣指數(shù),所以該算法的穩(wěn)定性不僅取決于積分求微計(jì)算即逆用DQ 原理本身,還與PTIM 的穩(wěn)定性有關(guān)。關(guān)于PTIM 穩(wěn)定性問(wèn)題已有相關(guān)研究[35],結(jié)論為PTIM 是有條件穩(wěn)定算法,但穩(wěn)定性條件極易滿(mǎn)足,對(duì)于常規(guī)工程結(jié)構(gòu)幾乎不會(huì)出現(xiàn)失穩(wěn)問(wèn)題。事實(shí)上,本文之所以只計(jì)算第一個(gè)時(shí)步的矩陣指數(shù),待求時(shí)段內(nèi)的其他時(shí)步的矩陣指數(shù)都是通過(guò)遞推的方式獲得,也從另一方面說(shuō)明了PTIM 具有非常好的穩(wěn)定性。

    綜上所述,執(zhí)行本文積分求微法的主要環(huán)節(jié)在于DQ 原理逆用和矩陣指數(shù)計(jì)算。通過(guò)引入矩陣指數(shù)及相應(yīng)算法,本文將結(jié)構(gòu)動(dòng)力時(shí)程分析的積分求微法推廣至多自由度情形,通過(guò)理論分析和數(shù)值試驗(yàn),可得到如下結(jié)論:

    (1) 結(jié)構(gòu)動(dòng)力時(shí)程分析的積分求微方法能夠適用于多自由度體系的直接積分,計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)是精確收斂的。該方法是一種高階時(shí)程分析方法,包含ρ 個(gè)時(shí)步的積分求微法理論上具有(ρ+1)階代數(shù)精度,與傳統(tǒng)的低階逐步積分方法相比有較高的計(jì)算精度。

    (2) 本文的積分求微法一般是有條件穩(wěn)定算法,其穩(wěn)定性與截?cái)囗?xiàng)數(shù)L 有關(guān)。當(dāng)L=4 時(shí),穩(wěn)定極限條件極易得到滿(mǎn)足,在較大的步距和較高的體系自振頻率參與下依然能得到穩(wěn)定的計(jì)算結(jié)果,因而具備較好的數(shù)值穩(wěn)定性。

    (3) 本文方法可視為顯式方法,無(wú)需進(jìn)行方程求解,動(dòng)力響應(yīng)分析表現(xiàn)為矩陣與向量的乘法運(yùn)算,一次計(jì)算能得到多個(gè)時(shí)刻的響應(yīng)值,且無(wú)需對(duì)計(jì)算結(jié)果進(jìn)行插值,因而計(jì)算效率較高,計(jì)算工作量要少于傳統(tǒng)的顯式中心差分法。

    (4) 建立本文方法的關(guān)鍵環(huán)節(jié)在于對(duì)DQ 原理的逆用,對(duì)于多自由度阻尼體系動(dòng)力響應(yīng)分析,計(jì)算動(dòng)力方程非齊次項(xiàng)時(shí)不再要求對(duì)定積分式進(jìn)行特殊處理,從而可以直接應(yīng)用精細(xì)積分法計(jì)算矩陣指數(shù),有利于充分發(fā)揮精細(xì)積分法的優(yōu)勢(shì)。

    (5) 本文方法屬于高階方法,在一個(gè)相對(duì)較長(zhǎng)的時(shí)段內(nèi)再劃分出多個(gè)時(shí)步,然后同時(shí)求解。這種模式為進(jìn)一步發(fā)展動(dòng)力響應(yīng)分析并行算法提供了有利條件。

    猜你喜歡
    時(shí)段計(jì)算結(jié)果動(dòng)力
    學(xué)習(xí)動(dòng)力不足如何自給自足
    不等高軟橫跨橫向承力索計(jì)算及計(jì)算結(jié)果判斷研究
    甘肅科技(2020年20期)2020-04-13 00:30:40
    四個(gè)養(yǎng)生黃金時(shí)段,你抓住了嗎
    胖胖一家和瘦瘦一家(10)
    動(dòng)力船
    傍晚是交通事故高發(fā)時(shí)段
    分時(shí)段預(yù)約在PICC門(mén)診維護(hù)中的應(yīng)用與探討
    超壓測(cè)試方法對(duì)炸藥TNT當(dāng)量計(jì)算結(jié)果的影響
    基于多動(dòng)力總成的六點(diǎn)懸置匹配計(jì)算
    噪聲對(duì)介質(zhì)損耗角正切計(jì)算結(jié)果的影響
    久久鲁丝午夜福利片| 操美女的视频在线观看| 黄色毛片三级朝国网站| 蜜桃国产av成人99| 男人添女人高潮全过程视频| 满18在线观看网站| 国产激情久久老熟女| 午夜福利影视在线免费观看| 无遮挡黄片免费观看| 国产欧美日韩精品亚洲av| 亚洲精品自拍成人| 爱豆传媒免费全集在线观看| 王馨瑶露胸无遮挡在线观看| 老汉色∧v一级毛片| 1024香蕉在线观看| 久久久亚洲精品成人影院| 亚洲欧美成人综合另类久久久| 国产视频首页在线观看| 成年人午夜在线观看视频| 男女午夜视频在线观看| 女人久久www免费人成看片| 亚洲国产欧美一区二区综合| 2018国产大陆天天弄谢| 我要看黄色一级片免费的| 亚洲av成人不卡在线观看播放网 | 99国产精品一区二区蜜桃av | 国产一区有黄有色的免费视频| 91精品三级在线观看| 亚洲av国产av综合av卡| 亚洲精品国产一区二区精华液| 亚洲中文日韩欧美视频| 97在线人人人人妻| 欧美日韩综合久久久久久| 男女边摸边吃奶| 可以免费在线观看a视频的电影网站| 啦啦啦视频在线资源免费观看| 久久久久国产精品人妻一区二区| 亚洲国产看品久久| 久久中文字幕一级| 91麻豆av在线| 国产精品99久久99久久久不卡| 日韩制服丝袜自拍偷拍| 你懂的网址亚洲精品在线观看| 欧美日韩一级在线毛片| 精品国产国语对白av| 亚洲av成人不卡在线观看播放网 | 婷婷色综合大香蕉| 国产av国产精品国产| 国产黄色视频一区二区在线观看| 十八禁高潮呻吟视频| 日韩 欧美 亚洲 中文字幕| 80岁老熟妇乱子伦牲交| 国产不卡av网站在线观看| 成年人免费黄色播放视频| 亚洲成色77777| 宅男免费午夜| 黄色一级大片看看| 精品免费久久久久久久清纯 | 欧美人与性动交α欧美精品济南到| 欧美日韩福利视频一区二区| 亚洲成国产人片在线观看| 国产高清国产精品国产三级| 亚洲 欧美一区二区三区| 亚洲欧美清纯卡通| 狂野欧美激情性bbbbbb| 少妇粗大呻吟视频| 在线观看免费视频网站a站| 国产免费现黄频在线看| 丝袜人妻中文字幕| 啦啦啦视频在线资源免费观看| 国产精品国产三级国产专区5o| 久久亚洲精品不卡| 满18在线观看网站| 熟女少妇亚洲综合色aaa.| 97人妻天天添夜夜摸| 精品人妻在线不人妻| 久久久久久久久久久久大奶| 久久综合国产亚洲精品| 麻豆乱淫一区二区| 99国产精品免费福利视频| 欧美精品一区二区免费开放| av国产久精品久网站免费入址| 丝瓜视频免费看黄片| 国产一区二区三区综合在线观看| 18禁黄网站禁片午夜丰满| 亚洲国产欧美网| 精品人妻一区二区三区麻豆| 国产有黄有色有爽视频| 日韩大片免费观看网站| 超碰成人久久| 人妻 亚洲 视频| 亚洲激情五月婷婷啪啪| 国产精品一区二区免费欧美 | 另类亚洲欧美激情| 如日韩欧美国产精品一区二区三区| 亚洲综合色网址| 国产精品久久久久久精品电影小说| 巨乳人妻的诱惑在线观看| 国产成人精品久久二区二区91| 一本色道久久久久久精品综合| 女人爽到高潮嗷嗷叫在线视频| av天堂在线播放| 欧美激情高清一区二区三区| 男女国产视频网站| 久久久久久免费高清国产稀缺| 精品人妻一区二区三区麻豆| 欧美+亚洲+日韩+国产| 亚洲专区中文字幕在线| 成年av动漫网址| 欧美精品高潮呻吟av久久| 热99久久久久精品小说推荐| 成在线人永久免费视频| 波多野结衣一区麻豆| av国产精品久久久久影院| 男人舔女人的私密视频| 国产免费福利视频在线观看| 欧美国产精品一级二级三级| 久久久久久久国产电影| 性色av乱码一区二区三区2| 2018国产大陆天天弄谢| 精品国产超薄肉色丝袜足j| 黄色视频不卡| 亚洲av男天堂| 亚洲自偷自拍图片 自拍| 久久狼人影院| 后天国语完整版免费观看| 免费一级毛片在线播放高清视频 | 99久久99久久久精品蜜桃| 天天影视国产精品| 亚洲精品日韩在线中文字幕| 国产片内射在线| 妹子高潮喷水视频| 欧美激情高清一区二区三区| 亚洲情色 制服丝袜| 欧美成人精品欧美一级黄| 99re6热这里在线精品视频| 美女主播在线视频| 精品少妇内射三级| 日韩电影二区| 欧美在线一区亚洲| 黄色毛片三级朝国网站| 精品一品国产午夜福利视频| 一区在线观看完整版| 亚洲国产欧美网| 丝袜喷水一区| 亚洲一区二区三区欧美精品| 久久精品国产综合久久久| 精品国产一区二区三区四区第35| 最新的欧美精品一区二区| 国产99久久九九免费精品| 久久精品国产a三级三级三级| 精品久久久久久久毛片微露脸 | 国产精品偷伦视频观看了| 亚洲欧美日韩高清在线视频 | a级毛片黄视频| 另类亚洲欧美激情| 别揉我奶头~嗯~啊~动态视频 | 亚洲中文字幕日韩| 色网站视频免费| 欧美日韩视频精品一区| 色播在线永久视频| 午夜福利乱码中文字幕| 婷婷色综合大香蕉| 国产老妇伦熟女老妇高清| 精品视频人人做人人爽| 一区二区三区精品91| 人人妻人人添人人爽欧美一区卜| 999精品在线视频| 国产av精品麻豆| 一级片'在线观看视频| 色婷婷久久久亚洲欧美| 免费看不卡的av| 大香蕉久久成人网| 9191精品国产免费久久| 欧美变态另类bdsm刘玥| 国产精品成人在线| 夜夜骑夜夜射夜夜干| 考比视频在线观看| 欧美日韩视频精品一区| 热99国产精品久久久久久7| 人妻 亚洲 视频| 欧美日韩视频精品一区| 老司机在亚洲福利影院| 少妇被粗大的猛进出69影院| 国产成人av教育| 在现免费观看毛片| 国产亚洲欧美在线一区二区| 亚洲欧洲精品一区二区精品久久久| 侵犯人妻中文字幕一二三四区| 国产精品久久久人人做人人爽| 亚洲国产日韩一区二区| 免费人妻精品一区二区三区视频| 啦啦啦在线观看免费高清www| 男人爽女人下面视频在线观看| 久久亚洲国产成人精品v| 免费观看a级毛片全部| 亚洲中文日韩欧美视频| 18禁国产床啪视频网站| 一区二区av电影网| 建设人人有责人人尽责人人享有的| 精品少妇黑人巨大在线播放| 可以免费在线观看a视频的电影网站| 久久久久国产精品人妻一区二区| www.av在线官网国产| 婷婷色综合大香蕉| 亚洲国产精品成人久久小说| 欧美激情极品国产一区二区三区| 大码成人一级视频| 亚洲专区中文字幕在线| 久久性视频一级片| 一边摸一边抽搐一进一出视频| 国产一卡二卡三卡精品| 自拍欧美九色日韩亚洲蝌蚪91| 天堂俺去俺来也www色官网| 天天添夜夜摸| 人人妻人人添人人爽欧美一区卜| 天天影视国产精品| 精品一区二区三区四区五区乱码 | 亚洲成色77777| 亚洲精品国产一区二区精华液| 老汉色∧v一级毛片| 久久99一区二区三区| 啦啦啦在线免费观看视频4| 日韩精品免费视频一区二区三区| 99香蕉大伊视频| 热99国产精品久久久久久7| 国产在线视频一区二区| 男女边吃奶边做爰视频| 在线观看一区二区三区激情| 亚洲欧洲国产日韩| 日韩 亚洲 欧美在线| 亚洲av日韩精品久久久久久密 | av线在线观看网站| 欧美在线一区亚洲| 最近手机中文字幕大全| 久久久欧美国产精品| 在线观看免费午夜福利视频| 亚洲中文字幕日韩| 久久免费观看电影| 2021少妇久久久久久久久久久| 国产不卡av网站在线观看| 国产精品亚洲av一区麻豆| 啦啦啦啦在线视频资源| 日韩av不卡免费在线播放| 日日摸夜夜添夜夜爱| 久久国产精品男人的天堂亚洲| 久久精品亚洲av国产电影网| 中文字幕高清在线视频| 亚洲欧美色中文字幕在线| 亚洲五月婷婷丁香| 久久ye,这里只有精品| 亚洲精品国产av蜜桃| 狠狠婷婷综合久久久久久88av| 久久午夜综合久久蜜桃| 久久精品熟女亚洲av麻豆精品| 一区在线观看完整版| a 毛片基地| 国产成人精品在线电影| 久久青草综合色| 午夜免费观看性视频| 精品熟女少妇八av免费久了| 亚洲av综合色区一区| 欧美精品人与动牲交sv欧美| av电影中文网址| 亚洲国产日韩一区二区| 丰满迷人的少妇在线观看| 欧美精品高潮呻吟av久久| 久久精品久久久久久噜噜老黄| 国产精品99久久99久久久不卡| videosex国产| 欧美日韩视频高清一区二区三区二| 丝袜喷水一区| 亚洲欧美成人综合另类久久久| 亚洲专区国产一区二区| 亚洲黑人精品在线| 欧美性长视频在线观看| 99精国产麻豆久久婷婷| 少妇猛男粗大的猛烈进出视频| 一区福利在线观看| 中文字幕人妻熟女乱码| 曰老女人黄片| 国产成人精品无人区| 精品少妇内射三级| 国产免费现黄频在线看| 国产免费视频播放在线视频| www.999成人在线观看| 免费久久久久久久精品成人欧美视频| 午夜两性在线视频| 日韩一本色道免费dvd| 超碰成人久久| 久久人妻熟女aⅴ| 国产成人a∨麻豆精品| av视频免费观看在线观看| 最新在线观看一区二区三区 | 天天操日日干夜夜撸| 99久久综合免费| 国产精品二区激情视频| 中国国产av一级| 91成人精品电影| 又大又黄又爽视频免费| 亚洲欧美激情在线| 亚洲激情五月婷婷啪啪| 我的亚洲天堂| 国产日韩欧美在线精品| 日韩一卡2卡3卡4卡2021年| 欧美国产精品va在线观看不卡| 日本91视频免费播放| 在线观看免费午夜福利视频| 中文字幕亚洲精品专区| 久久青草综合色| 久久久欧美国产精品| 久久久久视频综合| 无遮挡黄片免费观看| 久久久精品国产亚洲av高清涩受| 午夜视频精品福利| 亚洲人成电影免费在线| 少妇粗大呻吟视频| 精品国产乱码久久久久久小说| 亚洲精品第二区| 国产成人精品久久久久久| 亚洲中文日韩欧美视频| 午夜福利在线免费观看网站| 亚洲精品av麻豆狂野| 欧美激情 高清一区二区三区| 久久久精品国产亚洲av高清涩受| 久久久久久人人人人人| 欧美人与性动交α欧美精品济南到| 久久国产精品影院| 男女免费视频国产| 手机成人av网站| 亚洲欧美一区二区三区黑人| av不卡在线播放| 成人亚洲欧美一区二区av| 色视频在线一区二区三区| 亚洲精品一二三| 久久久久视频综合| 无遮挡黄片免费观看| 极品人妻少妇av视频| 伊人久久大香线蕉亚洲五| 69精品国产乱码久久久| 在线亚洲精品国产二区图片欧美| 婷婷色综合大香蕉| 男人舔女人的私密视频| 精品一区二区三区四区五区乱码 | 久久精品人人爽人人爽视色| 午夜福利视频精品| 九草在线视频观看| 高潮久久久久久久久久久不卡| 久久亚洲精品不卡| 亚洲国产中文字幕在线视频| 青春草视频在线免费观看| 亚洲伊人久久精品综合| 国产在线一区二区三区精| 波多野结衣一区麻豆| 久久久久久免费高清国产稀缺| 亚洲一码二码三码区别大吗| 97人妻天天添夜夜摸| 欧美精品人与动牲交sv欧美| 又紧又爽又黄一区二区| 交换朋友夫妻互换小说| tube8黄色片| 如日韩欧美国产精品一区二区三区| 久久久精品国产亚洲av高清涩受| 欧美日韩亚洲高清精品| 国产一区二区三区av在线| 美女国产高潮福利片在线看| 天天操日日干夜夜撸| 国产片特级美女逼逼视频| 丝袜美腿诱惑在线| 日本欧美国产在线视频| 母亲3免费完整高清在线观看| 少妇猛男粗大的猛烈进出视频| 精品少妇一区二区三区视频日本电影| 久久久久久久国产电影| 欧美成人精品欧美一级黄| 成年人免费黄色播放视频| 国产色视频综合| 九色亚洲精品在线播放| 香蕉丝袜av| 七月丁香在线播放| 大型av网站在线播放| 成年动漫av网址| 色精品久久人妻99蜜桃| 自拍欧美九色日韩亚洲蝌蚪91| 精品少妇久久久久久888优播| 午夜福利影视在线免费观看| 晚上一个人看的免费电影| 老汉色av国产亚洲站长工具| 亚洲精品久久久久久婷婷小说| 国产视频首页在线观看| 97精品久久久久久久久久精品| 亚洲精品一二三| 婷婷色麻豆天堂久久| 在线观看国产h片| www.自偷自拍.com| 日韩人妻精品一区2区三区| 久久热在线av| 午夜精品国产一区二区电影| 欧美日韩综合久久久久久| 日韩中文字幕视频在线看片| 午夜免费鲁丝| 女人久久www免费人成看片| 高潮久久久久久久久久久不卡| 成人午夜精彩视频在线观看| 午夜久久久在线观看| 日本欧美视频一区| 国产男女内射视频| 久久精品aⅴ一区二区三区四区| 性少妇av在线| 国产在线观看jvid| 午夜久久久在线观看| 欧美黑人精品巨大| 制服人妻中文乱码| 亚洲精品美女久久久久99蜜臀 | 亚洲国产看品久久| 亚洲精品国产av蜜桃| av欧美777| 极品人妻少妇av视频| 午夜两性在线视频| cao死你这个sao货| 色播在线永久视频| 午夜福利,免费看| 操美女的视频在线观看| 老司机靠b影院| 交换朋友夫妻互换小说| 美女福利国产在线| 久久亚洲国产成人精品v| 亚洲欧美成人综合另类久久久| 国产黄色视频一区二区在线观看| 岛国毛片在线播放| 少妇人妻 视频| 99热全是精品| 国产精品九九99| 日韩一卡2卡3卡4卡2021年| 一级a爱视频在线免费观看| 欧美日韩综合久久久久久| 在线观看免费高清a一片| 欧美激情极品国产一区二区三区| 国产成人精品久久二区二区免费| 国产有黄有色有爽视频| 男人爽女人下面视频在线观看| 久久99一区二区三区| a 毛片基地| 十八禁人妻一区二区| 亚洲精品国产区一区二| 久久这里只有精品19| 国产亚洲精品第一综合不卡| 又粗又硬又长又爽又黄的视频| 亚洲国产毛片av蜜桃av| √禁漫天堂资源中文www| 亚洲人成网站在线观看播放| 亚洲av电影在线观看一区二区三区| 欧美黄色片欧美黄色片| 美女福利国产在线| 高清黄色对白视频在线免费看| 亚洲精品日本国产第一区| 男女下面插进去视频免费观看| 我的亚洲天堂| 后天国语完整版免费观看| 日韩制服骚丝袜av| 免费不卡黄色视频| 国产老妇伦熟女老妇高清| 一区福利在线观看| 成人黄色视频免费在线看| www日本在线高清视频| 日韩,欧美,国产一区二区三区| 首页视频小说图片口味搜索 | 一区二区日韩欧美中文字幕| 国产麻豆69| 午夜视频精品福利| 母亲3免费完整高清在线观看| 91老司机精品| 午夜免费男女啪啪视频观看| 别揉我奶头~嗯~啊~动态视频 | 人人妻人人添人人爽欧美一区卜| 国产精品久久久av美女十八| 女警被强在线播放| 久久久精品区二区三区| 汤姆久久久久久久影院中文字幕| 欧美日韩黄片免| 日韩制服骚丝袜av| 国产精品麻豆人妻色哟哟久久| 久久久久国产一级毛片高清牌| 男女无遮挡免费网站观看| 亚洲精品国产色婷婷电影| 久久鲁丝午夜福利片| 啦啦啦在线免费观看视频4| 最近中文字幕2019免费版| 99国产精品一区二区三区| 1024视频免费在线观看| 亚洲 国产 在线| 少妇的丰满在线观看| 国产午夜精品一二区理论片| 脱女人内裤的视频| 中文字幕另类日韩欧美亚洲嫩草| 精品福利观看| 亚洲一卡2卡3卡4卡5卡精品中文| 午夜免费成人在线视频| 一级毛片电影观看| 久久久久精品人妻al黑| 亚洲欧美清纯卡通| 丝袜人妻中文字幕| 乱人伦中国视频| 美女扒开内裤让男人捅视频| 看免费成人av毛片| 久久久久网色| 91国产中文字幕| 男女午夜视频在线观看| svipshipincom国产片| 亚洲av美国av| 丰满少妇做爰视频| 老鸭窝网址在线观看| 欧美成人精品欧美一级黄| 女人精品久久久久毛片| 1024香蕉在线观看| 好男人视频免费观看在线| 欧美日韩综合久久久久久| 中文字幕精品免费在线观看视频| 亚洲国产欧美一区二区综合| 每晚都被弄得嗷嗷叫到高潮| 中文字幕色久视频| 欧美老熟妇乱子伦牲交| 咕卡用的链子| 国产成人影院久久av| av线在线观看网站| 日韩av免费高清视频| 中文字幕人妻丝袜一区二区| 午夜福利影视在线免费观看| 久久99精品国语久久久| 久久久久久久大尺度免费视频| 亚洲成国产人片在线观看| 欧美老熟妇乱子伦牲交| av电影中文网址| 大香蕉久久成人网| 国产亚洲精品第一综合不卡| 黑人欧美特级aaaaaa片| 国产熟女欧美一区二区| 另类精品久久| 色综合欧美亚洲国产小说| 日韩,欧美,国产一区二区三区| 夜夜骑夜夜射夜夜干| 国产精品国产三级专区第一集| 天天躁夜夜躁狠狠躁躁| 国产精品免费大片| 欧美激情极品国产一区二区三区| 国产在视频线精品| 日韩大码丰满熟妇| 男女午夜视频在线观看| 青青草视频在线视频观看| 99国产综合亚洲精品| 久久人人爽人人片av| 精品国产乱码久久久久久小说| 天天躁狠狠躁夜夜躁狠狠躁| 精品一区在线观看国产| 亚洲成色77777| 777久久人妻少妇嫩草av网站| 国产熟女午夜一区二区三区| 日韩 欧美 亚洲 中文字幕| 啦啦啦中文免费视频观看日本| 成人亚洲欧美一区二区av| 99国产精品一区二区三区| 日本午夜av视频| 少妇精品久久久久久久| 久久精品久久精品一区二区三区| 国产伦人伦偷精品视频| 国产一级毛片在线| 大码成人一级视频| 国产一区二区三区综合在线观看| 无限看片的www在线观看| 青草久久国产| 午夜影院在线不卡| 久久 成人 亚洲| 久久精品成人免费网站| 午夜免费成人在线视频| 人人澡人人妻人| 大型av网站在线播放| www.自偷自拍.com| 精品第一国产精品| 性色av乱码一区二区三区2| 国产成人欧美| 久久久国产一区二区| 久久精品熟女亚洲av麻豆精品| 看免费av毛片| 美女大奶头黄色视频| 自线自在国产av| 中文字幕最新亚洲高清| 丝袜美腿诱惑在线| 悠悠久久av| 国产激情久久老熟女| 久久国产精品大桥未久av| 亚洲成人免费电影在线观看 | 国产精品香港三级国产av潘金莲 | 最近最新中文字幕大全免费视频 | 精品欧美一区二区三区在线| 美女高潮到喷水免费观看| 久久人人爽人人片av| av片东京热男人的天堂| 高清视频免费观看一区二区| 美女国产高潮福利片在线看| 国产一区二区在线观看av| 欧美黑人欧美精品刺激| 我要看黄色一级片免费的| 成人黄色视频免费在线看| 成人午夜精彩视频在线观看| 好男人视频免费观看在线| 十八禁人妻一区二区| 日韩视频在线欧美| 黄色a级毛片大全视频| 人成视频在线观看免费观看| 日韩制服丝袜自拍偷拍| 国产日韩一区二区三区精品不卡| 色94色欧美一区二区| 欧美日韩黄片免| 欧美国产精品一级二级三级| 少妇裸体淫交视频免费看高清 |