国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

基于有限元的軸系回旋振動響應(yīng)計算及試驗

2021-03-06 03:15:06富琦晉周瑞平黃國兵高子航
造船技術(shù) 2021年1期
關(guān)鍵詞:渦動平衡力軸系

富琦晉,周瑞平,黃國兵,高子航

(武漢理工大學(xué) 能源與動力工程學(xué)院,湖北 武漢 430063)

0 引 言

船舶推進軸系回旋振動主要是由于軸系在運轉(zhuǎn)時受到螺旋槳旋轉(zhuǎn)時產(chǎn)生的不平衡力、不均勻伴流場流體力及螺旋槳偏心質(zhì)量的質(zhì)量力而產(chǎn)生的一種橫向振動。國內(nèi)外對于回旋振動響應(yīng)研究,大多在求解軸系臨界轉(zhuǎn)速的基礎(chǔ)上,基于傳遞矩陣法進行簡化處理,僅僅獲取軸系的振動位移,主要是從軸心軌跡上分析回旋振動的影響。關(guān)于回旋振動反映在軸承基座上點的振動加速度響應(yīng)的研究更是鮮有。周飛云[1]通過傳遞矩陣法計算某軸系的共振轉(zhuǎn)速,并進行實船測試。王蘇等[2]通過ANSYS計算艦船軸系的臨界轉(zhuǎn)速。馬召召等[3]通過命令流計算噴水推進軸系的固有頻率并繪制渦動振型。FIROUZI等[4]研究軸徑和剛度對回旋振動固有頻率的影響大小。LEE[5]計算50 000 t油船的回旋振動共振轉(zhuǎn)速與彎曲應(yīng)力,并通過應(yīng)變儀進行實船測試完成驗證。周凌波等[6-7]闡述回旋振動的激振力和影響因素及螺旋槳受橫向單位激勵力時回旋振動對艉部船體結(jié)構(gòu)的激勵特性。李小軍等[8]研究尾軸承剛度對回旋振動的影響,并計算在幅值一定的橫向力下螺旋槳與軸承節(jié)點處的位移響應(yīng),但其只研究不平衡力的橫向分量且未考慮隨轉(zhuǎn)速的幅值變化。

目前,由于推進軸系的轉(zhuǎn)速較低,在進行回旋振動固有特性的研究上往往忽略其與轉(zhuǎn)速之間的關(guān)系,但實際上推進軸系的回旋振動是軸系的渦動,為更好地利用系統(tǒng)仿真方法進行回旋振動特性及相應(yīng)研究,基于轉(zhuǎn)子動力學(xué)中的概念,以渦動頻率表述回旋振動頻率隨軸系轉(zhuǎn)速的變化關(guān)系。同時,部分學(xué)者關(guān)于回旋振動響應(yīng)的計算沒有考慮旋轉(zhuǎn)不平衡激勵力是存在兩個方向的分量且幅值隨轉(zhuǎn)速而變化的頻域力,將其簡化成單向且幅值一定的一般激勵力,且一般只研究至螺旋槳或軸承所在節(jié)點處的響應(yīng)值,而針對軸承及基座上點的響應(yīng)計算較少,也未見驗證其正確性的試驗。

以某推進軸系臺架為例,基于Workbench軟件提出完整的回旋振動渦動頻率及其響應(yīng)計算方法,將頻域響應(yīng)力加在軸承支承中心,計算在額定工況下軸承基座上點的加速度響應(yīng),并設(shè)計振動試驗,實測驗證該方法的正確性。

1 渦動頻率及振型

1.1 理論基礎(chǔ)

轉(zhuǎn)子動力學(xué)中考慮轉(zhuǎn)動效應(yīng)的結(jié)構(gòu)動力方程為

(1)

轉(zhuǎn)子動力學(xué)理論認為固有頻率是結(jié)構(gòu)的固有特性,只與結(jié)構(gòu)自身剛度、質(zhì)量分布有關(guān)。結(jié)構(gòu)在轉(zhuǎn)動時的振動頻率稱為渦動頻率。當轉(zhuǎn)速為零時渦動頻率與固有頻率相等。當轉(zhuǎn)速不為零時,因回轉(zhuǎn)效應(yīng)的影響,轉(zhuǎn)動結(jié)構(gòu)的渦動頻率隨轉(zhuǎn)速變化而改變,正進動 (Forward Whirling,F(xiàn)W) 頻率隨轉(zhuǎn)速的增加而增加,逆進動 (Backward Whirling,BW) 頻率則隨轉(zhuǎn)速的增加而減小。

