孫原理, 李兆俊

(海軍裝備研究院 北京,100161)

?

船用汽輪機(jī)轉(zhuǎn)子盤軸耦合振動特性研究*

孫原理, 李兆俊

(海軍裝備研究院 北京,100161)

基于小長徑比建模理論，建立了盤軸耦合轉(zhuǎn)子的有限元模型，并對采用工程實(shí)際的斜切方式和輪盤單元等不同建模方法下船用汽輪機(jī)轉(zhuǎn)子的振動特性進(jìn)行了分析,給出了較為系統(tǒng)的、適用于小長徑比、短粗轉(zhuǎn)子等結(jié)構(gòu)的建模方法對轉(zhuǎn)子系統(tǒng)振動特性的影響規(guī)律。不同的建模方法會對轉(zhuǎn)子盤軸耦合振動特性計(jì)算結(jié)果有顯著影響。輪盤單元會使得轉(zhuǎn)子軸承系統(tǒng)出現(xiàn)附加的振型，表現(xiàn)為輪盤與轉(zhuǎn)軸的耦合模態(tài)振型，輪盤單元建模的計(jì)算方法能準(zhǔn)確描述盤軸耦合轉(zhuǎn)子的振動特性。研究結(jié)果表明，對于厚度較大的輪盤，輪盤對轉(zhuǎn)子耦合振動的影響較小，工程常用的斜切建模會引入附加剛度，引起較大求解誤差，且這種誤差隨著輪盤厚度的減小而逐漸降低。

盤軸耦合； 有限元模型； 輪盤單元； 振動特性

引 言

艦船動力設(shè)備中的齒輪、汽輪機(jī)、燃?xì)廨啓C(jī)等結(jié)構(gòu)中存在小長徑比、短粗轉(zhuǎn)子，常規(guī)轉(zhuǎn)子動力學(xué)研究中轉(zhuǎn)子部件的梁假設(shè)在上述情況下不再適用，故對于這些小長徑比、短粗轉(zhuǎn)子的建模還存在問題[1-3]；同時(shí)在這些轉(zhuǎn)子中存在轉(zhuǎn)子輪盤與轉(zhuǎn)子轉(zhuǎn)軸之間的相互耦合，振動特性較為復(fù)雜。目前國內(nèi)外已有針對動力設(shè)備中轉(zhuǎn)子動力學(xué)特性的相關(guān)研究，并形成了許多可行計(jì)算分析方法和相關(guān)的軟件[4-5]，如傳遞矩陣法[6-8]、有限元法[9-10]、子結(jié)構(gòu)法[11]、結(jié)構(gòu)綜合法[12-13]和假設(shè)模態(tài)法[14]等。其中，Chun等[14]利用子結(jié)構(gòu)綜合法和假設(shè)模態(tài)法分析了軸和葉輪系統(tǒng)的耦合振動。總的來說，對于轉(zhuǎn)子盤軸耦合振動問題的耦合建模方法和分析方法研究還尚不成熟。因此，迫切需要針對轉(zhuǎn)子盤軸耦合振動的相關(guān)研究問題，研究考慮轉(zhuǎn)子盤軸耦合動力學(xué)的建模方法和分析方法，給出較為系統(tǒng)的、能完全適用于小長徑比、短粗轉(zhuǎn)子等結(jié)構(gòu)和工況特點(diǎn)的建模、計(jì)算分析方法，為動力設(shè)備中轉(zhuǎn)子盤軸耦合振動特性研究提供技術(shù)支撐。

基于小長徑比建模理論[15]，筆者將建立盤軸耦合轉(zhuǎn)子的有限元模型，并分析采用實(shí)際的工程斜切方式建模和輪盤單元等不同建模方式下船用汽輪機(jī)轉(zhuǎn)子的振動特性。

1 盤軸耦合轉(zhuǎn)子建模

輪盤是典型的小長徑比轉(zhuǎn)子，基于小長徑比建模理論，筆者研究了盤軸耦合轉(zhuǎn)子的建模計(jì)算方法，給出了盤軸耦合轉(zhuǎn)子的輪盤單元有限元模型和單元的運(yùn)動方程。

對于多自由度轉(zhuǎn)子運(yùn)動方程如下

