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

?

基于波數(shù)-頻率譜的腔體流動自持續(xù)振蕩頻率研究

2021-09-27 06:39:04倪小松肖新標(biāo)劉浩楠
機(jī)械 2021年9期
關(guān)鍵詞:波數(shù)后緣腔體

倪小松,肖新標(biāo),劉浩楠

基于波數(shù)-頻率譜的腔體流動自持續(xù)振蕩頻率研究

倪小松,肖新標(biāo)*,劉浩楠

(西南交通大學(xué) 牽引動力國家重點(diǎn)試驗(yàn)室,四川 成都 610031)

在目前普遍認(rèn)可及廣泛應(yīng)用的腔體流動自持續(xù)振蕩噪聲頻率預(yù)測方法Rossiter半經(jīng)驗(yàn)公式中,只給出了經(jīng)驗(yàn)常數(shù)的推薦值,并沒有給出具體的計(jì)算方法。為此,本文采用數(shù)值模擬的方法確定值。首先,采用LES湍流模型建立腔體自持續(xù)振蕩噪聲預(yù)測模型。進(jìn)而,分析腔體內(nèi)部監(jiān)測點(diǎn),得到前三階噪聲峰值頻率,并與取推薦值(0.57)的Rossiter公式預(yù)測頻率進(jìn)行對比,預(yù)測結(jié)果存在一定的誤差。然后,通過分析腔體后緣湍流脈動壓力的波數(shù)-頻率譜,得到值為0.63并代入Rossiter公式,得到前三階噪聲峰值頻率并與數(shù)值模擬結(jié)果吻合較好。結(jié)果表明,利用分析湍流脈動壓力的波數(shù)-頻率譜確定經(jīng)驗(yàn)常數(shù)的方法是可行的,并且能夠提高Rossiter公式預(yù)測頻率的準(zhǔn)確性。

腔體流致噪聲;峰值頻率;Rossiter公式;波數(shù)-頻率譜

腔體結(jié)構(gòu)廣泛存在于航空航天、航海、汽車、軌道交通等領(lǐng)域,如戰(zhàn)機(jī)上的武器艙、船舶上的流水孔、汽車的天窗以及高速列車的車廂連接處。當(dāng)腔體結(jié)構(gòu)高速運(yùn)行時,來流撞擊到腔體后緣,復(fù)雜的反饋過程會使空腔上方剪切層形成持續(xù)振蕩,從而導(dǎo)致空腔產(chǎn)生強(qiáng)烈的噪聲,屬于簡單的流體動力學(xué)作用。當(dāng)剪切層的自持續(xù)振蕩頻率與腔體結(jié)構(gòu)的聲學(xué)特性(Helmholtz共振或聲學(xué)駐波模態(tài))發(fā)生耦合時,腔體的噪聲峰值將進(jìn)一步加大,屬于流體聲學(xué)共振作用。當(dāng)剪切層的自持續(xù)振蕩頻率與腔體自身結(jié)構(gòu)的振動固有模態(tài)發(fā)生耦合時,會產(chǎn)生強(qiáng)烈的峰值噪聲且加速腔體結(jié)構(gòu)的疲勞,屬于流體彈性作用[1-2]。在流體聲學(xué)作用和流體彈性作用中,剪切層自持續(xù)振蕩都起到了一個“放大”的作用。為控制腔內(nèi)噪聲以及延長腔體結(jié)構(gòu)壽命,剪切層自持續(xù)振蕩頻率和機(jī)理一直是廣大學(xué)者研究的熱點(diǎn)。1964年,Rossiter[3]在對投彈艙噪聲進(jìn)行研究時推導(dǎo)出了預(yù)測剪切層自持續(xù)振蕩頻率的半經(jīng)驗(yàn)公式,并將經(jīng)驗(yàn)公式與試驗(yàn)數(shù)據(jù)對比,證明了公式的準(zhǔn)確性。

Rossiter公式自被提出以來,就受到相關(guān)學(xué)者的廣泛應(yīng)用和認(rèn)可,但對于公式中兩個經(jīng)驗(yàn)常數(shù)(相位延遲系數(shù))和(湍流運(yùn)動速度與自由來流速度之比)Rossiter并未明確給定。一般情況下取0.25[3-4];而的取值對預(yù)測結(jié)果有著重大影響,實(shí)際上是一個范圍值,它與腔體前緣湍流邊界層厚度和腔體開口長度有關(guān),一般需要通過試驗(yàn)測得。East等[4]通過試驗(yàn)測得該值的范圍為0.35~0.6;Ma等[5]研究以往的流致共振總結(jié)出的范圍則比East的小,為0.22~0.43。