1.2 臺架有限元模型建立

圖1為進行理論建模與試驗研究的某推進軸系模擬臺架,由模擬加載裝置、配重盤、后尾軸承、前尾軸承、推力軸承、高彈性聯(lián)軸器、減速齒輪箱及推進電機組成。

注:高彈性聯(lián)軸器、減速齒輪箱、推進電機在圖中省略圖1 軸系臺架布置

針對所研究的軸系臺架,采用Workbench對軸系臺架采取梁單元和實體單元兩種建模方法,在求解回旋振動渦動頻率時采用梁單元模型,在計算回轉(zhuǎn)不平衡力激勵對軸承基座的諧響應(yīng)時需要求解軸承基座上點的響應(yīng)值,因此通過SolidWorks對軸系臺架進行實體建模。

1.3 渦動頻率計算

轉(zhuǎn)動結(jié)構(gòu)部件的結(jié)構(gòu)動力學(xué)分析需要考慮慣性效應(yīng)的影響,即回旋振動的陀螺效應(yīng)。在計算時打開阻尼效應(yīng)和科里奧利效應(yīng)。在后尾軸承、前尾軸承、推力軸承所在節(jié)點建立彈簧單元并設(shè)置剛度為1×106N/mm、1×107N/mm、1×107N/mm;賦予軸系轉(zhuǎn)速并約束軸向平動及轉(zhuǎn)動自由度去除剛體模態(tài),在配重盤所在節(jié)點添加質(zhì)量點并賦予質(zhì)量4 800 kg、附水極慣性矩Jp為1 978 kg·m2、徑慣性矩Jd為989 kg·m2。將末端節(jié)點的邊界條件設(shè)置為簡支,代表模型末端簡化到高彈性聯(lián)軸器處,如圖2所示。

注: A為點質(zhì)量; B為轉(zhuǎn)速:0 r/min; C為位移; D為固定轉(zhuǎn)動: 0°; E為簡支4:0 mm圖2 邊界條件設(shè)置示例

在后處理模塊添加“Total Deformation”以查看每一階模態(tài)振型,插入“Campbell Diagram”以查看坎貝爾圖并求解渦動頻率和臨界轉(zhuǎn)速,如圖3所示??藏悹枅D是描述渦動頻率隨轉(zhuǎn)速不同而變化的圖形。圖3中斜率為1,代表1次回旋振動。當斜率為槳葉數(shù)時,代表葉頻次回旋振動。斜線與渦動頻率曲線的交點用三角標出,代表一、二、三階1次回旋的BW與FW的臨界轉(zhuǎn)速值??藏悹柷€中轉(zhuǎn)速為0時對應(yīng)的3個頻率點即為橫向振動的前三階固有頻率。每1個臨界轉(zhuǎn)速對應(yīng)1個渦動頻率,在數(shù)值上渦動頻率的值為其對應(yīng)的臨界轉(zhuǎn)速除以60,單位為Hz。在坎貝爾圖上的幾何含義為臨界轉(zhuǎn)速點對應(yīng)的y坐標值。

艦船推進軸系最高轉(zhuǎn)速一般在300 r/min以下,如果橫坐標按運行轉(zhuǎn)速范圍設(shè)置,坎貝爾圖有時會變得不易觀察。某階模態(tài)的回轉(zhuǎn)效應(yīng)影響比較顯著會使FW、BW頻率曲線分離較大,增加轉(zhuǎn)速可能會使其與其他階的渦動頻率曲線交叉從而影響結(jié)果。此時主要采取兩種解決辦法:(1)增加載荷步的數(shù)目,即增加求解坎貝爾圖的轉(zhuǎn)速點,在載荷步數(shù)增加時坎貝爾圖會更趨于一條光滑的曲線;

圖3 坎貝爾渦動曲線圖

(2)減小分析的轉(zhuǎn)速,在轉(zhuǎn)速減小時,正逆渦動頻率差值減小,不易產(chǎn)生頻率曲線交叉現(xiàn)象。