其中：x為描述轉(zhuǎn)子結(jié)構(gòu)的廣義坐標(biāo)向量；ω為自轉(zhuǎn)角速度；實(shí)對稱矩陣M，C，和K分別為整個(gè)系統(tǒng)的質(zhì)量、阻尼和剛度矩陣；Kω為由附加離心力造成的系統(tǒng)剛度矩陣；實(shí)反對稱陣G，Cr分別為陀螺矩陣和旋轉(zhuǎn)阻尼矩陣，其中旋轉(zhuǎn)阻尼矩陣主要與發(fā)生在該系統(tǒng)的旋轉(zhuǎn)部件的粘性耗散有關(guān)；Fn和Fr分別為由非旋轉(zhuǎn)力和不平衡引起的激勵(lì)力。

如果垂直于z軸的轉(zhuǎn)子在變形時(shí)保持水平，則彎曲運(yùn)動可以通過軸心坐標(biāo)X，Y轉(zhuǎn)換和旋轉(zhuǎn)角位移Φx,Φy進(jìn)行描述,如圖1所示。沿Z軸方向的位移和旋轉(zhuǎn)角位移Φz分別為軸向和扭轉(zhuǎn)運(yùn)動。

圖1 轉(zhuǎn)子坐標(biāo)系Fig.1 Coordinate system of rotor

假定轉(zhuǎn)子和輪盤為剛體，它們對轉(zhuǎn)子彎曲運(yùn)動的影響可以通過拉格朗日坐標(biāo)系中的X,Y方向位移和Φx,Φy旋轉(zhuǎn)角位移表示?？紤]到軸向和扭轉(zhuǎn)運(yùn)動，每個(gè)節(jié)點(diǎn)有3個(gè)平動自由度和3個(gè)轉(zhuǎn)動自由度。如果n是整個(gè)有限元離散的節(jié)點(diǎn)數(shù)，那么需要6n個(gè)自由度來描述整個(gè)轉(zhuǎn)子。

在動力設(shè)備中，有些輪盤的長徑比較小，那么剛體假設(shè)就無法成立，因此，分析小長徑比轉(zhuǎn)子輪盤結(jié)構(gòu)的變形時(shí)就必須增加一個(gè)額外的自由度。

由轉(zhuǎn)子軸向力和輪盤上熱拉伸引起的幾何附加剛度可以在剛度矩陣中K考慮，旋轉(zhuǎn)阻尼矩陣主要與發(fā)生在該系統(tǒng)的旋轉(zhuǎn)部件的黏性耗散有關(guān)。如果軸向和扭轉(zhuǎn)自由度與彎曲自由度能夠解耦，X,Y方向位移和Φx,Φy旋轉(zhuǎn)角位移可以被組合成復(fù)數(shù)位移Z=X+iY和復(fù)數(shù)角位移Φ=Φy-iΦx，這樣陀螺矩陣和旋轉(zhuǎn)阻尼矩陣可以用實(shí)對稱陣來表達(dá)。

1.1 輪盤單元有限元建模

在動力設(shè)備中典型結(jié)構(gòu)就是輪盤。圖1給出了坐標(biāo)系的變換，通過坐標(biāo)系的變換可以得到在慣性系和固定在輪盤單元上的隨動參考坐標(biāo)系之間點(diǎn)坐標(biāo)的關(guān)系式

(2)

其中:u,v,w分別為輪盤中間平面上任一點(diǎn)P0在徑向、切向和軸向位移分量；輪盤的極坐標(biāo)在Cξηz坐標(biāo)系里用r,β和z表示；P為P0變形后的點(diǎn)。

涉及的旋轉(zhuǎn)變換矩陣為

(3a)

(3b)

(3c)

(3d)

在方程(2)中忽略由幾何非線性所引起的非線性問題，并將相關(guān)的三角函數(shù)用泰勒級數(shù)展開，轉(zhuǎn)角有如下關(guān)系式：

(4)

P點(diǎn)在慣性參考坐標(biāo)系中變?yōu)?/p>

(5)

式(5)給出了P點(diǎn)在慣性參考系CX′Y′Z′中運(yùn)動的描述。未知量u,v,w是轉(zhuǎn)子在t時(shí)刻P0點(diǎn)位置c,?,z的函數(shù)，而轉(zhuǎn)角ΦX′,Φy,Φz和位移Xc,Yc,Zc僅和t有關(guān)。作為一個(gè)變形體，輪盤的動態(tài)描述可以用如下一組拉格朗日坐標(biāo)來表示。