通過風(fēng)洞試驗(yàn)的方法雖然能夠準(zhǔn)確計(jì)算出值,但腔體結(jié)構(gòu)樣件制作麻煩,且風(fēng)洞試驗(yàn)成本昂貴,會耗費(fèi)大量的財(cái)力物力。隨著計(jì)算機(jī)的發(fā)展,數(shù)值模擬能大幅度降低研發(fā)成本,提高產(chǎn)品質(zhì)量[6]。CFD(Computational Fluid Dynamics,計(jì)算流體動力學(xué))已經(jīng)逐漸成為研究腔體流動效應(yīng)的重要途徑。Daude等[7]采用大渦模擬(LES,Large Eddy Simulation)的方法對空腔流動中不穩(wěn)定性的被動控制進(jìn)行了研究并與試驗(yàn)數(shù)據(jù)對比,證明了LES方法的可靠性。楊國晶[8]采用LES方法對陷下式腔體的流動特性進(jìn)行了一系列的參數(shù)調(diào)查,并與試驗(yàn)結(jié)果相比較,證明了數(shù)值模擬的準(zhǔn)確性。Li等[9]采用隱式LES方法研究了矩形腔體在來流馬赫數(shù)大于1時渦脫落和聲激勵的關(guān)系,結(jié)果驗(yàn)證了聲激勵的反饋是腔體形成自持續(xù)振蕩的主要原因。鄧玉清等[10]采用LES方法對腔體壁面湍流脈動壓力進(jìn)行監(jiān)測并得到其波數(shù)-頻率譜,并從波數(shù)頻率譜中得到湍流運(yùn)動速度和能量分布等結(jié)果。

可見,數(shù)值模擬在研究腔體流動問題上得到了學(xué)者的認(rèn)可,且在試驗(yàn)中也得到了驗(yàn)證。因此,在腔體設(shè)計(jì)前期為了節(jié)約風(fēng)洞試驗(yàn)成本,本文提出采用LES方法,通過計(jì)算腔體后緣壁面湍流脈動壓力,得到湍流脈動壓力的波數(shù)-頻率譜,通過分析湍流脈動壓力的波數(shù)-頻率譜,得到湍流的運(yùn)動速度,從而確定Rossiter公式中的值,提高渦脫頻率預(yù)測的準(zhǔn)確性。

1 Rossiter公式與計(jì)算方法

1.1 Rossiter公式

渦撞擊到腔體后緣產(chǎn)生的聲渦將會向上游運(yùn)動,引發(fā)腔體前緣渦脫落,該渦運(yùn)動到腔體后緣時又會產(chǎn)生聲渦,如此往復(fù)循環(huán)形成腔體的自持續(xù)振蕩,如圖1所示。圖中:λ為聲波波長,λ為渦的波長,即兩個渦之間的間距。在初始時間=0時,向上傳播的聲波的相位假定為零,并且腔體跨越n個完整的波長λ,這時,距離腔體后緣γλγ為兩運(yùn)動之間的相位差)處有一個渦;當(dāng)=時,渦繼續(xù)以運(yùn)動速度kU向下游行進(jìn),向上游傳播的聲到達(dá)上游進(jìn)而誘發(fā)新渦的產(chǎn)生。

以上渦和聲的運(yùn)動,存在著以下的關(guān)系:

式中:為腔體自持續(xù)振蕩頻率,Hz;為聲速,m/s;U為渦的運(yùn)動速度,m/s;U聲渦的運(yùn)動速度,m/s;m、m分別為聲波和渦的波數(shù),皆為整數(shù)。

圖1 Rossiter渦脫分析示意圖[3]

結(jié)合式(1)~(3)可得:

令:

式(4)即可化為Rossiter公式:

式中:為自由來流速度,m/s;為腔體流向長度,m;為模態(tài)階數(shù);常數(shù)為相位延遲系數(shù);為湍流運(yùn)動速度與自由來流速度之比;為馬赫數(shù)。

