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

?

考慮附屬管的實尺寸鉆井隔水管系統(tǒng)渦激振動二維數(shù)值模擬研究

2021-01-29 05:57:34孔騰騰王嘉松吳文波許亮斌盛磊祥李朝瑋
振動與沖擊 2021年2期
關(guān)鍵詞:約化角為渦激

孔騰騰,王嘉松,吳文波,許亮斌,盛磊祥,李朝瑋

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院, 上海 200240; 2. 廣州大學(xué) 機(jī)械與電氣工程學(xué)院,廣州 510006;3. 中海油研究總院有限責(zé)任公司,北京 100028)

渦激振動(vortex-induced vibration, VIV)是一種流體與固體相互作用耦合的現(xiàn)象。當(dāng)流體流經(jīng)非流線型鈍體時會在鈍體后方發(fā)生分離,并產(chǎn)生交替泄放的漩渦。周期性的漩渦脫落會產(chǎn)生作用于鈍體的周期性變化的升阻力,誘發(fā)鈍體振動,即渦激振動。當(dāng)渦脫頻率和振動頻率與結(jié)構(gòu)的固有頻率相近時,會引發(fā)結(jié)構(gòu)的大幅值振動,即共振,導(dǎo)致結(jié)構(gòu)的疲勞損壞和結(jié)構(gòu)失效,嚴(yán)重威脅生產(chǎn)生活安全。

許多學(xué)者將目光投向了雙圓柱和多圓柱的渦激振動問題。Zhao等[1-2]研究了串列和并列剛性連接雙圓柱在低雷諾數(shù)下的渦激振動響應(yīng)。在串列圓柱中,太大或太小的間距將導(dǎo)致圓柱的鎖定區(qū)增大或減小,但圓柱系統(tǒng)的振幅均大于單圓柱的振幅。在并列圓柱中,圓柱系統(tǒng)的鎖定區(qū)比單圓柱的鎖定區(qū)更小,且渦激振動和馳振將同時出現(xiàn)在圓柱的響應(yīng)中。Zhao等[3]研究了2個不同直徑圓柱錯列排布時的渦激振動特征,討論了圓柱間隙和來流攻角對渦激振動鎖定區(qū)的影響。Zhao等[4]采用有限單元法研究了正方形排布的四圓柱雙自由度渦激振動,研究發(fā)現(xiàn)來流攻角對四圓柱系統(tǒng)的振動響應(yīng)影響很大。Zhu等[5-6]采用SST(shear stress transfer)k-ω湍流模型研究了控制桿對圓柱系統(tǒng)渦激振動的影響,討論了雷諾數(shù)、控制桿位置角、旋轉(zhuǎn)速度、數(shù)量、直徑比、間隙比對主圓柱上流體力和渦激振動響應(yīng)的影響。Song等[7]研究了來流攻角和間隙比對附有3根控制桿的多圓柱系統(tǒng)渦激振動響應(yīng)的影響。結(jié)果表明,在來流攻角為45°、間隙比為0.9時的抑制效果最好。田啟龍等[8]采用離散渦方法對均勻布置4~10根等直徑附屬管的隔水管系統(tǒng)的渦激振動進(jìn)行了數(shù)值模擬研究。結(jié)果表明,附加附屬管后隔水管與裸管相比受力系數(shù)與振幅均有不同程度的減小,對隔水管渦激振動的抑制效果明顯。不同數(shù)量的附屬管對來流攻角的敏感性也各不相同,其中附加10根附屬管對渦激振動的抑制效果最好。Wu等[9]采用二次開發(fā)的OpenFOAM程序?qū)χ鲌A柱后剛性連接一根小控制桿的隔水管系統(tǒng)在低雷諾數(shù)下的渦激振動進(jìn)行了數(shù)值模擬研究。結(jié)果表明,間隙比越小可以激發(fā)更大范圍的鎖定狀態(tài),并且可以觀察到更高階的次諧波共振??刂茥U在小間隙比時對隔水管系統(tǒng)的渦激振動產(chǎn)生很大影響甚至起決定性作用。隨著間隙比的增加,控制桿的作用逐漸變?nèi)?。Silva-ortega等[10]通過水中試驗研究了在主圓柱周圍環(huán)形均勻剛性連接2根、4根和8根控制桿的多圓柱系統(tǒng),討論直徑比和間隙比對多圓柱系統(tǒng)渦激振動的影響。結(jié)果表明,在多數(shù)布置下控制桿能抑制渦激振動,但在一些情況下會產(chǎn)生類似于馳振的響應(yīng),布置8根控制桿時抑制效果最好。Jiménez-gonzález等[11]對下游剛性連接一對控制桿的多圓柱系統(tǒng)進(jìn)行水中試驗。結(jié)果表明,在所有布置下多圓柱系統(tǒng)的振幅均有所減小,尤其在直徑比為0.4D時振幅的抑制效果達(dá)到90%以上。Gao等[12]對低雷諾數(shù)下彈性連接的正方形排布四圓柱系統(tǒng)進(jìn)行了數(shù)值模擬,研究了來流攻角為0°~45°、約化速度為3~14的渦激振動響應(yīng)和尾渦模式。結(jié)果表明,來流攻角對渦激振動幅值和力系數(shù)均有顯著影響,尤其是對下游圓柱影響更為明顯。Huera-huarte等[13]試驗研究了固定圓心間距為1.3D時,串列剛性連接單一控制桿的雙圓柱系統(tǒng)的流致振動響應(yīng),在小于等于1的范圍內(nèi)改變下游圓柱與上游圓柱的直徑比。結(jié)果表明,隨著直徑比的增大渦激振動響應(yīng)逐漸減弱。Munir等[14]利用直接數(shù)值模擬方法對剛性連接的2個并列等直徑圓柱的渦激振動響應(yīng)進(jìn)行研究,討論了間距比和雷諾數(shù)對雙圓柱系統(tǒng)渦激振動鎖定區(qū)和渦脫的影響。

