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

    載荷方向不確定條件下結(jié)構(gòu)動態(tài)穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)1)

    2023-12-16 11:48:40杜鼎新
    力學(xué)學(xué)報 2023年11期
    關(guān)鍵詞:穩(wěn)健性確定性標(biāo)準(zhǔn)差

    杜鼎新 王 棟

    (西北工業(yè)大學(xué)航空學(xué)院,西安 710072)

    引言

    傳統(tǒng)的結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)往往是在理想狀態(tài)下進(jìn)行的,即結(jié)構(gòu)所受的載荷環(huán)境、材料參數(shù)以及約束條件等均被認(rèn)為是確定性的.如此所獲得的結(jié)構(gòu)優(yōu)化設(shè)計(jì),也僅僅是針對某個(些)確定的名義設(shè)計(jì)狀態(tài)而言是最優(yōu)的[1-7].然而在實(shí)際工程結(jié)構(gòu)中,不確定性因素卻是普遍存在且無法避免的,如材料的物理性能、模型幾何尺寸、邊界約束條件和外界施加的載荷等[8-13].這些結(jié)構(gòu)設(shè)計(jì)參數(shù)會由于工業(yè)生產(chǎn)的誤差、工作環(huán)境的改變以及長時間服役引起性能的退化等因素而產(chǎn)生變化,導(dǎo)致結(jié)構(gòu)的實(shí)際工作狀態(tài)不再是理想和固定的,而是可變和不確定的.即使這些參數(shù)的變化范圍很小,無疑也會對結(jié)構(gòu)的承力性能產(chǎn)生巨大的作用,最終可能對結(jié)構(gòu)的響應(yīng)造成不利的影響.這就意味著,如果在結(jié)構(gòu)設(shè)計(jì)時僅考慮理想的確定性情形,那么所得到的系統(tǒng)就缺少必要的結(jié)構(gòu)備份或冗余,在實(shí)際應(yīng)用中缺乏足夠抵抗材料參數(shù)或載荷等不確定性因素變化的能力,其結(jié)構(gòu)功能甚至可能因?yàn)檫@些因素的變化而提前失效[13].因此,為了使所設(shè)計(jì)的結(jié)構(gòu)具有較強(qiáng)的抵抗不確定性因素的能力,設(shè)計(jì)者需要在其初步設(shè)計(jì)階段就主動地考慮它們的影響,而非僅關(guān)注少數(shù)幾個離散的理想設(shè)計(jì)狀態(tài),從而使結(jié)構(gòu)具有足夠的安全冗余,容忍工作狀態(tài)的可變性,且在實(shí)際應(yīng)用中能正常、穩(wěn)定地發(fā)揮預(yù)期的功能和效果[11-16].

    相較于結(jié)構(gòu)的材料屬性、邊界約束條件等因素,外載荷狀態(tài)對結(jié)構(gòu)拓?fù)錁?gòu)型設(shè)計(jì)的影響更大.實(shí)際工程中,由于外部環(huán)境的變化,作用于結(jié)構(gòu)上外載荷的大小、方向與位置等都可能具有一定范圍內(nèi)的隨機(jī)可變性[11,16-19].因此研究載荷不確定性對結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)的影響具有非?,F(xiàn)實(shí)的工程意義.另外,大多數(shù)工程結(jié)構(gòu),如航空航天器、船舶、橋梁等,其結(jié)構(gòu)在服役期內(nèi)會受到各種時變載荷的作用,而一個復(fù)雜的時變周期激勵通常能近似分解為幾個簡諧激勵的疊加.因此簡諧激勵是一種非?;A(chǔ)且典型的動載荷[6],研究簡諧激勵作用下結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)具有十分重要的理論和實(shí)際意義[17-18].

    研究外載荷不確定性條件下結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)可采用基于可靠性的拓?fù)鋬?yōu)化方法(reliability-based topology optimization,RBTO)[20-24],或基于穩(wěn)健性的拓?fù)鋬?yōu)化方法(robust topology optimization,RTO)[9,11,16-19,24-30].兩者均可使用概率或者非概率模型表示不確定變量的變化情況.不同之處在于前者旨在滿足可靠性約束條件,而后者聚焦于減少結(jié)構(gòu)性能對不確定性變量的敏感性.在穩(wěn)健性拓?fù)鋬?yōu)化的體系下,已經(jīng)有許多關(guān)于外載荷作用方向的不確定性的研究報告[11,16-17,25-28,30].通常,結(jié)構(gòu)的動柔順度響應(yīng)和外激勵作用方向的變化之間并沒有一個顯性的關(guān)系式[16,28],因此難以直接分析出外激勵方向擾動對結(jié)構(gòu)動柔順度的影響.本文旨在構(gòu)建一種簡單且高效的計(jì)算方法,可直觀、快速地獲得結(jié)構(gòu)動柔順度的概率特征指標(biāo).并基于這些特征指標(biāo),開展了考慮外載荷作用方向不確定性的穩(wěn)健性結(jié)構(gòu)動態(tài)拓?fù)鋬?yōu)化設(shè)計(jì).

    在實(shí)際工程中,外載荷作用狀態(tài)的小幅變化往往近似地符合正態(tài)的概率分布,即圍繞著預(yù)先設(shè)定的名義狀態(tài)有著一定概率性的變化規(guī)律,因此可使用正態(tài)隨機(jī)分布對外載荷作用狀態(tài)的可變性進(jìn)行描述[11,25-26].而且相應(yīng)的結(jié)構(gòu)動柔順度也不再是一個確定的值,但兩者之間也并非簡單的線性關(guān)系.為了高效地計(jì)算在外載荷作用方向變化情形下結(jié)構(gòu)的動柔順度響應(yīng),本文首先將外載荷沿坐標(biāo)軸進(jìn)行投影,然后對結(jié)構(gòu)的動柔順度響應(yīng)進(jìn)行二階泰勒展開,以便獲得動柔順度的近似表達(dá)式.并近一步推導(dǎo)出結(jié)構(gòu)動柔順度的期望和方差的計(jì)算表達(dá)式,以及它們對拓?fù)渥兞快`敏度的顯性計(jì)算公式.按照確定性的優(yōu)化策略,可完成結(jié)構(gòu)在動載荷作用下的穩(wěn)健性拓?fù)湓O(shè)計(jì).保證在外激勵作用方向出現(xiàn)擾動時,結(jié)構(gòu)的動柔順度響應(yīng)不發(fā)生較大的改變.用結(jié)構(gòu)動柔順度的變異系數(shù)(coefficient of variation,COV)[13]作為一種量化的穩(wěn)健度衡量指標(biāo),對比分析了不同設(shè)計(jì)策略結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)的穩(wěn)健程度.

    本文采用概率模型描述外載荷作用方向的不確定性,基于材料屬性的有理近似方法 (rational approximation of material properties,RAMP)建立了連續(xù)體結(jié)構(gòu)動態(tài)拓?fù)鋬?yōu)化模型.以結(jié)構(gòu)動柔順度的期望值及標(biāo)準(zhǔn)差(standard deviation,STD)的加權(quán)和最小化為設(shè)計(jì)目標(biāo)函數(shù)[11,16,25-27],以整體材料用量為約束條件,采用基于靈敏度分析的移動漸近方法(method of moving asymptotes,MMA)[31],研究簡諧激勵在作用方向不確定性條件下結(jié)構(gòu)受迫振動穩(wěn)態(tài)響應(yīng)的拓?fù)鋬?yōu)化問題.最后用兩個典型算例來證明本文所提方法的有效性,并與傳統(tǒng)的優(yōu)化結(jié)果進(jìn)行對比,以便驗(yàn)證穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)結(jié)果的可靠性.

    1 結(jié)構(gòu)穩(wěn)健性拓?fù)鋬?yōu)化模型

    1.1 線性振動系統(tǒng)穩(wěn)態(tài)響應(yīng)分析

    當(dāng)外激勵頻率遠(yuǎn)離結(jié)構(gòu)的固有頻率時,可以暫時忽略結(jié)構(gòu)內(nèi)的阻尼.采用有限元法對連續(xù)體結(jié)構(gòu)進(jìn)行離散化處理,則無阻尼結(jié)構(gòu)受迫振動微分方程可表示為[17]

    式中,K和M分別是振動系統(tǒng)N×N階的剛度和質(zhì)量矩陣,N表示結(jié)構(gòu)的自由度數(shù);u(t)和分別是結(jié)構(gòu)的位移和加速度響應(yīng)N維列陣;f(t)是系統(tǒng)所受外激勵列陣.假設(shè)結(jié)構(gòu)受到簡諧外力的作用

    其中,ω是外激勵頻率((°)/s),是一個給定的值;F(θ)是沿坐標(biāo)軸方向上的外激勵幅值N維列陣

    其中,奇(偶)數(shù)項(xiàng)表示在有限元結(jié)點(diǎn)處,外激勵沿x(y)軸的幅值分量;θ表示外激勵的作用方向,是一個隨機(jī)變量.假設(shè)其變化符合正態(tài)分布,即見圖1 所示.其中μθ代表外載荷作用的期望或名義方向,σθ是其標(biāo)準(zhǔn)差,代表著外載荷方向變化的概率統(tǒng)計(jì)偏差.嚴(yán)格地說,一個服從正態(tài)分布的隨機(jī)變量是無界的.然而按照工程實(shí)踐中廣泛使用的“3σ”法則,即對于特定的符合正態(tài)分布的隨機(jī)變量而言,它取值落在區(qū)間[μ-3σ,μ+3σ]之外的概率小于千分之三,因此便可以認(rèn)為區(qū)間[μ-3σ,μ+3σ]為該隨機(jī)變量實(shí)際可能的取值區(qū)間.依據(jù)該法則,σθ的值可以通過區(qū)間寬度基本確定.由于載荷作用方向θ具有可變性,因此外激勵幅值列陣F(θ)也具有一定的不確定性.

    圖1 結(jié)構(gòu)受具有方向不確定外載荷作用的示意圖Fig.1 Schematic of a load with the direction uncertainty imposed to a structure

    對方程式(1)兩邊同時進(jìn)行Fourier 變換,則可得振動系統(tǒng)響應(yīng)的特征平衡方程

    式中,D(ω)=K-ω2M是振動系統(tǒng)的動剛度矩陣.U(ω)是結(jié)構(gòu)位移響應(yīng)幅值列陣,在不考慮結(jié)構(gòu)阻尼的前提下是實(shí)數(shù),且與結(jié)構(gòu)的動態(tài)性能以及激振頻率密切相關(guān).本文采用直接方法精確計(jì)算結(jié)構(gòu)在簡諧激勵作用下的穩(wěn)態(tài)響應(yīng),不必考慮模態(tài)截尾帶來的計(jì)算誤差和虛假模態(tài)的干擾.求解方程(4)可得結(jié)構(gòu)響應(yīng)的幅值列陣U(ω),進(jìn)而由Fourier 逆變換可得時域內(nèi)結(jié)構(gòu)的穩(wěn)態(tài)動響應(yīng)u(t).

    1.2 結(jié)構(gòu)拓?fù)鋬?yōu)化數(shù)學(xué)模型

    實(shí)際工程中,由于結(jié)構(gòu)所處環(huán)境的變化,外載荷作用方向具有一定的可變性.因此結(jié)構(gòu)的響應(yīng)也具有一定的不確定性.本文以結(jié)構(gòu)有限元模型中單元的相對密度ρe作為設(shè)計(jì)變量,在給定的區(qū)域內(nèi)構(gòu)造結(jié)構(gòu)的最優(yōu)拓?fù)錁?gòu)型(材料布局).考慮到動柔順度的不確定性,將設(shè)計(jì)目標(biāo)函數(shù)設(shè)置為最小化結(jié)構(gòu)動柔順度均值和標(biāo)準(zhǔn)差的加權(quán)和[11,16,25-27],則結(jié)構(gòu)動態(tài)穩(wěn)健性拓?fù)鋬?yōu)化問題數(shù)學(xué)列式可表示為

    式中,Cd(θ)是結(jié)構(gòu)的動柔順度,它不僅與外激勵的大小有關(guān),還與激勵作用方向θ密切相關(guān),即隨著θ的變化而改變.因此,Cd(θ) 也是一個隨機(jī)變量.J為優(yōu)化設(shè)計(jì)的目標(biāo)函數(shù),是結(jié)構(gòu)動柔順度的期望,是結(jié)構(gòu)動柔順度的標(biāo)準(zhǔn)差,β是權(quán)重因子,其作用是調(diào)節(jié)結(jié)構(gòu)動柔順度標(biāo)準(zhǔn)差在目標(biāo)函數(shù)中的占比.ρ是n維拓?fù)渥兞苛嘘?n代表結(jié)構(gòu)的單元總數(shù).ρmin是為了避免結(jié)構(gòu)剛度矩陣奇異而設(shè)置的一個下限值(ρmin=0.001).ve是單元的體積,fv是給定的材料體分比,即材料用量的上限.

    以上模型與傳統(tǒng)的確定性拓?fù)鋬?yōu)化模型最大區(qū)別在于:結(jié)構(gòu)的動柔順度Cd不再是一個確切的量,而是隨機(jī)變量θ的函數(shù).且優(yōu)化的目標(biāo)函數(shù)不再僅僅針對名義設(shè)計(jì)狀態(tài),而是綜合考量在θ的隨機(jī)變化區(qū)間內(nèi)的整體設(shè)計(jì)狀態(tài),旨在讓最終的優(yōu)化設(shè)計(jì)結(jié)果具有能容忍外激勵方向變化的冗余性.這使得優(yōu)化設(shè)計(jì)過程更加復(fù)雜,但也更加真實(shí)地反映了工程結(jié)構(gòu)的實(shí)際受載情況.

    為了合理地描述結(jié)構(gòu)的動力學(xué)性能,根據(jù)以往的研究結(jié)果,本文采用RAMP 模型對單元質(zhì)量和剛度矩陣進(jìn)行插值,以避免虛假局部模態(tài)的產(chǎn)生[17]

    式中,q是插值模型的懲罰因子(q=5);m0和k0分別是實(shí)體材料單元的質(zhì)量和剛度矩陣.

    1.3 結(jié)構(gòu)動柔順度一階導(dǎo)數(shù)靈敏度的計(jì)算

    當(dāng)外激勵作用方向θ按照正態(tài)分布隨機(jī)擾動時,根據(jù)文獻(xiàn)[16,28]的定義,只需在外激勵的名義(理想)作用方向,即

    時,對振動系統(tǒng)進(jìn)行分析,計(jì)算結(jié)構(gòu)的動態(tài)響應(yīng).并在單元層面上計(jì)算結(jié)構(gòu)的動柔順度Cd0

    于是,其對單元的相對密度ρe的一階導(dǎo)數(shù)靈敏度可由下式計(jì)算得到

    其中,ue是單元的動響應(yīng)幅值列陣;de=ke-ω2me是單元的動剛度矩陣.由RAMP 模型的材料插值函數(shù)式(8)可得

    于是,通過有限元分析計(jì)算,即可獲得振動結(jié)構(gòu)的名義動柔順度Cd0及其對拓?fù)湓O(shè)計(jì)變量的一階導(dǎo)數(shù)靈敏度值.但當(dāng)外激勵的作用方向θ發(fā)生擾動時,如果依然通過傳統(tǒng)有限元方法進(jìn)行分析和計(jì)算,就需要再次求解特征平衡方程式(4).這樣的計(jì)算途徑工作效率很低,尤其考慮到外激勵的作用方向可以連續(xù)地隨機(jī)變化.因此需要一種能高效獲取在外激勵方向發(fā)生變化時結(jié)構(gòu)動柔順度的計(jì)算途徑.

    2 外激勵方向不確定性的影響規(guī)律

    2.1 載荷方向擾動對結(jié)構(gòu)動柔順度的影響

    假設(shè)外激勵的大小和位置保持不變,當(dāng)載荷作用方向發(fā)生改變時,其傳遞路徑亦將發(fā)生改變.即外激勵作用方向的變化必然引起結(jié)構(gòu)動柔順度的變化.利用二階泰勒級數(shù)展開公式,得到動柔順度Cd(θ)變化的顯性表達(dá)式,可避免振動系統(tǒng)動響應(yīng)方程的重復(fù)求解.

    在載荷位置固定的前提下,以載荷名義作用方向?yàn)橹行牡囊粋€較小的區(qū)域內(nèi),可得載荷方向發(fā)生改變后,結(jié)構(gòu)動柔順度Cd的近似表達(dá)式[16,28]

    按照式(3)對外激勵幅值F(θ)的表達(dá)形式,首先推導(dǎo)其對θ的一階和二階導(dǎo)數(shù)

    將結(jié)構(gòu)的名義動柔順度改寫成

    則動柔順度對載荷作用方向的一階導(dǎo)數(shù)可表示為

    據(jù)此,可進(jìn)一步計(jì)算結(jié)構(gòu)的動柔順度Cd0對載荷作用方向的二階導(dǎo)數(shù)

    根據(jù)以上推導(dǎo)結(jié)果,將式(10)、式(18) 和式(19)代入式(13),即可利用該表達(dá)式快速計(jì)算由于外激勵作用方向的擾動而引起結(jié)構(gòu)的動柔順度值.

    2.2 目標(biāo)函數(shù)及其靈敏度的計(jì)算

    雖然前文已經(jīng)假設(shè)外激勵作用方向的變化符合隨機(jī)正態(tài)分布,但結(jié)構(gòu)的動柔順度響應(yīng)變化一般不符合正態(tài)隨機(jī)分布.按照優(yōu)化設(shè)計(jì)目標(biāo)的要求,需要高效地計(jì)算結(jié)構(gòu)動柔順度響應(yīng)的均值和方差.對式(13)兩邊同時分別求期望和方差可得[16]

    根據(jù)θ的正態(tài)隨機(jī)分布規(guī)律,式(21)中的各協(xié)方差項(xiàng)均等于0[32].因此,結(jié)構(gòu)動柔順度響應(yīng)的均值和標(biāo)準(zhǔn)差表達(dá)式如下

    于是,結(jié)構(gòu)動柔順度響應(yīng)的變異系數(shù)為

    巖樣采用薄壁金剛石鉆頭沿垂直于巖層方向鉆取巖芯,經(jīng)過鋸、磨加工成直徑為50 mm,高為100 mm,試樣兩端面不平行度不大于0.05 mm,滿足《規(guī)程》要求,60組共計(jì)180個試樣,制備部分試樣,如圖2所示,其中頂板巖石用A編號,底板巖石用B編號。試驗(yàn)在RMT-150B型電液伺服巖石力學(xué)試驗(yàn)系統(tǒng)進(jìn)行單軸壓縮試驗(yàn),如圖3所示。軸向荷載采用1 000 kN力傳感器測量,軸向壓縮變形采用5.0 mm位移傳感器測量,變形精度為1.0×10-3 mm,采用位移控制方式,加載速率為0.002 mm/s,每組巖性重復(fù)進(jìn)行3次試驗(yàn)。煤層頂?shù)装鍘r石單軸壓縮試驗(yàn)結(jié)果見表1。

    顯然,式(22)和式(23)建立了結(jié)構(gòu)動柔順度響應(yīng)Cd(θ)的概率特征值顯式表達(dá)式.無需再通過其他額外的途徑,如Monte Carlo 模擬(Monte Carlo simulation,MCS)法等[11,25,30],即可獲得結(jié)構(gòu)動柔順度響應(yīng)的一、二階矩以及優(yōu)化設(shè)計(jì)的目標(biāo)函數(shù)

    這使得穩(wěn)健性結(jié)構(gòu)動態(tài)優(yōu)化設(shè)計(jì)效率大為提高.

    將式(25)對拓?fù)湓O(shè)計(jì)變量求導(dǎo),可得

    其中,式(26)右邊第一項(xiàng)由式(11)得到,其余各項(xiàng)可由式(11)分別對外載荷作用方向θ連續(xù)微分得到

    式中,單元位移列陣對θ的一階、二階導(dǎo)數(shù)可利用式(16)計(jì)算得到.另外,式(27)和式(28)同樣也可以分別利用式(18)和式(19)在單元層面上對ρe求一階導(dǎo)數(shù)獲得.根據(jù)以上推導(dǎo)結(jié)果,即可利用式(26)直接計(jì)算目標(biāo)函數(shù)對所有設(shè)計(jì)變量的一階靈敏度.在此基礎(chǔ)上,使用基于靈敏度的MMA 算法對設(shè)計(jì)變量進(jìn)行迭代更新,實(shí)現(xiàn)結(jié)構(gòu)拓?fù)錁?gòu)型的演化.此外,還要引入靈敏度過濾技術(shù),對動柔順度靈敏度進(jìn)行處理,以便消除優(yōu)化過程中的棋盤格現(xiàn)象,確保優(yōu)化結(jié)果的可行性和可靠性[16-19].

    3 算例分析

    下面我們用兩個算例來驗(yàn)證本文所提穩(wěn)健性拓?fù)鋬?yōu)化算法的可行性,并與確定性拓?fù)鋬?yōu)化(deterministic topology optimization,DTO)設(shè)計(jì)結(jié)果進(jìn)行對比分析,觀測結(jié)構(gòu)穩(wěn)健性構(gòu)型設(shè)計(jì)的變化情況,同時探尋結(jié)構(gòu)是如何通過增加額外構(gòu)件而有效提高其抵抗外載荷方向擾動能力的,進(jìn)而證明所得結(jié)果的合理性.假設(shè)材料的彈性模量E=200 GPa,密度ρ=7800 kg/m3,泊松比ν=0.3.有限元分析采用四結(jié)點(diǎn)平面應(yīng)力單元,薄板厚度為10 mm.優(yōu)化過程中取過濾半徑為單元尺寸的1.5 倍,收斂條件設(shè)置為相鄰兩次迭代各設(shè)計(jì)變量的最大改變量小于0.01[31].材料體積約束系數(shù)fv=0.35,初始設(shè)計(jì)時,將所給的材料均勻分布在整個設(shè)計(jì)區(qū)域內(nèi).

    3.1 簡支平板結(jié)構(gòu)

    簡支平板結(jié)構(gòu)設(shè)計(jì)區(qū)域尺寸及邊界約束如圖2所示.理論上,結(jié)構(gòu)在底邊中點(diǎn)承受一個集中簡諧激振力f(t)=10ejωtkN 的作用,頻率為150 Hz (ω=2π×150 (°)/s).但是由于加載環(huán)境的復(fù)雜性,外激振實(shí)際作用方向?qū)@著作用點(diǎn)產(chǎn)生一定的擾動.將設(shè)計(jì)區(qū)域均勻劃分成120×40 的有限元網(wǎng)格.初始設(shè)計(jì),結(jié)構(gòu)的第一階固有頻率為265.63 Hz.假設(shè)外激勵作用方向θ的標(biāo)準(zhǔn)差為 σθ=π/12.目標(biāo)函數(shù)中的權(quán)重因子β=2.

    圖2 簡支平面矩形板結(jié)構(gòu)Fig.2 A simply supported rectangular panel

    在開展結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)前,先對本文所提有關(guān)結(jié)構(gòu)動柔順度均值和方差估計(jì)公式的準(zhǔn)確性進(jìn)行驗(yàn)證.分別利用MCS (樣本數(shù)=105)以及式(22)和式(23)進(jìn)行計(jì)算,數(shù)值結(jié)果如表1 所示.可以看出,兩種方法的結(jié)果非常接近,尤其是均值的計(jì)算結(jié)果,這足以驗(yàn)證了式(22)和式(23)的有效性和準(zhǔn)確性.此外,用本文提出的方法計(jì)算結(jié)構(gòu)動柔順度均值和方差可以大幅降低計(jì)算機(jī)運(yùn)算所需的時間,顯著提高優(yōu)化設(shè)計(jì)的效率.

    表1 兩種方法所得結(jié)構(gòu)動柔順度概率特征指標(biāo)對比Table 1 Comparison of the obtained stochastic measurements of the structural dynamic compliance

    首先開展結(jié)構(gòu)的確定性拓?fù)鋬?yōu)化設(shè)計(jì),其最優(yōu)構(gòu)型如圖3(a)所示.可見,垂直方向的外載荷先通過中間的兩根斜桿(或構(gòu)件)傳遞到上緣,然后再通過外側(cè)的構(gòu)件傳遞到支撐點(diǎn)(基礎(chǔ)).由于載荷作用方向是固定的,沒有水平方向的分力,因此結(jié)構(gòu)設(shè)計(jì)區(qū)域的下緣沒有水平桿件將載荷作用點(diǎn)與支撐點(diǎn)直接連接在一起.

    圖3 確定性和穩(wěn)健性設(shè)計(jì)策略結(jié)構(gòu)拓?fù)錁?gòu)型對比Fig.3 Comparison of the topology optimizations on two different design strategies of DTO and RTO

    隨后,開展結(jié)構(gòu)穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì),所得最優(yōu)構(gòu)型設(shè)計(jì)如圖3(b)所示.兩種優(yōu)化策略得到的結(jié)構(gòu)構(gòu)型看似基本相同,但是穩(wěn)健性優(yōu)化設(shè)計(jì)在結(jié)構(gòu)的下緣增加了一根水平構(gòu)件,從而形成了額外的兩條傳力路徑.當(dāng)載荷作用方向發(fā)生擾動,不再是僅沿垂向激勵時,外激勵的水平分量將通過下緣桿件直接傳遞到支撐點(diǎn).由此可知,增加的水平桿件可有效抵抗外激勵方向的擾動,能顯著提高結(jié)構(gòu)在載荷作用方向不確定性情形下響應(yīng)的穩(wěn)健性.

    表2 比較了兩種設(shè)計(jì)策略所得結(jié)構(gòu)拓?fù)鋬?yōu)化數(shù)值結(jié)果.雖然在確定性拓?fù)鋬?yōu)化時沒有考慮外激勵作用方法的可變性,但仍可用第2 節(jié)中的方法計(jì)算其隨激勵方向隨機(jī)擾動時結(jié)構(gòu)動柔順度的概率特征值.從表2 中的數(shù)值對比可知,按照穩(wěn)健性拓?fù)鋬?yōu)化策略所得結(jié)構(gòu)動柔順度的期望、標(biāo)準(zhǔn)差以及變異系數(shù),比確定性拓?fù)鋬?yōu)化相應(yīng)結(jié)果均有明顯降低,尤其是標(biāo)準(zhǔn)差以及變異系數(shù)下降幅度均超過了95%.眾所周知,這兩個特征指標(biāo)能衡量不確定性數(shù)值的集中程度[32],其值越低,表明因外激勵方向變化所引起的結(jié)構(gòu)動柔順度響應(yīng)的變化幅度越小.即結(jié)構(gòu)對外激勵作用方向擾動不再敏感,結(jié)構(gòu)的穩(wěn)健度更高,材料布局設(shè)計(jì)更加合理、可靠.

    表2 兩種設(shè)計(jì)策略結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果對比Table 2 Comparison of the numerical results with the two optimization strategies

    圖4 示出了結(jié)構(gòu)的動柔順度隨外激勵作用方向的變化情況.值得注意的是,當(dāng)外載荷作用方向受到擾動以后,確定性設(shè)計(jì)結(jié)構(gòu)的動柔順度將快速增加,即結(jié)構(gòu)的動剛度明顯降低.這表明了該結(jié)構(gòu)抵抗外激勵方向擾動的能力較差.相比之下,穩(wěn)健性優(yōu)化設(shè)計(jì)結(jié)構(gòu)的動柔順度在激勵方向擾動狀態(tài)下能基本保持一致,且變化趨勢隨角度的增大而逐漸減小,這表明了穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)能避免結(jié)構(gòu)的動柔順度(剛度)發(fā)生劇烈變化,有效保持結(jié)構(gòu)動態(tài)性能的穩(wěn)定性.

    圖4 結(jié)構(gòu)動柔順度隨外激勵作用方向變化情況Fig.4 Variation of the structural dynamic compliance caused by the external loading direction disturbances

    然而,由圖4 還可以發(fā)現(xiàn),穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)的名義最小動柔順度Cd0比確定性相應(yīng)值稍大一些(15.30%),見表2 所示.進(jìn)一步比較圖3 可知,穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)中間兩根斜撐桿的截面積,較確定性的優(yōu)化設(shè)計(jì)要細(xì)一些,從而導(dǎo)致其名義最小動柔順度有所增大.這是因?yàn)闉榱颂岣呓Y(jié)構(gòu)抵抗外激勵方向的擾動能力,有一部分材料用來構(gòu)建額外的水平承力構(gòu)件,以便增加結(jié)構(gòu)對外載荷傳遞的可靠性.因此在相同材料體分比條件下,穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)結(jié)構(gòu)的名義動剛度會有一定的損失.但是這部分損失可以通過適當(dāng)?shù)卦黾咏Y(jié)構(gòu)材料得到有效補(bǔ)償[4].

    實(shí)際上,對于確定性拓?fù)鋬?yōu)化而言,其優(yōu)化設(shè)計(jì)目標(biāo)僅僅是名義加載狀態(tài)下結(jié)構(gòu)的動柔順度最小化,而忽略其變化情況,因此該優(yōu)化策略對應(yīng)的名義動柔順度值一般會更小.而穩(wěn)健性拓?fù)鋬?yōu)化要同時考慮在激勵方向隨機(jī)變化情形下,結(jié)構(gòu)的動柔順度值及其變化和分散情況.因此穩(wěn)健性設(shè)計(jì)所得動柔順度變化曲線比確定性設(shè)計(jì)對應(yīng)的曲線變化更加平緩,在載荷擾動狀態(tài)下的結(jié)構(gòu)動柔順度數(shù)值更加穩(wěn)定.由于受材料體積的制約,穩(wěn)健性設(shè)計(jì)的名義動柔順度值難免會有所上升,而優(yōu)化設(shè)計(jì)過程就是要在相互妥協(xié)中找到矛盾雙方最佳的解決方案[9].

    在本例中,只要將結(jié)構(gòu)材料體分比增大到fv=0.358(2.29%),就可以將穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)的名義最小動柔順度降到與確定性優(yōu)化設(shè)計(jì)相同的結(jié)果,并且結(jié)構(gòu)的動柔順度變化情況基本保持不變,見圖4 所示.即在少許增大結(jié)構(gòu)材料用量的基礎(chǔ)上,結(jié)構(gòu)的動柔順度響應(yīng)將全面優(yōu)(小)于確定性優(yōu)化設(shè)計(jì)的結(jié)果.這再次表明采用本文提出的穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)后,結(jié)構(gòu)抵抗外載荷作用方向擾動的能力有了顯著地增強(qiáng).

    比較結(jié)構(gòu)的前3 階固有頻率可以發(fā)現(xiàn),結(jié)構(gòu)的低階固有頻率,經(jīng)過拓?fù)鋬?yōu)化后都得到了明顯提高,如表3 所示.在當(dāng)前激勵情形下,結(jié)構(gòu)的總體動剛度較初始設(shè)計(jì)有了顯著增大,從而使結(jié)構(gòu)的動響應(yīng)幅值急劇減小.此外,在相同材料體分比約束條件下,穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)的結(jié)構(gòu)固有頻率要略小于確定性優(yōu)化的結(jié)果,這也是為了能夠有效提高結(jié)構(gòu)響應(yīng)的穩(wěn)健性,傳力路徑模式更新的結(jié)果.

    表3 簡支平板結(jié)構(gòu)前3 階固有頻率比較Table 3 Comparison of the first three natural frequencies of the rectangular panel structure

    此外,還有一點(diǎn)值得注意的是,隨著優(yōu)化過程的持續(xù)運(yùn)行,結(jié)構(gòu)的固有頻率也會不斷提升,這就使得優(yōu)化后結(jié)構(gòu)的第一階固有頻率和外激勵的頻率有了較大差距,從而使得動力學(xué)優(yōu)化設(shè)計(jì)問題越來越接近于靜力學(xué)優(yōu)化設(shè)計(jì)問題.為了驗(yàn)證該方法在高頻外激勵作用下的可行性,針對同一設(shè)計(jì)域,在外激振頻率為300 Hz,材料體積約束系數(shù)fv=0.55 (初始設(shè)計(jì)時,結(jié)構(gòu)的第一階固有頻率為262.72 Hz)的條件下,再次對結(jié)構(gòu)開展拓?fù)鋬?yōu)化設(shè)計(jì),結(jié)果分別如圖5 所示.可以看到,在高頻外激勵作用下,兩種不同優(yōu)化策略所得結(jié)構(gòu)之間的差異與低頻外激勵作用下的結(jié)果類似,穩(wěn)健性優(yōu)化設(shè)計(jì)仍然形成了額外兩條傳力路徑,可將外激勵的水平分量通過下緣桿件直接傳遞至支撐點(diǎn).除此之外,穩(wěn)健性設(shè)計(jì)的拱形傳力路徑中出現(xiàn)了空隙,這也是材料重新分布,額外傳力路徑消耗材料的結(jié)果.

    圖5 兩種設(shè)計(jì)策略結(jié)構(gòu)拓?fù)鋬?yōu)化在高頻外激勵條件下結(jié)果對比Fig.5 Comparison of the topology optimizations on two different design strategies under an external load of a higher frequency

    表4 列出了在高頻外激勵作用下兩種拓?fù)鋬?yōu)化設(shè)計(jì)所得數(shù)值結(jié)果.可以看出,穩(wěn)健性優(yōu)化設(shè)計(jì)仍以略微增大結(jié)構(gòu)名義動柔順度的代價,大幅降低了結(jié)構(gòu)動柔順度的期望、標(biāo)準(zhǔn)差以及變異系數(shù),使得結(jié)構(gòu)對高頻外激勵方向的變化不再敏感.

    表4 兩種設(shè)計(jì)策略結(jié)構(gòu)拓?fù)湓诟哳l外激勵條件下優(yōu)化結(jié)果對比Table 4 Comparison of the optimized results under an external load of a higher frequency

    3.2 MBB 梁結(jié)構(gòu)

    圖6 所示為MBB 梁結(jié)構(gòu)設(shè)計(jì)區(qū)域及尺寸,左、右兩端固定.理想狀況下,一個集中簡諧力f(t)=10ejωtkN 施加在梁上邊緣的中點(diǎn)處,外激振頻率是200 Hz.但實(shí)際上外激振實(shí)際作用方向會繞著理想方向產(chǎn)生一定的擾動.將結(jié)構(gòu)設(shè)計(jì)區(qū)域均勻劃分成150×25 的有限元網(wǎng)格.假設(shè)外激勵作用方向θ的標(biāo)準(zhǔn)差為 σθ=π/54.

    圖6 MBB 梁設(shè)計(jì)區(qū)域及外激勵Fig.6 Design domain of the MBB beam and the external force

    圖7 分別示出了結(jié)構(gòu)確定性和穩(wěn)健性拓?fù)鋬?yōu)化得到的構(gòu)型設(shè)計(jì).可以看到,兩種設(shè)計(jì)策略得到的結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果仍存在著較大的差異.確定性拓?fù)鋬?yōu)化設(shè)計(jì)中間兩根斜桿張成90°的夾角.表明外激勵先沿著相互正交的兩個方向傳遞到下弦,然后再傳遞到固定基礎(chǔ)上[9].而穩(wěn)健性優(yōu)化設(shè)計(jì)在外激勵作用點(diǎn)兩側(cè)又構(gòu)建了額外的兩根斜桿,產(chǎn)生兩條輔助的傳力路徑.且內(nèi)側(cè)兩根斜桿和下弦桿的截面尺寸均有明顯減小,以滿足對材料體積的約束.由于外激勵的作用方向具有不確定性,施加在結(jié)構(gòu)上的外激勵會存在一定的水平分量,而增加的外側(cè)斜桿件能有效地將這一部分載荷傳遞至約束邊界,以抵抗外激勵方向的擾動,提高結(jié)構(gòu)性能的穩(wěn)健性.

    圖7 兩種設(shè)計(jì)策略結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果對比Fig.7 Comparison of the topology optimizations on two different design strategies

    對比圖7(b)和圖7(c)可以看出,當(dāng)增大目標(biāo)函數(shù)中標(biāo)準(zhǔn)差的權(quán)重因子β以后,所得穩(wěn)健性優(yōu)化設(shè)計(jì)中外側(cè)傾斜構(gòu)件的截面積有所減小,底部路徑截面積有所增大.這使得在外激勵作用方向擾動情形下,結(jié)構(gòu)動柔順度的響應(yīng)值更加集中,即標(biāo)準(zhǔn)差更小.此外,文獻(xiàn)[9]開展了考慮外激勵作用位置不確定性的結(jié)構(gòu)動態(tài)穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì),采用區(qū)間變量表示外激勵位置的擾動,并成功獲得了該算例的穩(wěn)健性最優(yōu)構(gòu)型設(shè)計(jì),如圖7(d)所示.可以看到,穩(wěn)健性設(shè)計(jì)只在激勵點(diǎn)周圍構(gòu)建了作用平臺,明顯擴(kuò)寬了激勵作用區(qū)域.但是,由于外載荷仍然沿豎直方向作用在結(jié)構(gòu)上,因此缺少額外的傾斜傳力構(gòu)件.

    表5 列出了兩種設(shè)計(jì)策略結(jié)構(gòu)拓?fù)鋬?yōu)化數(shù)值結(jié)果.經(jīng)過比較可知,考慮了外激勵方向不確定性的穩(wěn)健性拓?fù)鋬?yōu)化的動柔順度標(biāo)準(zhǔn)差和變異系數(shù),均明顯小于確定性優(yōu)化結(jié)果,使得結(jié)構(gòu)抵抗外載荷作用方向擾動的能力更強(qiáng).當(dāng)權(quán)重因子β增大時,更加強(qiáng)調(diào)動柔順度標(biāo)準(zhǔn)差在目標(biāo)函數(shù)中的占比,優(yōu)化設(shè)計(jì)所得結(jié)構(gòu)的動柔順度標(biāo)準(zhǔn)差會進(jìn)一步降低,使得結(jié)構(gòu)動柔順度的變異系數(shù)進(jìn)一步減小,即結(jié)構(gòu)性能的穩(wěn)健度更高.

    表5 結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果對比Table 5 Comparison of the optimized results

    如圖8 所示,當(dāng)外載荷不再沿著理想方向作用時,確定性結(jié)構(gòu)設(shè)計(jì)的動柔順度將快速增大.相比之下,穩(wěn)健性結(jié)構(gòu)設(shè)計(jì)的動柔順度變化幅度相對較小,而且隨著擾動角度的增大在逐漸減小.這說明穩(wěn)健性設(shè)計(jì)的結(jié)構(gòu)在外激勵名義作用方向時動柔順度最大,隨著外激勵作用方向的變化,相應(yīng)的動柔順度不會超過名義激勵作用方向時的值.與前例類似,由于挪用了一定的結(jié)構(gòu)材料構(gòu)建額外的傳力路徑,穩(wěn)健性設(shè)計(jì)結(jié)果的名義動柔順度Cd0比確定性優(yōu)化相應(yīng)結(jié)果稍大一些.而對比兩條穩(wěn)健性設(shè)計(jì)的曲線還可以發(fā)現(xiàn),當(dāng)權(quán)重因子β增大以后,相應(yīng)的名義動柔順度值會有所減小,而外激勵擾動狀態(tài)下的動柔順度將有所增加,總體上曲線更加平緩,即動柔順度值的離散程度以及標(biāo)準(zhǔn)差值更小,這與表4 中的數(shù)值結(jié)果相一致.總體而言,結(jié)構(gòu)的動柔順度對激勵方向變化更加不敏感.

    圖8 結(jié)構(gòu)動柔順度隨外激勵作用方向變化情況Fig.8 Variation of the structural dynamic compliance caused by the external loading direction disturbances

    圖9 分別繪出了兩種不同優(yōu)化策略的名義結(jié)構(gòu)動柔順度迭代收斂曲線,同時還示出了迭代的中間構(gòu)型和最終優(yōu)化設(shè)計(jì)結(jié)果.在穩(wěn)健性設(shè)計(jì)中,額外的傾斜構(gòu)件從開始即存在,為外激勵作用方向的變化提供可靠的傳力路徑.

    圖9 兩種設(shè)計(jì)策略拓?fù)鋬?yōu)化收斂過程Fig.9 Convergence curves of the dynamic compliance on the two design strategies

    4 結(jié)論

    本文研究了在外激勵作用方向不確定條件下,連續(xù)體結(jié)構(gòu)動態(tài)穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)問題.假設(shè)外激勵作用方向的不確定性變化遵循正態(tài)隨機(jī)分布,借助于二階泰勒展開公式,推導(dǎo)了結(jié)構(gòu)的動柔順度期望和標(biāo)準(zhǔn)差及其對拓?fù)渥兞康囊浑A靈敏度顯式表達(dá)公式.以結(jié)構(gòu)動柔順度均值和標(biāo)準(zhǔn)差的加權(quán)和為目標(biāo)函數(shù),充分考慮了載荷方向的不確定性對結(jié)構(gòu)動響應(yīng)的影響,且利用動柔順度的變異系數(shù)作為結(jié)構(gòu)穩(wěn)健量的衡量指標(biāo).本文用兩個典型算例,驗(yàn)證了所提結(jié)構(gòu)穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)方法的可行性.分析了結(jié)構(gòu)構(gòu)型優(yōu)化設(shè)計(jì)變化的規(guī)律,并比較了結(jié)構(gòu)的動柔順度對外激勵方向擾動的抵抗能力.研究結(jié)果如下.

    (1) 考慮載荷作用方向不確定性條件下得到的穩(wěn)健性結(jié)構(gòu)拓?fù)錁?gòu)型設(shè)計(jì),明顯異與確定性條件下的拓?fù)鋬?yōu)化結(jié)果.并且具有抵抗外激勵作用方向可變性的必要構(gòu)件,在載荷方向隨機(jī)擾動情況下因而具有更可靠的傳力路徑.

    (2) 當(dāng)外載荷作用方向發(fā)生擾動時,穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)的結(jié)構(gòu)動柔順度變化幅度和變異系數(shù)遠(yuǎn)低于相應(yīng)的確定性優(yōu)化結(jié)果.充分表明了穩(wěn)健性優(yōu)化設(shè)計(jì)所得結(jié)構(gòu)抵抗外激勵方向擾動的能力更強(qiáng),動響應(yīng)變化更小,結(jié)構(gòu)的性能更加穩(wěn)健.

    (3) 在相同體積約束情況下,穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)的名義動柔順度稍大于相應(yīng)的確定性拓?fù)鋬?yōu)化結(jié)果.這是因?yàn)樾枰靡徊糠值牟牧蠘?gòu)建抵抗外載荷作用方向擾動的必要構(gòu)件,形成額外的傳力路徑.因而結(jié)構(gòu)穩(wěn)健性拓?fù)鋬?yōu)化設(shè)計(jì)會損失一部分名義動剛度,以換取其抵抗載荷方向擾動能力的提高.但這部分損失的結(jié)構(gòu)名義動剛度,可以通過增加少許結(jié)構(gòu)材料而得到有效彌補(bǔ).

    猜你喜歡
    穩(wěn)健性確定性標(biāo)準(zhǔn)差
    論中國訓(xùn)詁學(xué)與經(jīng)典闡釋的確定性
    論法律解釋的確定性
    法律方法(2022年1期)2022-07-21 09:18:56
    含混還是明證:梅洛-龐蒂論確定性
    用Pro-Kin Line平衡反饋訓(xùn)練儀對早期帕金森病患者進(jìn)行治療對其動態(tài)平衡功能的影響
    會計(jì)穩(wěn)健性的定義和計(jì)量
    商情(2019年3期)2019-03-29 12:04:52
    會計(jì)穩(wěn)健性的文獻(xiàn)綜述
    財訊(2018年22期)2018-05-14 08:55:57
    法律確定性的統(tǒng)合理性根據(jù)與法治實(shí)施
    貨幣政策、會計(jì)穩(wěn)健性與銀行信貸關(guān)系探析
    Effect of Ammonia on the Performance of Catalysts for Selective Hydrogenation of 1-Methylnaphthalene
    對于平均差與標(biāo)準(zhǔn)差的數(shù)學(xué)關(guān)系和應(yīng)用價值比較研究
    日本撒尿小便嘘嘘汇集6| 男人的好看免费观看在线视频 | 亚洲色图 男人天堂 中文字幕| 亚洲av第一区精品v没综合| 最近视频中文字幕2019在线8| 亚洲男人天堂网一区| 欧美一区二区精品小视频在线| 亚洲av中文字字幕乱码综合| 在线播放国产精品三级| 狂野欧美激情性xxxx| 女人被狂操c到高潮| x7x7x7水蜜桃| 毛片女人毛片| 欧美成人性av电影在线观看| 99在线人妻在线中文字幕| 久久久久亚洲av毛片大全| 搞女人的毛片| 欧美中文日本在线观看视频| av福利片在线| 别揉我奶头~嗯~啊~动态视频| 久久精品成人免费网站| 我要搜黄色片| 99久久无色码亚洲精品果冻| svipshipincom国产片| 在线十欧美十亚洲十日本专区| 久久人人精品亚洲av| a级毛片a级免费在线| 日韩有码中文字幕| 免费观看精品视频网站| 少妇粗大呻吟视频| 国产精品精品国产色婷婷| 久久九九热精品免费| 亚洲人成电影免费在线| 母亲3免费完整高清在线观看| 欧美日韩中文字幕国产精品一区二区三区| 五月伊人婷婷丁香| 亚洲精品一区av在线观看| av中文乱码字幕在线| av天堂在线播放| 搡老岳熟女国产| 日韩免费av在线播放| 男人舔女人的私密视频| 精品人妻1区二区| 国产在线观看jvid| 国产私拍福利视频在线观看| 免费无遮挡裸体视频| 国产99白浆流出| 91麻豆精品激情在线观看国产| 亚洲国产欧美一区二区综合| 中文字幕最新亚洲高清| 亚洲精品中文字幕在线视频| 视频区欧美日本亚洲| 99精品久久久久人妻精品| 高潮久久久久久久久久久不卡| 日本黄色视频三级网站网址| 国产一区二区在线av高清观看| 欧美色视频一区免费| 亚洲在线自拍视频| 此物有八面人人有两片| 亚洲精品中文字幕一二三四区| 首页视频小说图片口味搜索| 亚洲一区中文字幕在线| 91麻豆精品激情在线观看国产| 99riav亚洲国产免费| 熟女少妇亚洲综合色aaa.| 成人高潮视频无遮挡免费网站| 日韩成人在线观看一区二区三区| 欧美黑人巨大hd| 欧美中文综合在线视频| 全区人妻精品视频| 亚洲中文av在线| 女警被强在线播放| 搞女人的毛片| 国产91精品成人一区二区三区| 日韩高清综合在线| 精品第一国产精品| 亚洲全国av大片| 琪琪午夜伦伦电影理论片6080| 狂野欧美白嫩少妇大欣赏| 免费一级毛片在线播放高清视频| 午夜免费成人在线视频| 国产成人精品久久二区二区91| 老司机福利观看| 国产精品国产高清国产av| 午夜福利欧美成人| 两性夫妻黄色片| 午夜精品一区二区三区免费看| 老司机靠b影院| 久久香蕉国产精品| www国产在线视频色| 757午夜福利合集在线观看| 一a级毛片在线观看| 亚洲黑人精品在线| 国产在线精品亚洲第一网站| 久久精品人妻少妇| 国产高清有码在线观看视频 | 999精品在线视频| 久久久久久国产a免费观看| 非洲黑人性xxxx精品又粗又长| 欧美日本视频| 国产精品av视频在线免费观看| 18禁裸乳无遮挡免费网站照片| 国产成人av激情在线播放| 三级国产精品欧美在线观看 | 欧美色视频一区免费| 天堂av国产一区二区熟女人妻 | 香蕉国产在线看| 国产激情欧美一区二区| av国产免费在线观看| 亚洲国产欧美一区二区综合| 丁香六月欧美| 国产在线观看jvid| 毛片女人毛片| 香蕉久久夜色| 丰满人妻一区二区三区视频av | 麻豆国产97在线/欧美 | 国产男靠女视频免费网站| 少妇裸体淫交视频免费看高清 | 色在线成人网| 国产成人影院久久av| 九色国产91popny在线| 麻豆久久精品国产亚洲av| 欧美日韩中文字幕国产精品一区二区三区| 色av中文字幕| 12—13女人毛片做爰片一| 国产精品爽爽va在线观看网站| 国产精品亚洲av一区麻豆| 婷婷亚洲欧美| 床上黄色一级片| 99热6这里只有精品| 国产伦一二天堂av在线观看| 免费电影在线观看免费观看| aaaaa片日本免费| 欧美绝顶高潮抽搐喷水| 久久久久久久久久黄片| 免费在线观看黄色视频的| 香蕉久久夜色| 男女床上黄色一级片免费看| 久久久国产成人精品二区| 亚洲一区中文字幕在线| 国产精品九九99| 桃色一区二区三区在线观看| 国产精品亚洲一级av第二区| 成在线人永久免费视频| 亚洲avbb在线观看| 国产高清视频在线播放一区| 岛国在线免费视频观看| 看免费av毛片| 国产一区在线观看成人免费| 欧美色视频一区免费| 91九色精品人成在线观看| 中文字幕人成人乱码亚洲影| 午夜精品在线福利| 99国产精品一区二区三区| 黄色女人牲交| 非洲黑人性xxxx精品又粗又长| 欧美性长视频在线观看| 久久精品影院6| 精品无人区乱码1区二区| 欧美日韩瑟瑟在线播放| 长腿黑丝高跟| 国产成人啪精品午夜网站| 亚洲一区二区三区色噜噜| 免费高清视频大片| 久久久久九九精品影院| 日本 欧美在线| 国产精品美女特级片免费视频播放器 | av视频在线观看入口| 精品久久久久久久久久久久久| 一进一出抽搐gif免费好疼| 日韩欧美免费精品| 黄色丝袜av网址大全| 久久久久久九九精品二区国产 | 女同久久另类99精品国产91| 国产精品自产拍在线观看55亚洲| 老司机福利观看| www日本在线高清视频| 国产在线观看jvid| 国产伦人伦偷精品视频| 亚洲中文字幕一区二区三区有码在线看 | 久久香蕉国产精品| 校园春色视频在线观看| 天天躁夜夜躁狠狠躁躁| 手机成人av网站| 免费高清视频大片| 欧美中文综合在线视频| 成人av一区二区三区在线看| 母亲3免费完整高清在线观看| 啦啦啦韩国在线观看视频| 亚洲乱码一区二区免费版| 最近最新免费中文字幕在线| 国产熟女xx| 两个人免费观看高清视频| 精品熟女少妇八av免费久了| 超碰成人久久| 九九热线精品视视频播放| 精品久久久久久成人av| 亚洲精品一区av在线观看| 一边摸一边抽搐一进一小说| 久久精品亚洲精品国产色婷小说| 国产一区二区在线av高清观看| 成熟少妇高潮喷水视频| 久久久久亚洲av毛片大全| 久久精品91无色码中文字幕| 亚洲自偷自拍图片 自拍| 丰满人妻一区二区三区视频av | 免费观看人在逋| 亚洲天堂国产精品一区在线| 亚洲成人中文字幕在线播放| 别揉我奶头~嗯~啊~动态视频| 国内精品一区二区在线观看| 国产亚洲欧美在线一区二区| 日韩欧美在线二视频| 亚洲 国产 在线| 日本 欧美在线| 午夜日韩欧美国产| 久久久久久久久久黄片| 蜜桃久久精品国产亚洲av| 日本免费一区二区三区高清不卡| 高清毛片免费观看视频网站| 国产精品乱码一区二三区的特点| 波多野结衣高清作品| 国产v大片淫在线免费观看| 日本黄大片高清| 岛国在线免费视频观看| 中文资源天堂在线| 美女黄网站色视频| 午夜福利成人在线免费观看| 女人高潮潮喷娇喘18禁视频| 久99久视频精品免费| 少妇熟女aⅴ在线视频| 国产成年人精品一区二区| 精品久久蜜臀av无| 午夜免费成人在线视频| 欧美日韩一级在线毛片| 一边摸一边做爽爽视频免费| 色播亚洲综合网| 国产激情偷乱视频一区二区| 91麻豆精品激情在线观看国产| 久久久久久亚洲精品国产蜜桃av| 欧美一区二区精品小视频在线| 夜夜爽天天搞| 日韩国内少妇激情av| 哪里可以看免费的av片| 天堂影院成人在线观看| 美女免费视频网站| 夜夜夜夜夜久久久久| 不卡一级毛片| 一边摸一边抽搐一进一小说| 亚洲午夜精品一区,二区,三区| 久久久国产精品麻豆| 国内久久婷婷六月综合欲色啪| 亚洲真实伦在线观看| 成人手机av| 亚洲国产日韩欧美精品在线观看 | 午夜a级毛片| 久久久精品国产亚洲av高清涩受| 亚洲一区二区三区色噜噜| 免费人成视频x8x8入口观看| 超碰成人久久| 中文字幕熟女人妻在线| 国产精品野战在线观看| 亚洲国产欧洲综合997久久,| 欧美大码av| 国产在线观看jvid| 久久久久久亚洲精品国产蜜桃av| 国产精品久久久久久亚洲av鲁大| 欧美色视频一区免费| 婷婷精品国产亚洲av| 香蕉丝袜av| 夜夜躁狠狠躁天天躁| avwww免费| 亚洲人成伊人成综合网2020| 18禁观看日本| 别揉我奶头~嗯~啊~动态视频| 亚洲专区国产一区二区| 亚洲aⅴ乱码一区二区在线播放 | 久久九九热精品免费| 国产乱人伦免费视频| 少妇裸体淫交视频免费看高清 | 欧美黄色片欧美黄色片| 中文字幕最新亚洲高清| 国产精品久久久人人做人人爽| 成年免费大片在线观看| 婷婷精品国产亚洲av在线| 欧美成狂野欧美在线观看| 又黄又爽又免费观看的视频| 国产精品精品国产色婷婷| 欧美最黄视频在线播放免费| 五月玫瑰六月丁香| 久久午夜综合久久蜜桃| 久久亚洲精品不卡| 亚洲欧美精品综合久久99| 久久精品人妻少妇| 日本一区二区免费在线视频| 三级国产精品欧美在线观看 | 中文字幕最新亚洲高清| 色综合亚洲欧美另类图片| 欧美一级毛片孕妇| 午夜影院日韩av| 欧美日韩一级在线毛片| 久久久久久人人人人人| 亚洲色图 男人天堂 中文字幕| 变态另类成人亚洲欧美熟女| 亚洲黑人精品在线| 成人高潮视频无遮挡免费网站| 亚洲精品美女久久久久99蜜臀| svipshipincom国产片| 啦啦啦免费观看视频1| 黄色片一级片一级黄色片| 精品国产超薄肉色丝袜足j| 日本 欧美在线| 欧美精品啪啪一区二区三区| 午夜亚洲福利在线播放| 午夜日韩欧美国产| 亚洲精品粉嫩美女一区| 国产激情欧美一区二区| 亚洲精品美女久久av网站| svipshipincom国产片| 美女免费视频网站| 国内精品久久久久久久电影| 日本黄色视频三级网站网址| 欧美性猛交黑人性爽| 国产成人一区二区三区免费视频网站| 国产高清videossex| 久久久久久九九精品二区国产 | 欧美+亚洲+日韩+国产| 香蕉国产在线看| 午夜两性在线视频| a在线观看视频网站| 亚洲一区中文字幕在线| 欧美成人性av电影在线观看| 国产精品乱码一区二三区的特点| 日本免费a在线| 人人妻人人看人人澡| 日韩欧美一区二区三区在线观看| 亚洲国产精品久久男人天堂| 亚洲美女视频黄频| 午夜成年电影在线免费观看| 国产精品香港三级国产av潘金莲| 很黄的视频免费| 午夜精品在线福利| 真人一进一出gif抽搐免费| 欧美日韩瑟瑟在线播放| 国产高清有码在线观看视频 | 88av欧美| 欧美一级a爱片免费观看看 | 不卡av一区二区三区| 91av网站免费观看| 欧美 亚洲 国产 日韩一| 一个人免费在线观看的高清视频| 一级毛片高清免费大全| 亚洲狠狠婷婷综合久久图片| 中文字幕av在线有码专区| 亚洲熟妇熟女久久| 精品国产乱子伦一区二区三区| 精品乱码久久久久久99久播| 级片在线观看| 美女高潮喷水抽搐中文字幕| 最新美女视频免费是黄的| 国产私拍福利视频在线观看| svipshipincom国产片| 国产高清激情床上av| 成年女人毛片免费观看观看9| 国产三级中文精品| x7x7x7水蜜桃| 亚洲av成人不卡在线观看播放网| 亚洲人与动物交配视频| 日本五十路高清| 非洲黑人性xxxx精品又粗又长| 又粗又爽又猛毛片免费看| 99在线视频只有这里精品首页| 在线观看免费视频日本深夜| 亚洲成av人片免费观看| 免费在线观看影片大全网站| 久9热在线精品视频| 成人av一区二区三区在线看| 美女 人体艺术 gogo| 亚洲专区中文字幕在线| xxxwww97欧美| 国产成人系列免费观看| 校园春色视频在线观看| 久久天堂一区二区三区四区| 两个人看的免费小视频| 床上黄色一级片| 国产激情偷乱视频一区二区| 欧美性猛交╳xxx乱大交人| 99久久99久久久精品蜜桃| 精品久久久久久久毛片微露脸| 成人精品一区二区免费| 日韩欧美免费精品| 激情在线观看视频在线高清| 18美女黄网站色大片免费观看| 亚洲专区国产一区二区| 午夜精品在线福利| 欧美+亚洲+日韩+国产| 老司机午夜十八禁免费视频| 久久国产乱子伦精品免费另类| www国产在线视频色| 啦啦啦免费观看视频1| 校园春色视频在线观看| 亚洲激情在线av| 一本精品99久久精品77| 波多野结衣高清作品| 日日夜夜操网爽| 91成年电影在线观看| 国产亚洲精品一区二区www| 日本免费a在线| 又黄又爽又免费观看的视频| 丝袜人妻中文字幕| 亚洲一区二区三区色噜噜| 欧美日韩一级在线毛片| 国产野战对白在线观看| 熟女少妇亚洲综合色aaa.| 国产免费av片在线观看野外av| 午夜久久久久精精品| 亚洲成人国产一区在线观看| 国产激情欧美一区二区| 免费看美女性在线毛片视频| 手机成人av网站| 这个男人来自地球电影免费观看| 悠悠久久av| 亚洲人成电影免费在线| 国产av在哪里看| 女人被狂操c到高潮| 一a级毛片在线观看| 在线播放国产精品三级| 老汉色av国产亚洲站长工具| 黑人操中国人逼视频| 精品免费久久久久久久清纯| 777久久人妻少妇嫩草av网站| 日韩大码丰满熟妇| 黄色 视频免费看| 午夜精品久久久久久毛片777| 神马国产精品三级电影在线观看 | 午夜免费观看网址| 免费在线观看日本一区| or卡值多少钱| 婷婷亚洲欧美| 毛片女人毛片| 亚洲精品在线观看二区| 亚洲精品国产一区二区精华液| 国产精品久久久久久精品电影| 国产久久久一区二区三区| 国产私拍福利视频在线观看| 国产激情偷乱视频一区二区| 搡老妇女老女人老熟妇| 成年版毛片免费区| 一本一本综合久久| 亚洲av成人不卡在线观看播放网| 精品日产1卡2卡| 久久天堂一区二区三区四区| 免费一级毛片在线播放高清视频| 精品免费久久久久久久清纯| 亚洲熟女毛片儿| 成人av在线播放网站| 欧美+亚洲+日韩+国产| 国产成人欧美在线观看| 一进一出好大好爽视频| 看片在线看免费视频| 制服丝袜大香蕉在线| 欧美不卡视频在线免费观看 | 日本五十路高清| 午夜福利成人在线免费观看| 男女床上黄色一级片免费看| 在线观看舔阴道视频| 欧美成人性av电影在线观看| 动漫黄色视频在线观看| xxxwww97欧美| 国产精品久久久人人做人人爽| 精品国产乱子伦一区二区三区| 露出奶头的视频| 亚洲国产精品合色在线| 91国产中文字幕| 久久欧美精品欧美久久欧美| 精品午夜福利视频在线观看一区| 99热只有精品国产| 中文字幕高清在线视频| 免费人成视频x8x8入口观看| 黄色 视频免费看| 久久久久久亚洲精品国产蜜桃av| 午夜福利欧美成人| 欧美一级a爱片免费观看看 | 精品少妇一区二区三区视频日本电影| 欧美色视频一区免费| 人人妻人人澡欧美一区二区| 国产精品永久免费网站| 亚洲中文av在线| www日本在线高清视频| 日本一二三区视频观看| 此物有八面人人有两片| 日韩欧美在线二视频| 久久精品影院6| 一本综合久久免费| 十八禁人妻一区二区| 12—13女人毛片做爰片一| 女人被狂操c到高潮| 国产私拍福利视频在线观看| 国产亚洲精品综合一区在线观看 | 黄色 视频免费看| 成人国语在线视频| 丰满人妻一区二区三区视频av | 首页视频小说图片口味搜索| 免费电影在线观看免费观看| 黄色毛片三级朝国网站| 国产成人精品无人区| 亚洲精品中文字幕在线视频| 国产成人欧美在线观看| 久久久久久久久久黄片| 一本综合久久免费| 男人的好看免费观看在线视频 | 亚洲欧美激情综合另类| 国产精品久久久久久亚洲av鲁大| 亚洲一码二码三码区别大吗| 免费在线观看日本一区| 人人妻,人人澡人人爽秒播| 可以在线观看的亚洲视频| 女人高潮潮喷娇喘18禁视频| 日韩欧美国产在线观看| videosex国产| 变态另类成人亚洲欧美熟女| 午夜激情福利司机影院| 99久久精品热视频| 人人妻人人看人人澡| 69av精品久久久久久| 成人特级黄色片久久久久久久| 久久久久性生活片| 久久99热这里只有精品18| 看片在线看免费视频| 99久久无色码亚洲精品果冻| 成人av一区二区三区在线看| 一卡2卡三卡四卡精品乱码亚洲| 男女床上黄色一级片免费看| 特级一级黄色大片| 成人手机av| 亚洲国产精品sss在线观看| 制服人妻中文乱码| 精品久久久久久久毛片微露脸| 变态另类成人亚洲欧美熟女| 在线看三级毛片| 日韩精品免费视频一区二区三区| 免费一级毛片在线播放高清视频| 一级片免费观看大全| or卡值多少钱| 国产亚洲av高清不卡| av片东京热男人的天堂| 丁香欧美五月| 精品午夜福利视频在线观看一区| 哪里可以看免费的av片| www日本在线高清视频| 午夜精品久久久久久毛片777| 色噜噜av男人的天堂激情| 看免费av毛片| 啪啪无遮挡十八禁网站| 国产爱豆传媒在线观看 | 亚洲av美国av| 在线十欧美十亚洲十日本专区| 欧美人与性动交α欧美精品济南到| 精华霜和精华液先用哪个| 免费在线观看日本一区| 99久久无色码亚洲精品果冻| 此物有八面人人有两片| 国产又黄又爽又无遮挡在线| 日本免费a在线| 久久婷婷成人综合色麻豆| 欧美极品一区二区三区四区| 黄色视频不卡| 欧美精品亚洲一区二区| 午夜激情福利司机影院| 宅男免费午夜| 超碰成人久久| tocl精华| 欧美黑人精品巨大| 中文资源天堂在线| 亚洲五月婷婷丁香| 无遮挡黄片免费观看| 91老司机精品| 免费人成视频x8x8入口观看| 神马国产精品三级电影在线观看 | 俄罗斯特黄特色一大片| 不卡一级毛片| 久久草成人影院| 亚洲中文av在线| 国产精品亚洲一级av第二区| 狂野欧美激情性xxxx| 国产一区二区在线av高清观看| 亚洲va日本ⅴa欧美va伊人久久| 成年人黄色毛片网站| 国产精品久久久av美女十八| 日日夜夜操网爽| 免费看十八禁软件| 在线视频色国产色| 天堂av国产一区二区熟女人妻 | 成人国产综合亚洲| 热99re8久久精品国产| 亚洲激情在线av| 99国产精品一区二区三区| 免费在线观看黄色视频的| 香蕉av资源在线| 成人永久免费在线观看视频| 国产真实乱freesex| 一进一出好大好爽视频| 日韩国内少妇激情av| 国产麻豆成人av免费视频| 美女黄网站色视频| 成人国产综合亚洲| aaaaa片日本免费| 一二三四社区在线视频社区8| 国产精品一区二区精品视频观看| aaaaa片日本免费| АⅤ资源中文在线天堂| av福利片在线观看| 最新在线观看一区二区三区| 国产一区二区在线av高清观看| 久久香蕉精品热|