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

?

考慮受壓屈曲的圓鋼管桿單元等效彈塑性滯回模型

2012-02-13 11:56:36謝道清
振動與沖擊 2012年6期
關(guān)鍵詞:網(wǎng)殼彈塑性桿件

謝道清,沈 金,鄧 華,張 瑞

(1.浙江大學(xué) 空間結(jié)構(gòu)研究中心,杭州 310058;2.浙江大學(xué) 建筑設(shè)計研究院,杭州 310027)

網(wǎng)架、網(wǎng)殼等網(wǎng)格結(jié)構(gòu)的彈塑性地震響應(yīng)分析受到越來越多的重視,主要有兩個方面的原因。首先,作為應(yīng)用最廣的一類大跨度結(jié)構(gòu)形式,網(wǎng)格結(jié)構(gòu)被普遍認(rèn)為有相當(dāng)好的抗震性能。然而近些年的震害情況表明[1],網(wǎng)格結(jié)構(gòu)發(fā)生各種程度破壞的情況并不少見,而定量分析這些震害發(fā)生的原因(包括設(shè)計的不合理性)就不能回避彈塑性分析的問題。其次,對于按非地震作用或僅考慮小震作用來設(shè)計桿件截面的網(wǎng)格結(jié)構(gòu),在中震或大震條件下會呈現(xiàn)怎樣的動力響應(yīng)特點(diǎn)以及破壞形態(tài),這也必須進(jìn)行彈塑性分析。

關(guān)于網(wǎng)格結(jié)構(gòu)彈塑性地震響應(yīng)方面的研究總體上還比較薄弱,這也導(dǎo)致了我國相關(guān)抗震設(shè)計規(guī)范對網(wǎng)格結(jié)構(gòu)在中震和大震下的抗震措施基本上缺乏詳細(xì)規(guī)定。此外在具體進(jìn)行彈塑性分析時,目前的文獻(xiàn)中桿單元普遍采用的是理想彈塑性本構(gòu)模型(圖2),即拉壓極限狀態(tài)均為軸向屈服[2-4]。然而在中震、大震等強(qiáng)震作用下,網(wǎng)格結(jié)構(gòu)中桿件在受壓狀態(tài)下顯然發(fā)生的是屈曲破壞,因此理想彈塑性模型并不能客觀地反映壓桿的實(shí)際工作性態(tài)。正是出于對此基本問題的考慮,本文將對網(wǎng)格結(jié)構(gòu)中常用規(guī)格的圓鋼管桿件進(jìn)行不同長細(xì)比條件下的非線性屈曲分析,重點(diǎn)考察這些桿單元的受壓極限承載力、平衡路徑以及屈曲后卸載路徑,并繪制拉壓往復(fù)作用下桿件軸力P和兩端節(jié)點(diǎn)相對位移Δ的滯回曲線。進(jìn)一步統(tǒng)計這些不同規(guī)格鋼管桿單元的滯回曲線控制點(diǎn)與桿件長細(xì)比λ之間的關(guān)系,將P-Δ曲線轉(zhuǎn)換成統(tǒng)一描述的應(yīng)力-應(yīng)變(σε)曲線,以此提出一個能夠同時考慮受拉屈服和受壓屈曲的圓鋼管桿單元的等效彈塑性滯回模型。文中以一個球面網(wǎng)殼作為算例,應(yīng)用該等效彈塑性滯回模型進(jìn)行了網(wǎng)殼在罕遇地震作用下的彈塑性時程響應(yīng)計算。通過與理想彈塑性模型計算結(jié)果的比較,分析了桿件應(yīng)力應(yīng)變反應(yīng)和結(jié)構(gòu)薄弱區(qū)域分布上的差異,以此來考察本文提出的等效彈塑性滯回模型的有效性。

1 桿件有限元分析

1.1 考察的鋼管規(guī)格

根據(jù)網(wǎng)格結(jié)構(gòu)常用鋼管規(guī)格,選擇表1所列的46根不同截面和長度的圓鋼管桿件作為考察拉壓滯回規(guī)律的樣本桿件。鋼管材質(zhì)為Q235鋼,屈服強(qiáng)度fy=235 MPa,彈性模量E=2.06e5 MPa,泊松比為 0.3。表1中樣本桿件的長細(xì)比范圍在40~180之間。