輪盤剛體運(yùn)動

(6)

輪盤變形

(7)

定義形函數(shù)來描述P點(diǎn)位置r,?，z的函數(shù)u(r,?,z,t)、w(r,?,z,t)、v(r,?,z,t)，輪盤被分成如圖2所示的環(huán)狀單元，位移場u,w,v為

(8)

圖2 輪盤單元示意圖Fig.2 Schematic diagram of disk element

1.2 單元運(yùn)動方程

輪盤單元如圖2所示。圖中：h0,hi分別為半徑r0,ri的厚度。在獲得單元形函數(shù)之后，通過計(jì)算單元?jiǎng)幽芎蛣菽芗纯傻玫絾卧倪\(yùn)動方程，單元?jiǎng)幽芎蛣菽艿挠?jì)算如下。

(9)

其中：ρ為質(zhì)量密度。

這樣就可以將位移場方程帶入方程(9)，從而獲得單元在拉格朗日坐標(biāo)系下的動能。

單元的勢能一般包含兩部分：材料的彈性變形引起的勢能Ue和涉及到“幾何剛化效應(yīng)”的勢能Ug。

(10)

廣義的應(yīng)變向量ε可以從輪盤內(nèi)的位移場中獲得

(11)

矩陣D通常表示單元的彈性性能，是輪盤厚度和材料特性的函數(shù)

(12)

2 盤軸耦合轉(zhuǎn)子動力學(xué)特性分析

由于轉(zhuǎn)子輪盤是典型的小長徑比轉(zhuǎn)子，且不能用常規(guī)的梁單元進(jìn)行處理，振動過程中容易與轉(zhuǎn)子轉(zhuǎn)軸發(fā)生耦合作用，使得振動特性較為復(fù)雜。針對這一問題，本研究首先以簡單的單輪盤轉(zhuǎn)子為例，建模中采用實(shí)際的工程斜切方式建模，并將動力學(xué)特性計(jì)算結(jié)果和輪盤單元建模的計(jì)算結(jié)果進(jìn)行了對比。在此基礎(chǔ)上，針對實(shí)際汽輪機(jī)轉(zhuǎn)子模型提出了轉(zhuǎn)子輪盤的3種不同建模方法，并對這3種不同建模方法下的實(shí)際汽輪機(jī)轉(zhuǎn)子動力學(xué)特性進(jìn)行了分析。

筆者提出的3種轉(zhuǎn)子輪盤建模方法分別為直切建模、斜切建模和輪盤單元建模。其中，直切建模是將輪盤部分全部作為集中質(zhì)量和轉(zhuǎn)動慣量加在轉(zhuǎn)子-軸承系統(tǒng)相應(yīng)的節(jié)點(diǎn)上；工程中常用的45°斜切建模是將圖3中虛線內(nèi)部分作為梁單元建模，其余部分作為集中質(zhì)量和轉(zhuǎn)動慣量加于轉(zhuǎn)子-軸承系統(tǒng)相應(yīng)節(jié)點(diǎn)上。而第3種輪盤單元建模是基于1.1節(jié)的輪盤單元有限元建模理論。

圖3 輪盤單元斜切建模示意圖Fig.3 Schematic diagram of disk oblique cutting model

2.1 輪盤單元與斜切建模對單輪盤轉(zhuǎn)子動力學(xué)特性的影響

以單輪盤轉(zhuǎn)子模型為例，本節(jié)采用實(shí)際的工程斜切方式建模，并將其動力學(xué)特性計(jì)算結(jié)果和輪盤單元建模的計(jì)算結(jié)果進(jìn)行了對比。圖4和圖5給出了斜切有限元建模和輪盤單元法有限元建模示意圖。

圖6給出了斜切建模軸系的前3階振型圖；圖7給出了輪盤厚度為5 mm時(shí)輪盤建模時(shí)的前3階振型圖；圖8給出了輪盤厚度為100 mm時(shí)輪盤建模時(shí)的前3階振型圖。

圖4 工程有限元建模Fig.4 Finite element model of rotor in engineering

圖5 輪盤單元有限元建模Fig.5 Finite element model of rotor in disk unit

圖6 輪盤單元斜切建模的振型圖Fig.6 Modal graph of disk oblique cutting model