海洋鉆井平臺是深海油氣資源勘探開發(fā)的主要工具,而隔水管是連接深海鉆井平臺與海底井口的“咽喉”。在海洋石油鉆井操作中,常常需要在鉆井隔水管附近布置一些附屬管線,如壓井管線、泥漿增壓管線、化學(xué)管線等。在實際工程中為了安裝操作方便及成本考慮,附屬管通常直徑大小不同,且不均勻地緊密布置在隔水管主管周圍。隔水管附屬管主要從功能上考慮其布局,而不是從流動控制或渦激振動控制角度考慮。Wu等[15]從流動控制角度開展數(shù)值模擬研究,討論了來流攻角對多圓柱系統(tǒng)流體力和渦脫模式的影響。結(jié)果表明,附屬管線對隔水管主管有很好的流動控制效果,但受來流攻角的顯著影響,個別時候可能增大曳力,多數(shù)情況下可顯著減小升阻力。然而,實際附屬管線對鉆井隔水管渦激振動有怎樣的影響尚未有報道,但值得關(guān)注。本文以正服役于南海油氣開發(fā)的某實尺寸鉆井隔水管系統(tǒng)為模型,基于RANS(Reynolds-averaged Navier-Stokes)方程及k-ω湍流模型,采用二次開發(fā)的嵌入多圓柱運動求解模塊的OpenFOAM求解器,分析了實尺寸鉆井隔水管系統(tǒng)(含有錯列圓柱、不同間隙比、不同直徑比等特性)在不同來流攻角和雷諾數(shù)下的二維渦激振動響應(yīng)特征。

1 數(shù)學(xué)模型

本文采用的流動控制方程為二維Navier-Stokes方程結(jié)合k-ω湍流模型,有限體積法進(jìn)行離散求解。采用PIMPLE算法進(jìn)行壓力速度解耦。笛卡爾坐標(biāo)系中黏性不可壓縮流體非定常流動的控制方程的解析形式可以表示為

(1)

(2)

式中:U為流體的速度向量;x為空間坐標(biāo)向量;μ為分子動態(tài)黏度系數(shù);p為流體壓強(qiáng);ρ為流體密度;τ為雷諾應(yīng)力張量。這里使用k-ω湍流模型,用Boussinesq表達(dá)式在2個方程模型中計算雷諾應(yīng)力

(3)

(4)

式中:μt為湍流黏度系數(shù);δ為克羅內(nèi)克爾張量算子;k為湍動能;ω為湍流耗散率。湍流的k和ω的輸運方程可表示為

(5)

(6)

(7)

(8)

式中:S為流體的平均應(yīng)變率張量;σk=2.0;σω=2.0;γ1=0.553;β1=0.075;β*=0.09; 常量取值可參見文獻(xiàn)[16]。