表1 樣本桿件的長細(xì)比λ及其對應(yīng)的截面和長度Tab.1 Slenderness ratios of steel circular-tube specimens and their corresponding cross sections and lengths

1.2 有限元模型

鋼管桿件的計算簡圖如圖1所示,兩端鉸接,受軸向荷載P作用。采用LS-DYNA軟件建立鋼管桿件的有限元模型,管壁采用Shell163殼單元模擬。殼單元分別沿截面周圈和桿長方向按16和80等分劃分,單元總數(shù)共計1 280個。殼單元采用理想彈塑性本構(gòu)模型(圖2),符合Von Mises屈服準(zhǔn)則。受壓屈曲分析時計入桿件初彎曲的影響[5],初彎曲形狀按該殼單元有限元模型的一階特征值屈曲模態(tài)引入,最大撓度取v0=L/1 000[6]。

1.3 桿件加載和卸載

使用LS-DYNA軟件進(jìn)行加卸載計算。加卸載按位移控制,速度為10-3m/s。首先計算桿件受壓屈曲前后的平衡路徑,并獲得受壓臨界荷載Pcr。受壓屈曲后的最大加載位移為8Δcr,其中Δcr為Pcr對應(yīng)的桿端相對位移。然后再考察桿件在屈曲后平衡路徑上2Δcr、4Δcr、8Δcr位移點(diǎn)處的卸載情況,跟蹤卸載路徑直至受拉屈服。最后再從受拉屈服段反向卸載并重新加壓至桿件屈曲,以得到完整的拉壓P-Δ滯回曲線。

1.4 拉壓滯回特點(diǎn)

在46根樣本桿件中選擇三根典型長細(xì)比桿件來闡述圓鋼管桿單元的加載和卸載滯回曲線特點(diǎn)。三根桿件的規(guī)格分別為:Ф180×8(L=2.5 m,λ =41)、Ф140×4.5(L=5.0 m ,λ =104)和 Ф75.5 ×3.75(L=4.5 m,λ=177)。圖3給出了三根桿件在不同加、卸載情況下的P/Py-Δ/Δcr曲線,其中Py=-Afy,A為桿件截面積。為方便對比,故圖3中的縱橫坐標(biāo)均用相對值來表達(dá),如橫坐標(biāo)Δcr/Δ=1對應(yīng)的縱坐標(biāo)值為Pcr/Py。

圖1 鋼管桿件的計算模型Fig.1 The computational model of circular-tube members

圖2 理想彈塑性模型Fig.2 The ideal elasto-plastic model

(1)單向受壓加載

圖3(a)為三根桿件僅在壓力加載下的P/Py-Δ/Δcr曲線??梢钥闯?,小長細(xì)比(λ=41)桿件的屈曲臨界荷載與極限屈服荷載之比(Pcr/Py)高達(dá)0.988 7,到達(dá)臨界點(diǎn)前基本呈線彈性反應(yīng),但屈服后承載力較快下降。對于中等長細(xì)比(λ=104)桿件,Pcr/Py下降至0.641 6,屈曲前還基本呈線彈性反應(yīng),屈曲后承載力也較快下降,但后段的下降速度相對趨緩。對于λ高達(dá)177的桿件,Pcr/Py下降到0.278 6,屈曲前后曲線過渡平緩,臨界點(diǎn)并不明顯,且屈曲后路徑下降平緩。

(2)2Δcr處卸載

圖3 三根典型桿件在加卸載作用下的P/Py-Δ/Δcr曲線Fig.3 The P/Py-Δ/Δcrcurves of three circular-tube specimens under loading and unloading

在三根桿件受壓屈曲后變形至2Δcr時,進(jìn)行卸載并拉至屈服,然后再反向卸載和加壓以形成一個拉壓循環(huán),其P/Py-Δ/Δcr滯回曲線如圖3(b)所示??梢钥闯?,屈曲后路徑上的卸載曲線并不像理想彈塑性模型那樣呈直線變化。對于小長細(xì)比(λ=41)桿件,卸載初期基本還按直線卸載,曲線斜率接近彈性加載斜率,但是在受拉狀態(tài)后期至屈服前,曲線斜率快速減小。而大長細(xì)比(λ=177)桿件的卸載特點(diǎn)卻相反,卸載曲線前期斜率小,在進(jìn)入受拉區(qū)后斜率增大并總體呈線彈性變化。對于中等長細(xì)比(λ=104)桿件,卸載曲線在到達(dá)受拉屈服前斜率雖有所變化,但幅度很小。