1.2 大渦模擬方法(LES)

LES是對湍流脈動的一種空間平均,通過采用某種濾波函數(shù)將大尺度的渦和小尺度的渦分開,對大尺度的渦直接進(jìn)行Navier-Stokes方程的求解,小尺度的渦用模型來封閉,從而減小計(jì)算成本。對于過濾掉的小尺度渦,通過引入亞格子(SubGrid-Scale,SGS)模型附加應(yīng)力項(xiàng)來體現(xiàn)其對大尺度渦的影響。

則:

對Navier-Stokes方程進(jìn)行濾波,則有:

根據(jù)Smagorinsky的基本SGS模型,假定SGS應(yīng)力具有下面的形式:

2 簡易腔體的建模

研究表明[11],腔體長深比小于1時,腔體剪切層自持續(xù)振蕩很容易與腔體自身的聲學(xué)模態(tài)發(fā)生耦合作用。為避免發(fā)生流體聲學(xué)作用而僅考慮流體動力學(xué)作用,本文腔體計(jì)算模型長寬比1.5:1。如圖2所示,腔體長900 mm、高600 mm。在腔體內(nèi)部設(shè)置了1個監(jiān)測點(diǎn)用來監(jiān)測腔體內(nèi)部脈動壓力的變化。為得到湍流渦運(yùn)動速度,在腔體后緣中心線位置設(shè)置了一個湍流脈動壓力監(jiān)測陣列用以進(jìn)行波數(shù)-頻率譜分析,陣列由48個監(jiān)測點(diǎn)組成,相鄰兩個監(jiān)測點(diǎn)的間距為4.22 mm,如圖3所示。

圖2 腔體幾何模型(單位:mm)

圖3 計(jì)算域與脈動壓力監(jiān)測點(diǎn)(單位:mm)

如圖4所示,采用ICEM進(jìn)行網(wǎng)格劃分,整個計(jì)算域采用非結(jié)構(gòu)網(wǎng)格劃分,最大網(wǎng)格尺寸為70 mm,為了準(zhǔn)確捕捉到腔體內(nèi)部和附近的流場特征,對腔體附近流場區(qū)域采用網(wǎng)格加密處理,加密網(wǎng)格尺寸大小為15 mm,總網(wǎng)格數(shù)量在300萬左右。

圖4 腔體數(shù)值模擬網(wǎng)格劃分

進(jìn)行數(shù)值計(jì)算時,考慮氣體的可壓縮性,將氣體屬性設(shè)置為理想氣體。計(jì)算域入口設(shè)置為質(zhì)量氣流入口,出口設(shè)置為壓力出口。設(shè)置為先采用可實(shí)現(xiàn)的-兩方程模型進(jìn)行穩(wěn)態(tài)流場計(jì)算,計(jì)算結(jié)果作為瞬態(tài)計(jì)算的初始值。然后采用LES方法計(jì)算瞬態(tài)流場。瞬態(tài)時間步長可以直接確定空氣動力學(xué)噪聲的最大分析頻率,因此時間步長設(shè)置必須充分考慮研究的頻率范圍。根據(jù)Nyquist的采樣定理:要獲得原始信號而不會從采樣信號中丟失信號,采樣率必須至少是最大分析頻率的兩倍。由于腔體流致噪聲以中低頻為主,本文計(jì)算分析的最高頻率設(shè)定為=1000 Hz,為獲得更準(zhǔn)確的計(jì)算結(jié)果,采樣率設(shè)置為2000 Hz。時間步長Δ=0.0005 s,迭代步數(shù)為20步,總采樣時間為1 s。

3 計(jì)算結(jié)果分析

3.1 腔體內(nèi)部聲壓級

當(dāng)自由來流速度為200 km/h,首先得到腔體內(nèi)部監(jiān)測點(diǎn)脈動壓力的時域數(shù)據(jù),然后對其進(jìn)行FFT(Fast Fourier Transform,快速傅立葉變換)變換得到腔體內(nèi)部監(jiān)測點(diǎn)的聲壓級頻譜,如圖5所示。