考慮附屬管與主管整體一起運動,結(jié)構(gòu)控制方程采用單一彈簧振子模型對多圓柱同步運動系統(tǒng)進(jìn)行模擬

(9)

在圓柱渦激振動中,無量綱橫流向均方根振幅(Yrms·D-1)的表達(dá)式為

(10)

式中:D為實尺寸鉆井隔水管系統(tǒng)中主圓柱直徑,取D=0.533 4 m;N為穩(wěn)定后至少20個周期內(nèi)的數(shù)據(jù)個數(shù)。

本文采用的實尺寸鉆井隔水管系統(tǒng)的橫截面示意圖,如圖1所示。圖1表示當(dāng)來流攻角為0°時的布置情況。圖1中:θ為主管上的周向角;U∞為來流速度。參考實際工況中兩端鉸支的鉆井隔水管的第一階固有頻率約為fn=0.024 4。實尺寸鉆井隔水管系統(tǒng)中附屬管的直徑和幾何配置參數(shù),如表1所示。

圖1 實尺寸鉆井隔水管系統(tǒng)模型示意圖Fig.1 Sketch of the model for a real drilling riser system

表1 實尺寸鉆井隔水管系統(tǒng)中附屬管的直徑和位置參數(shù)

如圖2所示,α為來流方向與x軸正方向的夾角,即來流攻角。本文所采用的實尺寸鉆井隔水管系統(tǒng)模型可分別在順流向和橫流向振動。系統(tǒng)中的所有管線之間剛性連接,即隔水管與附屬管線在任意時刻同步運動,且圓柱系統(tǒng)在兩方向上的剛度系數(shù)ks和阻尼系數(shù)cs相同。海水中隔水管的質(zhì)量比較小,本文取質(zhì)量比為m*=3,阻尼比為ζ=5×10-4,剛度系數(shù)根據(jù)約化速度計算得到。

圖2 實尺寸鉆井隔水管系統(tǒng)渦激振動模型示意圖Fig.2 Sketch of the model for VIV of a real drilling riser system

本文計算域選取65D×30D的矩形域,主圓柱中心距入口邊界及上、下邊界的距離皆為15D,距出流邊界的距離為50D。采用非結(jié)構(gòu)網(wǎng)格離散計算域,在各圓柱近表面采用貼體網(wǎng)格處理,且第一層網(wǎng)格厚度滿足y+≈1,以準(zhǔn)確模擬圓柱邊界層中的流動。主圓柱周向網(wǎng)格數(shù)為200,各附屬管為140。入口邊界初始湍流強(qiáng)度設(shè)為1%,渦黏性系數(shù)比設(shè)為10。實尺寸鉆井隔水管系統(tǒng)的局部網(wǎng)格,如圖3所示。

圖3 實尺寸鉆井隔水管系統(tǒng)局部網(wǎng)格Fig.3 The mesh around the cylinders of a real drilling riser system

2 結(jié)果與討論

2.1 單圓柱渦激振動數(shù)值模擬結(jié)果

為驗證模型和計算方法,本文首先對Khalak等[17]的試驗?zāi)P瓦M(jìn)行了數(shù)值模擬,并將數(shù)值計算結(jié)果與試驗結(jié)果進(jìn)行了對比,Guilmineau等[18-19]的數(shù)值工作中也采用了相同的立管響應(yīng)模型參數(shù)。根據(jù)試驗條件:圓柱的直徑為D=0.05 m;固有頻率為fn=0.5;質(zhì)量比為m*=2.4;質(zhì)量阻尼比為m*ζ=0.013;f為圓柱在橫流方向上的振動頻率。通過改變雷諾數(shù)來改變來流的約化速度,從而觀測圓柱系統(tǒng)渦激振動隨約化速度的改變。采用k-ω湍流模型進(jìn)行數(shù)值模擬,入流邊界上的初始湍流強(qiáng)度為5%。圖4給出本文計算和其他數(shù)值模擬得到的渦激振動響應(yīng)與試驗結(jié)果的比較。

圖4 高雷諾數(shù)下單圓柱渦激振動的幅值和頻率Fig.4 The VIV amplitude and frequency of single circular cylinder at high Reynolds number