當(dāng)卸載至受拉屈服后,再反向卸載,此時受拉區(qū)的P/Py-Δ/Δcr曲線基本上和理想彈塑性模型一樣呈線性變化。繼續(xù)反向卸載進(jìn)入受壓狀態(tài)時,屈曲前后曲線形狀也基本和初始的加壓曲線相同,無非是由于受拉屈服使得曲線沿橫坐標(biāo)向左平移了一個相對伸長值。

(3)4Δcr和 8Δcr處卸載

圖3(c)和(d)分別為三根桿件在受壓屈曲后路徑上的4Δcr和8Δcr位移點(diǎn)進(jìn)行卸載,并反向加載的P/Py-Δ/Δcr滯回曲線圖。總體上看,滯回曲線的變化規(guī)律和2Δcr處卸載的滯回曲線相差不大。小長細(xì)比桿件卸載曲線在受壓區(qū)斜率較大,在受拉區(qū)后期斜率變小。而大長細(xì)比桿件在卸載初期曲線斜率較小,進(jìn)入受拉區(qū)后斜率增大。中等長細(xì)比桿件的卸載曲線斜率總體變化依然相對較小。此外,8Δcr位移點(diǎn)的卸載曲線斜率相對4Δcr位移點(diǎn)的卸載曲線變化平緩些。當(dāng)桿件由受拉屈服段反向卸載至受壓屈曲時,P/Py-Δ/Δcr曲線變化規(guī)律和2Δcr處卸載的滯回曲線基本一致。

2 等效彈塑性滯回模型

2.1 基本思路

從圖3中可以看出,考慮受壓屈曲的圓管桿單元P-Δ滯回曲線與理想彈塑性模型的主要不同之處是受壓區(qū)的臨界點(diǎn)、屈曲后路徑以及屈曲后的卸載路徑,而且以上三方面內(nèi)容與桿件長細(xì)比λ密切相關(guān)。而在受拉區(qū),P-Δ曲線是基本符合理想彈塑性模型的。因此為便于計算機(jī)處理,可采用如圖4所示(以拉為正、壓為負(fù))的分段線性化辦法來描述受壓屈曲后路徑和卸載路徑。首先,將桿件屈曲后路徑用 Δcr、2Δcr、4Δcr、8Δcr對應(yīng)的坐標(biāo)點(diǎn)來分段模擬,即圖4中的3~6點(diǎn)。其次,考慮到不同長細(xì)比桿件屈曲后卸載曲線的斜率變化情況,將 2Δcr、4Δcr和 8Δcr處的卸載曲線按兩個直線段來模擬,即通過增加如圖4所示的7~9三個過渡點(diǎn)來定義。如果屈曲后的卸載點(diǎn)在3~6點(diǎn)之間,可相應(yīng)在2、7~9點(diǎn)間按等比例線性插值來確定該卸載曲線的過渡點(diǎn),并最終定義其卸載路徑。至此,考慮受壓屈曲的圓鋼管桿單元滯回模型實(shí)際上就可通過圖4中的1~9號點(diǎn)來描述。

圖4 等效彈塑性滯回模型Fig.4 The equivalent elasto-plastic hysteretic model

進(jìn)一步將表1所列46根樣本桿件的P-Δ曲線按圖4形狀進(jìn)行無量綱歸一化處理,即縱坐標(biāo)表示桿件的名義應(yīng)力σ=(P/A)和fy的比值,橫坐標(biāo)表示桿件的名義應(yīng)變ε=(Δ/L)和εy的比值,其中εy=fy/E。通過歸一化,易知九個控制點(diǎn)中,點(diǎn)1和點(diǎn)2的坐標(biāo)分別為(1,1)、(0,0)。設(shè)臨界點(diǎn) 3 的坐標(biāo)為(γ1,α1),則點(diǎn)4 ~6的坐標(biāo)可表示為(2γ1,α2)、(4γ1,α3)和(8γ1,α4)。根據(jù)對46根桿件的歸一化滯回曲線分析,確定卸載路徑上7~9過渡點(diǎn)的相對應(yīng)變值取(1-2α1)/2、(1-8α1)/3、(1-24α1)/4比較恰當(dāng),而此三點(diǎn)的相對應(yīng)力值分別定義為 β2、β3、β4,如圖 4。