由于本軸系通過配重盤模擬螺旋槳,只存在軸頻而不存在葉頻次,因此,在仿真計算中只考慮1次正回旋(h=1)與1次逆回旋(h=-1)。將有限元回旋振動計算所得渦動頻率與傳統(tǒng)傳遞矩陣法計算所得固有頻率結(jié)果進行對比,如表1所示。

表1 固有頻率結(jié)果對比

由表1可看出:有限元法計算所得頻率與傳統(tǒng)傳遞矩陣法所得固有頻率誤差小于1%,表明邊界條件設(shè)置及有限元計算是準確的。

得到渦動頻率后選擇“User Defined Result”對渦動振型進行繪制后處理。目前以傳遞矩陣法為基礎(chǔ)進行振型的推算時,由于累加誤差的影響,難以獲取滿意的振動形態(tài)。如果考慮實際軸承支承的各向異性,傳遞矩陣法振型繪制會更加麻煩,且準確性較差。通過有限元法可在“Tabular Data”中獲得軸上每個節(jié)點的振動幅值,得到較為精確的軸系撓度曲線,即正逆回旋的各階振型。圖4為正回旋振型。由圖4可知:不同階回旋振動振型不同;由于回轉(zhuǎn)效應(yīng),正回旋與逆回旋的渦動頻率不同,但同一階正回旋和逆回旋的振型與幅值都完全相同,只是方向不同而已。

圖4 正回旋振型

2 回旋不平衡力響應(yīng)及軸承座響應(yīng)計算

2.1 激勵力

關(guān)于推進軸系回旋振動的激勵力,特別是對于潛艇推進軸系,其激勵主要來自螺旋槳水動力及旋轉(zhuǎn)過程中產(chǎn)生的不平衡力等,關(guān)于螺旋槳激勵的計算國內(nèi)外文獻已有充分的研究。為模擬螺旋槳水動力的影響,試驗臺架一般設(shè)計有軸向及橫向加載裝置。在試驗軸系臺架正常運轉(zhuǎn)工況下,電機激振力及偏心重力引起的二次激振力的影響較小,可將其忽略。由于試驗條件所限,研究中只考慮模擬螺旋槳的配重盤旋轉(zhuǎn)時產(chǎn)生的不平衡離心力。配重盤質(zhì)量較大且懸臂支承,因此在旋轉(zhuǎn)時由于偏心質(zhì)量的存在產(chǎn)生旋轉(zhuǎn)不平衡力F0。針對單一激勵力,研究回旋振動不平衡力響應(yīng)的計算方法:

F0=meω2=meΩ2

(2)

式中:m為偏心質(zhì)量,kg;e為偏心質(zhì)量中心到轉(zhuǎn)軸的距離,m;ω為結(jié)構(gòu)轉(zhuǎn)速,rad/s;Ω為激勵圓頻率,rad/s。

由式(2)可知:旋轉(zhuǎn)不平衡激勵力區(qū)別于大小不隨轉(zhuǎn)速而變化的一般激勵力,其幅值與轉(zhuǎn)速有關(guān)。在激勵范圍內(nèi),不平衡力是一個幅值隨激勵頻率(轉(zhuǎn)速)變化的頻域力。

對于推進軸系來說,其激勵與轉(zhuǎn)軸的結(jié)構(gòu)轉(zhuǎn)速不一定是同步的。在船舶運行時,流體的激勵頻率為葉頻,螺旋槳偏心產(chǎn)生的重力激勵頻率為兩倍軸頻,質(zhì)量偏心旋轉(zhuǎn)產(chǎn)生的不平衡激勵頻率為軸頻。因此,在每個激勵頻率下采用結(jié)構(gòu)轉(zhuǎn)速來計算回轉(zhuǎn)效應(yīng)。結(jié)構(gòu)轉(zhuǎn)速ω計算式為

ω=Ω/s=2πf/s

(3)

式中:f為激勵頻率,Hz;s為激勵圓頻率Ω與結(jié)構(gòu)轉(zhuǎn)速ω之比。在Workbench的“Rotating Force”中,“Synchronous Ratio”即為s。

在軸系的轉(zhuǎn)速方向定義為軸向(x向)時,其激勵力F的相位滯后角為α,作用在垂直轉(zhuǎn)軸的平面內(nèi)。旋轉(zhuǎn)不平衡激勵力的大小與方向都隨轉(zhuǎn)動而變化,但其方向始終為運動軌跡的切線方向,且幅值是一定的。因此,旋轉(zhuǎn)不平衡激勵力可分解成幅值相等的橫向、垂向兩個簡諧力,即