圖7 輪盤厚度為5 mm時(shí)的振型圖Fig.7 Modal graph when the disc thickness is 5 mm

圖8 輪盤厚度為100 mm時(shí)振型圖Fig.8 Modal graph when the disc thickness is 100 mm

從圖6～8可以看出，當(dāng)輪盤厚度(100 mm)較大時(shí)，輪盤剛度較大，彈性變形小，進(jìn)而對轉(zhuǎn)子耦合振動的影響較??；當(dāng)輪盤厚度(5 mm)較小時(shí)，輪盤對轉(zhuǎn)子振動的影響較大，轉(zhuǎn)子的振動為輪盤的振型與轉(zhuǎn)軸的振型耦合。

表1給出了隨著輪盤厚度的變化斜切建模1階臨界轉(zhuǎn)速和2階臨界轉(zhuǎn)速的偏差，其中轉(zhuǎn)子的直徑為400 mm。從表1中可以看出，隨著輪盤厚度的減小，輪盤單元建模與斜切建模計(jì)算得到的1，2階臨界轉(zhuǎn)速結(jié)果偏差也逐漸降低。當(dāng)輪盤厚度較大時(shí)斜切建模的誤差較大，這是由于輪盤厚度較大時(shí)，采用斜切法建模中選取的直徑大于理論的等效直徑從而引入了一部分附加剛度造成臨界轉(zhuǎn)速偏高。

表1 不同輪盤厚度時(shí)兩種建模計(jì)算結(jié)果的偏差

Tab.1 Difference between two modeling methods in different disk thickness

直徑/mm厚度/mm1階臨界轉(zhuǎn)速/(r/min)2階臨界轉(zhuǎn)速/(r/min)1階臨界轉(zhuǎn)速/(r/min)2階臨界轉(zhuǎn)速/(r/min)1階臨界轉(zhuǎn)速偏差/%2階臨界轉(zhuǎn)速偏差/%0．2501005822225206299231188．192．660．225906033225176462230347．112．300．200806268225206646229576．031．940．175706532225326859228925．011．600．150606830225597102228453．981．270．125507173226097388228263．000．960．100407570226947732228502．140．690．075308036228388141229391．310．440．050208594230758649231340．640．260．025109270234689286235070．170．170．01259663237619669238000．060．16

從表1還可以看出，當(dāng)輪盤厚度直徑比值為0.175時(shí)一階臨界轉(zhuǎn)速的誤差會達(dá)到5%，這類結(jié)構(gòu)在實(shí)際工程中大量存在，如后續(xù)2.2節(jié)中汽輪機(jī)雙列復(fù)速級輪盤的厚度直接比為0.16，因此，輪盤單元的建模計(jì)算方法對于實(shí)際工程有很大意義。

2.2 不同建模方法下汽輪機(jī)轉(zhuǎn)子動力學(xué)特性分析

研究的實(shí)際汽輪機(jī)轉(zhuǎn)子模型，由一個(gè)雙列復(fù)速級和7個(gè)壓力級，兩個(gè)軸承組成。由于輪盤單元的不同建模方法對盤軸耦合轉(zhuǎn)子動力學(xué)特性的影響，因此，兩個(gè)軸承均采用線性化的8參數(shù)模型，且不同工況下軸承的8參數(shù)保持不變。

1) 直切建模。整個(gè)轉(zhuǎn)子總長2.2 m，劃分60個(gè)單元，建模中將雙列復(fù)速級和7個(gè)壓力級的輪盤部分全部作為集中質(zhì)量和轉(zhuǎn)動慣量加在相應(yīng)的節(jié)點(diǎn)上，兩個(gè)軸承采用8參數(shù)模型。汽輪機(jī)轉(zhuǎn)子的模型如圖9所示。

圖9 轉(zhuǎn)子有限元模型示意圖Fig.9 Schematic diagram of rotor finite element model

2) 斜切建模。整個(gè)轉(zhuǎn)子總長2.2 m，劃分60個(gè)單元，建模中將雙列復(fù)速級和7個(gè)壓力級的輪盤部分采用45 °斜切，其余部分作為集中質(zhì)量和轉(zhuǎn)動慣量加在相應(yīng)的節(jié)點(diǎn)上，兩個(gè)軸承采用8參數(shù)模型。汽輪機(jī)轉(zhuǎn)子的模型如圖10所示。