2.2 控制參數(shù)的擬合

在圖4定義的圓鋼管桿單元等效彈塑性滯回模型中,有八個控制參數(shù)需要確定,即 γ1、α1、α2、α3、α4、β2、β3、β4。于是在46根樣本桿件的歸一化滯回曲線上,確定該八個控制參數(shù)并統(tǒng)計其與長細(xì)比λ的關(guān)系,見圖5(a)~(h)所示。可以看出,有些控制參數(shù)和長細(xì)比λ的關(guān)系為單調(diào)遞減,而有些為先遞減再遞增,因此統(tǒng)一采用簡單的多項(xiàng)式形式進(jìn)行擬合,并發(fā)現(xiàn)三次曲線就能達(dá)到較好的擬合精度。于是,進(jìn)一步可得出以上八個參數(shù)和長細(xì)比λ∈[40,180]的近似關(guān)系為:

圖5 八個控制參數(shù)和長細(xì)比λ的關(guān)系Fig.5 The relationships between 8 governing parameters and slenderness ratios λ

2.3 等效彈塑性模型的校驗(yàn)

用表1所列樣本桿件以外的多根圓鋼管桿件來考察等效彈塑性滯回模型的有效性。首先采用殼單元有限元法計算出這些桿件的拉壓滯回曲線,然后將計算結(jié)果和等效彈塑性滯回模型的擬合結(jié)果進(jìn)行比較。結(jié)果表明,等效彈塑性模型的計算結(jié)果和有限元法的計算結(jié)果非常吻合,特別在控制點(diǎn)處??紤]到篇幅限制,本文僅給出了一根校驗(yàn)桿件Ф127×4.5(L=3.25 m,λ=75)的計算結(jié)果,見圖6。

其次,臨界點(diǎn)控制參數(shù)α1實(shí)際上與我國《鋼結(jié)構(gòu)設(shè)計規(guī)范》(GB50017-2003)中的壓桿穩(wěn)定系數(shù)φ有非常相近的含義,盡管后者是根據(jù)邊緣屈服準(zhǔn)則確定的。由于本文桿件的初始缺陷僅考慮了初彎曲,接近GB50017-2003規(guī)定的A類截面[7],故將 α1和A類截面的穩(wěn)定系數(shù)φ進(jìn)行比較,如圖7??梢园l(fā)現(xiàn),α1和φ在長細(xì)比λ∈[40,180]區(qū)域的誤差在5%之內(nèi),吻合很好。

3 罕遇地震作用下的球面網(wǎng)殼計算

3.1 網(wǎng)殼模型

圖8所示為一個Kiewitt型雙層球面網(wǎng)殼屋蓋。網(wǎng)殼跨度為60 m,矢高6 m,厚度2 m。屋蓋均勻支承在周邊12根10 m高的C30混凝土柱上,柱截面為1.2 m×0.8 m。網(wǎng)殼承受的荷載標(biāo)準(zhǔn)值為:靜荷載 0.5 kN/m2(不包括結(jié)構(gòu)自重);活荷載為1.0 kN/m2;溫度差為±20℃;風(fēng)荷載0.5 kN/m2??紤]以上非地震荷載的工況組合,采用Mstcad軟件對該網(wǎng)殼進(jìn)行滿應(yīng)力設(shè)計以確定桿件的截面,其中滿應(yīng)力設(shè)計時所選擇的鋼管截面規(guī)格如表1。

3.2 地震波

參照我國《建筑抗震設(shè)計規(guī)范》(GB50011-2001),對結(jié)構(gòu)進(jìn)行8度罕遇地震下的彈塑性時程分析。地震加速度時程采用Big-bear波[8],按三向輸入(X為主方向,如圖9)。截取了該波的前38 s記錄(共1 900個點(diǎn)),X向加速度最大值按8度罕遇地震取400 cm/s2,三向加速度峰值的比例調(diào)整系數(shù)為1∶0.85∶0.65,結(jié)構(gòu)阻尼比取0.035。在輸入Big-bear波之前的2 s內(nèi),先以漸增豎向加速度的方式施加重力荷載代表值(1.0靜+0.5活),然后再疊加上Big-bear波的加速度時程,以此同時考慮結(jié)構(gòu)的重力荷載代表值和地震動的作用。

3.3 結(jié)構(gòu)響應(yīng)計算