由圖5可看到,腔體內(nèi)部測點(diǎn)位置存在三階峰值噪聲頻率,分別為27.3 Hz、62.5 Hz和99.6 Hz。故猜測最大渦脫落頻率為99.6 Hz,根據(jù)頻率與周期間的關(guān)系=1/,計(jì)算得到渦脫落最小周期為0.01 s?;赒準(zhǔn)則作出一個周期的渦圖,如圖6所示。

圖5 腔體內(nèi)部監(jiān)測點(diǎn)聲壓級頻譜

圖6 一個周期內(nèi)的渦量圖

由圖6可看出,在整個周期內(nèi)共有三個主渦:當(dāng)=0時,第三個渦剛好撞擊到腔體后緣,腔體前緣新渦還沒有生成;當(dāng)=/4時,第三個渦往腔體后緣運(yùn)動,腔體前緣新渦逐漸生成;當(dāng)=/2時,新渦已經(jīng)生成并脫離腔體前緣,此時腔體上方重新有三個渦;當(dāng)=3/4時,三個渦繼續(xù)向前運(yùn)動,新的第三個渦接近腔體后緣;當(dāng)=時,新的第三個渦撞擊到腔體后緣,新的一個周期開始。對于這一現(xiàn)象Rossiter認(rèn)為流體流過腔體開口時的周期性渦脫落是邊緣音現(xiàn)象的特征,并且渦脫落和聲音輻射之間存在某些關(guān)系。聲輻射會導(dǎo)致渦流脫落,而渦向后運(yùn)動撞擊到后緣會導(dǎo)致聲輻射,從而形成一個聲渦反饋現(xiàn)象。正是基于以上假設(shè),他通過將自由剪切層簡化為離散渦模型推導(dǎo)出以后一直被廣泛采用的Rossiter半經(jīng)驗(yàn)預(yù)測公式。

由于在一個周期內(nèi)腔體開口上方共有三個主渦,則Rossiter周期性渦脫落公式中=1、2、3時,代入腔體幾何數(shù)據(jù)和自由來流速度,?。?.57、=0.25計(jì)算得渦脫頻率,與大渦模擬結(jié)果對比如表1所示。

一般認(rèn)為,相對誤差小于5%在工程可接受的范圍內(nèi),由表1可看出,預(yù)測頻率與數(shù)值模擬結(jié)果還有一定誤差,需更精確地確定值。

表1 峰值模態(tài)頻率對比

3.2 湍流脈動壓力的波數(shù)-頻率譜分析

波數(shù)-頻率譜分析法是波動分析法中的一種新方法,這種方法主要是把空間時間域信號變換到頻率波數(shù)域,進(jìn)而通過對波數(shù)頻率譜的分析得到波傳播的信息。近年來,該方法在地震工程、電子通訊等方面有廣泛的應(yīng)用。目前在振動噪聲方面的應(yīng)用也逐漸增多。

表面壓力波動的波數(shù)-頻率圖可以從時空表面壓力的二維傅里葉變換獲得。它顯示了波數(shù)-頻域中特定波分量的大小在波數(shù)和頻率上的分布,進(jìn)而得到該波動的傳播速度和傳播方向。二維傅立葉表達(dá)式為:

在實(shí)際中,由于測量的空間和時間有限,為了得到湍流脈動壓力作用下聲壓的波數(shù)頻率譜,一般通過對脈動壓力離散的時空信號進(jìn)行離散的二維傅里葉變換,表達(dá)式為:

式中:為時間結(jié)點(diǎn)數(shù);為空間脈動壓力監(jiān)測點(diǎn)數(shù);(x)和(t)均為窗函數(shù),本文采用漢寧窗;為監(jiān)測點(diǎn)的脈動壓力,Pa。

下面通過分析湍流脈動壓力的波數(shù)頻率譜來確定值。通過前面對Rossiter公式的推導(dǎo)知道,只需確定湍流運(yùn)動速度即可確定值。

對腔體后緣中心線的湍流脈動壓力監(jiān)測陣列進(jìn)行二維傅里葉變換,得到湍流脈動壓力作用下聲壓的波數(shù)-頻率譜,如圖7所示。基于波的傳播速度不同,波數(shù)-頻率譜可將脈動壓力分為可壓縮部分和不可壓縮部分:當(dāng)波的傳播速度小于聲速時,可認(rèn)為是由于湍流渦直接作用到壁面產(chǎn)生的脈動壓力;當(dāng)波數(shù)大于等于聲速時,可認(rèn)為是由于聲傳遞時聲壓引起的。