圖10 斜切建模轉(zhuǎn)子有限元模型示意圖Fig.10 Finite element model of rotor oblique cutting model

3) 輪盤建模。本節(jié)中采用直接對輪盤單元進(jìn)行建模，將葉片部分作為集中質(zhì)量和轉(zhuǎn)動慣量附加于相應(yīng)的節(jié)點(diǎn)，軸承采用8參數(shù)模型，整個(gè)軸系總長2.2 m共80個(gè)單元，其中梁單元60個(gè)。圖11給出了模型示意圖。

圖11 輪盤建模轉(zhuǎn)子有限元模型示意圖Fig.11 Finite element model of rotor in disk unit

4) 3種建模結(jié)果對比。表2給出了3種建模方式下轉(zhuǎn)子盤軸耦合振動特性的計(jì)算結(jié)果。從表2中可以看出，但是斜切與輪盤單元建模結(jié)果相差不大。在直切建模中沒有把輪盤單元的剛度計(jì)算在內(nèi)，因此系統(tǒng)的總剛度偏小，使得系統(tǒng)臨界轉(zhuǎn)速值偏小；相反，把整個(gè)的輪盤單元均作為梁單元處理必然會使得附加剛度過大，在45°斜切建模中，整個(gè)輪盤單元只有一部分的剛度對系統(tǒng)有貢獻(xiàn)，45°斜切相交線的直徑外輪盤單元的剛度對總的剛度沒有貢獻(xiàn)。

同時(shí)由于汽機(jī)轉(zhuǎn)子的軸承剛度較小，1階模態(tài)振型表現(xiàn)轉(zhuǎn)子的平動、2階模態(tài)振型表現(xiàn)為轉(zhuǎn)子的錐動；此時(shí)，支承軸承剛度對臨界轉(zhuǎn)速的影響大于轉(zhuǎn)子結(jié)構(gòu)剛度對其的影響。由于汽輪機(jī)各壓力級的輪盤厚度較小，因此，不同建模之間的誤差較小，這與2.1節(jié)單輪盤轉(zhuǎn)子的計(jì)算結(jié)果一致。

表2 3種建模臨界轉(zhuǎn)速結(jié)果

3 結(jié)束語

盤類零件本身的模態(tài)振型會引起輪盤和轉(zhuǎn)子轉(zhuǎn)軸發(fā)生耦合振動，使得對于典型的輪盤單元出現(xiàn)一些附加的振型，主要表現(xiàn)為輪盤軸系耦合振型。筆者提出的基于輪盤單元建模計(jì)算方法能準(zhǔn)確描述盤軸耦合轉(zhuǎn)子的振動特性。當(dāng)輪盤厚度較大時(shí)，輪盤對轉(zhuǎn)子耦合振動的影響較??；當(dāng)輪盤厚度較小時(shí)，輪盤對轉(zhuǎn)子耦合振動的影響較大，轉(zhuǎn)子振動為輪盤振型與轉(zhuǎn)子振型的耦合。輪盤單元建模與斜切建模計(jì)算得到的臨界轉(zhuǎn)速結(jié)果偏差隨著輪盤厚度的減小而逐漸降低。

[1] 孟光. 轉(zhuǎn)子動力學(xué)研究的回顧與展望[J]. 振動工程學(xué)報(bào), 2002,15(1): 1-9.

Meng Guang. Retrospect and prospect to the research on rotordynamics[J]. Journal of Vibration Engineering, 2002,15(1): 1-9.(in Chinese)

[2] 虞烈, 劉恒. 轉(zhuǎn)子-軸承系統(tǒng)動力學(xué)[M]. 西安: 西安交通大學(xué)出版社, 2001: 155-158.

[3] 聞邦椿, 顧家柳, 夏松波，等. 高等轉(zhuǎn)子動力學(xué)-理論、技術(shù)與應(yīng)用[M]. 北京: 機(jī)械工業(yè)出版社, 1999:24-30.

[4] 鐘一諤, 何衍宗, 王正，等. 轉(zhuǎn)子動力學(xué)[M]. 北京: 清華大學(xué)出版社, 1987: 243-247.

[5] Newmark N M. A Method of computation for structural dynamics[J]. Journal of Engineering Mechanics Division,1959, 85(1): 67-94.