采用Ansys軟件進(jìn)行結(jié)構(gòu)在Big-bear波作用下的彈塑性時程計算。網(wǎng)殼桿件考慮兩種材料模型:模型一為圖2所示的理想彈塑性模型;模型二為圖4所示的考慮壓桿屈曲的等效彈塑性模型。在采用模型二計算時,利用Ansys軟件的APDL語言控制計算流程,即不斷記錄每一個荷載步所求得的桿件變形量,然后通過重啟動(restart)命令并利用前一荷載步的桿件變形量來修改材料本構(gòu)模型,以此實(shí)現(xiàn)等效彈塑性模型加、卸載迭代求解。

(1)失效桿件分布和薄弱區(qū)域分析

從理想彈塑性模型的計算結(jié)果發(fā)現(xiàn),在罕遇地震作用下網(wǎng)殼桿件始終處于線彈性階段,即最大拉壓應(yīng)力絕對值均沒有超過屈服應(yīng)力fy。而等效彈塑性模型的計算結(jié)果中,有110根桿件出現(xiàn)了受壓屈曲,并最終存在殘余塑性應(yīng)變(指名義應(yīng)變Δ/L)。發(fā)生受壓屈曲的桿件主要分布在由外至內(nèi)的2~4環(huán)的上下弦層內(nèi),尤其是下弦層以及部分支座附近(如圖9),其中最大名義塑性應(yīng)變(總名義應(yīng)變扣除彈性應(yīng)變)為-1.015×103με。進(jìn)一步考察這些桿件的非地震作用控制內(nèi)力,發(fā)現(xiàn)其中大多數(shù)桿件位于殼體邊緣彎矩效應(yīng)和薄膜內(nèi)力效應(yīng)相互抵消的區(qū)域,桿件(特別是弦桿)軸力較小,因此在滿應(yīng)力設(shè)計時其截面配置主要由長細(xì)比控制,截面普遍較小。

(2)桿件的應(yīng)力和應(yīng)變

圖10和11分別給出了名義塑性應(yīng)變最大的兩根桿件1 191、1 201(如圖9)的應(yīng)力、應(yīng)變時程曲線??梢钥闯?,采用等效彈塑性模型所求得的桿件名義壓應(yīng)變要大于理想彈塑性模型的計算結(jié)果,兩者最大值相差2倍以上。而采用等效彈塑性模型所求得的名義壓應(yīng)力要小于理想彈塑性模型的結(jié)果。究其原因,主要是因?yàn)榈刃椝苄阅P椭挟?dāng)桿件發(fā)生受壓屈曲后,其承載力和軸向剛度將大大低于不考慮屈曲的理想彈塑性模型。此外,理想彈塑性模型中構(gòu)件的時程響應(yīng)始終處于彈性狀態(tài),而等效彈塑性模型的結(jié)果是構(gòu)件最終存在明顯的殘余應(yīng)變和應(yīng)力。

4 結(jié)論

要研究網(wǎng)格結(jié)構(gòu)對應(yīng)于中震、大震設(shè)防水準(zhǔn)的抗震措施,就不能回避對結(jié)構(gòu)進(jìn)行動力彈塑性性能分析。而現(xiàn)代數(shù)值計算技術(shù)的快速發(fā)展,為網(wǎng)格結(jié)構(gòu)動力彈塑性響應(yīng)的精細(xì)化分析提供了有力支持。本文的目的就是試圖將強(qiáng)震作用下網(wǎng)格結(jié)構(gòu)中桿件受壓屈曲這一基本問題盡可能細(xì)致考慮,以克服理想彈塑性模型高估了桿件受壓承載力的缺點(diǎn)。較為有利的是,無論是受拉屈服還是受壓屈曲,桿單元實(shí)際上僅需通過桿件軸力和伸長量這兩個簡單參數(shù)來定義其受力狀態(tài),因此形式上就可以將屈服和屈曲問題進(jìn)行統(tǒng)一描述,而且是在構(gòu)件層面上。此外,計算分析也發(fā)現(xiàn)盡管考察的桿件規(guī)格不同,但它們的受壓屈曲路徑和屈曲后卸載路徑從大的趨勢上看是相同的,P-Δ曲線的控制參數(shù)也主要跟長細(xì)比相關(guān),這也為建立能同時考慮受拉屈服和受壓屈曲的圓鋼管桿單元等效彈塑性滯回模型提供了條件。

