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

    輪盤概念設(shè)計中拓撲和形狀同時優(yōu)化方法

    2015-12-20 05:30:56范俊尹澤勇王建軍米棟閆成
    北京航空航天大學學報 2015年3期
    關(guān)鍵詞:輪盤形狀約束

    范俊,尹澤勇,王建軍,米棟,3,閆成

    (1.北京航空航天大學 能源與動力工程學院,北京100191;2.陸航研究所,北京101121;3.中航工業(yè) 航空動力機械研究所,株洲412002)

    結(jié)構(gòu)優(yōu)化包括拓撲、形狀和尺寸優(yōu)化3個階段,對于前2個優(yōu)化階段,通常的優(yōu)化順序是,先進行拓撲優(yōu)化,在拓撲優(yōu)化結(jié)果的基礎(chǔ)上,再進行形狀優(yōu)化.在這種分步進行的優(yōu)化過程中,需要人為設(shè)定一個大小一定的不可調(diào)的初始設(shè)計區(qū)域用于拓撲優(yōu)化;并且在基于變密度法(即SIMP法)的拓撲優(yōu)化過程結(jié)束后,還需要人為設(shè)定一個用于刪除中間密度單元的相對密度門檻值.然而,假如這兩種人為設(shè)定參數(shù)選擇不恰當,就有可能使得優(yōu)化結(jié)果并不是最優(yōu)解.

    為了減小上述兩類人為選擇參數(shù)對優(yōu)化結(jié)果的影響,本文提出了一種拓撲和形狀同時優(yōu)化(STSO)的方法.國外學者 Deaton等[1]認為從拓撲優(yōu)化結(jié)果到形狀優(yōu)化初始設(shè)計模型的轉(zhuǎn)化,是未來連續(xù)體多學科拓撲優(yōu)化研究重點之一.Ansola等[2]論述了一個形狀和拓撲組合優(yōu)化的方法,在該方法中,交替進行形狀和拓撲優(yōu)化步驟,而在每一個步驟中,先進行形狀優(yōu)化再進行拓撲優(yōu)化;這些步驟不斷重復,直到得到一個收斂的全局解.Hassani等[3]論述了一個板殼結(jié)構(gòu)形狀和拓撲同時優(yōu)化的方法,在該方法中,形狀和拓撲優(yōu)化不是分開進行的,而是并行的.上述兩篇文獻除了優(yōu)化流程不一致之外,在優(yōu)化算法上也是不一樣的,Ansola采取的是單純形法,Hassani采取的是移動漸進線法.國內(nèi)學者張衛(wèi)紅等[4]為了解決壓力載荷作用下的結(jié)構(gòu)輕量化設(shè)計問題,直接采用B樣條曲線描述壓力加載面,通過拓撲和形狀變量的聯(lián)合優(yōu)化滿足了工程實際對結(jié)構(gòu)輕量化與邊界的功能性與光滑性設(shè)計要求.雖然本文中的拓撲和形狀優(yōu)化也是并行的,但與上述文獻只考慮應(yīng)變能和應(yīng)力優(yōu)化響應(yīng)不同,本文的優(yōu)化模型中還考慮了以頻率為約束的三維實體結(jié)構(gòu)(輪盤)動力學優(yōu)化,采用的優(yōu)化算法是序列二次規(guī)劃優(yōu)化算法(SQP法).

    輪盤是航空發(fā)動機中關(guān)鍵部件,已有文獻對輪盤進行了拓撲優(yōu)化設(shè)計研究.文獻[5]以體積剛度比為目標,應(yīng)力、破裂轉(zhuǎn)速、低循環(huán)疲勞壽命為約束建立拓撲優(yōu)化數(shù)學模型,采用基于隨機抽樣敏度分析的雙向漸進結(jié)構(gòu)拓撲優(yōu)化方法,得到了航空發(fā)動機雙輻板形式的渦輪盤結(jié)構(gòu).文獻[6]以體積剛度比為目標、應(yīng)力為約束建立了拓撲優(yōu)化數(shù)學模型,采用基于自適應(yīng)隨機抽樣策略的周期結(jié)構(gòu)拓撲優(yōu)化方法得到了多輻板風扇盤.文獻[7]以最小體積為目標、應(yīng)力為約束建立了典型風扇盤子午面結(jié)構(gòu)優(yōu)化設(shè)計模型,得到了一個多輻板盤結(jié)構(gòu).從國外先進的發(fā)動機來看,許多高性能發(fā)動機寬弦風扇盤經(jīng)過拓撲優(yōu)化設(shè)計,選用了多輻板的盤轂混合式輪盤結(jié)構(gòu)[6],如CFM56-7,RB211-535E4等采用了雙輻板輪盤結(jié)構(gòu).但是,上述文獻一般只考慮了靜力學目標和約束,且拓撲和形狀優(yōu)化都是分步進行的.為了設(shè)計出既輕且剛性適當又安全可靠的輪盤,有必要研究輪盤的模態(tài)特性,進行與實體輪盤頻率等有關(guān)的動力學優(yōu)化設(shè)計.

    本文首先建立了基于SIMP法的拓撲和形狀同時優(yōu)化(STSO)數(shù)學模型,分析了相應(yīng)的靈敏度,使用SQP法進行求解.然后,以板殼結(jié)構(gòu)為例,進行了拓撲和形狀同時優(yōu)化,將其結(jié)果與分步優(yōu)化進行比對,探討了相對密度門檻值的選取和兩種優(yōu)化變量之間的相互作用對最終優(yōu)化結(jié)果的影響,說明拓撲和形狀同時優(yōu)化的優(yōu)點.最后,本文將同時優(yōu)化的方法應(yīng)用于輪盤概念設(shè)計中,對比分析了不同振型所對應(yīng)的頻率約束時實體輪盤拓撲和形狀同時優(yōu)化方法與單獨拓撲優(yōu)化方法的不同結(jié)果.

    1 STSO方法數(shù)學模型

    SIMP法是拓撲優(yōu)化方法的一種,由于其數(shù)學論證嚴密、便于實施等優(yōu)點而得到廣泛應(yīng)用[8-9].本文的STSO方法將基于SIMP法來展開研究.

    SIMP法建立在結(jié)構(gòu)有限元模型的基礎(chǔ)上,設(shè)計變量是各單元的相對密度,是在[0,1]區(qū)間的連續(xù)變量,第 i個單元的相對密度用 ρi來表示(i=1,2,…,N),N指單元總數(shù).設(shè)單元的原始密度是 ρ,優(yōu)化后密度是 ρe,則存在關(guān)系式 ρi= ρe/ρ;SIMP法中,單元材料屬性是單元相對密度的指數(shù)函數(shù),設(shè)第i個單元初始和優(yōu)化后的彈性模量分別是E0,Ei(ρi),第i個單元初始和優(yōu)化后的剛度矩陣分別是K0i,Ki(ρi),則存在關(guān)系式:

    式中P為懲罰因子,其作用是對中間相對密度進行懲罰,使結(jié)構(gòu)單元相對密度盡可能趨于0或者1.用于航空發(fā)動機輪盤拓撲和形狀優(yōu)化的數(shù)學模型是以體積、應(yīng)力和頻率為約束,以最小化結(jié)構(gòu)應(yīng)變能為目標.

    式中,C為應(yīng)變能;F為載荷;U為位移;σeq為當量應(yīng)力;σ0.1為 0.1% 屈服強度;V0i為第i個單元的初始體積;V為優(yōu)化之后結(jié)構(gòu)體積;V0為結(jié)構(gòu)原始體積;V0iρi為優(yōu)化后單元的等效體積;g(f)為頻率f的函數(shù);K為全局剛度陣;M為全局質(zhì)量陣;λ為特征值;φ為對應(yīng)的特征向量;Di為第i個節(jié)點的位移矢量;D為一個位移常量;j和為形狀優(yōu)化變量yj的上下限.

    2 靈敏度分析

    靈敏度分析包括拓撲優(yōu)化變量和形狀優(yōu)化變量靈敏度分析兩個部分,本文著重分析目標和約束對拓撲優(yōu)化變量的靈敏度.

    2.1 頻率靈敏度

    頻率約束的靈敏度,是通過多自由度無阻尼系統(tǒng)自由振動的運動方程,利用特征值法求得.對式(2)中的運動方程,通過對特征向量用質(zhì)量矩陣進行規(guī)格化[10-11],則特征值對設(shè)計變量的導數(shù)可表示為

    式中,λj為第j階特征值;φj為第j階特征值對應(yīng)的特征向量.由于

    式中,M0i為第i個單元的初始質(zhì)量矩陣;Mi為第i個單元優(yōu)化后的質(zhì)量矩陣,所以

    將式(5)代入式(3),可得

    由 λj=f2j,可得

    2.2 應(yīng)變能、體積和位移靈敏度

    由于式(2)中外力F是一個常量,則有

    由于結(jié)構(gòu)剛度矩陣的對稱性,并將式(8)代入可得應(yīng)變能的靈敏度為

    容易得到,體積約束的靈敏度為

    由式(8)可得位移靈敏度為

    2.3 應(yīng)力靈敏度

    現(xiàn)階段應(yīng)力相關(guān)的優(yōu)化問題面臨著諸多的難題[12],如優(yōu)化過程中由于應(yīng)力約束不連續(xù)導致的應(yīng)力奇異現(xiàn)象、局部應(yīng)力約束過多導致的超大計算規(guī)模等計算難題.隋允康等[13]等利用Mises強度理論,提出了應(yīng)力約束全局化策略,將局部的應(yīng)力約束問題轉(zhuǎn)化為結(jié)構(gòu)整體的應(yīng)變能約束問題.榮見華等[14]在優(yōu)化迭代循環(huán)的每一輪子循環(huán)迭代求解開始時,通過形成和引進新的位移和應(yīng)力約束限,自動構(gòu)建設(shè)計變量移動限,將結(jié)構(gòu)應(yīng)力約束歸并為幾個最可能的有效應(yīng)力約束,從而減少應(yīng)力靈敏度的分析量.París等[15]總結(jié)了連續(xù)體拓撲優(yōu)化中應(yīng)力約束靈敏度分析的一般方法.

    式(2)中第i號單元的VonMises應(yīng)力近似地可以表示為

    式中,Bhi為第i號單元的第h個高斯積分點的幾何矩陣;ni為第i號單元的高斯積分點數(shù);Ui為僅與第i號單元所有自由度相應(yīng)的U中元素構(gòu)成的單元位移矢量.則應(yīng)力靈敏度為

    由式(11)可以求得?Ui/?ρi為?U/?ρi中與第i號單元所有自由度相應(yīng)元素構(gòu)成的矢量,將其代入式(13),則可以求得應(yīng)力靈敏度.

    對于形狀優(yōu)化變量的靈敏度,并不能十分明確地列出[3].在本文中,采用一種半解析的辦法,這樣,形狀優(yōu)化變量的靈敏度一部分由解析解推出來,一部分由有限差分近似推出來,具體的推導見文獻[16-18].

    在求得目標和約束函數(shù)的靈敏度基礎(chǔ)上,運用SQP法進行求解,SQP算法具體原理可參考文獻[19-20].

    3 分步優(yōu)化和同時優(yōu)化對比分析

    以圖1中的板殼結(jié)構(gòu)為例來進行對比分析,待優(yōu)化的板長度為100mm,寬度為10mm,厚度為1 mm,劃分成1000個1×1的單元.

    板的左端受全約束,板的右端上頂點處施加100 N的集中力.材料的彈性模量是 2.1×105MPa,泊松比是0.3.優(yōu)化的設(shè)計區(qū)域是整個黃色區(qū)域.為了對比分步優(yōu)化和同時優(yōu)化結(jié)果,相比于式(2),減少了頻率和應(yīng)力約束,將優(yōu)化目標改為體積最小,約束是板的右端下頂點的位移小于3 mm,建立數(shù)學模型見式(14).

    式中d101為第101個節(jié)點的垂直位移,該節(jié)點位于板的右端下頂點.

    圖1 結(jié)構(gòu)優(yōu)化的初始有限元模型Fig.1 Original finite element model for structural optimization

    3.1 先拓撲后形狀分步優(yōu)化

    圖2(a)描述了通常一般優(yōu)化順序:先進行拓撲優(yōu)化,在拓撲優(yōu)化結(jié)果的基礎(chǔ)上,再進行形狀優(yōu)化.其中,在基于SIMP法的拓撲優(yōu)化過程中,并不刪除單元,而是通過改變單元的相對密度值來滿足約束條件.在拓撲優(yōu)化過程結(jié)束后,需要人為設(shè)定相對密度門檻值,刪除低于該值的單元后,得到最終拓撲優(yōu)化結(jié)果,用于形狀優(yōu)化.

    圖2 結(jié)構(gòu)優(yōu)化流程Fig.2 Flow of structure optimization

    按照這個優(yōu)化順序,先對圖1中的板殼結(jié)構(gòu)進行拓撲優(yōu)化,迭代25步收斂,取相對密度門檻值為0.16,刪除相對密度小于該值的單元,得到的拓撲優(yōu)化結(jié)果如圖3(a)所示.然后,定義該結(jié)果中的形狀優(yōu)化變量.本文出于演示論證的目的,只定義了一個形狀優(yōu)化變量.設(shè)形狀變量為圖形上邊界中點處節(jié)點的位置,該節(jié)點可上下移動.同樣采用SQP算法,迭代3步后收斂,得到的結(jié)果如圖3(b)所示.從結(jié)果看出,在拓撲結(jié)構(gòu)一定的情況下,進行形狀優(yōu)化,其拓撲構(gòu)型并沒有改變,只是形狀的邊界尺寸發(fā)生了變化,優(yōu)化的效果不明顯.

    圖3 拓撲優(yōu)化結(jié)果及基于拓撲優(yōu)化的形狀優(yōu)化結(jié)果Fig.3 Topology optimization result and shape optimization result based on topology optimization

    因此,拓撲優(yōu)化結(jié)果對于后續(xù)的優(yōu)化有著重要的作用,而人為選定的相對密度門檻值對拓撲優(yōu)化結(jié)果有著重要的影響.為了減小這種主觀選擇參數(shù)對結(jié)果的影響,有必要拓撲和形狀同時優(yōu)化.

    3.2 先形狀后拓撲分步優(yōu)化

    不同的初始設(shè)計模型會產(chǎn)生不同的拓撲優(yōu)化結(jié)果,為了得到較好的拓撲優(yōu)化結(jié)果,文獻[2]認為,應(yīng)該先通過形狀優(yōu)化來確定初始設(shè)計模型,再進行拓撲優(yōu)化.本文將以算例表明,先形狀后拓撲的優(yōu)化流程也并不一定能得到最優(yōu)解.

    對于圖1所示模型,先進行形狀優(yōu)化,以確定用于拓撲優(yōu)化的初始設(shè)計模型.定義兩個形狀優(yōu)化變量:一是模型上邊界右端頂點處節(jié)點的位置,該節(jié)點可以上下移動;二是模型上邊界中間點處節(jié)點的位置,該節(jié)點可以上下移動.優(yōu)化在迭代5步后收斂,得到的形狀優(yōu)化結(jié)果如圖4(a)所示.將該形狀優(yōu)化結(jié)果作為初始設(shè)計模型,再進行拓撲優(yōu)化,迭代5步后收斂,得到的拓撲優(yōu)化結(jié)果如圖4(b)所示.從圖中可以看出,在形狀優(yōu)化結(jié)果基礎(chǔ)上再拓撲優(yōu)化,發(fā)生的改變并不多.這也表明,經(jīng)過形狀優(yōu)化確定的初始設(shè)計模型,不一定是最合適的用于拓撲優(yōu)化的初始設(shè)計模型;同時從圖3(a)和圖4(b)的比對看出,不同的初始設(shè)計模型進行拓撲優(yōu)化會產(chǎn)生不同的優(yōu)化解.

    圖4 形狀優(yōu)化結(jié)果及基于形狀優(yōu)化的拓撲優(yōu)化結(jié)果Fig.4 Shape optimization result and topology optimization result based on shape optimization

    3.3 拓撲和形狀同時優(yōu)化

    在分步優(yōu)化過程中,人為設(shè)定的初始設(shè)計區(qū)域和相對密度門檻值,對最終優(yōu)化結(jié)果有著重要的影響,為了減小這種人為設(shè)定參數(shù)的影響,本節(jié)使用拓撲和形狀同時優(yōu)化的方法對該板殼結(jié)構(gòu)進行優(yōu)化.

    對圖1中的初始設(shè)計模型,進行拓撲和形狀同時優(yōu)化,其中,形狀優(yōu)化變量與3.2節(jié)保持一致.優(yōu)化的流程如圖2(b)所示,迭代59步后得到STSO優(yōu)化結(jié)果如圖5(a)所示,優(yōu)化結(jié)果類似于桁架結(jié)構(gòu).該結(jié)果的重量已經(jīng)較小,可直接進入詳細設(shè)計階段對應(yīng)的尺寸優(yōu)化.

    圖5 形狀和拓撲同時優(yōu)化結(jié)果Fig.5 Simultaneous shape and topology optimization result

    因為3.1節(jié)和3.2節(jié)所述分步優(yōu)化過程中,第2步優(yōu)化的效果相比于第1步要小得多,圖5(b)中描述了兩種分步優(yōu)化方法的第1步優(yōu)化結(jié)果(體積)和同時優(yōu)化結(jié)果隨優(yōu)化步數(shù)變化趨勢.從圖中可以看出,相比于分步優(yōu)化,同時優(yōu)化后的結(jié)構(gòu)體積要小得多,這說明同時優(yōu)化方法的效果是最好的.

    此外,從上述的優(yōu)化結(jié)果對比分析看,形狀和拓撲優(yōu)化是一個整體,是相互作用的兩個過程,割裂地進行優(yōu)化有可能得不到全局最優(yōu)解.

    4 STSO方法在輪盤概念設(shè)計階段的應(yīng)用

    4.1 初始設(shè)計模型

    使用STSO方法優(yōu)化的輪盤是一個全三維模型,該模型所參考的航空發(fā)動機渦輪盤子午面如圖6(a)所示[4].其中,中心孔半徑 Rb=65 mm,輪緣半徑Rr=278.6 mm,輪緣寬度Wr=44 mm,輪緣高度Hr=11.4mm,輪轂寬度Wb=91mm;輪緣均布載荷Pr=169.39 MPa,工作轉(zhuǎn)速 N=13333 r/min;輪盤所用的材料為 GH4169,文中計算采用了600°C時該材料的性能參數(shù).

    圖6(b)是用于結(jié)構(gòu)優(yōu)化的初始設(shè)計模型,其中藍色非設(shè)計區(qū)域與參考盤的輪緣部分是一致的;綠色設(shè)計區(qū)域?qū)挾葏⒖驾嗇瀸挾?,將其暫定為Wb=160 mm,綠色設(shè)計區(qū)域高度為Rr-Rb-Hr=202.2 mm.

    STSO方法在輪盤應(yīng)用中的優(yōu)化數(shù)學模型除了不考慮位移約束外,與式(2)所示相同.需要說明的是,本文目的主要是闡述STSO方法和概念,如圖2(c)所示,該方法是應(yīng)用于輪盤概念設(shè)計階段,其優(yōu)化結(jié)果距實際應(yīng)用還有較長距離.因此,雖然輪盤不同的地方應(yīng)有不同的應(yīng)力約束條件,但是為了突出動力學優(yōu)化以及對比STSO方法和單獨拓撲優(yōu)化方法的優(yōu)劣,本文將各個地方的靜力學要求簡化為0.1%屈服極限的6/10.

    另外,在STSO方法中,形狀優(yōu)化設(shè)計變量為圖6(b)中4個角點(1,2,3,4 點)水平方向的位置,為了獲取最佳的初始設(shè)計區(qū)域,將其變化范圍設(shè)為[-80,80];拓撲優(yōu)化設(shè)計變量為綠色區(qū)域單元的相對密度值.

    4.2 STSO方法在輪盤概念設(shè)計階段的應(yīng)用

    圖6中初始設(shè)計模型的第4階和第5階頻率是679 Hz和1060 Hz,其振型如圖7所示,分別是節(jié)圓和節(jié)徑振動.從第3節(jié)可知,分步優(yōu)化中第2步優(yōu)化的效果相比于第1步要小得多,并且同時優(yōu)化方法所得的優(yōu)化結(jié)果如有必要也可以用形狀和尺寸優(yōu)化方法進行進一步的細節(jié)設(shè)計,因此,本節(jié)將運用STSO方法和單獨拓撲優(yōu)化方法對該輪盤進行結(jié)構(gòu)優(yōu)化,對比分析兩種不同方法下不同振型對應(yīng)的頻率約束對輪盤優(yōu)化結(jié)果結(jié)構(gòu)形式的影響.

    圖6 初始設(shè)計模型Fig.6 Initial design model

    圖7 初始設(shè)計模型的第4和第5階振型Fig.7 Fourth and fifth mode of vibration of initial design model

    4.2.1 降低節(jié)徑振動頻率

    將式(2)中的頻率約束g(f)≤0確定為降低節(jié)徑振動頻率約束為f5≤900.單獨進行拓撲優(yōu)化,得到的優(yōu)化結(jié)果見圖8(a)所示;拓撲和形狀同時優(yōu)化的結(jié)果如圖8(b)所示.

    圖8 降低節(jié)徑振動頻率的優(yōu)化結(jié)果對比Fig.8 Comparison of optimization results with lower frequency of nodal diameter vibration

    從結(jié)果看,兩種方法得到優(yōu)化結(jié)果拓撲形式是不一樣的,同時優(yōu)化結(jié)果拓撲形式較單獨拓撲優(yōu)化結(jié)果更為簡潔.

    圖9分別列出了拓撲和形狀同時優(yōu)化,以及單獨拓撲優(yōu)化的應(yīng)變能、頻率、體積分數(shù)(優(yōu)化后的體積與初始體積的比值)隨迭代步數(shù)變化情況.從圖中可以看出,兩種優(yōu)化的頻率和體積約束在最后迭代步均分別達到相應(yīng)上限值,而優(yōu)化的目標——應(yīng)變能則不同,拓撲和形狀同時優(yōu)化結(jié)果的應(yīng)變能較單獨拓撲優(yōu)化的應(yīng)變能更小,減小了約2.4%,并且其收斂速度較快.

    圖9 降低節(jié)徑振動頻率的優(yōu)化結(jié)果隨迭代步數(shù)變化趨勢對比Fig.9 Variation tendency comparison of optimization results with lower frequency of nodal diameter vibration in different iteration steps

    4.2.2 提高節(jié)徑振動頻率

    明確式(2)中的頻率約束,提高節(jié)徑振動頻率約束為f5≥1100,優(yōu)化結(jié)果如圖10所示.圖中給出的是優(yōu)化后相對密度的分布情況(刪除了較小相對密度單元的優(yōu)化結(jié)果),紅色表示優(yōu)化后相對密度較高的單元顏色,藍色表示優(yōu)化后相對密度較低的單元顏色.

    在圖10(a)和圖10(b)中,當使用單獨拓撲優(yōu)化方法時,優(yōu)化結(jié)果出現(xiàn)了開孔現(xiàn)象.而在圖10(c)和圖10(d)中,當使用STSO方法時,優(yōu)化結(jié)果與參考盤類似,是一個單幅板的輪盤,并未出現(xiàn)開孔現(xiàn)象.

    相較于單獨拓撲優(yōu)化方法的結(jié)果,使用STSO方法得到的結(jié)果,更符合現(xiàn)役發(fā)動機少有開孔輪盤的現(xiàn)狀.

    圖11也分別列出了拓撲和形狀同時優(yōu)化,以及單獨拓撲優(yōu)化的應(yīng)變能、頻率、體積隨迭代步數(shù)變化情況.從圖中可以看出,兩種優(yōu)化的體積約束在最后迭代步均達到相應(yīng)上限值,而優(yōu)化的目標——應(yīng)變能則不同,拓撲和形狀同時優(yōu)化結(jié)果的應(yīng)變能較單獨拓撲優(yōu)化的應(yīng)變能更小,減小了約8.5%,并且其收斂速度較快.

    圖10 提高節(jié)徑振動頻率的優(yōu)化結(jié)果對比Fig.10 Comparison of optimization results with higher frequency of nodal diameter vibration

    圖11 提高節(jié)徑振動頻率的優(yōu)化結(jié)果隨迭代步數(shù)變化趨勢對比Fig.11 Variation tendency comparison of optimization results with higher frequency of nodal diameter vibration in different iteration steps

    4.2.3 提高節(jié)圓振動頻率

    明確式(2)中的頻率約束,提高節(jié)圓振動頻率約束為f4≥700,優(yōu)化結(jié)果如圖12所示.從圖中可以看出,與單獨拓撲優(yōu)化結(jié)果類似,拓撲和形狀同時優(yōu)化得到了類似于雙幅板形式的結(jié)構(gòu).二者不同的是,同時優(yōu)化結(jié)果的幅板之間的間隔更小.

    圖12 提高節(jié)圓振動頻率的優(yōu)化結(jié)果對比Fig.12 Comparison of optimization results with higher frequency of umbrella vibration

    圖12(d)所示的雙幅板渦輪盤拓撲形式,經(jīng)過完善約束條件,以及在該優(yōu)化結(jié)果基礎(chǔ)上進一步的形狀和尺寸優(yōu)化,是能夠應(yīng)用于工程實踐的.美國空軍 CRRT計劃[5]中采用了類似于圖12(d)(如圖12)的雙幅板渦輪盤而取得了明顯效果.

    圖13分別列出了拓撲和形狀同時拓撲優(yōu)化,以及單獨拓撲優(yōu)化的應(yīng)變能、頻率、體積隨迭代步數(shù)變化情況.從圖中可以看出,兩種優(yōu)化的體積約束在最后迭代步均達到相應(yīng)上限值,而優(yōu)化的目標——應(yīng)變能則不同,拓撲和形狀同時優(yōu)化結(jié)果的應(yīng)變能較單獨拓撲優(yōu)化的應(yīng)變能更小,減小了約28.3%,并且其收斂速度較快.

    4.2.4 降低節(jié)圓振動頻率

    明確式(2)中的頻率約束,降低節(jié)圓振動頻率約束為f4≤500,優(yōu)化結(jié)果如圖14所示.

    圖13 提高節(jié)圓振動頻率的優(yōu)化結(jié)果隨迭代步數(shù)變化趨勢對比Fig.13 Variation tendency comparison of optimization results with higher frequency of umbrella vibration in different iteration steps

    圖14 降低節(jié)圓振動頻率的拓撲和形狀同時拓撲優(yōu)化結(jié)果Fig.14 Results of simultaneous topology and shape optimization with lower frequency of umbrella vibration

    形狀優(yōu)化變量使得初始設(shè)計模型中輪盤中心孔處的厚度不斷變化,但是如果厚度過小,如圖14(b)所示,會使該處單元的邊長比過大,而導致在優(yōu)化過程中由于網(wǎng)格畸變過大,使得優(yōu)化過程提前終止,這樣的優(yōu)化結(jié)果是無效的.圖15描述了此時中心孔處單元邊長比達到117,超過了上限值(該值等于100).

    圖15 邊長比過大的單元Fig.15 Elements with overlarge aspect ratio

    為了得到合適的優(yōu)化結(jié)果,需修改形狀優(yōu)化設(shè)計變量的取值范圍,由前述的[-80,80]改為[-25,25],使用同時優(yōu)化方法得到的結(jié)果如圖16(c)、圖16(d)所示.與單獨拓撲優(yōu)化方法所得結(jié)果類似,均為單幅板形式的拓撲.

    圖16 降低節(jié)圓振動頻率的優(yōu)化結(jié)果對比Fig.16 Comparison of optimization results with lower frequency of umbrella vibration

    使用同時優(yōu)化方法所得結(jié)構(gòu)的應(yīng)變能較單獨拓撲優(yōu)化方法所得結(jié)構(gòu),減小3.25%,如圖17所示.

    圖17 降低節(jié)圓振動頻率優(yōu)化結(jié)果的應(yīng)變能隨迭代步數(shù)變化趨勢對比Fig.17 Variation tendency comparison of strain energy of optimization results with lower frequency of umbrella vibration in different iteration steps

    總的來看,STSO方法的優(yōu)點是系統(tǒng)、準確、收斂快.它的局限性除了易引起網(wǎng)格畸變過大外,還在于將STSO方法應(yīng)用于輪盤優(yōu)化中時,并不能直接將優(yōu)化結(jié)果應(yīng)用于工程實踐中.這是由于STSO方法只能包含初始設(shè)計模型中的形狀優(yōu)化變量,對于STSO方法的優(yōu)化結(jié)果中新出現(xiàn)的拓撲構(gòu)型,依然需要設(shè)置相應(yīng)的形狀優(yōu)化變量,進行進一步的形狀和尺寸優(yōu)化等細化設(shè)計.因此,從這個方面看,STSO方法只是對應(yīng)于結(jié)構(gòu)的概念設(shè)計階段,如圖2(c)所示.

    從概念設(shè)計階段的結(jié)果到工程實際應(yīng)用的輪盤,需要再考慮的因素較多,比如滿足其他靜力學約束,降低某階次頻率同時、提高另外階次頻率,熱應(yīng)力約束等.本文僅僅是運用STSO方法考慮了振型對應(yīng)的頻率約束對同時優(yōu)化方法所得結(jié)果的影響,得到了幾種概念設(shè)計方案,在本文所用STSO方法得到的優(yōu)化結(jié)果的基礎(chǔ)上,還需要進一步的形狀和尺寸優(yōu)化.

    5 結(jié)論

    本文首先建立了STSO方法的數(shù)學模型,并推導了目標和約束函數(shù)的靈敏度,運用SQP算法進行求解.用板殼結(jié)構(gòu)優(yōu)化的例子對比分析了同時優(yōu)化和分步優(yōu)化結(jié)果.將同時優(yōu)化的方法用于輪盤動力學優(yōu)化,分析了不同振型時頻率約束對同時優(yōu)化方法所得結(jié)果的影響,并將之與單獨拓撲優(yōu)化結(jié)果進行比較,結(jié)果表明:

    1)拓撲和形狀同時優(yōu)化方法較分步優(yōu)化減少了人為參與的環(huán)節(jié),從優(yōu)化過程看,更為客觀科學.

    2)拓撲和形狀同時優(yōu)化考慮了形狀優(yōu)化變量和拓撲優(yōu)化變量之間的相互作用,回避了文獻中所提出的拓撲和形狀優(yōu)化先后順序問題,從整體看,更為系統(tǒng)化.

    3)將拓撲和形狀同時優(yōu)化的方法應(yīng)用于板殼結(jié)構(gòu)優(yōu)化和實體輪盤結(jié)構(gòu)動力學優(yōu)化時,結(jié)果表明,相較于分步拓撲優(yōu)化和單獨拓撲優(yōu)化,同時優(yōu)化方法的結(jié)果較分步優(yōu)化效果更好,且收斂更快.

    4)相較于輪盤單獨拓撲優(yōu)化方法結(jié)果,使用STSO方法得到的結(jié)果開孔較少甚至不開孔,更加符合現(xiàn)役發(fā)動機少有開孔輪盤的現(xiàn)狀.

    5)不同振型對應(yīng)的頻率約束不同時,得到了不同的優(yōu)化結(jié)果.降低節(jié)徑振動頻率,均得到了輻條式拓撲;提高節(jié)圓振動頻率,均得到了雙幅板式拓撲;降低節(jié)圓振動頻率,均得到了單幅板式拓撲.

    6)由于形狀優(yōu)化變量取值范圍選取不當,可能導致網(wǎng)格畸變過大而得到無效優(yōu)化解,此時需要重新選擇形狀優(yōu)化變量的取值范圍.

    References)

    [1] Deaton J D,Grandhi R V.A survey of structural and multidisciplinary continuum topology optimization:post 2000[J].Struct Multidisc Optim,2014,49(1):1-38.

    [2] Ansola R,Canales J,Tarrago J A,et al.An integrated approach for shape and topology optimization of shell structures[J].Computer & Structures.2002,80(5):449-458.

    [3] Hassani B,Tavakkoli S M,Ghasemnejad H.Simultaneous shape and topology optimizationof shell structures[J].Struct Multidisc Optim.2013,48(1):221-233.

    [4] 張衛(wèi)紅,楊軍剛,朱繼宏.壓力載荷下的結(jié)構(gòu)拓撲-形狀協(xié)同優(yōu)化[J].航空學報,2009,30(12):2335-2341.Zhang W H,Yang J G,Zhu J H.Simultaneous topology and shape optimization of pressure loaded structures[J].Acta Aeronautica et Astronautica Sinica,2009,30(12):2335-2341(in Chinese).

    [5] 劉超.渦輪盤結(jié)構(gòu)拓撲與形狀優(yōu)化方法研究[D].南京:南京航空航天大學,2010.Liu C.Topology and shape optimization method for turbine disk[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2010(in Chinese).

    [6] 宋健,溫衛(wèi)東,崔海濤,等.航空發(fā)動機多輻板風扇盤拓撲優(yōu)化與形狀優(yōu)化設(shè)計技術(shù)[J].推進技術(shù),2013,34(9):1188-1196.Song J,Wen W D,Cui H T,et al.Topology and shape optimization method for multi-web fan disk in aero-engine[J].Journal of Propulsion Technology,2013,34(9):1188-1196(in Chinese).

    [7] 李倫未,陸山.基于ANSYS的多輻板風扇盤結(jié)構(gòu)優(yōu)化設(shè)計技術(shù)[J].航空動力學報,2011,26(10):2245-2250.Li L W,Lu S.Structure optimum design techniques for multi-web fan disk based on ANSYS[J].Journal of Aerospace Power,2011,26(10):2245-2250(in Chinese).

    [8] 榮見華,唐國軍,羅銀燕,等.考慮位移要求的大型三維連續(xù)體結(jié)構(gòu)拓撲優(yōu)化數(shù)值技術(shù)研究[J].工程力學,2007,24(3):20-27.Rong J H,Tang G J,Luo Y Y,et al.A research on the numerical topology optimization technology of large three-dimensional continuum structures considering displacement requirements[J].Engineering Mechanics,2007,24(3):20-27(in Chinese).

    [9] 左孔天,陳立平,鐘毅芳,等.基于人工材料密度的新型拓撲優(yōu)化理論和算法研究[J].機械工程學報,2004,40(12):31-37.Zuo K T,Chen L P,Zhong Y F,et al.New theory and algorithm research about topology optimization based on artificial material density[J].Chinese Journal of Mechanical Engineering,2004,40(12):31-37(in Chinese).

    [10] 鄒春江,左孔天,向宇,等.基于SIMP方法微電容加速度計結(jié)構(gòu)固有頻率拓撲優(yōu)化[J].科學技術(shù)與工程,2011,11(29):7086-7091.Zhou C J,Zuo K T,Xiang Y,et al.Natural frequency topology optimization for structure of micro-capacitive accelerometer based on SIMP method[J].Science Technology and Engineering,2011,11(29):7086-7091(in Chinese).

    [11] 葉紅玲,沈靜嫻,隋允康.頻率約束的三維連續(xù)體結(jié)構(gòu)動力拓撲優(yōu)化設(shè)計[J].力學學報,2012,44(6):1037-1045.Ye H L,Shen J X,Sui Y K.Dynamic topological optimal design of three-dimensional continuum structures with frequencies constraints[J].Chinese Journal of Theoretical and Applied Mechanics,2012,44(6):1037-1045(in Chinese).

    [12] 張維聲.基于水平集方法的應(yīng)力相關(guān)拓撲優(yōu)化問題研究[D].大連:大連理工大學,2013.Zhang W S.Studies on stress-related topology optimization problem via level set approach[D].Dalian:Dalian University of Technology,2013(in Chinese).

    [13] 隋允康,葉紅玲,彭細榮.應(yīng)力約束全局化策略下的連續(xù)體結(jié)構(gòu)拓撲優(yōu)化[J].力學學報,2006,38(3):364-370.Sui Y K,Ye H L,Peng X R.Topological optimization of continuum structure under the strategy of globalization of stress constraints[J].Chinese Journal of Theoretical and Applied Mechanics,2006,38(3):364-370(in Chinese).

    [14] 榮見華,葛森,鄧果,等.基于位移和應(yīng)力靈敏度的結(jié)構(gòu)拓撲優(yōu)化設(shè)計[J].力學學報,2009,41(4):518-529.Rong J H,Ge S,Deng G,et al.Structure topological optimization based on displacement and stress sensitivity analysis[J].Chinese Journal of Theoretical and Applied Mechanics,2009,41(4):518-529(in Chinese).

    [15] Paris J,Navarria F,Colominas I,et al.Stress constraints sensitivity analysis in structural topology optimization[J].Computer Methodsin Applied Mechanicsand Engineering,2010,199(133):2110-2122.

    [16] Kai-Uwe B,Matthias F,F(xiàn)ernass D.Approximation of derivatives in semi-analytical structural optimization[J].Computers &Structures,2008,86(13):1404-1416.

    [17] De Boer H,Van Keulen F.Refined semi-analytical design sensitivities[J].International Journal of Solids and Structures,2000,37(46):6961-6980.

    [18] Van Keulen F,Haftka R T,Kim N H.Review of options for structural design sensitivity analysis.Part 1 linear systems[J].Computer Methods in Applied Mechanics and Engineering,2005,194(30):3213-3243.

    [19] Fletche R.Practical methods of optimization[M].New York:John Wiley and Sons,1987:304-319.

    [20] Nocedal J,Wright.S J.Numerical optimization[M].New York:Springer Series in Operations Research,1999:529-561.

    猜你喜歡
    輪盤形狀約束
    挖藕 假如悲傷有形狀……
    “碳中和”約束下的路徑選擇
    某型航空發(fā)動機鈦合金輪盤模擬疲勞試驗件設(shè)計
    約束離散KP方程族的完全Virasoro對稱
    你的形狀
    基于ANSYS的輪盤轉(zhuǎn)子模態(tài)影響因素分析
    看到的是什么形狀
    適當放手能讓孩子更好地自我約束
    人生十六七(2015年6期)2015-02-28 13:08:38
    不等式約束下AXA*=B的Hermite最小二乘解
    玩玩算算
    讀寫算(上)(2012年7期)2012-02-03 01:22:16
    免费av毛片视频| 精品第一国产精品| 亚洲欧洲精品一区二区精品久久久| 美女扒开内裤让男人捅视频| 免费在线观看视频国产中文字幕亚洲| 国产精品一区二区精品视频观看| 黑人巨大精品欧美一区二区蜜桃| 亚洲专区中文字幕在线| 亚洲人成网站在线播放欧美日韩| 国产精品久久久av美女十八| 免费在线观看完整版高清| 国产av精品麻豆| 色精品久久人妻99蜜桃| 精品日产1卡2卡| 9色porny在线观看| 欧美成人午夜精品| 国产野战对白在线观看| 成人国产一区最新在线观看| 丰满的人妻完整版| videosex国产| 桃红色精品国产亚洲av| √禁漫天堂资源中文www| 日日夜夜操网爽| 看黄色毛片网站| 中文亚洲av片在线观看爽| 日韩精品青青久久久久久| 久久国产精品男人的天堂亚洲| 日韩欧美三级三区| 18禁美女被吸乳视频| 国产熟女午夜一区二区三区| 欧洲精品卡2卡3卡4卡5卡区| 99久久99久久久精品蜜桃| 久久久久久免费高清国产稀缺| av天堂在线播放| 午夜福利,免费看| 亚洲性夜色夜夜综合| 精品福利永久在线观看| 色综合欧美亚洲国产小说| 国产人伦9x9x在线观看| 777久久人妻少妇嫩草av网站| 在线观看66精品国产| 国产亚洲精品第一综合不卡| 亚洲精品美女久久久久99蜜臀| 日韩人妻精品一区2区三区| av在线天堂中文字幕 | 日韩av在线大香蕉| 午夜免费观看网址| 日韩欧美一区视频在线观看| av视频免费观看在线观看| 巨乳人妻的诱惑在线观看| 人人妻,人人澡人人爽秒播| 黄网站色视频无遮挡免费观看| 老司机靠b影院| 麻豆成人av在线观看| 水蜜桃什么品种好| 老司机亚洲免费影院| 国产精品98久久久久久宅男小说| 岛国在线观看网站| 又黄又爽又免费观看的视频| 久久午夜综合久久蜜桃| 久久久久国产一级毛片高清牌| 免费在线观看亚洲国产| 成人黄色视频免费在线看| 99riav亚洲国产免费| 日韩欧美一区视频在线观看| 女人爽到高潮嗷嗷叫在线视频| 精品日产1卡2卡| 熟女少妇亚洲综合色aaa.| 久久伊人香网站| 精品第一国产精品| 性色av乱码一区二区三区2| a级片在线免费高清观看视频| 99re在线观看精品视频| 91成人精品电影| 亚洲,欧美精品.| 久久久精品欧美日韩精品| 超色免费av| 欧美日韩福利视频一区二区| 久久中文字幕人妻熟女| 视频在线观看一区二区三区| 国产精品成人在线| 久久亚洲真实| 精品无人区乱码1区二区| 亚洲第一欧美日韩一区二区三区| 亚洲男人的天堂狠狠| 免费久久久久久久精品成人欧美视频| 国产激情久久老熟女| 欧美大码av| 99国产精品99久久久久| 国产片内射在线| www.熟女人妻精品国产| 欧洲精品卡2卡3卡4卡5卡区| 人人妻人人添人人爽欧美一区卜| 91老司机精品| 国产熟女xx| 午夜精品国产一区二区电影| 国产又色又爽无遮挡免费看| 国产精品爽爽va在线观看网站 | 乱人伦中国视频| 一级a爱片免费观看的视频| 亚洲久久久国产精品| 视频区欧美日本亚洲| 免费人成视频x8x8入口观看| 欧美激情高清一区二区三区| 一级片免费观看大全| 久久精品91无色码中文字幕| 男人的好看免费观看在线视频 | 午夜91福利影院| 日日摸夜夜添夜夜添小说| 色老头精品视频在线观看| 变态另类成人亚洲欧美熟女 | 91老司机精品| 美女扒开内裤让男人捅视频| 搡老乐熟女国产| 亚洲午夜理论影院| 99国产精品99久久久久| 欧美精品亚洲一区二区| 夜夜看夜夜爽夜夜摸 | 国产一区在线观看成人免费| 国产三级黄色录像| 在线看a的网站| 啪啪无遮挡十八禁网站| 高清毛片免费观看视频网站 | 黄色片一级片一级黄色片| 天堂影院成人在线观看| 91精品国产国语对白视频| 欧美日韩瑟瑟在线播放| 777久久人妻少妇嫩草av网站| 热re99久久精品国产66热6| av欧美777| 男女床上黄色一级片免费看| 十八禁网站免费在线| 日韩高清综合在线| 1024视频免费在线观看| tocl精华| 中出人妻视频一区二区| 午夜精品久久久久久毛片777| 狂野欧美激情性xxxx| 国产激情欧美一区二区| 亚洲av成人一区二区三| 中文亚洲av片在线观看爽| 亚洲av成人不卡在线观看播放网| 新久久久久国产一级毛片| 亚洲av美国av| 成人免费观看视频高清| 精品久久久久久久久久免费视频 | 深夜精品福利| 99国产精品免费福利视频| 51午夜福利影视在线观看| 最近最新免费中文字幕在线| 中文欧美无线码| 中文字幕高清在线视频| 国产精品香港三级国产av潘金莲| 日韩欧美国产一区二区入口| 色哟哟哟哟哟哟| 午夜福利影视在线免费观看| 一本综合久久免费| 精品久久久久久久久久免费视频 | 黄色女人牲交| 美女高潮到喷水免费观看| 久久性视频一级片| xxxhd国产人妻xxx| 国产又色又爽无遮挡免费看| 亚洲精华国产精华精| 亚洲情色 制服丝袜| 变态另类成人亚洲欧美熟女 | 国产亚洲精品久久久久5区| 欧美日韩亚洲国产一区二区在线观看| 久久久久国产精品人妻aⅴ院| 亚洲精品一二三| 国产精品成人在线| 热re99久久国产66热| 不卡一级毛片| 日本a在线网址| 自拍欧美九色日韩亚洲蝌蚪91| 国产精品免费视频内射| 19禁男女啪啪无遮挡网站| 国产单亲对白刺激| 亚洲成人免费电影在线观看| 精品免费久久久久久久清纯| 如日韩欧美国产精品一区二区三区| 搡老乐熟女国产| 国产精品一区二区在线不卡| 国产亚洲精品久久久久5区| 国产精品偷伦视频观看了| 一级作爱视频免费观看| 久久精品国产99精品国产亚洲性色 | 丁香六月欧美| 国产精品免费视频内射| 乱人伦中国视频| 精品乱码久久久久久99久播| 免费人成视频x8x8入口观看| 成年女人毛片免费观看观看9| 十分钟在线观看高清视频www| 脱女人内裤的视频| 999久久久精品免费观看国产| 午夜福利一区二区在线看| av网站免费在线观看视频| 婷婷丁香在线五月| 色播在线永久视频| 99久久人妻综合| 欧美激情极品国产一区二区三区| 国产av又大| 黄片播放在线免费| 国产亚洲精品久久久久久毛片| 日本 av在线| 中文字幕另类日韩欧美亚洲嫩草| 精品久久久久久电影网| 99久久99久久久精品蜜桃| 国产精品影院久久| 99国产综合亚洲精品| 桃色一区二区三区在线观看| 老司机深夜福利视频在线观看| 精品国产超薄肉色丝袜足j| av在线天堂中文字幕 | 女人被躁到高潮嗷嗷叫费观| 午夜亚洲福利在线播放| 在线观看免费高清a一片| 久久精品aⅴ一区二区三区四区| 亚洲色图 男人天堂 中文字幕| 国产精品一区二区免费欧美| 午夜影院日韩av| 久久性视频一级片| 久久久久国产精品人妻aⅴ院| 国产一区二区三区综合在线观看| 大码成人一级视频| 一区二区三区激情视频| 久久亚洲精品不卡| 757午夜福利合集在线观看| 精品久久久久久久久久免费视频 | 两个人看的免费小视频| 9热在线视频观看99| 午夜免费成人在线视频| 成在线人永久免费视频| 亚洲国产欧美一区二区综合| 在线看a的网站| 久9热在线精品视频| 成人三级做爰电影| 精品高清国产在线一区| 亚洲国产毛片av蜜桃av| 1024视频免费在线观看| 91九色精品人成在线观看| 日韩欧美一区二区三区在线观看| 天堂中文最新版在线下载| 如日韩欧美国产精品一区二区三区| 丁香六月欧美| 欧美日韩黄片免| 身体一侧抽搐| 久久精品91无色码中文字幕| 两个人免费观看高清视频| 成在线人永久免费视频| 一个人观看的视频www高清免费观看 | 国产熟女xx| 欧美激情 高清一区二区三区| 国产欧美日韩综合在线一区二区| 午夜激情av网站| 亚洲成av片中文字幕在线观看| 校园春色视频在线观看| 日韩中文字幕欧美一区二区| 丝袜美腿诱惑在线| 国产精品久久电影中文字幕| 亚洲全国av大片| 满18在线观看网站| 精品国产国语对白av| 亚洲片人在线观看| 精品久久蜜臀av无| 国产免费男女视频| 国产精品爽爽va在线观看网站 | 少妇被粗大的猛进出69影院| 日本撒尿小便嘘嘘汇集6| 国产日韩一区二区三区精品不卡| 久久久国产欧美日韩av| 久9热在线精品视频| 亚洲自偷自拍图片 自拍| 日韩三级视频一区二区三区| 免费久久久久久久精品成人欧美视频| 久久久国产欧美日韩av| 一个人观看的视频www高清免费观看 | 日本五十路高清| 久久精品亚洲av国产电影网| 精品久久久精品久久久| 国产野战对白在线观看| 在线观看免费日韩欧美大片| 国产精品免费视频内射| 99精国产麻豆久久婷婷| 夜夜躁狠狠躁天天躁| 亚洲国产精品sss在线观看 | 日日爽夜夜爽网站| 日韩欧美一区视频在线观看| 国产精品亚洲av一区麻豆| 黄色 视频免费看| 一边摸一边抽搐一进一小说| 搡老熟女国产l中国老女人| 天堂动漫精品| 天堂影院成人在线观看| 另类亚洲欧美激情| 中出人妻视频一区二区| 久久国产乱子伦精品免费另类| 99香蕉大伊视频| 欧美日韩精品网址| 亚洲一码二码三码区别大吗| 搡老熟女国产l中国老女人| 黄色视频不卡| 高潮久久久久久久久久久不卡| 国产精品香港三级国产av潘金莲| 欧美成人免费av一区二区三区| 999久久久国产精品视频| 国产精品成人在线| 精品国产亚洲在线| 嫁个100分男人电影在线观看| 免费在线观看视频国产中文字幕亚洲| 99精国产麻豆久久婷婷| 久久精品国产综合久久久| www.999成人在线观看| 国产精品久久电影中文字幕| 男人舔女人下体高潮全视频| 国产单亲对白刺激| 亚洲自拍偷在线| 亚洲人成电影免费在线| 亚洲国产欧美日韩在线播放| 国产精品偷伦视频观看了| 成年版毛片免费区| 男女下面插进去视频免费观看| 国产熟女xx| 五月开心婷婷网| 亚洲情色 制服丝袜| 婷婷六月久久综合丁香| 成熟少妇高潮喷水视频| 99re在线观看精品视频| 久久草成人影院| 99国产极品粉嫩在线观看| 麻豆国产av国片精品| 午夜两性在线视频| 电影成人av| 757午夜福利合集在线观看| 国产精品自产拍在线观看55亚洲| 一级毛片精品| av有码第一页| 亚洲色图综合在线观看| 日韩欧美一区视频在线观看| 午夜久久久在线观看| 男人舔女人的私密视频| bbb黄色大片| 热re99久久精品国产66热6| 久久天躁狠狠躁夜夜2o2o| 又黄又粗又硬又大视频| 黄色视频不卡| 欧美最黄视频在线播放免费 | 亚洲精品国产一区二区精华液| 亚洲成人免费av在线播放| 很黄的视频免费| 91字幕亚洲| 真人一进一出gif抽搐免费| 亚洲人成网站在线播放欧美日韩| 免费人成视频x8x8入口观看| 欧美国产精品va在线观看不卡| 99国产综合亚洲精品| 亚洲国产欧美网| 成人免费观看视频高清| 99国产精品99久久久久| 国产精品日韩av在线免费观看 | 精品国产超薄肉色丝袜足j| 中国美女看黄片| 亚洲色图av天堂| 国产精品久久久久久人妻精品电影| 久久精品亚洲精品国产色婷小说| 成人亚洲精品一区在线观看| 国产激情久久老熟女| 水蜜桃什么品种好| 在线观看舔阴道视频| 在线视频色国产色| 精品免费久久久久久久清纯| 在线播放国产精品三级| 一边摸一边做爽爽视频免费| 一a级毛片在线观看| 亚洲色图av天堂| 亚洲伊人色综图| 国产有黄有色有爽视频| 精品久久久久久久久久免费视频 | 精品一品国产午夜福利视频| 久久中文字幕一级| 国产精品 欧美亚洲| 好男人电影高清在线观看| 久久亚洲真实| 一区二区三区国产精品乱码| 精品人妻在线不人妻| 国产av又大| 欧美国产精品va在线观看不卡| 黄片小视频在线播放| 国产成人av教育| 99国产精品一区二区三区| 美女高潮喷水抽搐中文字幕| 十分钟在线观看高清视频www| 亚洲欧美激情在线| 久久中文字幕人妻熟女| 黄色a级毛片大全视频| 国产不卡一卡二| 国产极品粉嫩免费观看在线| av免费在线观看网站| 欧美日韩中文字幕国产精品一区二区三区 | 亚洲 国产 在线| 日韩人妻精品一区2区三区| 午夜日韩欧美国产| 日本精品一区二区三区蜜桃| 亚洲免费av在线视频| 日韩中文字幕欧美一区二区| 高清黄色对白视频在线免费看| 免费观看精品视频网站| 天堂影院成人在线观看| 欧美精品啪啪一区二区三区| 少妇被粗大的猛进出69影院| 又大又爽又粗| 免费在线观看视频国产中文字幕亚洲| 日韩欧美一区视频在线观看| 18禁裸乳无遮挡免费网站照片 | 国产伦人伦偷精品视频| 狂野欧美激情性xxxx| 欧美成人午夜精品| 国产黄a三级三级三级人| e午夜精品久久久久久久| 精品少妇一区二区三区视频日本电影| 日韩视频一区二区在线观看| 国产片内射在线| 精品电影一区二区在线| av在线天堂中文字幕 | 欧美国产精品va在线观看不卡| 国产欧美日韩精品亚洲av| 国产亚洲欧美98| 精品国产亚洲在线| 国产国语露脸激情在线看| 精品国产超薄肉色丝袜足j| 午夜精品久久久久久毛片777| 大型黄色视频在线免费观看| 亚洲免费av在线视频| 亚洲欧美日韩无卡精品| 亚洲欧美精品综合久久99| av欧美777| www.熟女人妻精品国产| 久久中文看片网| 国产成人av激情在线播放| 免费高清视频大片| 久久久久久久久中文| 亚洲精品美女久久av网站| 国产片内射在线| 男人舔女人的私密视频| 乱人伦中国视频| 国产成人精品无人区| 一级毛片高清免费大全| 国产成年人精品一区二区 | 在线观看66精品国产| 国产一区二区三区在线臀色熟女 | 亚洲国产精品sss在线观看 | 99国产精品99久久久久| 深夜精品福利| 男女高潮啪啪啪动态图| 男女午夜视频在线观看| 中文字幕人妻熟女乱码| 麻豆成人av在线观看| 51午夜福利影视在线观看| 男男h啪啪无遮挡| 国产精品免费一区二区三区在线| 亚洲一区高清亚洲精品| 超碰97精品在线观看| 亚洲精品中文字幕一二三四区| 9色porny在线观看| av欧美777| 在线av久久热| 国产深夜福利视频在线观看| 国产成人影院久久av| 丁香六月欧美| 国产精品免费一区二区三区在线| 美女扒开内裤让男人捅视频| 嫩草影视91久久| 麻豆av在线久日| 色老头精品视频在线观看| 日韩欧美免费精品| 一级毛片高清免费大全| 人人妻人人爽人人添夜夜欢视频| 精品久久久久久成人av| 99国产极品粉嫩在线观看| 99久久国产精品久久久| 亚洲av成人不卡在线观看播放网| 色精品久久人妻99蜜桃| 一边摸一边做爽爽视频免费| 日本一区二区免费在线视频| 欧美在线黄色| www.精华液| 村上凉子中文字幕在线| 十分钟在线观看高清视频www| 91在线观看av| 91字幕亚洲| 后天国语完整版免费观看| 亚洲三区欧美一区| 免费搜索国产男女视频| 激情在线观看视频在线高清| 狠狠狠狠99中文字幕| 国产激情欧美一区二区| 深夜精品福利| 亚洲精品国产色婷婷电影| 精品久久久久久,| 国产麻豆69| 精品少妇一区二区三区视频日本电影| 国产精品 国内视频| 免费观看人在逋| 99国产精品一区二区三区| 性色av乱码一区二区三区2| 美国免费a级毛片| av在线天堂中文字幕 | 免费在线观看亚洲国产| 色综合婷婷激情| 亚洲精品国产精品久久久不卡| 高潮久久久久久久久久久不卡| 亚洲精品美女久久久久99蜜臀| 亚洲欧美精品综合久久99| 天天躁夜夜躁狠狠躁躁| 两性夫妻黄色片| 国产视频一区二区在线看| 在线看a的网站| 两性午夜刺激爽爽歪歪视频在线观看 | 国产精品av久久久久免费| 黄色丝袜av网址大全| 国产一区在线观看成人免费| 十分钟在线观看高清视频www| 夜夜看夜夜爽夜夜摸 | 香蕉国产在线看| 亚洲成人免费电影在线观看| 一级,二级,三级黄色视频| 欧美亚洲日本最大视频资源| 波多野结衣一区麻豆| 999精品在线视频| av福利片在线| 男人舔女人下体高潮全视频| 黄色怎么调成土黄色| 91精品三级在线观看| 高潮久久久久久久久久久不卡| 视频区欧美日本亚洲| 9热在线视频观看99| 午夜精品久久久久久毛片777| 国产精品久久久av美女十八| 我的亚洲天堂| 午夜亚洲福利在线播放| 热re99久久精品国产66热6| 桃色一区二区三区在线观看| 久久天堂一区二区三区四区| 精品国产国语对白av| 麻豆av在线久日| 神马国产精品三级电影在线观看 | 国产av一区二区精品久久| 国产成人免费无遮挡视频| 日日干狠狠操夜夜爽| 999精品在线视频| 99久久综合精品五月天人人| 国产精品免费一区二区三区在线| 女人爽到高潮嗷嗷叫在线视频| 国产成人影院久久av| 黑人巨大精品欧美一区二区蜜桃| 日韩欧美在线二视频| 国产免费现黄频在线看| 国产精品99久久99久久久不卡| 男人操女人黄网站| 在线观看一区二区三区激情| 脱女人内裤的视频| 国产成人精品久久二区二区91| 女人被躁到高潮嗷嗷叫费观| 精品日产1卡2卡| 国产一区二区激情短视频| 老熟妇乱子伦视频在线观看| bbb黄色大片| 男女做爰动态图高潮gif福利片 | 精品高清国产在线一区| 午夜福利在线观看吧| 波多野结衣高清无吗| 搡老乐熟女国产| av在线天堂中文字幕 | 交换朋友夫妻互换小说| 国产精品综合久久久久久久免费 | 黄色a级毛片大全视频| 日韩视频一区二区在线观看| 亚洲熟女毛片儿| 操美女的视频在线观看| 黄片小视频在线播放| 无限看片的www在线观看| 欧美色视频一区免费| av在线播放免费不卡| 99国产精品99久久久久| a级毛片在线看网站| 午夜免费成人在线视频| 欧美久久黑人一区二区| 中文字幕另类日韩欧美亚洲嫩草| 国产无遮挡羞羞视频在线观看| 欧美成人免费av一区二区三区| 19禁男女啪啪无遮挡网站| 国产精品99久久99久久久不卡| 国产精品美女特级片免费视频播放器 | 免费看a级黄色片| 久久久国产成人精品二区 | 久久这里只有精品19| 精品欧美一区二区三区在线| 欧美av亚洲av综合av国产av| www.熟女人妻精品国产| 宅男免费午夜| 欧美大码av| 美女高潮喷水抽搐中文字幕| 国产精品一区二区在线不卡| 国产一卡二卡三卡精品| 久久午夜综合久久蜜桃| 国产一区二区三区综合在线观看| 日韩大尺度精品在线看网址 | 男女下面进入的视频免费午夜 | 黄片大片在线免费观看| 999久久久精品免费观看国产| 午夜精品在线福利| 日韩免费av在线播放| 免费不卡黄色视频| 日韩精品免费视频一区二区三区|