亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        金屬粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法研究

        2025-08-07 00:00:00李海東史曉鳴張橫
        關(guān)鍵詞:阻尼器粉末顆粒

        中圖分類號(hào):TJ760.3 文獻(xiàn)標(biāo)志碼:A

        Topology optimization design method of metal powder particle damping composite rudder structure

        LI Haidong1,SHI Xiaoming1,ZHANG Heng2 (1.Shanghai Electro-Mechanical Engineering Institute,Shanghai2Ol1O9,China; 2.SchoolofMechanicalEngineering, University of Shanghai for Science and Technology, Shanghai 2Ooo93, China)

        Abstract: The mechanical characteristics of the rudder structure have a critical impact on the flight performance of the aircraft, and these characteristics are governed by the structure’ s geometric configuration. A topology optimization method for designing the metal powder particle damping composite rudder structure was proposed. The load-bearing function was satisfied by optimizing the stiffners of the rudder structure and the vibration suppression was achieved by optimizing the position distribution of the powder particle dampers. The internal stiffener distribution of the rudder structure was determined via the adaptive growth method, while the equivalent performance of the powder particle damper was evaluated using the homogenization method. Based on the homogenized parameters, the position distribution of the powder particle dampers within the rudder was optimized to realize a collaborative design that concurrently addressed both load-bearing and vibration suppression functions. The simulation results show that, compared to the original rudder structure, the static mechanical performance of the optimized metal powder particle damping composite rudder structure is improved by 10.35% , the first natural frequency is increased by 22.10% , the frequency response at the first natural frequency is reduced by 22.88% , and the flutter speed is increased by more than 10% when the Mach number is 4, 5, 6. The proposed design method provides a new solution for the design of lightweight, high-rigidity structures with vibration suppression functions.

        Keywords: topology optimization; stiffener distribution design; powder particle damper; rudder structure

        舵面結(jié)構(gòu)是飛行器飛行過程中實(shí)現(xiàn)高速飛行和機(jī)動(dòng)控制的核心部件之一,不僅要求其重量輕,還要求其具有較好的承載能力和抗振性能[1]。目前的舵面結(jié)構(gòu)多采用蒙皮-骨架結(jié)構(gòu)形式,蒙皮與骨架通過鉚釘、螺釘或釬焊連接在一起,骨架一般是輻射梁式或網(wǎng)格式等相對(duì)規(guī)則的經(jīng)驗(yàn)設(shè)計(jì)構(gòu)型[2]。然而,在日益嚴(yán)苛的飛行載荷作用下,這種經(jīng)驗(yàn)設(shè)計(jì)的結(jié)構(gòu)構(gòu)型無法滿足高速飛行器的高剛度設(shè)計(jì)要求。與此同時(shí),如何有效提升飛行器的顫振邊界也成為亟待攻克的技術(shù)難題[3]

        輕量化設(shè)計(jì)方面,近年來,隨著結(jié)構(gòu)拓?fù)鋬?yōu)化理論的不斷發(fā)展,國內(nèi)外學(xué)者積極開展該理論在飛行器領(lǐng)域的研究與應(yīng)用,在飛行器結(jié)構(gòu)設(shè)計(jì)制造一體化方面取得了眾多具有代表性的研究成果,推動(dòng)了飛行器領(lǐng)域相關(guān)技術(shù)的飛速發(fā)展。王昕江等[4提出了一種基于拓?fù)鋬?yōu)化的舵面結(jié)構(gòu)輕量化設(shè)計(jì)方法,并驗(yàn)證了靜氣動(dòng)彈性變形與強(qiáng)度。朱繼宏研究團(tuán)隊(duì)[5將舵面作為獨(dú)立部件,采用拓?fù)鋬?yōu)化獲取主傳力路徑,提取結(jié)構(gòu)特征,通過尺寸參數(shù)優(yōu)化對(duì)重構(gòu)模型進(jìn)行詳細(xì)設(shè)計(jì),實(shí)現(xiàn)結(jié)構(gòu)減重。丁曉紅研究團(tuán)隊(duì)[7-9]結(jié)合仿生優(yōu)化思想,通過研究自然界分支系統(tǒng)的形成機(jī)理,提出了一種結(jié)構(gòu)加筋分布優(yōu)化設(shè)計(jì)自適應(yīng)成長(zhǎng)法,并應(yīng)用于高速飛行器舵面結(jié)構(gòu)設(shè)計(jì)中,得到樹狀分布的骨架宏觀構(gòu)型,在實(shí)現(xiàn)1階固有頻率提升的同時(shí),達(dá)到輕量化設(shè)計(jì)的目的[10]。但是,目前針對(duì)舵類結(jié)構(gòu)優(yōu)化設(shè)計(jì)的研究多以剛度為目標(biāo),未考慮結(jié)構(gòu)振動(dòng)設(shè)計(jì)問題

        針對(duì)提升顫振性能的設(shè)計(jì),常用的手段有增加結(jié)構(gòu)剛度、施加配重、增加阻尼等方式,其核心思想是通過控制空氣動(dòng)力、彈性力和慣性力的耦合作用,減小結(jié)構(gòu)發(fā)生自激振動(dòng)的振幅,提高自激振動(dòng)趨于發(fā)散的顫振臨界速度,達(dá)到提高顫振邊界的目的。然而,施加配重會(huì)導(dǎo)致結(jié)構(gòu)質(zhì)量增加,與輕量化設(shè)計(jì)背道而馳。因此,對(duì)于提升顫振性能而言,如何在提高輕質(zhì)結(jié)構(gòu)剛度的同時(shí)提高結(jié)構(gòu)的阻尼,顯得尤為關(guān)鍵。

        顆粒阻尼器通過顆粒與顆?;蝾w粒與壁面之間的碰撞和摩擦耗能實(shí)現(xiàn)結(jié)構(gòu)振動(dòng)控制,具有耐久性好、對(duì)環(huán)境溫度不敏感、可在多個(gè)方向和寬頻域內(nèi)使用等優(yōu)點(diǎn)。因此,顆粒阻尼器在機(jī)械工程、建筑工程和高端裝備制造等領(lǐng)域得到了廣泛應(yīng)用[1。近年來,將顆粒阻尼器與金屬增材制造工藝相結(jié)合形成了一種創(chuàng)新技術(shù),即在增材制造過程中特意保留結(jié)構(gòu)內(nèi)部未熔融的粉末,由此得到一種新型的金屬粉末顆粒阻尼器(以下簡(jiǎn)稱粉末顆粒阻尼器)[12]。該阻尼器不僅兼具高剛度和高阻尼特性,還能實(shí)現(xiàn)與主結(jié)構(gòu)的一體化成型。Scott-Emuakpor等[13]在懸梁臂結(jié)構(gòu)中設(shè)置空腔,在增材制造過程中,將金屬粉末留置于空腔內(nèi),形成粉末顆粒阻尼器。研究表明,含有粉末顆粒阻尼器的懸臂梁振動(dòng)性能與完全熔融的懸臂梁相比,在僅有 4% 未熔融金屬粉末的情況下,結(jié)構(gòu)振動(dòng)減小了 88%~95% 。Scott-Emuakpor團(tuán)隊(duì)還對(duì)粉末顆粒阻尼器進(jìn)行了更為深入的研究,他們認(rèn)為:增材制造方式以及掃描路徑等因素對(duì)粉末顆粒阻尼器的阻尼性能影響較小[14];彎曲阻尼器上的表面結(jié)構(gòu)比平整阻尼器上的表面結(jié)構(gòu)具有更高的阻尼性能[15];粉末顆粒阻尼器的阻尼性能與粉末空腔體積之間有明顯的相關(guān)關(guān)系[等。Niedermeyer 等[17]對(duì)粉末顆粒阻尼器在渦輪葉片上的應(yīng)用進(jìn)行了研究,結(jié)果表明,添加粉末顆粒阻尼器的葉片相較原始葉片具有更高的阻尼性能。Ehlers等[18]研究了激光打印方向、腔體尺寸等對(duì)結(jié)構(gòu)阻尼性能的影響。Celli等[19]通過實(shí)驗(yàn)證明添加粉末顆粒阻尼器的結(jié)構(gòu)在實(shí)現(xiàn)有效減振的同時(shí)不會(huì)顯著降低結(jié)構(gòu)的疲勞壽命?,F(xiàn)有研究多是對(duì)粉末顆粒阻尼器的阻尼性能進(jìn)行分析,并未涉及粉末顆粒阻尼器在結(jié)構(gòu)中分布位置的設(shè)計(jì)。

        為了提升舵面結(jié)構(gòu)的承載和顫振性能,本文將舵面結(jié)構(gòu)骨架分布設(shè)計(jì)與粉末顆粒阻尼器分布位置設(shè)計(jì)相結(jié)合。通過基于自適應(yīng)成長(zhǎng)法的高剛度骨架結(jié)構(gòu)設(shè)計(jì)方法,在實(shí)現(xiàn)結(jié)構(gòu)輕量化目標(biāo)的同時(shí),有效提高結(jié)構(gòu)的承載能力。在此基礎(chǔ)上,對(duì)粉末顆粒阻尼器分布位置進(jìn)行設(shè)計(jì),提高舵面結(jié)構(gòu)阻尼性能,實(shí)現(xiàn)抑振效果,提升舵面結(jié)構(gòu)的顫振邊界。

        1 粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)設(shè)計(jì)流程

        針對(duì)舵面結(jié)構(gòu)的承載和抑振協(xié)同設(shè)計(jì)問題,提出一種粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法。通過對(duì)舵面結(jié)構(gòu)內(nèi)部骨架分布和粉末顆粒阻尼器位置的設(shè)計(jì),實(shí)現(xiàn)輕質(zhì)高承載且具有高顫振邊界的舵面結(jié)構(gòu)設(shè)計(jì)目標(biāo),設(shè)計(jì)流程如圖1所示。具體設(shè)計(jì)步驟如下:

        a.基于自適應(yīng)成長(zhǎng)法的骨架結(jié)構(gòu)設(shè)計(jì)

        舵面結(jié)構(gòu)的內(nèi)部骨架分布直接影響結(jié)構(gòu)的整體性能?;谧赃m應(yīng)成長(zhǎng)法,以剛度最大為目標(biāo),以重心位置為約束,對(duì)舵面結(jié)構(gòu)的內(nèi)部骨架分布進(jìn)行優(yōu)化,以實(shí)現(xiàn)輕質(zhì)高承載骨架結(jié)構(gòu)的設(shè)計(jì)

        b.等效力學(xué)性能計(jì)算

        由于粉末顆粒阻尼器中金屬粉末顆粒直徑非常小,其平均直徑一般為 30μm ,直接通過離散元方法對(duì)整體結(jié)構(gòu)進(jìn)行設(shè)計(jì)的計(jì)算代價(jià)極高。因此,基于均勻化方法,計(jì)算粉末顆粒阻尼器等效力學(xué)性能,包括剛度、質(zhì)量和阻尼性能參數(shù),并在此基礎(chǔ)上進(jìn)行后續(xù)的阻尼器分布優(yōu)化設(shè)計(jì)。

        c.最優(yōu)分布設(shè)計(jì)

        采用密度法,以粉末顆粒阻尼器的偽密度為設(shè)計(jì)變量,對(duì)其等效力學(xué)性能進(jìn)行插值計(jì)算,構(gòu)建基于密度法的拓?fù)鋬?yōu)化模型。通過優(yōu)化求解,實(shí)現(xiàn)粉末顆粒阻尼器在舵面結(jié)構(gòu)中的最優(yōu)分布,進(jìn)一步提升結(jié)構(gòu)阻尼性能。

        d.工程化設(shè)計(jì)

        基于優(yōu)化結(jié)果,對(duì)舵面結(jié)構(gòu)進(jìn)行幾何重構(gòu),建立粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)三維幾何模型

        e.性能驗(yàn)證

        對(duì)優(yōu)化后舵面結(jié)構(gòu)靜力學(xué)、模態(tài)、幅頻響應(yīng)和顫振等性能進(jìn)行仿真分析,驗(yàn)證優(yōu)化后粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)的設(shè)計(jì)效果。

        圖1粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)設(shè)計(jì)流程Fig.1Design process of the powder particle damping composite rudder structure

        2 粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)優(yōu)化設(shè)計(jì)算例

        2.1 基于自適應(yīng)成長(zhǎng)法的舵面結(jié)構(gòu)骨架設(shè)計(jì)

        經(jīng)過億萬年的進(jìn)化,植物等自然生物分支系統(tǒng)發(fā)展出了復(fù)雜的自適應(yīng)機(jī)制,大多數(shù)分支系統(tǒng)都表現(xiàn)出最優(yōu)的分支結(jié)構(gòu)特征。這些特征使它們能夠適應(yīng)外部環(huán)境載荷,包括自重、附著物載荷(積雪、附冰、果實(shí)等)以及風(fēng)載等。丁曉紅研究團(tuán)隊(duì)[7,20]通過研究自然生物分支系統(tǒng)中的生長(zhǎng)機(jī)制,提出了自適應(yīng)成長(zhǎng)法。類比植物根系的生長(zhǎng)過程,主根總是從種子出發(fā),沿著能使其整體功能最優(yōu)的方向成長(zhǎng)。基于自適應(yīng)成長(zhǎng)法的結(jié)構(gòu)加筋分布優(yōu)化設(shè)計(jì),采用基結(jié)構(gòu)法建立優(yōu)化幾何模型,如圖2(a)所示。加強(qiáng)筋從選定的種子點(diǎn)開始根據(jù)一定的規(guī)則進(jìn)行成長(zhǎng)—分歧—退化。

        a.成長(zhǎng)

        生物分支系統(tǒng)的生長(zhǎng)受到諸如水分或陽光等物質(zhì)的引導(dǎo),類比這種機(jī)制,自適應(yīng)成長(zhǎng)法中,加強(qiáng)筋從種子點(diǎn)開始,沿著對(duì)優(yōu)化目標(biāo)靈敏度最大的方向進(jìn)行成長(zhǎng),如圖2(b)所示。由于靈敏度不同,加強(qiáng)筋沿不同方向成長(zhǎng)出具有不同厚度的新加強(qiáng)筋,如圖2(b)中1、2、3所示的3個(gè)加強(qiáng)筋。

        b.分歧

        當(dāng)加強(qiáng)筋成長(zhǎng)到一定厚度時(shí),加強(qiáng)筋具有分歧能力,如圖2(b)中加強(qiáng)筋1和2具有分歧能力,加強(qiáng)筋繼續(xù)分歧成長(zhǎng),圖2(c)為加強(qiáng)筋1和2的分歧成長(zhǎng)。

        圖2自適應(yīng)成長(zhǎng)過程Fig.2Processof adaptive growth

        c.退化

        在分歧成長(zhǎng)過程中,如果已生成的加強(qiáng)筋對(duì)設(shè)計(jì)目標(biāo)的靈敏度變小,加強(qiáng)筋會(huì)退化甚至消失。如圖2(c),加強(qiáng)筋3由于對(duì)設(shè)計(jì)目標(biāo)貢獻(xiàn)小,最終退化消失。

        上述成長(zhǎng)—分歧—退化過程反復(fù)進(jìn)行,直到滿足收斂條件,得到最優(yōu)的加強(qiáng)筋分布形態(tài)。

        舵面結(jié)構(gòu)為一典型的薄壁三明治結(jié)構(gòu),本研究通過平板殼單元模擬蒙皮和內(nèi)部加強(qiáng)筋骨架,采用八節(jié)點(diǎn)六面體線性完全積分單元模擬粉末顆粒阻尼器。平板殼單元由四邊形膜元與剪應(yīng)變混合插值板元耦合構(gòu)造而成,詳見文獻(xiàn)[21]。原型舵面結(jié)構(gòu)的外形尺寸如圖3(a)所示。為降低網(wǎng)格質(zhì)量對(duì)自適應(yīng)成長(zhǎng)法結(jié)果的影響,長(zhǎng)寬方向上采用長(zhǎng)寬比相近的四邊形單元,邊緣位置采用三角形單元,最終建立用于拓?fù)鋬?yōu)化的舵面基結(jié)構(gòu)模型,如圖3(b)所示。采用自適應(yīng)成長(zhǎng)法對(duì)舵面結(jié)構(gòu)內(nèi)部骨架分布進(jìn)行設(shè)計(jì),種子點(diǎn)選取為舵軸位置,如圖3(b)所示。在迭代中以筋板厚度為設(shè)計(jì)變量,通過厚度變化表示成長(zhǎng)、分歧和退化,實(shí)現(xiàn)骨架分布設(shè)計(jì)。

        舵面內(nèi)部骨架的分布形態(tài),不僅對(duì)其結(jié)構(gòu)的剛度起著關(guān)鍵作用,還直接影響著舵面結(jié)構(gòu)的重心位置,而重心位置又與舵面結(jié)構(gòu)的操控性和抗顫振能力密切相關(guān)。通常來說,舵面結(jié)構(gòu)的重心越靠近舵軸,結(jié)構(gòu)的操控性越好,同時(shí)抗顫振能力也越強(qiáng)。因此,本研究在舵面的骨架分布設(shè)計(jì)中考慮重心位置對(duì)設(shè)計(jì)結(jié)果的影響,最終建立以靜柔度為目標(biāo)、以重心位置為約束的舵面骨架分布優(yōu)化設(shè)計(jì)模型:

        式中: F 為結(jié)構(gòu)的靜柔度; U 為位移矩陣; K 為總剛度矩陣; Vs 為加強(qiáng)筋總體積; fs 為加強(qiáng)筋體積分?jǐn)?shù); ′ 為設(shè)計(jì)變量; V0 為舵面結(jié)構(gòu)的總體積;Ccog 為舵面結(jié)構(gòu)的重心位置; Xcog 為重心約束位置坐標(biāo); ti 表示第 i 個(gè)加強(qiáng)筋厚度; tmax,tmin 分別為加強(qiáng)筋厚度的最大值和最小值; n 為加強(qiáng)筋總數(shù);下標(biāo)s表示加強(qiáng)筋。

        設(shè)計(jì)目標(biāo)靜柔度對(duì)設(shè)計(jì)變量 ti 的靈敏度為

        Fig.3Original rudderstructureand basic structureof rudder

        式中: ks,i 為第 i 個(gè)加強(qiáng)筋的剛度矩陣; us,i 為第i 個(gè)加強(qiáng)筋的位移場(chǎng)向量。

        基于自適應(yīng)成長(zhǎng)法和式(1)所示的優(yōu)化模型,以移動(dòng)漸近線法(methodofmovingasymptotes,MMA)為尋優(yōu)算法,通過迭代更新準(zhǔn)則對(duì)結(jié)構(gòu)分布進(jìn)行尋優(yōu),具體流程如下:

        步驟1建立基結(jié)構(gòu),設(shè)置設(shè)計(jì)域。

        步驟2選取種子點(diǎn),設(shè)置優(yōu)化參數(shù)。參考結(jié)構(gòu)的載荷和支撐情況,選取基結(jié)構(gòu)上若干節(jié)點(diǎn)作為種子點(diǎn),給定加強(qiáng)筋的初始厚度、分歧臨界值、體積約束因子、收斂容差以及最大迭代次數(shù)。

        步驟3靈敏度分析。對(duì)結(jié)構(gòu)進(jìn)行有限元分析,通過解析法計(jì)算活動(dòng)加強(qiáng)筋的靈敏度。

        圖3原型舵面結(jié)構(gòu)及基結(jié)構(gòu)模型"

        步驟4加強(qiáng)筋成長(zhǎng)、分歧和退化。如果滿足加強(qiáng)筋成長(zhǎng)條件,通過MMA算法更新加強(qiáng)筋厚度;如果厚度達(dá)到分歧臨界值,則表明該加強(qiáng)筋具有分歧能力;如果滿足退化條件,則該加強(qiáng)筋退化。

        步驟5收斂判斷。計(jì)算目標(biāo)函數(shù),如果滿足收斂條件,則優(yōu)化結(jié)束,否則返回步驟3,繼續(xù)迭代,直至收斂條件得到滿足。

        基于圖3所示的基結(jié)構(gòu),對(duì)舵面結(jié)構(gòu)骨架的分布形態(tài)進(jìn)行優(yōu)化設(shè)計(jì)。相關(guān)參數(shù)取值如下:蒙皮厚度為 0.5mm ,加強(qiáng)筋厚度最大為 9mm ,最小為 0.1mm ,分歧厚度為 0.3mm ,過濾半徑為7mm 。舵軸處作為種子點(diǎn),加強(qiáng)筋結(jié)構(gòu)體積約束為 30% 。采用自適應(yīng)成長(zhǎng)法優(yōu)化設(shè)計(jì)后的舵面結(jié)構(gòu)骨架如圖4所示,呈現(xiàn)為典型的樹狀分支形態(tài)。該形態(tài)使得舵面結(jié)構(gòu)具有較好的承載性能,實(shí)現(xiàn)了輕質(zhì)高剛骨架結(jié)構(gòu)設(shè)計(jì)的目標(biāo)。

        2.2 基于均勻化法的粉末顆粒阻尼器等效性能計(jì)算

        如圖5(a)所示,在增材制造過程中,通過在結(jié)構(gòu)內(nèi)部設(shè)置空腔,在每一次鋪粉過程中將未熔融的金屬粉末逐層堆積于空腔結(jié)構(gòu)內(nèi)部,如此往復(fù),可在結(jié)構(gòu)內(nèi)部直接嵌入生成一種具備高填充率且顆粒材料與主結(jié)構(gòu)一致的顆粒阻尼器——金屬粉末顆粒阻尼器。該阻尼器的阻尼機(jī)制與宏觀的顆粒阻尼器類似,圖5(b)為顆粒阻尼器的工作原理,結(jié)構(gòu)體在發(fā)生振動(dòng)時(shí),通過顆粒與顆粒之間以及顆粒與壁面之間的碰撞和摩擦耗能實(shí)現(xiàn)減振。

        為了實(shí)現(xiàn)粉末顆粒阻尼器位置的優(yōu)化設(shè)計(jì),基于均勻化法,對(duì)阻尼器單元的等效性能進(jìn)行計(jì)算,然后基于密度法對(duì)等效性能進(jìn)行插值,通過對(duì)阻尼器偽密度的0-1分布設(shè)計(jì),實(shí)現(xiàn)阻尼器位置分布的優(yōu)化。圖6為粉末顆粒阻尼器等效性能計(jì)算過程示意圖,包括六面體單元?jiǎng)偠染仃嚒①|(zhì)量矩陣以及阻尼矩陣的等效構(gòu)建過程。

        a.等效剛度矩陣的計(jì)算

        基于均勻化法求解得到粉末顆粒阻尼器單元結(jié)構(gòu)的等效彈性本構(gòu)矩陣 DH ,其第 r 行和第 s 列的 DrsH

        圖5增材制造方式與顆粒阻尼器減振原理

        式中:上標(biāo)H表示經(jīng)過均勻化法得到的等效性能; 為粉末顆粒阻尼器單元域; kp,e 為粉末顆粒阻尼器單元?jiǎng)澐志W(wǎng)格后對(duì)應(yīng)的第 e 個(gè)單元的剛度矩陣; χe0 為通過施加單位應(yīng)變后計(jì)算得到的節(jié)點(diǎn)位移場(chǎng)向量; χe 為通過施加全局單位應(yīng)變后計(jì)算得到的節(jié)點(diǎn)位移場(chǎng)向量; r 和 s 表示本構(gòu)矩陣D中第r行和第s列,r=1,2,3,4,5,6,s=1,2,3,4,5,6;下標(biāo)p表示粉末顆粒阻尼器

        圖6金屬粉末顆粒阻尼器等效性能計(jì)算過程Fig.6Equivalent performance calculation process for the matal powder particle dampers

        對(duì)粉末顆粒阻尼器單元結(jié)構(gòu)進(jìn)行網(wǎng)格劃分,以宏觀體積為 10mm×10mm×10mm 、壁厚為 1mm 的單元顆粒阻尼器為例,如果以 1mm 大小對(duì)其進(jìn)行網(wǎng)格劃分,將阻尼器劃分為 10×10×10 的空間體素矩陣。體素矩陣中1元素代表實(shí)體,0元素代表空腔中的粉末,體素矩陣沿任意方向的每一層信息如圖7所示。

        經(jīng)過均勻化法求解后得到粉末顆粒阻尼器等效彈性本構(gòu)矩陣 DH ,則粉末顆粒阻尼器單元的等效剛度矩陣 kp

        圖7體素矩陣信息Fig.7Voxel matrix information

        式中: B 為全局坐標(biāo)系下的應(yīng)力-應(yīng)變矩陣; B* 為局部坐標(biāo)系 5-17-5 下的應(yīng)力-應(yīng)變矩陣; J 為坐標(biāo)變換的雅克比矩陣。

        b.等效質(zhì)量矩陣的計(jì)算

        粉末顆粒阻尼器單元的等效質(zhì)量矩陣 mp 采用 協(xié)調(diào)質(zhì)量矩陣來計(jì)算,即

        式中: N 為單元形函數(shù)矩陣; ρH 為粉末顆粒阻尼器單元的等效密度; Re 為圖7中體素矩陣的元素總數(shù); R1 為體素矩陣中1元素的個(gè)數(shù),1元素表示鈦合金實(shí)體,對(duì)應(yīng)的材料密度為 4.44g/cm3 ; R0 為體素矩陣中0元素的個(gè)數(shù),0元素表示鈦合金粉末,對(duì)應(yīng)的材料密度為 2.62g/cm3

        c.等效阻尼矩陣的計(jì)算

        由于粉末顆粒阻尼器的特殊性,其阻尼能力由碰撞產(chǎn)生,阻尼能力主要在節(jié)點(diǎn)位移方向上體現(xiàn),在剪切方向上不體現(xiàn)。本研究通過構(gòu)建集中質(zhì)量矩陣的思想來構(gòu)造單元阻尼矩陣,即阻尼矩陣對(duì)角線元素不為0,非對(duì)角線元素為0。粉末顆粒阻尼器單元的等效阻尼矩陣 cp 中的元素值為

        式中: cjH 為第 j 個(gè)粉末顆粒阻尼器單元的阻尼系數(shù); ck 為參考粉末顆粒阻尼器的阻尼系數(shù); Vk 為空腔體積; Vj 為第 j 個(gè)粉末顆粒阻尼器單元內(nèi)部粉末的體積; 、 和分別為單元在 x,y 和 z 這3個(gè)方向上的長(zhǎng)度。

        d.非規(guī)則六面體單元的等效性能計(jì)算

        如圖6所示,實(shí)際的舵面結(jié)構(gòu)中,由于其前緣、后緣和尖弦等部位較薄,而在舵軸處較厚,并不能用均勻的立方體單元進(jìn)行網(wǎng)格劃分。對(duì)于形狀各異的六面體單元模型,需按上述方法對(duì)具有特殊形狀的六面體單元進(jìn)行等效性能計(jì)算。對(duì)每個(gè)單元的平均長(zhǎng)度進(jìn)行分析后發(fā)現(xiàn),舵面結(jié)構(gòu)模型等參單元具有以下特點(diǎn):

        (a)六面體單元在 z 軸投影下近似為標(biāo)準(zhǔn)的正方形;(b)六面體單元關(guān)于 xy 平面對(duì)稱;(c)六面體單元在 z 方向上尺寸差異較小,可近似為橫置的扁平長(zhǎng)方體。

        舵面蒙皮及粉末顆粒阻尼器壁厚設(shè)置為0.5mm ,構(gòu)建體素矩陣,體素矩陣 x 、y向的單元長(zhǎng)度設(shè)置為 0.5mm ,單層體素矩陣為方陣。體素矩陣 z 向的層數(shù) Nz 定義為 z 向平均長(zhǎng)度與 0.5mm 相除的商,即

        式中, lt;gt; 表示求整運(yùn)算。

        根據(jù)式(3),采用均勻化法確定粉末顆粒阻尼器單元的等效彈性本構(gòu)矩陣,進(jìn)而根據(jù)式(4)確定單元的等效剛度矩陣。

        文獻(xiàn)[22]理論測(cè)試得到的邊長(zhǎng)為 10mm 、壁厚為 1mm 的正方體單元顆粒阻尼器的阻尼系數(shù) ck 為 0.007N?s/mm ,其內(nèi)部粉末空腔體積 Vk 為512mm3 ?;诖耍_定不同空腔體積的粉末顆粒阻尼器對(duì)應(yīng)的等效阻尼系數(shù)后,由式 (7)~(9) 確定粉末顆粒阻尼器單元等效阻尼矩陣。最終,獲得舵面不同位置處對(duì)應(yīng)的粉末顆粒阻尼器單元的等效剛度、質(zhì)量和阻尼矩陣。

        2.3 粉末顆粒阻尼器分布位置設(shè)計(jì)

        減小振動(dòng)位移可以有效提高結(jié)構(gòu)的顫振性能,而粉末顆粒阻尼器的分布位置對(duì)結(jié)構(gòu)的振動(dòng)位移影響很大。分布位置優(yōu)化設(shè)計(jì)以動(dòng)柔度為優(yōu)化目標(biāo),以粉末顆粒阻尼器單元偽密度為設(shè)計(jì)變量,以粉末顆粒阻尼器體積分?jǐn)?shù)為約束條件,對(duì)顆粒阻尼器在舵面內(nèi)部的分布位置進(jìn)行優(yōu)化設(shè)

        計(jì),優(yōu)化模型為

        式中: 為結(jié)構(gòu)的動(dòng)柔度; Vp 為粉末顆粒阻尼器單元的總體積; V0 為舵面結(jié)構(gòu)的總體積; fp 為粉末顆粒阻尼器體積分?jǐn)?shù); ρ 為設(shè)計(jì)變量, ρj 為第j 個(gè)粉末顆粒阻尼器單元的偽密度; n 為粉末顆粒阻尼器單元總個(gè)數(shù)。

        由于阻尼矩陣的存在,結(jié)構(gòu)的響應(yīng)為復(fù)數(shù),相應(yīng)的動(dòng)柔度也為復(fù)數(shù)形式,表示為

        式中,Re 分別表示設(shè)定頻域內(nèi)動(dòng)柔度的實(shí)部和虛部。

        ρj 關(guān)于結(jié)構(gòu)動(dòng)柔度的靈敏度為

        式中, ω 為角頻率。

        式中: α ! β 為瑞利阻尼系數(shù); cp,j 為第 j 個(gè)粉末顆 粒阻尼器單元結(jié)構(gòu)阻尼矩陣; 為瑞利阻尼矩 陣; p 為懲罰因子;i為虛數(shù); up,j 為第 j 個(gè)粉末顆 粒阻尼器的位移場(chǎng)向量; kp,j 為第 j 個(gè)粉末顆粒阻 尼器的剛度矩陣; mp,j 為第 j 個(gè)粉末顆粒阻尼器單 元的質(zhì)量矩陣。

        在粉末顆粒阻尼器分布設(shè)計(jì)中,選擇粉末質(zhì)量占總質(zhì)量的 1.5% 作為約束條件。對(duì)于圖4中的舵面結(jié)構(gòu),粉末顆粒阻尼器的質(zhì)量占結(jié)構(gòu)總質(zhì)量的比例為

        式中: M 為舵面結(jié)構(gòu)質(zhì)量。解得 fp 為 3.4% ,即優(yōu)化模型式(11)中體積分?jǐn)?shù) fp 設(shè)置為 3.4% 0

        對(duì)粉末顆粒阻尼器進(jìn)行分布優(yōu)化設(shè)計(jì),設(shè)計(jì)結(jié)果如圖8(a)所示,圖中帶編號(hào)的綠色單元為粉末顆粒阻尼器單元,編號(hào)對(duì)應(yīng)阻尼器的阻尼系數(shù)如表1所示。由圖可知,粉末顆粒阻尼器位置多分布于振動(dòng)響應(yīng)比較大的邊緣區(qū)域,相對(duì)比較集中。阻尼器較為集中的分布形態(tài)可以在短時(shí)間內(nèi)提供較大的耗能效果,對(duì)振動(dòng)有著比較強(qiáng)的抑制作用,能夠有效提高舵面的顫振性能。粉末顆粒阻尼器位置優(yōu)化過程迭代曲線如圖8(b)所示。優(yōu)化迭代過程中,結(jié)構(gòu)動(dòng)柔度 與結(jié)構(gòu)初始動(dòng)柔度 的比值在前10次迭代中下降較快,隨后,迭代結(jié)果逐漸趨于平緩,直至第58次時(shí)實(shí)現(xiàn)收斂。

        2.4 工程化設(shè)計(jì)

        對(duì)設(shè)計(jì)結(jié)果進(jìn)行工程化,原型舵面結(jié)構(gòu)如圖3(a)所示,優(yōu)化后并進(jìn)行工程化設(shè)計(jì)的舵面結(jié)構(gòu)如圖9

        表1各單元粉末顆粒阻尼器阻尼系數(shù)Tab.1 Damping coefficients of powder particle dampers for each unit

        所示。通過三維模型計(jì)算可得優(yōu)化后舵面結(jié)構(gòu)的質(zhì)量相較于原型舵面結(jié)構(gòu)的減少了 22%

        圖9優(yōu)化后工程化設(shè)計(jì)的舵面結(jié)構(gòu)

        3 舵面力學(xué)性能評(píng)估

        為驗(yàn)證優(yōu)化設(shè)計(jì)后粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)的性能,對(duì)原型和優(yōu)化后舵面結(jié)構(gòu)的靜力學(xué)、模態(tài)、幅頻響應(yīng)和顫振等性能進(jìn)行分析,通過這4個(gè)性能的對(duì)比,驗(yàn)證優(yōu)化設(shè)計(jì)效果

        3.1 靜力學(xué)性能

        將工程化舵面結(jié)構(gòu)導(dǎo)入有限元分析軟件中進(jìn)行建模,在舵軸末端施加固定約束,將氣動(dòng)載荷等效為一個(gè)垂直于表面的節(jié)點(diǎn)力,大小為 800N 載荷位置為展向距根弦 74.3mm ,弦向距后緣154.2mm ,如圖10所示點(diǎn)1位置。舵面結(jié)構(gòu)的材料為鈦合金,材料密度為 4.4g/cm3 ,彈性模量為110GPa ,泊松比為0.33。網(wǎng)格劃分為四面體網(wǎng)格,原型舵面和優(yōu)化舵面結(jié)構(gòu)網(wǎng)格數(shù)分別為221403和203706。求解得到原型舵面與優(yōu)化舵面的位移云圖如圖11所示,最大位移分別為 19.62mm 和17.59mm 。相較于原型舵面結(jié)構(gòu),優(yōu)化后的舵面結(jié)構(gòu)最大位移值減少了 10.35% ,顯示出更好的靜力學(xué)性能。

        3.2 模態(tài)性能

        模態(tài)頻率對(duì)舵面綜合性能有較大的影響,模態(tài)分析中結(jié)構(gòu)的邊界約束與3.1一致。經(jīng)有限元分析求解,原型舵面和優(yōu)化舵面結(jié)構(gòu)前4階模態(tài)振型如圖12所示,具體模態(tài)頻率結(jié)果如表2所示。由圖和表可知:優(yōu)化前后舵面結(jié)構(gòu)各階模態(tài)振型基本相同,優(yōu)化后舵面結(jié)構(gòu)的1階固有頻率為50.50Hz ,較原型舵面結(jié)構(gòu)的1階固有頻率 41.36Hz 提高了 22.10% , 2~4 階模態(tài)頻率也都有不同程度的提高。

        Fig.9Optimized rudder structure for engineering design圖10載荷位置示意圖Fig.10Schematic diagram of load position
        圖11優(yōu)化前后舵面靜力學(xué)分析結(jié)果
        圖12優(yōu)化前后舵面結(jié)構(gòu)前4階模態(tài)振型Fig.12First four-order modal vibration modes of the rudder structure before and after optimization

        3.3 幅頻響應(yīng)性能

        在靜力學(xué)以及模態(tài)分析中一般不考慮阻尼項(xiàng),阻尼對(duì)靜力學(xué)和模態(tài)性能影響較小。粉末顆粒阻尼器的阻尼特性采用COMBIN14彈簧阻尼單元進(jìn)行簡(jiǎn)化。彈簧阻尼單元中設(shè)置彈簧剛度為0,阻尼值為表1中對(duì)應(yīng)的阻尼系數(shù)值。在舵面舵軸位置處施加固定約束。設(shè)置分析頻率區(qū)間為0~500Hz ,力值為 10N ,采樣間隔為 1Hz ,振動(dòng)測(cè)量點(diǎn)為后緣尖弦角點(diǎn)。幅頻響應(yīng)分析結(jié)果如圖13所示,縱坐標(biāo)位移值為對(duì)數(shù)坐標(biāo),取自然對(duì)數(shù)為底。結(jié)果顯示,原型舵面幅頻響應(yīng)最大位移為15.3mm ,優(yōu)化舵面最大位移為 11.8mm ,幅頻響應(yīng)最大位移減少了 22.88% 。因此,粉末顆粒阻尼器的存在有效提升了舵面的動(dòng)力學(xué)性能。

        表2優(yōu)化前后舵面結(jié)構(gòu)前4階固有頻率對(duì)比

        3.4 顫振性能

        采用PK法進(jìn)行顫振邊界搜索,顫振分析時(shí),設(shè)置大氣密度、飛行速度以及馬赫數(shù),在粉末顆粒阻尼器處添加表1中對(duì)應(yīng)阻尼系數(shù)的COMBIN14彈簧阻尼單元。PK法通過求解系統(tǒng)特征值的實(shí)部(阻尼比)和虛部(瀕率),判斷能量狀態(tài)。優(yōu)化前后,舵面在馬赫數(shù)分別為4、5、6時(shí)的顫振速度如表3所示。由表3可見,優(yōu)化舵面在不同馬赫數(shù)下顫振速度提高均大于 10% ,優(yōu)化舵面結(jié)構(gòu)顯示出更好的顫振性能。舵面結(jié)構(gòu)的顫振速度隨著馬赫數(shù)的增加而增大,優(yōu)化舵面相較于原型舵面在高速飛行時(shí)擁有更加良好的綜合力學(xué)性能。取前4階模態(tài)的速度與阻尼繪制成速度-阻尼圖,取前4階模態(tài)的速度與頻率繪制成速度-頻率圖。不同模態(tài)的固有頻率(如彎曲模態(tài)和扭轉(zhuǎn)模態(tài))隨速度變化的曲線若發(fā)生交匯和分叉,則預(yù)示顫振發(fā)生,對(duì)于速度-阻尼圖,阻尼為系統(tǒng)總阻尼(結(jié)構(gòu)阻尼 + 氣動(dòng)阻尼)。當(dāng)阻尼為0時(shí),系統(tǒng)處于臨界穩(wěn)定狀態(tài),此時(shí)對(duì)應(yīng)的速度為顫振速度的起始點(diǎn)。當(dāng)阻尼從負(fù)(能量耗散)轉(zhuǎn)為正(能量累積)時(shí),對(duì)應(yīng)的速度即為顫振臨界速度。在不同馬赫數(shù)下對(duì)舵面進(jìn)行顫振仿真,圖14為在馬赫數(shù)為5時(shí)原型舵面和優(yōu)化舵面的速度-阻尼圖和速度-頻率圖。由圖可見,在所有速度下,舵面結(jié)構(gòu)發(fā)生顫振的模態(tài)為兩個(gè)彎扭模態(tài)。

        圖13優(yōu)化前后舵面幅頻響應(yīng)分析結(jié)果Fig.13 Rudder amplitude-frequency response analysis results before and after optimization
        表3優(yōu)化前后舵面顫振性能對(duì)比Tab.3Comparison of the rudder flutter performance before and after optimization

        4結(jié)論

        提出一種金屬粉末顆粒阻尼復(fù)合舵面結(jié)構(gòu)優(yōu)化設(shè)計(jì)方法。首先,基于自適應(yīng)成長(zhǎng)法對(duì)舵面結(jié)構(gòu)的內(nèi)部骨架分布進(jìn)行設(shè)計(jì),滿足高承載要求。其次,基于均勻化法進(jìn)行粉末顆粒阻尼器剛度、質(zhì)量、阻尼等效性能計(jì)算。在此基礎(chǔ)上,基于密度法,以粉末顆粒阻尼器的偽密度為設(shè)計(jì)變量,對(duì)得到的等效計(jì)算結(jié)果進(jìn)行插值,實(shí)現(xiàn)粉末顆粒阻尼器最優(yōu)位置設(shè)計(jì)。最后,對(duì)設(shè)計(jì)結(jié)果進(jìn)行工程化建模,并對(duì)其靜力學(xué)、模態(tài)、幅頻響應(yīng)和顫振等性能進(jìn)行分析,驗(yàn)證設(shè)計(jì)效果。相關(guān)研究結(jié)果表明:優(yōu)化后的粉末顆粒阻尼器主要分布在結(jié)構(gòu)振動(dòng)位移較大的位置,這種分布形式能夠最大程度地激發(fā)粉末顆粒的運(yùn)動(dòng)狀態(tài),通過顆粒間碰撞耗能實(shí)現(xiàn)減振,從而實(shí)現(xiàn)提升舵面顫振性能的目的。通過性能對(duì)比可知,相較于原型舵面,優(yōu)化舵面質(zhì)量降低 22% ,最大位移減小 10.35% 1階固有頻率提升 22.10% ,馬赫數(shù)在4、5、6時(shí),顫振速度提升均在 10% 以上。

        圖14馬赫數(shù)為5時(shí)優(yōu)化前后舵面顫振結(jié)果Fig.14Rudder flutter results before and after optimization at Mach number 5

        點(diǎn)陣填充既可以提高結(jié)構(gòu)的可制造性,還可以進(jìn)一步增強(qiáng)結(jié)構(gòu)的剛度,在未來的研究中,可進(jìn)一步對(duì)共形點(diǎn)陣填充設(shè)計(jì)方法進(jìn)行研究。此外,熱失效是舵面結(jié)構(gòu)的重要失效形式。因此,如何在考慮結(jié)構(gòu)熱防護(hù)的同時(shí),實(shí)現(xiàn)舵面結(jié)構(gòu)的承載、抑振、熱控協(xié)同設(shè)計(jì),這也是未來重要的研究方向。

        參考文獻(xiàn):

        [1]鈕耀斌.高超聲速飛行器機(jī)翼非線性顫振研究[D].長(zhǎng)沙:國防科學(xué)技術(shù)大學(xué),2013

        [2]池沛,楊艷艷,陳宗基.空天飛行器建模分析與自主重構(gòu)[J].航空學(xué)報(bào),2008,29(S1):S163-S169.

        [3]劉燕斌.高超聲速飛行器建模及其先進(jìn)飛行控制機(jī)理的研究[D].南京:南京航空航天大學(xué),2007.

        [4]王昕江,郭力,金朋,等.一種基于拓?fù)鋬?yōu)化的舵面仿生多級(jí)分叉結(jié)構(gòu)設(shè)計(jì)[J].氣體物理,2020,5(6):45-51.

        [5]WANG C, ZHUJH, WUMQ,et al. Multi-scale designandoptimization for solid-lattice hybrid structures and theirapplication to aerospace vehicle components[J].ChineseJournal ofAeronautics, 2021, 34(5): 386-398.

        [6]ZHUJH, ZHAO YB, ZHANG WH, et al. Bio-inspiredfeature-driven topology optimization for rudder structure

        design[J].Engineered Science,2019, 5: 46-55.

        [7]DING X H, YAMAZAKI K. Stiffener layout design for plate structures by growing and branching tree model (application to vibration-proof design)[J]. Structural and Multidisciplinary Optimization, 2004, 26(1): 99-110.

        [8]HU T N,DING X H, ZHANG H, et al. Geometry and sizeoptimization of stiffener layout for three-dimensional boxstructures with maximization of natural frequencies[J].Chinese Journal of Aeronautics, 2023,36(1): 324-341.

        [9]張德慧,丁曉紅,胡天男,等.基于改進(jìn)自適應(yīng)成長(zhǎng)法的薄壁結(jié)構(gòu)頻率優(yōu)化設(shè)計(jì)[J].航空學(xué)報(bào),2023,44(19):228378.

        [10]鄭昌隆,丁曉紅,沈洪,等.基于自適應(yīng)成長(zhǎng)法的舵面結(jié)構(gòu)動(dòng)力學(xué)拓?fù)鋬?yōu)化設(shè)計(jì)方法研究[J].空天防御,2021,4(2): 7-12.

        [11]魯正,呂西林,閆維明.顆粒阻尼技術(shù)研究綜述[J].振動(dòng)與沖擊,2013,32(7): 1-7.

        [12]張學(xué)軍,唐思熠,肇恒躍,等.3D打印技術(shù)研究現(xiàn)狀和關(guān)鍵技術(shù)[J].材料工程,2016,44(2):122-128.

        [13]SCOTT-EMUAKPORO,GEORGET,RUNYONB,etal. Investigating damping performance of laser powder bed fused components with unique internal structures[C]//Proceedings of the ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition.

        Oslo:ASME,2018.

        [14]SCOTT-EMUAKPOR O,GEORGE T,RUNYONB, et al. Assessingmanufacturingrepeatabilityofinherently damped nickel alloy components via forced-response testing[C]//Proceedings of the ASME Turbo Expo 2019: Turbomachinery Technical Conference and Exposition. Phoenix:ASME,2019.

        [15]SCOTT-EMUAKPOR O, SCHOENING A,GOLDINA, et al. Internal geometry effectson inherent damping performance of additively manufactured components[J]. AIAA Journal,2021, 59(1): 379-385.

        [16]SCOTT-EMUAKPORO,BECKJ,RUNYONB,etal. Determining unfused powder threshold for optimal inherent dampingwith additivemanufacturing[J].Additive Manufacturing,2021, 38:101739.

        [17]NIEDERMEYERJ,EHLERST,LACHMAYERR Potential of additively manufactured particle damped compressor blades:a literature review[J].Procedia CIRP, 2023,119:570-575.

        [18]EHLERS T, TATZKO S, WALLASCHEK J, et al. Design of particle dampers for additive manufacturing[J]. Additive Manufacturing,2021,38:101752.

        [19]CELLI D A,JANCZEWSKI T,SHERIDAN L,et al. Dampingandfatigueperformanceofadditive manufactured particle damper infused instrumentation rake[C]//Proceedings of the ASME Turbo Expo 2023: Turbomachinery Technical Conference and Exposition. Boston:ASME,2023.

        [20] DING X H, YAMAZAKI K. Adaptive growth technique of stiffener layout pattern for plate and shell structures to achieve minimum compliance[J]. Engineering Optimization,2005,37(3):259-276.

        [21]王謙,丁曉紅,張橫.厚薄通用四邊形平板殼元在薄壁結(jié) 構(gòu)加筋布局優(yōu)化中的應(yīng)用[J].空天防御,2023,6(2): 55-61.

        [22]鮑澤源,丁曉紅,李海東,等.粉末顆粒阻尼器阻尼性能 的影響因素分析[J].機(jī)械設(shè)計(jì)與研究,2024,40(6): 110-115.

        (編輯:丁紅藝)

        猜你喜歡
        阻尼器粉末顆粒
        Li7La3Zr2O12 固態(tài)電解質(zhì)致密度提升策略研究進(jìn)展
        遼寧化工(2025年7期)2025-08-18 00:00:00
        超細(xì)鎳粉的制取工藝
        遼寧化工(2025年7期)2025-08-18 00:00:00
        基于EDEM的攪拌式糧食烘干筒倉攪拌均勻性研究
        補(bǔ)藻模式對(duì)菌藻顆粒污泥形成及微生物群落結(jié)構(gòu)的影響
        新型兩階段屈曲約束支撐滯回性能試驗(yàn)研究
        離心壓縮機(jī)葉輪沖蝕失效的可靠性研究
        日本视频一中文有码中文| 久久99精品久久久久久齐齐百度| 爱v天堂在线观看| 国产专区亚洲专区久久| 大肉大捧一进一出好爽视频动漫| 日本熟女精品一区二区三区| 亚洲精品乱码久久久久久中文字幕| 中文无码制服丝袜人妻av| 91中文人妻丝袜乱一区三区| 一本色道88久久加勒比精品| 亚洲av无码无限在线观看| 中日av乱码一区二区三区乱码| 亚洲成aⅴ人片在线观看天堂无码| 丝袜美腿诱惑区在线播放| 人人鲁人人莫人人爱精品| 最新国产av无码专区亚洲| 探花国产精品三级在线播放| 午夜视频一区二区三区在线观看| 激情综合色五月丁香六月欧美 | 国产精品乱码一区二区三区 | 91精品国产91久久久无码95| av资源在线播放网站| 在线观看免费日韩精品| 无码一区二区三区在线| 亚洲日韩欧美一区二区三区 | 涩涩鲁精品亚洲一区二区| 7m精品福利视频导航| 亚洲欧美另类自拍| 亚洲日本一区二区在线观看| 亚洲国产精品亚洲一区二区三区| 国产高清一区二区三区视频| 国产福利小视频在线观看| 性感的小蜜桃在线观看| 亚洲av日韩综合一区久热 | 国产小屁孩cao大人| 偷拍熟女露出喷水在线91| 成人国产一区二区三区| 国产在线无码一区二区三区| 无码精品一区二区三区免费16 | 人妻中文字幕一区二区二区| 开心久久婷婷综合中文字幕|