本文提出的等效彈塑性滯回模型與理想彈塑性模型的差別主要反映在桿件受壓時的加、卸載曲線上,特別是屈曲后的卸載路徑隨桿件長細(xì)比不同存在一些差異,即小長細(xì)比桿件的初期卸載剛度(曲線斜率)大而后期小,而大長細(xì)比桿件正好相反,因此在等效彈塑性模型中引入了八個控制參數(shù)來加以描述,形式上略顯復(fù)雜。但是,由于這八個參數(shù)僅與長細(xì)比相關(guān),在實(shí)際計算(包括插值計算)時并不復(fù)雜,用較短的程序段就可以實(shí)現(xiàn)。更為重要的是,只要在常規(guī)的結(jié)構(gòu)彈塑性分析中對材料的本構(gòu)模型進(jìn)行修改就能將桿件屈曲問題加以考慮,總體上應(yīng)該是方便的。此外,與理想彈塑性模型相比,本文算例結(jié)果表明了采用等效彈塑性滯回模型確實(shí)能夠有效地反映網(wǎng)格結(jié)構(gòu)在地震作用下的受壓屈曲桿件分布情況,并有助于分析結(jié)構(gòu)的薄弱區(qū)域。

[1]中國建筑科學(xué)研究院編.2008年汶川地震建筑震害圖片集[M].北京:中國建筑工業(yè)出版社,2008.

[2]范 峰.空間網(wǎng)殼結(jié)構(gòu)彈塑性地震響應(yīng)及抗震性能分析[J].哈爾濱建筑大學(xué)學(xué)報,1999,32(1):132-37.

[3]薛素鐸,王建寧,曹 資,等.鋼網(wǎng)殼彈塑性地震反應(yīng)分析[J].北京工業(yè)大學(xué)學(xué)報,2001,27(1):50-53.

[4]沈世釗,支旭東.球面網(wǎng)殼結(jié)構(gòu)在強(qiáng)震下的失效機(jī)理[J].土木工程學(xué)報,2005,38(1):11-20.

[5] Kala Z. Sensitivity assessmentofsteelmembersunder compression[J].Engineering Structures,2009,31:1344-1348.

[6]陳 驥.鋼結(jié)構(gòu)穩(wěn)定-理論與設(shè)計[M].北京:科學(xué)出版社,2001.

[7]夏志斌,姚 諫.鋼結(jié)構(gòu)——原理與設(shè)計[M].北京:中國建筑工業(yè)出版社,2010.

[8] PEER Strongmotion database[EB OL]. http://peer.berkeley.edu/products/strong_ground_motion_db.html.

猜你喜歡
網(wǎng)殼彈塑性桿件
基于臨時支撐結(jié)構(gòu)的桿件初彎曲對其軸壓性能的影響
四川建筑(2021年1期)2021-03-31 01:01:46
塔式起重機(jī)拼裝式超長附著桿設(shè)計與應(yīng)用
矮塔斜拉橋彈塑性地震響應(yīng)分析
基于CFD模擬的球面網(wǎng)殼風(fēng)壓分布分析
彈塑性分析在超高層結(jié)構(gòu)設(shè)計中的應(yīng)用研究
江西建材(2018年4期)2018-04-10 12:36:52
大型拱頂儲罐三角形板式節(jié)點(diǎn)網(wǎng)殼正裝施工工藝
KD379:便攜折疊式衣架
某網(wǎng)架桿件彎曲的原因分析及處理
地震動斜入射對樁-土-網(wǎng)殼結(jié)構(gòu)地震響應(yīng)影響
動載荷作用下冪硬化彈塑性彎曲裂紋塑性區(qū)
独山县| 锡林浩特市| 边坝县| 天柱县| 余姚市| 新巴尔虎右旗| 黄山市| 秦安县| 泰宁县| 子洲县| 洮南市| 丹棱县| 普兰店市| 岳普湖县| 广丰县| 陆河县| 安福县| 林口县| 临泉县| 鸡西市| 廊坊市| 磐安县| 苏尼特右旗| 红安县| 汤阴县| 扬州市| 塘沽区| 洛浦县| 津市市| 观塘区| 田林县| 三明市| 周至县| 上虞市| 清远市| 读书| 姚安县| 金湖县| 泗水县| 灵山县| 凤庆县|