Fy=F0cos(Ωt-α)

(4)

Fz=F0sin(Ωt-α)

(5)

式中:Fy、Fz分別為橫向、垂向的簡諧激勵力;Ωt、α分別為激勵力的相位角與相位滯后角。

圖5為一般激振力分解示例。

圖5 一般激振力分解示例

對于質(zhì)量偏心旋轉(zhuǎn)產(chǎn)生的不平衡激勵力而言,其激勵頻率與結(jié)構(gòu)轉(zhuǎn)速相等,相位滯后角α為0,將激勵用復(fù)數(shù)形式表示可得

Fy=F0ej Ω t

(6)

Fz=-jF0ej Ω t

(7)

即y向、z向激勵力幅值相等且相位角相差90°。

2.2 回旋振動響應(yīng)計算

在諧響應(yīng)分析中需要考慮系統(tǒng)阻尼。根據(jù)實測和相關(guān)資料,鋼結(jié)構(gòu)的阻尼比一般在0.01~0.02,將結(jié)構(gòu)阻尼系數(shù)設(shè)置為0.02。滑動軸承阻尼通常為103~105N,因此在軸承對應(yīng)的彈簧單元中設(shè)置軸承阻尼值為105N。

在Workbench的“Rotating Force”中,“Synchronous Ratio”定義激勵頻率與結(jié)構(gòu)轉(zhuǎn)速的比值,即諧響應(yīng)激勵頻率與結(jié)構(gòu)轉(zhuǎn)動頻率是否同步。其可對諧響應(yīng)的每個頻率步通過激勵頻率更新轉(zhuǎn)速,在該項缺省時則定義為不平衡激勵,否則為不隨轉(zhuǎn)速變化的一般轉(zhuǎn)動諧波激勵力。據(jù)此,將“Ratio”項設(shè)置為空,并輸入質(zhì)量m、不平衡質(zhì)量旋轉(zhuǎn)半徑r,即不平衡力F0=mr。ANSYS會自動在每一頻率步將轉(zhuǎn)速的平方項相乘進行運算求解。采用這種定義方式可避免激勵力手動計算及表格加載。在配重盤所在節(jié)點設(shè)置m=45 kg,r=0.03 m,可得旋轉(zhuǎn)不平衡力F0=1.35 kg。

在實測中軸承所在節(jié)點位移響應(yīng)無法測量,因此在后尾軸承前端10 mm的軸上設(shè)置1個節(jié)點,并計算其位移響應(yīng),方便與試驗數(shù)據(jù)進行對比。在“Analysis Setting”中將最小頻率設(shè)置為1.33 Hz,最大頻率設(shè)置為3.33 Hz,迭代次數(shù)設(shè)置為13。求解測點在80~200 r/min穩(wěn)定工作轉(zhuǎn)速下的垂直和水平方向振動位移響應(yīng)值。

由響應(yīng)仿真計算可得:在激勵頻率為3.33 Hz時,即軸系臺架在200 r/min額定轉(zhuǎn)速工況下運轉(zhuǎn)時,軸上測點的橫向與垂向振動位移響應(yīng)值均為0.010 02 mm。由于激勵力是y、z方向幅值相等,且為相位角相差90°的旋轉(zhuǎn)不平衡力,因此橫向、垂向位移響應(yīng)幅值相等且橫向響應(yīng)相位角比垂向滯后90°。

2.3 軸承基座響應(yīng)計算

求得軸上測點響應(yīng)后,計算在額定工況下不平衡力激勵時軸承座上測點的加速度響應(yīng)。通過插入探針獲取后尾軸承處彈簧單元的響應(yīng)變形量和頻域響應(yīng)力,即為后尾軸承受到的頻域力。在后尾軸承支承中心施加此頻域力并再次進行諧響應(yīng)計算,求解軸承基座處測點的加速度響應(yīng)。新建兩個參考坐標系,切分出支承中心所在截面和下軸瓦的受力面,并在軸承頂端、軸承兩層基座上建立與試驗相同的3個測點。在振動試驗中后尾軸承基座的測點布置如圖6所示。