[6] Muthy V R. Dynamic characteristic of rotor blades[J]. Journal o f Sound and Vibration, 1978, 49(4):483-500.

[7] 費(fèi)鐘秀. 復(fù)雜轉(zhuǎn)子耦合系統(tǒng)有限元建模及其動力特性研究[D]. 杭州：浙江大學(xué)，2013.

[8] 商大中, 李宏亮. 葉輪與軸耦合振動的傳遞矩陣法分析及仿真計(jì)算[J]. 華東船舶工業(yè)學(xué)院學(xué)報(bào)：自然科學(xué)版,2002,16(6): 37-42.

Shang Dazhong, Li Hongliang. Anlaysis and emulation of coupled vibration of bladed disc and shaft by means of transfer matrix[J]. Journal of East China Shipbuilding Institute：Natural Science Edition, 2002, 16(6): 37-42. (in Chinese)

[9] Soares C A M, Petyt M, Salama A M. Finite element analysis of bladed disks, structural dynamic aspects of bladed disk assemblies[C]∥ American Society of Mechianical Engineers Winter Annual Meeting. New York: The American Society of Mechanical Engineers, 1976: 73 -91.

[10]Soares C A M, Petyt M. Finite element dynamic analysis of practical bladed discs[J]. Journal of Sound and Vibration, 1978, 61(4): 561-570.

[11]Michimura S, Nagamatsu A, Ikevchi T, et al. Vibration analysis on bladed disk assemblies of axial-flow turbines[J]. Bulletin of the Japan Society of Mechanical Engineering, 1981, 24(197): 1988-1993.

[12]張文, 陳運(yùn)西. 葉片-盤-軸整體轉(zhuǎn)子振動特性分析[J]. 航空動力學(xué)報(bào), 1994, 9(2): 120-124.

Zhang Wen, Chen Yunxi. Vibration anlaysis of the whole rotor with blade-disk and shaft[J]. Journal of Aerospace Power, 1994, 9(2): 120-124. (in Chinese)

[13]張文. 轉(zhuǎn)子動力學(xué)理論基礎(chǔ)[M].北京: 科學(xué)出版社, 1991: 57-62.

[14]Chun S B, Lee C W. Vibration Analysis of shaft-bladed disk system by using substructure synthesis and assumed modes method[J]. Journal of Sound and Vibration, 1996, 189(5): 587-608.

[15]Genta G. Dynamics of rotating systems[M]. New York: Springer US, 2005: 170-173.

歡迎訂閱《振動、測試與診斷》

《振動、測試與診斷》由工業(yè)和信息化部主管，南京航空航天大學(xué)和全國高校機(jī)械工程測試技術(shù)研究會聯(lián)合主辦，是反映振動、動態(tài)測試及故障診斷學(xué)科領(lǐng)域的科研成果及其應(yīng)用情況的技術(shù)性刊物。主要刊登國內(nèi)外以振動測試與故障診斷為中心的動態(tài)測試?yán)碚?、方法和手段的研究及?yīng)用方面的技術(shù)文獻(xiàn)，包括實(shí)驗(yàn)測試技術(shù)、測試儀器的研制、方法和系統(tǒng)組成、信號分析、數(shù)據(jù)處理、參數(shù)識別與故障診斷以及有關(guān)裝置的設(shè)計(jì)、使用、控制、標(biāo)定和校準(zhǔn)等，不拘泥于行業(yè)和測試項(xiàng)目。

本刊為EICompendex數(shù)據(jù)庫收錄期刊和中文核心期刊，雙月刊，每逢雙月末出版，每本定價(jià)30元，全年180元。歡迎訂閱和投稿，歡迎在本刊刊登各類廣告和科技信息。

編輯部地址：南京市御道街29號 郵政編碼：210016 電話：(025)84893332

傳真：(025)84893332 E-mail:qchen@nuaa.edu.cn 網(wǎng)址：http:∥zdcs.nuaa.edu.cn

*國家自然科學(xué)基金資助項(xiàng)目(51106035)

2015-03-09；

2015-05-22

10.16450/j.cnki.issn.1004-6801.2017.02.017

TK14; TH113.1

孫原理，男，1983年10月生，工程師。主要研究方向?yàn)楹藙恿こ獭?/p>

E-mail:sy1850122@126.com