由圖4(a)可看出:本文的模型和計算方法能夠得到單圓柱渦激振動中的不同分支(初始支、上支和下支),能夠模擬初始支和上支之間的突變現(xiàn)象;本文結(jié)果比Khalak等研究的結(jié)果偏?。籊uilmineau等和Pan等研究的數(shù)值模型未能在Ur=5附近捕捉到大位移振幅,當(dāng)前數(shù)值模型的模擬結(jié)果優(yōu)于這些模型的結(jié)果。由圖4(b)中可看出,本文的渦激振動頻率與Khala等研究的結(jié)果和Guilmineau等及Pan等研究的結(jié)果在4

2.2 雙圓柱渦激振動數(shù)值模擬結(jié)果

為進(jìn)一步驗證模型和計算方法對多圓柱系統(tǒng)的適用性,本文對Huera-huarte等的試驗?zāi)P瓦M(jìn)行了數(shù)值模擬。單一控制桿位于主圓柱的下游,且與主圓柱剛性連接(兩圓柱同步運動)。主圓柱直徑為D,控制桿直徑為d,且圓柱之間的間距比為G/D=1.3。本文固定圓柱直徑比為d/D=0.2,固定來流速度為U∞。主圓柱直徑為D=0.05 m,基于主圓柱直徑的雷諾數(shù)為Re=13 000。圓柱系統(tǒng)的質(zhì)量比為m*=2.57,阻尼比為ζ=0.003 5。采用k-ω湍流模型,入口邊界上的初始湍流強(qiáng)度為3%。通過改變系統(tǒng)的剛度系數(shù)ks改變系統(tǒng)的固有頻率,進(jìn)而模擬不同約化速度下的渦激振動。采用的計算區(qū)域為55D×30D,主圓柱的中心位于原單圓柱中心位置。采用非結(jié)構(gòu)網(wǎng)格對計算域進(jìn)行離散,第一層網(wǎng)格厚度滿足y+≈1。

表2和表3為本文計算結(jié)果與試驗結(jié)果的比較。由表2和表3可知,本文結(jié)果與Huera-huarte等研究的結(jié)果契合度較高,即本文中的數(shù)值模型和計算方法能夠可靠地模擬高雷諾數(shù)下雙圓柱系統(tǒng)的渦激振動。

表2 串列圓柱系統(tǒng)無量綱渦激振動幅值

表3 串列圓柱系統(tǒng)無量綱渦激振動頻率

2.3 實尺寸鉆井隔水管系統(tǒng)的渦激振動數(shù)值模擬結(jié)果

為揭示實尺寸鉆井隔水管系統(tǒng)在海洋環(huán)境中的渦激振動特性,本文對正在南海服役中的某實尺寸鉆井隔水管系統(tǒng)進(jìn)行模型簡化。由于實尺寸鉆井隔水管系統(tǒng)結(jié)構(gòu)復(fù)雜、所處環(huán)境變化大、計算量大、耗時長,因此,只對其在不同來流攻角(0°,30°,60°,90°,…,330°)下的典型約化速度(Ur=3,5,7,9,12,即典型雷諾數(shù)相當(dāng)于21 000,35 000,49 000,63 000,83 000)時的渦激振動進(jìn)行二維數(shù)值模擬研究,討論來流攻角和約化速度對其振幅、振動頻率和在水平面上的運動軌跡的影響。

2.3.1 渦激振動幅值

圖5為實尺寸鉆井隔水管系統(tǒng)在不同來流攻角下,橫流向無量綱渦激振動幅值隨約化速度的變化。以下所討論的渦激振動強(qiáng)弱均指隔水管渦激振動幅值的大小。為分析比較,圖中實線為與主管直徑相同的單圓柱的橫流向渦激振動幅值。由圖5可知,由于附屬管分布的多管、多尺寸和不規(guī)則排布等特性,實尺寸鉆井隔水管系統(tǒng)的渦激振動幅值變化規(guī)律差異明顯。在約化速度3

圖5 實尺寸鉆井隔水管系統(tǒng)渦激振動幅值隨約化速度的變化Fig.5 The VIV amplitude of the drilling riser system with reduced velocity

當(dāng)來流攻角為0°和30°時,隔水管系統(tǒng)渦激振動幅值先隨約化速度的增大而增大到Ur=9,然后在9