在進行后尾軸承的諧響應(yīng)分析時,激勵為計算回旋振動響應(yīng)所得到的頻域力。其中,頻域激勵力的y、z向幅值是相等的,且z向的相位角比y向滯后90°。考慮到軸承真實受力情況僅為下軸瓦約30°~45°部分受力,而非整個軸承軸瓦面均受力,因此選取受力面時只取下軸瓦面的一部分。

圖6 后尾軸承基座測點布置

軸承的真實受力位置,即支承中心所在位置,并不在軸向中心1/2處。針對后尾軸承,通常軸向約1/3處才是其真實支承中心所在截面,因此軸承激勵力不能直接通過“Force”進行簡單的加載,而是需要通過“Remote Force”指定其真實受力點確定支承中心所在位置,再將激勵力耦合至下軸瓦面,并通過力的平移定理將力和力矩傳遞到軸承的下軸瓦面上。將遠端力的施力點選在軸向1/3處所在截面的中心。激勵的施力點與下軸瓦受力面的耦合、頻域激勵力加載情況如圖7所示。

此時激勵力的加載,不僅考慮軸承真實受力面情況,而且考慮軸承支承點位置選取對結(jié)果的影響,更符合實際情況。由于實體單元的響應(yīng)計算較復(fù)雜,網(wǎng)格密度較小,所需計算時間較長,并且考慮到測試工況的穩(wěn)定程度會對結(jié)果產(chǎn)生影響,因此選取軸系試驗臺架200 r/min的額定轉(zhuǎn)速工況進行仿真計算與對比。將計算頻率設(shè)置為3.33 Hz,在后處理模塊查看3個測點y向、z向的響應(yīng)加速度。通過仿真計算所得,在200 r/min轉(zhuǎn)速工況下測點1的垂向、橫向加速度響應(yīng)值分別為0.129 20 m/s2、0.213 89 m/s2,測點2分別為0.236 31 m/s2、0.195 45 m/s2,測點3分別為0.107 18 m/s2、0.073 26 m/s2。

3 試驗與仿真數(shù)據(jù)分析

3.1 軸系振動測試系統(tǒng)

在臺架軸系回旋振動測試時,采用國內(nèi)外回旋振動測量通用的非接觸式法,采用多功能數(shù)據(jù)采集儀結(jié)合電渦流傳感器進行測量。振動測點布置在尾承前端10 cm處。在軸承及其支座振動測試中,采

圖7 遠端力耦合及頻域力加載情況

用加速度傳感器測量軸承及其支座的振動加速度信號,通過信號采集儀進行加速度信號采集,并連接至測試電腦。振動測試傳感器型號及參數(shù)如表2所示。振動測試系統(tǒng)如圖8所示。

表2 傳感器型號及參數(shù)

圖8 振動測試系統(tǒng)

3.2 位移響應(yīng)測試結(jié)果及對比分析

在臺架軸系位移響應(yīng)測試時,將電機轉(zhuǎn)速從80 r/min上升至額定轉(zhuǎn)速200 r/min,通過電渦流傳感器記錄軸上測點的振動位移,將測試結(jié)果與第2.2節(jié)中仿真計算所得同一點的振動位移隨轉(zhuǎn)速變化結(jié)果進行繪制,如圖9所示。圖9中虛線為軸上測點的振動位移,實線為有限元仿真計算結(jié)果。

由圖9可看出:軸的振動位移仿真分析結(jié)果與實測的趨勢是一致的,但實測軸的振動位移隨轉(zhuǎn)速的變化存在一定的波動,而仿真計算相對比較光滑。這主要是由于電渦流傳感器的測試結(jié)果是相對值,且軸運轉(zhuǎn)時的不穩(wěn)定性會直接導(dǎo)致軸的振動位移存在波動。同時,在穩(wěn)定運轉(zhuǎn)工況200 r/min下,后尾軸承前端10 mm軸上測點的振動位移仿真計算結(jié)果為0.010 02 mm,實測幅值為0.011 07 mm,二者相差較小,說明仿真計算模型、邊界條件的設(shè)置及加載方法是正確的。

3.3 軸承加速度響應(yīng)測試結(jié)果及對比分析

在軸承座加速度響應(yīng)測試時,將電機轉(zhuǎn)速設(shè)置為額定轉(zhuǎn)速200 r/min,通過軸承及基座上的加速

圖9 測點橫向振動響應(yīng)曲線圖