圖7 腔體后緣壁面脈動壓力波數(shù)-頻率譜

由圖7(a)可看出,存在一條脊,該脊即為波數(shù)-頻率譜的湍流遷移脊,是由于湍流運(yùn)動過程直接作用到后緣表面形成的。另外可以看到,湍流大部分能量都集中在遷移脊上,且主要在400 Hz以下,在100 Hz以下最為顯著,這符合湍流噪聲以低頻為主的特征。對脈動壓力監(jiān)測陣列時間序列進(jìn)行FFT變換,可得到各個監(jiān)測點(diǎn)脈動壓力在頻域內(nèi)的信息;對脈動壓力監(jiān)測點(diǎn)陣列各個時刻空間序列進(jìn)行FFT變換,可得到各個監(jiān)測點(diǎn)脈動壓力在波數(shù)域的信息。結(jié)合波數(shù)和頻率的關(guān)系即可得到湍流運(yùn)動時的傳播速度。

由圖7(b)可算出遷移脊的斜率Δ/Δ=5.56,根據(jù):

自由來流速度為200 km/h,即為55.6 m/s,則湍流渦運(yùn)動速度與自由來流運(yùn)動速度之比為=0.63。將代入式(6),計(jì)算得前3階峰值噪聲頻率,與數(shù)值模擬結(jié)果對比如表2所示。

表2 峰值模態(tài)頻率對比

對比表2與表1可看出,當(dāng)取0.63時,相比于取經(jīng)典值0.57代入Rossiter公式預(yù)測效果更好。因?yàn)橥牧鬟\(yùn)動速度受到湍流邊界層厚和腔體開口長度影響,即不同工況下湍流運(yùn)動速度與自由來流速度之比不同。因此可以通過分析湍流脈動壓力的波數(shù)-頻率譜來計(jì)算湍流運(yùn)動速度,確定Rossiter公式中的值,從而提高Rossiter公式的預(yù)測準(zhǔn)確性。

4 結(jié)論

本文采用大渦模擬的方法,對腔體流動現(xiàn)象進(jìn)行了數(shù)值模擬分析。首先確定了腔體流動自持續(xù)振蕩噪聲峰值頻率,并對腔體上方流動現(xiàn)象進(jìn)行了詳細(xì)對比分析,然后構(gòu)建并分析了腔體后緣脈動壓力作用下聲壓的波數(shù)-頻率譜,基于波數(shù)和頻率的關(guān)系得到了湍流運(yùn)動速度,進(jìn)一步得到了湍流運(yùn)動速度和自由來流速度之比。主要結(jié)論如下:

(1)腔體流動過程存在自持續(xù)振蕩現(xiàn)象,并會產(chǎn)生特定頻率的噪聲峰值。

(2)可以通過分析腔體后緣湍流脈動壓力的波數(shù)-頻率譜得到湍流的運(yùn)動速度,進(jìn)一步得到湍流運(yùn)動速度與自由來流速度之比,以提高Rossiter公式的預(yù)測準(zhǔn)確性。

[1]Rockweel D,Naudascher E. Review—self-sustaining oscillations of flow past cavities[J]. Journal of Fluids Engineering,1978(100):152-164.

[2]李輝,肖新標(biāo),朱旻昊,等. 高速列車車廂連接處氣動噪聲特性初探[J]. 振動與沖擊,2016,35(6):109-114.

[3]Rossiter J. E. Wind-tunnel experiments on the ?ow over rectangular cavities at subsonic and transonic speeds[R]. Aeronautical Research Council Reports and Memoranda,October,1964:1-32.

[4]EAST L F. Aerodynamically induced resonance in rectangular cavities[J]. Journal of Sound and Vibration,1966,3(3):277-287.

[5]MA R,SLABOCH P E,MORRIS S C. Fluid mechanics of the flow-excited Helmholtz resonator[J]. Journal of Fluid Mechanics,2009(623):1-26.

[6]董小杏. 基于CFD除霜風(fēng)道結(jié)構(gòu)分析及優(yōu)化設(shè)計(jì)[J]. 機(jī)械,2018,45(12):26-28,32.

[7]Comte P,Daude F,Mary I. Simulation of the reduction of unsteadiness in a passively controlled transonic cavity ?ow[J]. Fluids Struct,2008,24(8):1252-1261.