圖6為實尺寸鉆井隔水管系統(tǒng)在不同約化速度下,橫流向無量綱渦激振動幅值隨來流攻角的變化。由圖6可知,在約化速度Ur=3時,鉆井隔水管系統(tǒng)在任意來流攻角下的渦激振動幅值都很小,來流攻角為330°時最大約為0.1D。隨約化速度增大,渦激振動幅值整體也隨之增大,在約化速度為Ur=7、來流攻角為150°時達(dá)到最大,但仍小于單圓柱的振幅。隨約化速度繼續(xù)增大,渦激振動幅值又整體減小。當(dāng)來流攻角為210°和330°時,任意約化速度下渦激振動幅值均很小,說明在此來流攻角下附屬管對鉆井隔水管的渦激振動抑制效果最好。

圖6 實尺寸鉆井隔水管系統(tǒng)渦激振動幅值隨來流攻角的變化Fig.6 The VIV amplitude of the drilling riser system with angle of attack

2.3.2 渦激振動頻率

圖7為實尺寸鉆井隔水管系統(tǒng)在不同來流攻角下的橫流向渦激振動位移功率譜。由圖7中可看出:當(dāng)來流攻角為0°,120°,180°和270°、約化速度為5,7和9時,功率譜有峰值較大的唯一主頻,且主頻大小接近于隔水管系統(tǒng)的固有頻率,表明此時振動較為劇烈,相應(yīng)的渦激振動幅值較大;當(dāng)來流攻角為30°,90°和150°、約化速度為5,7,9和12時,隔水管系統(tǒng)有峰值較大的唯一主頻,且主頻大小接近于隔水管系統(tǒng)的固有頻率,表明此時渦激振動較劇烈;當(dāng)來流攻角為60°、約化速度為5和7時,功率譜有峰值較大的唯一主頻,且主頻大小接近于隔水管系統(tǒng)的固有頻率,表明此時振動較為劇烈振動幅值較大,而在約化速度為9時功率譜有2個較為明顯的峰值,表明此時為多頻振動;當(dāng)來流攻角為210°和330°、功率譜只在約化速度為3時,有一個峰值較小的主頻,其余約化速度下未見明顯主頻,表明此時渦激振動整體較為微弱;當(dāng)來流攻角為240°、功率譜在約化速度為5和12時,各有一個峰值較大的主頻,且主頻大小接近于隔水管系統(tǒng)的固有頻率,表明此時渦激振動較為劇烈,相應(yīng)的振動幅值較大;當(dāng)來流攻角為300°、功率譜只在約化速度為5時,有一個峰值較大的主頻,且主頻大小接近于隔水管系統(tǒng)的固有頻率,表明此時渦激振動較劇烈。

圖7 實尺寸鉆井隔水管系統(tǒng)在不同來流攻角下的渦激振動頻譜Fig.7 The VIV frequency power spectrum of the drilling riser system at different angle of attack

2.3.3 渦激振動軌跡

圖8為實尺寸鉆井隔水管系統(tǒng)渦激振動穩(wěn)定后在水平面上的運動軌跡。由圖8可知,在任意來流攻角下隨約化速度的增大運動軌跡整體向右移動。這表明隨著來流速度的增大隔水管系統(tǒng)的順流向平均位置逐漸右移。由于附屬管的多尺寸和非均勻布置,隔水管系統(tǒng)渦激振動穩(wěn)定后的運動情況也較復(fù)雜,相應(yīng)的運動軌跡也復(fù)雜多變,呈現(xiàn)出扁平橢圓形、橢圓形、倒立雨滴形、雨滴形、傾斜“8”字形及雜亂不規(guī)則等形狀。