度傳感器記錄軸承振動加速度值。將第2.3節(jié)中仿真計算所得后尾軸承上3個振動測點的加速度響應(yīng)值與實測值進行對比,如表3所示。

表3 軸承振動結(jié)果對比表 m·s-2

由表3可看出:

(1) 軸承座橫向振動響應(yīng)加速度為軸承頂部最大(測點1),軸承座第1層次之(測點2),且軸承頂部與第1層基座橫向響應(yīng)值比較接近,均大于軸承座底部(測點3)。此相對關(guān)系與實測值一致。

(2) 軸承座垂向振動響應(yīng)加速度為軸承座第1層最大(測點2),軸承頂部次之(測點1),軸承座第2層最小(測點3),但三者在同一數(shù)量級,其中測點1、測點3仿真值均略大于實測值。對于測點2,其實測與仿真的振動響應(yīng)加速度均最大。

(3) 軸承座同一測點的橫向及垂向加速度響應(yīng)仿真與實測完全一致。軸承頂部的橫向響應(yīng)值大于垂向響應(yīng)值(測點1),而軸承基座第1層和第2層的垂向響應(yīng)值均大于橫向響應(yīng)值(測點2、測點3)。

(4) 對比仿真結(jié)果與試驗結(jié)果可知:仿真與實測在數(shù)值上存在一定的差異,但其結(jié)果都在同一數(shù)量級,且仿真與實測除垂向的個別測點大小關(guān)系存在差異外,軸承振動響應(yīng)加速度的仿真與實測結(jié)果在大小及相對關(guān)系上都是匹配的。分析原因,可能是仿真中的激勵力理想化為螺旋槳模擬配重盤的不平衡力,忽略其他激勵因素的影響,如齒輪嚙合、軸與軸承的摩擦、電機運轉(zhuǎn)的不平穩(wěn)性等。

4 結(jié) 論

(1) 通過有限元法仿真所得渦動頻率與傳統(tǒng)傳遞矩陣法所得固有頻率相對偏差小于1%,且較傳統(tǒng)傳遞矩陣法易于獲取軸系的橫向振動模態(tài),可較好地滿足推進軸系回旋振動模態(tài)的計算。

(2) 由于回轉(zhuǎn)效應(yīng)的影響,軸系回旋振動FW與BW的渦動頻率雖不同,但同階FW與BW的振型與幅值完全相同,只是方向不同。

(3) 通過所給邊界條件的設(shè)置方法、載荷加載處理等,方便地實現(xiàn)軸系回旋振動響應(yīng)計算,且通過實測驗證響應(yīng)分析的正確性,其方法可為軸系設(shè)計中回旋振動的響應(yīng)預(yù)報提供理論支撐。

(4) 通過對軸承座三維及有限元的精確建模,將軸系振動響應(yīng)等效為作用于軸承座上的激勵力,對軸承座振動響應(yīng)進行分析,實測表明:軸承基座上振動響應(yīng)的仿真計算與實測結(jié)果趨勢一致。

猜你喜歡
渦動平衡力軸系
你能區(qū)分平衡力與相互作用力嗎
臥式異步電機軸系支撐載荷研究
防爆電機(2022年3期)2022-06-17 01:41:24
平衡力與相互作用力辨析
平衡力與相互作用力辨
雙機、雙槳軸系下水前的安裝工藝
平衡力好,可以保命
女子世界(2017年3期)2017-03-13 00:38:31
BTA鉆桿渦動數(shù)學(xué)建模及實驗研究
軸系校中參數(shù)與軸系振動特性相關(guān)性仿真研究
基于ANSYS的高速艇艉軸架軸系振動響應(yīng)分析
船海工程(2015年4期)2016-01-05 15:53:26
理想條件下BTA鉆鉆桿的渦動分析
马山县| 麻阳| 平阴县| 花莲市| 宁海县| 隆德县| 江门市| 遂溪县| 苗栗县| 淮安市| 监利县| 尖扎县| 康保县| 丰台区| 澄城县| 桦南县| 八宿县| 阳城县| 陈巴尔虎旗| 息烽县| 甘洛县| 儋州市| 灵璧县| 项城市| 栾川县| 海阳市| 新建县| 嵊州市| 陵水| 萍乡市| 尼木县| 仪征市| 鲁甸县| 高要市| 乌什县| 金山区| 安化县| 墨竹工卡县| 吉木乃县| 五常市| 忻州市|