[8]楊國晶. 陷落式腔體水動力特性研究[D]. 哈爾濱:哈爾濱工程大學(xué),2009.

[9]Li W,Nonomura T,F(xiàn)uji K. On the feedback mechanism in supersonic cavity flows[J]. Physics of Fluids,2013,25(5):2118-2128.

[10]鄧玉清,張楠. 孔腔脈動壓力及其波數(shù)-頻率譜的大渦模擬研究[J]. 船舶力學(xué),2017,21(10):1199-1209.

[11]East L F. Aerodynamically induced resonance in rectangular cavities[J]. Journal of Sound and Vibration,1966,3(3):277-287.

The Frequency of Self-Sustained Oscillation of Cavity Flow Based on Wavenumber-Frequency Spectrum

NI Xiaosong,XIAO Xinbiao,LIU Haonan

( State Key Laboratory of Traction Power, Southwest Jiaotong University, Chengdu 610031, China )

In the semi-empirical Rossiter formula, which is currently generally recognized and widely used in self-sustained oscillation of cavity flow noise frequency prediction, only the recommended value of the empirical constantis given, and no specific calculation method is given. Therefore, this paper uses numerical simulation to determine the value of. First, the LES turbulence model is used to establish the cavity self-sustained oscillation noise prediction model. Furthermore, the first three-order noise peak frequency is obtained by analyzing the monitoring points inside the cavity and compared with the predicted frequency of Rossiter formula with the recommended value of(0.57). The predicted results have certain deviation. Then, by analyzing the wavenumber-frequency spectrum of the turbulent pulsating pressure at the trailing edge of the cavity, thevalue of 0.63 is obtained and substituted into the Rossiter formula. The first three-order noise peak frequency is obtained, which is in good agreement with the numerical simulation results. The results show that it is feasible to determine the empirical constantby analyzing the wavenumber-frequency spectrum of the turbulent pulsating pressure, and it can improve the accuracy of the Rossiter formula to predict the frequency.

cavity flow induced noise;peak frequency;Rossiter formula;wavenumber-frequency spectrum

TB123

A

10.3969/j.issn.1006-0316.2021.09.001

1006-0316 (2021) 09-0001-07

2021-01-21

國家自然科學(xué)基金(U1934203)

倪小松(1996-),男,湖南衡陽人,碩士研究生,主要研究方向?yàn)檐壍澜煌p振降噪,E-mail:1512015956@qq.com。

通訊作者:肖新標(biāo)(1978-),男,廣東陽春人,博士,副研究員,主要研究方向?yàn)檐壍澜煌p振降噪,E-mail:xinbiaoxiao@163.com。

猜你喜歡
波數(shù)后緣腔體
聲場波數(shù)積分截?cái)嗖〝?shù)自適應(yīng)選取方法
一種基于SOM神經(jīng)網(wǎng)絡(luò)中藥材分類識別系統(tǒng)
電子測試(2022年16期)2022-10-17 09:32:26
高鐵復(fù)雜腔體鑄造數(shù)值仿真及控制技術(shù)研究
高鐵制動系統(tǒng)復(fù)雜腔體鑄造成形數(shù)值模擬
橡膠擠出裝置
機(jī)翼后緣連續(xù)變彎度對客機(jī)氣動特性影響
柔性后緣可變形機(jī)翼氣動特性分析
TNF-α和PGP9.5在椎體后緣離斷癥軟骨終板的表達(dá)及意義
開孔金屬腔體場強(qiáng)增強(qiáng)效應(yīng)分析
重磁異常解釋的歸一化局部波數(shù)法
安图县| 峨山| 乌鲁木齐县| 庆元县| 广河县| 海阳市| 合水县| 营口市| 吉安市| 蓬安县| 文安县| 临漳县| 固安县| 隆尧县| 灌南县| 新宁县| 尉氏县| 霞浦县| 界首市| 东丰县| 文登市| 盘山县| 屏东县| 滕州市| 读书| 湖南省| 凉城县| 道孚县| 绥化市| 宁南县| 孝义市| 昌乐县| 阳山县| 南溪县| 正定县| 泗阳县| 义乌市| 海盐县| 崇礼县| 静安区| 深州市|