當(dāng)來流攻角為0°時,可以觀察到四種不同類型的運動軌跡:約化速度為3和12時雜亂無章;約化速度為5時呈近似橢圓形;約化速度為7時呈傾斜的倒立雨滴形;約化速度為9時呈略微傾斜的非對稱“8”字形。當(dāng)來流攻角為30°時,可以觀察到兩種運動軌跡:約化速度為3,7和9時雜亂無章;約化速度為5和12時呈傾斜的倒立雨滴形。當(dāng)來流攻角為60°時,可以觀察到三種不同類型的運動軌跡:約化速度為3和12時雜亂無章;約化速度為5和7時呈近似于一條斜線的扁平橢圓形;約化速度為9時呈扁平傾斜的倒立雨滴形。當(dāng)來流攻角為90°時,運動軌跡在所有約化速度下都雜亂無章。當(dāng)來流攻角為120°約化速度為5時,運動軌跡呈傾斜橢圓形,其他約化速度下雜亂無章。當(dāng)來流攻角為150°時,可以觀察到:約化速度為3,5,7和9時運動軌跡雜亂無章;約化速度為12時呈扁平的傾斜雨滴形。當(dāng)來流攻角為180°、約化速度為5,7和9時,運動軌跡呈傾斜的雨滴形,且約化速度為9時的運動軌跡處于約化速度為7時的包絡(luò)線內(nèi)。當(dāng)來流攻角為210°和330°時,渦激振動在所有約化速度下都很微弱,運動軌跡雜亂。當(dāng)來流攻角為240°約化速度為5時,運動軌跡呈傾斜的扁平橢圓形,其他約化速度下雜亂無章。當(dāng)來流攻角為270°時,可以觀察到三種不同的運動軌跡:約化速度為3和12時雜亂無章;約化速度為5時呈橢圓形;約化速度為7和9時呈雨滴形。當(dāng)來流攻角為300°、約化速度為5時,運動軌跡呈傾斜的扁平橢圓形,其他約化速度下雜亂無章。

圖8 實尺寸鉆井隔水管系統(tǒng)在不同來流攻角下的渦激振動軌跡Fig.8 The VIV trajectory of the drilling riser system at different angle of attack

3 結(jié) 論

本文基于RANS方程,采用k-ω湍流模型,用嵌入多圓柱運動求解模塊的OpenFOAM求解器,對不可壓縮流體進(jìn)行求解,再耦合結(jié)構(gòu)動力學(xué)模型進(jìn)行流固耦合求解。研究了亞臨界區(qū)雷諾數(shù)為20 000~85 000、典型約化速度為3~12、來流攻角為0°~330°的實尺寸鉆井隔水管系統(tǒng)的渦激振動問題。對不同來流攻角和約化速度下實尺寸鉆井隔水管系統(tǒng)的渦激振動幅值、頻率和運動軌跡進(jìn)行了分析研究。結(jié)果表明:

(1) 附屬管對鉆井隔水管的渦激振動總體上具有抑制作用,抑制效果與來流攻角和約化速度密切相關(guān)。當(dāng)來流攻角為210°和330°時,隔水管系統(tǒng)在所有約化速度下的渦激振動都很微弱,此時的抑制效果最好。

(2) 鉆井隔水管系統(tǒng)的橫流向渦激振動較劇烈時功率譜有峰值較大的唯一主頻。此時主頻大小接近于隔水管系統(tǒng)的固有頻率。

(3) 鉆井隔水管系統(tǒng)的渦激振動運動軌跡有多種形狀,如扁平橢圓形、橢圓形、倒立雨滴形、雨滴形、傾斜“8”字形及雜亂不規(guī)則等。

猜你喜歡
約化角為渦激
《空間點、直線、平面之間的位置關(guān)系》專題訓(xùn)練
不同間距比下串聯(lián)圓柱渦激振動數(shù)值模擬研究
約化的(3+1)維Hirota方程的呼吸波解、lump解和半有理解
渦激振動發(fā)電裝置及其關(guān)鍵技術(shù)
盤球立管結(jié)構(gòu)抑制渦激振動的數(shù)值分析方法研究
電子制作(2018年14期)2018-08-21 01:38:42
柔性圓管在渦激振動下的模態(tài)響應(yīng)分析
M-強(qiáng)對稱環(huán)
立體幾何測試卷(B卷)
一種以雙環(huán)鉑為有效成分的藥物檢測方法
(3+1)-維廣義Kadomtsev-Petviashvili方程的對稱約化與精確解
太仓市| 红安县| 田林县| 南木林县| 舟曲县| 莒南县| 迁安市| 奉新县| 墨竹工卡县| 常山县| 琼中| 太仓市| 余姚市| 霍邱县| 澄迈县| 巫溪县| 习水县| 松桃| 吐鲁番市| 舒城县| 青阳县| 广宁县| 长丰县| 张家港市| 大连市| 拜城县| 封开县| 嘉善县| 女性| 会泽县| 井冈山市| 溧阳市| 乐至县| 淮北市| 冕宁县| 射阳县| 蕲春县| 浦县| 永胜县| 麻城市| 广宁县|