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

?

串列布置三圓柱渦激振動頻譜特性研究1)

2021-11-09 08:46:10涂佳黃譚瀟玲梁經(jīng)群
力學(xué)學(xué)報 2021年6期
關(guān)鍵詞:雷諾數(shù)升力振幅

涂佳黃 胡 剛 譚瀟玲 梁經(jīng)群 張 平

(湘潭大學(xué)土木工程與力學(xué)學(xué)院,湖南湘潭 411105)

引言

渦激振動現(xiàn)象廣泛存在于實(shí)際工程中,當(dāng)柱體結(jié)構(gòu)處于流場時,交替的渦脫落在結(jié)構(gòu)表面產(chǎn)生周期性的流體力,導(dǎo)致結(jié)構(gòu)產(chǎn)生振動.當(dāng)海洋立管[1-2]、橋梁結(jié)構(gòu)[3-4]、高層建筑[5-6]結(jié)構(gòu)振動頻率接近于自身固有頻率時,結(jié)構(gòu)會發(fā)生“鎖定”現(xiàn)象[7],誘發(fā)結(jié)構(gòu)產(chǎn)生較大的振動幅度,帶來災(zāi)難性后果.另一方面,強(qiáng)烈的結(jié)構(gòu)振動可以進(jìn)行能源轉(zhuǎn)換和利用[8-9].因此,有關(guān)渦激振動問題一直是研究者們所關(guān)注的熱點(diǎn)之一.

目前,學(xué)者們對單圓柱渦激振動問題進(jìn)行了相關(guān)研究,取得了一系列成果[10-13].與單圓柱渦激振動問題相比,串列多圓柱渦激振動情況更加復(fù)雜,對實(shí)際工程的參考意義更大.Papaioannou 等[14]研究了不同間距比工況下串列雙圓柱渦激振動問題,發(fā)現(xiàn)間距比較小時,上游圓柱的鎖定區(qū)間范圍較單圓柱工況明顯擴(kuò)大且最大振幅增大,超過臨界間距比后,下游圓柱對上游圓柱振幅的影響較小.Prasanth 等[15]對Re=100 和L=5.5D工況串列雙圓柱進(jìn)行了數(shù)值模擬研究,發(fā)現(xiàn)下游圓柱的鎖定區(qū)間受到了上游圓柱的影響,明顯大于單圓柱渦激振動的現(xiàn)象.及春寧等[16]對低雷諾數(shù)下串列雙圓柱渦激振動中下游圓柱大振幅響應(yīng)的耦合機(jī)制進(jìn)行了全面分析,發(fā)現(xiàn)上游圓柱脫落漩渦激勵下游圓柱大振幅響應(yīng)的機(jī)理來自于脫落漩渦所誘發(fā)的位于下游圓柱的運(yùn)動方向上低壓區(qū).Mysa 等[17]研究了L=(1.0~4.0)D,Re=100下串列雙圓柱渦激振動中對響應(yīng)起關(guān)鍵作用的因素,發(fā)現(xiàn)上游圓柱尾流與下游圓柱之間的耦合作用的強(qiáng)弱取決于橫流向流體力與圓柱位移之間的相位差,與下游圓柱運(yùn)動同相的部分對其響應(yīng)起到了促進(jìn)作用.Mittal 等[18]通過數(shù)值模擬發(fā)現(xiàn)Re=100,L=5.5D時上游圓柱的振動響應(yīng)趨近于單圓柱渦激振動,而下游圓柱經(jīng)歷了大振幅振動.

關(guān)于串列三圓柱渦激振動問題,圓柱體之間相互作用對振動影響很大,潛在的機(jī)理也值得更多的努力去研究.Igarashi 和Suzuki[19]通過試驗(yàn)研究了L=(1.0~4.0)D,Re=10 900~39 200 的串列三圓柱繞流特性,根據(jù)下游圓柱上從上游圓柱分離的具有動態(tài)效應(yīng)的剪切層,將6 種尾流模態(tài)進(jìn)行了分類,并得到了發(fā)生雙穩(wěn)態(tài)現(xiàn)象的區(qū)域.Yu 等[20]通過研究發(fā)現(xiàn)三圓柱系統(tǒng)的振動響應(yīng)劇烈,橫流向最大振幅提高了25%.同時,結(jié)構(gòu)運(yùn)動軌跡大多呈現(xiàn)為不規(guī)則形狀.Chen 等[21]對Re=100,L=(1.2~5.0)D工況下串列三圓柱進(jìn)行數(shù)值模擬研究,發(fā)現(xiàn)L=1.2D時,串列三圓柱均表現(xiàn)為馳振現(xiàn)象.指出馳振現(xiàn)象的3 個關(guān)鍵因素是平衡位置偏移、低頻振動以及漩渦脫落與圓柱之間的時機(jī).Mahmoud 和Atef[22]對串列布置三圓柱系統(tǒng)的流致振動問題進(jìn)行了研究,著重分析了上中游兩靜止圓柱體間距比的變化對下游圓柱體動力響應(yīng)及其流場特性的影響.張志猛等[23]對Re=100 時上游圓柱固定串列三圓柱渦激振動進(jìn)行數(shù)值模擬研究,發(fā)現(xiàn)當(dāng)振幅較小時,上游圓柱的剪切層將三圓柱包裹在一起,尾流表現(xiàn)為單鈍體脫渦模式.當(dāng)振幅較大時,上游圓柱的旋渦重附著于下游圓柱上,尾流呈現(xiàn)出兩列并排的漩渦.Behara 等[24]對串列布置三圓柱進(jìn)行數(shù)值模擬分析,對于下游兩圓柱根據(jù)其振幅和頻率的響應(yīng),將約化速度劃分為3 個區(qū)間:初始激勵區(qū)、上鎖定區(qū)、下鎖定區(qū).涂佳黃等[25]對Re=100下剪切來流作用下串列三圓柱體雙自由度流致振動問題進(jìn)行了數(shù)值計(jì)算研究,隨固有頻率比的增大,上游圓柱順流向鎖定區(qū)間范圍會減小,而中下游圓柱雙向鎖定區(qū)間會擴(kuò)大.

綜上所述,頻譜特性分析在多柱體結(jié)構(gòu)渦激振動問題中的研究討論較少,本文基于四步半隱式特征線分裂算子有限元方法,對均勻來流作用下串列布置三圓柱系統(tǒng)雙自由度渦激振動問題進(jìn)行了數(shù)值模擬,重點(diǎn)分析了圓柱體功率譜密度隨不同雷諾數(shù)、固有頻率比、約化速度的變化規(guī)律,對不同工況下功率譜密度的形態(tài)進(jìn)行對比,深入研究了功率譜密度對結(jié)構(gòu)動力響應(yīng)的影響并揭示了其內(nèi)在機(jī)理.為實(shí)際工程應(yīng)用提供參考依據(jù).

1 流固耦合數(shù)值方法

1.1 流體控制方程

基于ALE 方法下的不可壓縮黏性流體的N-S 方程的無量綱形式表達(dá)如下

式中,xi為空間坐標(biāo);ui為流體速度;cj=uj-wj,cj為流體相對于網(wǎng)格移動速度的對流速度,wj為網(wǎng)格移動速度;p為流體壓力;t為時間;Re=U∞D(zhuǎn)/ν,Re為雷諾數(shù),D與U∞分別為特征長度尺寸與特征速度,ν為動力黏性系數(shù).

1.2 固體控制方程

彈性支撐結(jié)果運(yùn)動體系可假設(shè)為之質(zhì)量-阻尼-彈簧系統(tǒng),考慮雙自由度運(yùn)動的控制方程量綱歸一化的形式如下

式中,X與Y分別為結(jié)構(gòu)順流向和橫流向無量綱位移;ξ 為結(jié)構(gòu)阻尼系數(shù),為了得到最大結(jié)構(gòu)位移響應(yīng),取ξ=0;Ur,x=U∞/(fn,xD) 和Ur,y=U∞/(fn,yD)和分別為x方向(順流向) 與y方向(橫流向) 約化速度,其中fn,x和fn,y分別為彈性支撐圓柱體結(jié)構(gòu)的順流向和橫流向自然頻率;和分別為阻力系數(shù)和升力系數(shù);mr=m/(ρD2)為結(jié)構(gòu)折合質(zhì)量,其中m為單位長度結(jié)構(gòu)質(zhì)量,ρ 為流體密度.本文采用Newmark-β 時間積分法求解結(jié)構(gòu)動力方程.

1.3 網(wǎng)格更新

為了避免網(wǎng)格失效問題,本文采用改進(jìn)Laplace方程的邊值問題方法對網(wǎng)格坐標(biāo)進(jìn)行更新[26]

式中,Si為網(wǎng)格節(jié)點(diǎn)i方向位移;Γm和Γf分別為網(wǎng)格動邊界和固定邊界;gi是運(yùn)動邊界上的節(jié)點(diǎn)位移;τ 是網(wǎng)格形變控制參數(shù),表達(dá)式如下

式中,Δe為計(jì)算網(wǎng)格單元的面積(或體積);Δmin和Δmax分別為網(wǎng)格單元中最小與最大的面積(或體積).本文采用Galerkin 有限元方法求解該Laplace 方程的邊值問題.

1.4 計(jì)算流程

本文采用分區(qū)迭代方法求解鈍體結(jié)構(gòu)渦激振動問題,該算法的計(jì)算流程如下:

(1) 求解流體控制方程.運(yùn)用基于四步半隱式CBS 穩(wěn)定化有限元方法求解流體控制方程式(1) 和式(2),獲得t(n+1)時刻流場速度與壓力,從而得到流體作用于結(jié)構(gòu)上的流體力.

(2) 求解固體控制方程.將流體力施加到結(jié)構(gòu),以Newmark-β 時間積分方法求解結(jié)構(gòu)運(yùn)動控制方程式(3),得到t(n+1)時刻結(jié)構(gòu)的動力響應(yīng).

(3) 網(wǎng)格更新.采用Galerkin 有限元方法求解Laplace 方程,獲得網(wǎng)格節(jié)點(diǎn)位移及網(wǎng)格速度,并更新網(wǎng)格坐標(biāo).

(4)返回第一步開始計(jì)算t(n+2)時刻,并循環(huán)至系統(tǒng)達(dá)到穩(wěn)定為止.

本文的數(shù)值方法已運(yùn)用于渦激振動相關(guān)問題求解過程中[25,27-28],并能得到較好的數(shù)值結(jié)果,驗(yàn)證其正確性與可靠性.

2 問題描述

2.1 計(jì)算模型

為了研究上游結(jié)構(gòu)尾流充分發(fā)展后對下游結(jié)構(gòu)振動響應(yīng)的影響,選取圓柱體結(jié)構(gòu)間距比Lx=5.5D.其他參數(shù)為:雷諾數(shù)Re=80,100,160,固有頻率比r=fn,x/fn,y=1.0,1.5,2.0,約化速度Ur=Ur,x=3~21,質(zhì)量比mr=2.0.計(jì)算域尺寸為[-30D,60D]×[-20D,20D],上游圓柱(upstream cylinder,UC)、中游圓柱(midstream cylinder,MC) 和下游圓柱(downstream cylinder,DC)的圓心位置分別為(-5.5D,0)、(0,0) 與(5.5D,0),本文計(jì)算域堵塞率B=5%,如圖1 所示.邊界條件設(shè)置如下,入口處采用速度入口:ux=U∞=1.0,uy=0;出口處采用壓力出口:p=0;上下邊界均為自由滑移邊界:?ux/?y=0,uy=0;圓柱表面均采用無滑移邊界:ux=uy=0.為了獲得較大的動力響應(yīng),不考慮結(jié)構(gòu)阻尼的影響.對于彈性支撐圓柱體結(jié)構(gòu),可簡化為質(zhì)量-彈簧系統(tǒng)模型,圓柱表面速度為

圖1 計(jì)算模型和邊界條件Fig.1 Computational model and boundary conditions

2.2 網(wǎng)格劃分

流場計(jì)算域采用非結(jié)構(gòu)化三角形網(wǎng)格進(jìn)行劃分,其中圓柱附近及尾流區(qū)域進(jìn)行網(wǎng)格加密處理.表1為不同網(wǎng)格參數(shù)下串列三圓柱結(jié)構(gòu)渦激振動動力響應(yīng)特性的統(tǒng)計(jì)結(jié)果對比.由表1 可知:與粗網(wǎng)格GI相比,網(wǎng)格GII 的Xmax/D與Ymax/D的計(jì)算結(jié)果最大變化率分別為5.26%和4.00%;與密網(wǎng)格GIII 相比,其最大變化率分別下降為1.25%和0.85%.網(wǎng)格GII已滿足數(shù)值模擬結(jié)果網(wǎng)格收斂性要求和計(jì)算精度要求.綜上所述,本文所有算例采用的網(wǎng)格密度和分布情況與GII 類似,無量綱計(jì)算時間步長Δt=0.002.

表1 網(wǎng)格獨(dú)立性驗(yàn)證:串列布置三圓柱體流致振動計(jì)算結(jié)果(Re=100,r=1.0,Ur=6.0)Table 1 Grid independence test:The results for the three tandem circular cylinders at Re=100,r=1.0 and Ur=6.0

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

本文分析了不同雷諾數(shù)、固有頻率比與約化速度工況下串列三圓柱體結(jié)構(gòu)體系雙自由度流致振動振幅特性、頻譜特性、流體力系數(shù)的變化情況.本文根據(jù)UC 的振動特性對約化速度進(jìn)行如下劃分:3 ≤Ur≤5(區(qū)間A);5 <Ur≤9(區(qū)間B);9 <Ur≤15(區(qū)間C1),Ur>15(區(qū)間C2).

3.1 振幅特性

由圖2 可知,在區(qū)間A 內(nèi),UC 在x和y兩個方向的振幅隨約化速度的增大而增大.在大振幅區(qū)間內(nèi),頻率比對UC 振幅的影響較大,具體表現(xiàn)為:當(dāng)r≤1.5 時,頻率比的變化對UC 的振動響應(yīng)的影響較小.各工況下橫流向振幅十分接近,在Ur=5 工況下達(dá)到最大振幅值后逐漸下降.頻率比較小時結(jié)構(gòu)剛度較大,圓柱順流向振動響應(yīng)微弱,x方向的振幅趨于0.當(dāng)r=2.0,Ur=5 時,UC 橫流向振動頻率fs,y=0.174(Re=80,100)和0.178(Re=160)接近其固有頻率,同時順流向振動頻率與固有頻率的比值在1 附近.此時UC 發(fā)生雙頻共振,其振動響應(yīng)顯著增強(qiáng),順流向振幅遠(yuǎn)大于其他頻率比工況下的振幅,橫流向振幅也達(dá)到較大值.值得注意地是,當(dāng)Re=80時,在雙頻共振作用下的尾渦模式為2S 模式,漩渦強(qiáng)度逐漸減弱并在下游不遠(yuǎn)處會消失,如圖3(a)所示;Re=100 時的尾渦模式也呈2S 模式,但下游的漩渦重新卷起并脫落形成新的渦街;不同地是,Re=160時尾渦模式呈現(xiàn)P+S 模式,在較遠(yuǎn)的下游能保持較高的強(qiáng)度.雷諾數(shù)對UC 的振動響應(yīng)的影響較小,各工況下UC 的振幅變化規(guī)律類似.進(jìn)入C1區(qū)間后,UC 在x和y兩個方向的振幅逐漸趨于穩(wěn)定,x方向的振動響應(yīng)極其微弱.

圖2 不同雷諾數(shù)與頻率比工況下,UC 在x 和y 兩個方向最大振幅隨約化速度的變化Fig.2 Variation of the maximum vibrating amplitudes of upstream cylinder with reduced velocity under different Reynolds number and frequency ratios

圖3 Ur=5,r=2.0 時,流場瞬時渦量圖Fig.3 Instantaneous vorticity diagram of flow field at Ur=5 and r=2.0

與UC 相比,雷諾數(shù)和頻率比對MC 的振幅的影響更加顯著.在區(qū)間A 內(nèi),MC 受到“弱鎖定”作用,如圖4 所示,Ur=4 時,x和y兩個方向的振幅取得第一個極大值,當(dāng)r=1.5 時順流向振幅接近MC 的順流向振幅最大值.特別地,當(dāng)r=2.0 時,MC 在x和y兩個方向的運(yùn)動在區(qū)間A 內(nèi)的振幅很小,這種“屏蔽”現(xiàn)象可以用流場進(jìn)行解釋.如圖5 所示,當(dāng)Ur=4時,UC 在靠近MC 的方向產(chǎn)生漩渦,而漩渦脫落的位置與UC 的距離非常近,連續(xù)渦結(jié)構(gòu)的垂直距離較大,說明尾流發(fā)展過程中沒有渦撞擊的發(fā)生,漩渦直接從MC 和DC 的兩側(cè)通過,在UC 的下游不遠(yuǎn)處形成穩(wěn)定的兩列渦街.因此,MC 的漩渦脫落被顯著抑制,進(jìn)一步導(dǎo)致橫向動力響應(yīng)降低.在Ur=5 處,UC在P+S 尾渦模式中表現(xiàn)出較強(qiáng)的振動響應(yīng),渦結(jié)構(gòu)的垂直距離被拉大,從而對振動產(chǎn)生實(shí)質(zhì)性的抑制.Bao 等[29]也發(fā)現(xiàn)了類似的現(xiàn)象.在區(qū)間B 內(nèi),各工況下MC 的橫流向運(yùn)動在Ur=7 附近發(fā)生頻率鎖定現(xiàn)象,此時y方向的振幅最大,如圖4 所示,且鎖定區(qū)間的范圍會隨著雷諾數(shù)的增大而擴(kuò)大.頻率比對y方向的振幅影響較小,各工況下的振幅比較接近.特別的,當(dāng)Re=160 時,在Ur=6 處的橫流向振幅,r=2.0工況Ymax/D=0.322 遠(yuǎn)小于1.073 (r=1.0) 和1.025(r=1.5).頻率比和雷諾數(shù)對x方向的振幅影響較大,Re<160 時,x方向的振幅隨頻率比的增大而增大,特別的是,當(dāng)Re=100,Ur=9 時,x方向的振幅與頻率比成反比關(guān)系.當(dāng)Re=160 時,r=1.0 工況下MC 在x方向上的振幅是其他頻率比工況下振幅的2~3 倍,但在Ur=7 時Xmax/D=0.073 與其他工況下的振幅接近.在區(qū)間C 內(nèi)(包括C1和C2),MC 逐漸退出鎖定區(qū)間,圓柱振動響應(yīng)逐漸減弱,因此MC在x和y方向的振幅逐漸減小并趨于穩(wěn)定.值得注意的是,當(dāng)r=1.0 時,在C1區(qū)間內(nèi)各雷諾數(shù)工況下MC 在x方向上的振動仍會保持相當(dāng)大的振幅,達(dá)到該工況下順流向振幅的最大值后逐漸下降趨于平穩(wěn).總的來說,MC 的動力響應(yīng)受UC 尾流振蕩的影響較大,相反,頻率比的影響較小.

圖4 不同雷諾數(shù)與頻率比工況下,MC 在x 和y 兩個方向最大振幅隨約化速度的變化Fig.4 Variation of the maximum vibrating amplitudes of midstream cylinder with reduced velocity under different Reynolds numbers and frequency ratios

圖5 Re=160,r=2.0 時,流場瞬時渦量圖Fig.5 Instantaneous vorticity diagram of flow field at Re=160 and r=2.0

DC 振幅的變化趨勢與MC 類似,但DC 的“鎖定”區(qū)間的范圍比MC 更寬,且會隨著雷諾數(shù)的增大而明顯擴(kuò)大.在區(qū)間A 內(nèi)DC 在x和y兩個方向的振幅較小,Ur=4 時出現(xiàn)與MC 類似的“弱鎖定”現(xiàn)象,此時振幅出現(xiàn)小幅度的“上升-下降”過程,如圖6 所示.特別的是,當(dāng)r=2.0 工況DC 在區(qū)間A 內(nèi)的振幅波動很小.由于UC 和MC 尾流的共同“屏蔽”作用,當(dāng)Ur≤6 時,各工況下DC 的兩個方向的振幅均小于UC 和MC 的振幅;但當(dāng)Ur>6 時,DC 在兩個方向的振幅均大于UC 和MC 的振幅.由此可見,本文的臨界約化速度在6~7 之間,且受雷諾數(shù)與頻率比的影響較小[16,30].在區(qū)間B 內(nèi),DC 在y方向上突然劇烈振動,隨Ur增大,橫流向振幅取得最大值后下降,但仍保持較大的振幅值.但Re=160 時,在整個B 區(qū)間內(nèi)橫流向振幅均不斷增大.值得注意的是,與其他頻率比工況相比,r=2.0 時DC 進(jìn)入大振幅區(qū)間的約化速度Ur=7.在區(qū)間C1內(nèi),隨著Ur的增大,不同雷諾數(shù)和頻率比工況下DC 的橫流向振幅在取得最大值后逐漸減小.當(dāng)Re=160 時,r=1.5 與其他工況完全不同,此時橫流向振幅在整個區(qū)間C1內(nèi)均保持較大的數(shù)值.值得注意地是,當(dāng)r=2.0 時,橫流向振幅在Ur=14 處突然減小,Ymax/D=0.849(r=2.0)約為r=1.5 工況的53%(Ymax/D=1.583).在區(qū)間C 內(nèi)DC順流向振幅受頻率比的影響比橫流向振幅更敏感,具體地說:當(dāng)r=1.0 時,隨著Ur的增大,DC 的順流向運(yùn)動在Ur=11 附近達(dá)到該工況下的最大振幅值后逐漸下降,此時DC 在C1區(qū)間內(nèi)的順流向振幅遠(yuǎn)大于其他頻率比工況下的振幅.當(dāng)r>1 時,不同雷諾數(shù)下DC 順流向運(yùn)動在區(qū)間C1內(nèi)較弱,振幅較小且呈現(xiàn)逐漸下降的趨勢.特殊的是,當(dāng)Re=160,r=1.5時,DC 順流向振幅在Ur>12 后出現(xiàn)明顯的躍升,并達(dá)到該工況下順流向振幅最大值Xmax/D=1.12.C2區(qū)間內(nèi)DC 在兩個方向的振幅進(jìn)一步減小,逐漸趨于穩(wěn)定.值得注意地是,當(dāng)Re<160,r=1.0 時,DC 的橫流向振幅在C2區(qū)間內(nèi)出現(xiàn)增大的趨勢.同時,MC也具有相同的變化規(guī)律.

圖6 不同雷諾數(shù)與頻率比工況下,DC 在x 和y 兩個方向最大振幅隨約化速度的變化Fig.6 Variation of the maximum vibrating amplitudes of downstream cylinder with reduced velocity under different Reynolds numbers and frequency ratios

3.2 流體力系數(shù)分析

Lx=5.5D工況下,MC 對UC 的影響極其微小,UC 的尾流能夠得到充分發(fā)展,導(dǎo)致其振動響應(yīng)類似于單圓柱工況[18].另一方面,受到UC 尾流和自身渦脫落的共同作用,使得MC 和DC 振動響應(yīng)發(fā)生顯著變化[15].

由圖7(a) 可知,隨約化速度的增大,UC 在各工況下阻力系數(shù)的平均值(CD-mean)逐漸增大,頻率比對CD-mean 的影響很小,雷諾數(shù)越大CD-mean 越大.特別的,在r=2.0 工況下,當(dāng)Re=160,Ur=4 時CD-mean=2.405 遠(yuǎn)大于其他雷諾數(shù)工況,此時UC 順流向振幅Xmax/D=0.192 3 分別是0.015 8 (Re=80)和0.022 9(Re=100)的12 倍和8 倍.UC 的CD-mean在Ur=5 處達(dá)到最大值隨后逐漸降低并最終保持穩(wěn)定,CD-mean 隨頻率比的增大而增大,但與雷諾數(shù)成反比.區(qū)間A 內(nèi)MC 受UC 尾流的屏蔽作用,MC的CD-mean 幅值較小且小幅下降,此時頻率比小的工況下CD-mean 反而較大.隨著MC 進(jìn)入大振幅區(qū)間,CD-mean 顯著增大并在Ur=7 時達(dá)到最大值并逐漸減小趨于穩(wěn)定,此時CD-mean 受頻率比的影響很小.值得注意地是,MC 的CD-mean 穩(wěn)定時的幅值大約是UC 的一半.DC 的CD-mean 隨約化速度的變化十分復(fù)雜,同時雷諾數(shù)與頻率比對CD-mean 的影響也較為顯著.CD-mean 在區(qū)間A 和區(qū)間B 內(nèi)的波動較大,分別在Ur=4 和Ur=8 時達(dá)到兩個極值.在區(qū)間C 內(nèi),雷諾數(shù)越大DC 的CD-mean 越大,隨Ur增大,CD-mean 逐漸減小并趨于穩(wěn)定.

圖7 不同雷諾數(shù)工況下,圓柱阻力系數(shù)平均值的變化Fig.7 Variation of mean value of drag coefficient under different Reynolds numbers

UC 在各工況下升力系數(shù)的均方根值(CL-rms)隨約化速度的變化趨勢與CD-mean 變化類似,從圖8(a)可以看出,各工況下UC 的CL-rms 大致呈現(xiàn)先增大后降低直到平穩(wěn)的變化趨勢,并在Ur=4 處達(dá)到最大值.隨雷諾數(shù)的增大,UC 的CL-rms 變大,此時在y方向的振動響應(yīng)變化劇烈程度變強(qiáng).但當(dāng)r≤1.5,Ur=5 時,雷諾數(shù)與UC 的CL-rms 卻成反比.MC 的CL-rms 在Ur=3 處就達(dá)到最大值,隨著Ur增大,CL-rms 迅速減小,頻率比和雷諾數(shù)對CL-rms的影響很小.在區(qū)間B 內(nèi),MC 的CL-rms 曲線變化的波動性明顯增強(qiáng).當(dāng)Re=160 時,不同頻率比工況下CL-rms 曲線出現(xiàn)多個突變點(diǎn).值得注意地是,當(dāng)Ur=6 時,MC 在r=2.0 工況下的CL-rms=0.081 遠(yuǎn)小于0.565(r=1.0)和0.577(r=1.5),導(dǎo)致此MC 的在y方向的振幅遠(yuǎn)小于其他頻率比工況下的振幅.在區(qū)間C 內(nèi),CL-rms 變化很小呈平穩(wěn)趨勢.“弱鎖定”作用使DC 的CL-rms 在Ur=4 處突然增大,隨后逐漸減小.當(dāng)5 <Ur<15 時,雷諾數(shù)和頻率比對DC 的CL-rms 的影響很強(qiáng),DC 的CL-rms 曲線出現(xiàn)很強(qiáng)的波動.特別的,在Re=160 工況下,當(dāng)r=1.5,Ur=15時,MC 脫落的漩渦會與DC 上、下側(cè)的剪切層融合形成更強(qiáng)的漩渦,加劇DC 的振動響應(yīng),導(dǎo)致其升力系數(shù)的波動性會增強(qiáng),如圖9 所示.i 時刻,DC 處于上方最大位移處,其尾流上方的漩渦基本已脫離,同時MC 尾流已傳播至其下方,使得DC 下表面受到較大吸力.在ii 時刻,DC 的尾流區(qū)會出現(xiàn)一個大漩渦,同時迎流面會受到MC 尾流的沖擊作用.隨后,DC 運(yùn)動至下方最大位移處,其尾流上方的漩渦與MC 尾流已融合成一個整體,使得DC 上表面受到較大吸力.

圖8 不同雷諾數(shù)工況下,圓柱升力系數(shù)均方根值的變化Fig.8 Variation of root mean square value of lift coefficient under different Reynolds numbers

圖8 不同雷諾數(shù)工況下,圓柱升力系數(shù)均方根值的變化(續(xù))Fig.8 Variation of root mean square value of lift coefficient under different Reynolds numbers(continued)

圖9 Re=160,r=1.5,Ur=15 時,DC 運(yùn)動位移及升力時程曲線圖及對應(yīng)不同時刻渦量場云圖Fig.9 Time-history curve of displacement and lift coefficient of downstream cylinders and instantaneous pressure contours at different times when Re=160,r=1.5 and Ur=15

3.3 流體力功率譜密度分析

通過對圓柱體流體力系數(shù)功率譜密度(PSD)對比分析,發(fā)現(xiàn)圓柱流體力系數(shù)PSD 曲線在不同雷諾數(shù)工況下呈現(xiàn)不同的主峰幅值、頻譜成分及波動性,導(dǎo)致圓柱動力響應(yīng)特性發(fā)生較大變化.另外,研究發(fā)現(xiàn)r與Ur的變化對PSD 曲線有顯著影響,導(dǎo)致圓柱動力學(xué)響應(yīng)呈現(xiàn)明顯的差異性.本節(jié)針對r=1.0 和2.0 工況下三圓柱體流體力系數(shù)PSD 曲線的變化及其對圓柱的振動響應(yīng)的影響進(jìn)行詳細(xì)分析.

在區(qū)間A 內(nèi),當(dāng)r=1.0 且Ur=3 時,各工況下UC 的流體力系數(shù)PSD 曲線出現(xiàn)兩階頻率峰值且曲線均比較光滑,隨著約化速度和雷諾數(shù)的增大,PSD曲線的波動性增強(qiáng),如圖10(a) 所示,在一階頻率聚集了更多的能量,增強(qiáng)了UC 的振動響應(yīng).Ur增大到5 時,各雷諾數(shù)工況下升力系數(shù)PSD 曲線波動性減弱,曲線包含的頻率峰值增大至三階,有效增大了升力系數(shù)PSD 曲線包圍的面積,使得UC 獲得更多動能,在y方向上產(chǎn)生較大的振幅,具體分析見第3.4節(jié).在區(qū)間B 內(nèi),隨Ur增大,升力系數(shù)PSD 曲線波動性增強(qiáng),但阻力系數(shù)PSD 曲線受約化速度的影響較小.各工況下升力系數(shù)PSD 曲線包含了主峰幅值與次峰幅值同等量級的多階頻率峰值,此時能量在次峰的分布比Ur=5 工況更加豐富,而主頻能量值更低,流體流向結(jié)構(gòu)的能量值減少,導(dǎo)致UC 的振動響應(yīng)比Ur=5 時有所減弱.隨Ur進(jìn)一步增大,升、阻力系數(shù)PSD 曲線趨于光滑,頻率峰值逐漸減少至二階,主峰幅值逐漸穩(wěn)定,次峰幅值逐漸減小直至可以忽略,UC 在區(qū)間C 內(nèi)的振動響應(yīng)逐漸減弱,在兩個方向的振幅趨于穩(wěn)定,最終PSD 曲線逐漸形成光滑的雙峰譜形態(tài).

圖10 不同工況下,上游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)Fig.10 PSD of lift coefficient(red)and drag coefficient(black)of upstream cylinder under different case

當(dāng)r增大至2.0 時,在區(qū)間A 內(nèi),如圖10(e) 和圖10(f)所示,各工況下流體力系數(shù)PSD 曲線的波動性很弱,隨著約化速度和雷諾數(shù)的增大,PSD 曲線的頻率峰值成分增大,UC 的運(yùn)動變得復(fù)雜,其他特征與r=1.0 時類似.在區(qū)間B 內(nèi),各工況下圓柱流體力PSD 曲線均只出現(xiàn)三階頻率峰值,在Ur>7 后,PSD曲線的波動性逐漸減弱,主導(dǎo)頻率的能量值越來越低,UC 的動力學(xué)響應(yīng)也越來越弱,即圓柱振動響應(yīng)的強(qiáng)弱與主頻能量值成正比.需要說明的是,Ur=8時,雷諾數(shù)對升力系數(shù)PSD 曲線有顯著的影響.當(dāng)Re=160 時,升力系數(shù)PSD 曲線波動性與其他工況相比明顯增強(qiáng),次峰及雜頻占據(jù)的能量值增多,此時各工況下升力系數(shù)PSD 曲線的主頻能量值幾乎相同.升力系數(shù)PSD 曲線的波動性對UC 運(yùn)動的產(chǎn)生較大影響,Re=160 時UC 橫流向振幅Ymax/D=0.131 相較于Re=100 時的Ymax/D=0.366 下降了64%,同時比Ur=7 時的Ymax/D=0.705 下降了81%.研究發(fā)現(xiàn),當(dāng)r=1.0 時同樣具有類似規(guī)律.進(jìn)入?yún)^(qū)間C后,PSD 曲線隨約化速度的增大而變得光滑,雷諾數(shù)越小曲線的波動也越來越弱,最終形成光滑的兩階頻率峰值形態(tài).

與UC 相比,MC 由于受到上游圓柱的尾流影響,流體力系數(shù)PSD 曲線隨雷諾數(shù)和約化速度的變化更加復(fù)雜.當(dāng)r=1.0 時,在區(qū)間A 內(nèi),不同工況下的升、阻力系數(shù)PSD 曲線均比較光滑,分別包含三階和兩階頻率峰值,但主頻能量值較小,圓柱振動響微弱.特別地,當(dāng)Ur=4 時,PSD 曲線的波動性隨雷諾數(shù)的增大而增強(qiáng),如圖11(a)~圖11(c) 所示,PSD 曲線出現(xiàn)多階同等量級的頻率峰值,表明此時的泄渦頻率不位于鎖定區(qū)間.在區(qū)間B 內(nèi),當(dāng)Re<160,升力系數(shù)PSD 曲線的頻率峰值數(shù)量隨Ur的增大而增大,圓柱運(yùn)動趨于復(fù)雜.Re=160 時,升力系數(shù)PSD 曲線波動性增強(qiáng),并包含多階頻率峰值,主頻維持較高能量值,故此時MC 在y方向的運(yùn)動能以較大的振幅持續(xù)到區(qū)間B 的末端.隨著Ur增大,阻力系數(shù)PSD 曲線變化較小,僅表現(xiàn)為波動性的差異,所以在區(qū)間B內(nèi)順流向振動相對平穩(wěn),振幅較區(qū)間A 有小幅度增大.特別地,當(dāng)Re=160 時,如圖11(c)所示,Ur=6 時阻力系數(shù)PSD 曲線主、次頻峰值的能量值接近,且能量在各雜頻的分布幅值也較多,圓柱順流向振動響應(yīng)較強(qiáng).隨著Ur增大,阻力系數(shù)PSD 曲線的頻率幅值減少到三階(Ur=7),主頻能量值與Ur=6 工況相比差了一個數(shù)量級,順流向振幅Xmax/D=0.073 與Ur=6 工況比降低了三倍.流體力系數(shù)的PSD 曲線在區(qū)間C 后波動逐漸減弱,并逐漸形成光滑的三階頻率峰值曲線.

圖11 不同工況下,中游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)Fig.11 PSD of lift coefficient(red)and drag coefficient(black)of midstream cylinder under different cases

圖11 不同工況下,中游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)(續(xù))Fig.11 PSD of lift coefficient(red)and drag coefficient(black)of midstream cylinder under different cases(continued)

當(dāng)r=2.0 時,在區(qū)間A 內(nèi)PSD 曲線的波動性隨雷諾數(shù)的增大而減弱,如圖11(f)所示,Ur=4 時MC的“弱鎖定”作用逐漸消失,流體經(jīng)過MC 時傳遞給圓柱的能量減少,導(dǎo)致圓柱的動力學(xué)響應(yīng)十分微弱.當(dāng)Ur=5 時,雷諾數(shù)對升、阻力系數(shù)PSD 的影響十分顯著.當(dāng)Re=100 時,阻力系數(shù)PSD 的主峰值頻率fscd=0.352 是升力系數(shù)PSD 的兩倍,從UC 脫落的漩渦直接從MC 上、下側(cè)通過,如圖12 所示,導(dǎo)致MC 的運(yùn)動軌跡呈現(xiàn)為對稱的“8”字型,與Prasanth等[31]研究的結(jié)論一致.但當(dāng)Re=160 時,升、阻力系數(shù)PSD 曲線的各階峰值頻率完全重疊,UC 在一個周期內(nèi)脫落的孤立漩渦A 和一對反向的漩渦B+C,分別從MC 的上、下側(cè)向下游發(fā)展,由于上下側(cè)渦量不對稱,導(dǎo)致MC 的運(yùn)動軌跡由“8”字形變成不對稱形狀.值得注意地是,上述工況下MC 受到UC 的“屏蔽”作用較大,且此時升、阻力曲線與位移的時程曲線的幾乎反向,導(dǎo)致MC 在Ur=5 時的振動響應(yīng)很弱.在區(qū)間B 內(nèi),MC 進(jìn)入鎖定區(qū)間,當(dāng)Ur接近8 時,PSD 曲線中會出現(xiàn)高階的頻率峰值,曲線波動性也大大減弱.此時MC 被完全鎖定,更多的能量由流體傳遞給圓柱,MC 獲得更多的動能而劇烈振動,x和y方向的振幅達(dá)到最大值.進(jìn)入?yún)^(qū)間C 后,Re<160 時升力系數(shù)的PSD 曲線的變化規(guī)律與r=1.0 時大體相同.但當(dāng)Re=160 時,在C1區(qū)間內(nèi),升力系數(shù)PSD曲線中出現(xiàn)了多階頻率峰值,包含的總能量值高,使得MC 產(chǎn)生較大的橫流向振幅.在C2區(qū)間中,升、阻力系數(shù)PSD 曲線波動性逐漸減弱,逐漸形成光滑的三階頻率峰值曲線.

圖12 Re=100 與160 時中游圓柱流體力系數(shù)與位移時程曲線、運(yùn)動軌跡及流場圖(r=2.0,Ur=5)Fig.12 Time history of the fluid force coefficient and displacement,the trajectory and the flow field of midstream cylinder at Re=100 and 160(r=2.0,Ur=5)

DC 受到上中游兩圓柱的尾流共同作用,雷諾數(shù)、頻率比和約化速度對升、阻力系數(shù)PSD 的影響更加顯著.當(dāng)r=1.0 時,與MC 類似,PSD 曲線在Ur=4 時出現(xiàn)數(shù)量眾多的頻率峰值,主頻能量值較高,PSD 曲線的波動性隨雷諾數(shù)的增大而增強(qiáng),如圖13(a)~圖13(c)所示,導(dǎo)致DC 在A 區(qū)間內(nèi)的運(yùn)動在Ur=4 時最強(qiáng)烈.在區(qū)間B 內(nèi),隨Ur增大PSD 曲線中出現(xiàn)了更多數(shù)量的的頻率峰值,同時曲線的波動性也明顯增強(qiáng),導(dǎo)致DC 的運(yùn)動越趨復(fù)雜.特殊地,當(dāng)Re=100,Ur=7 時,升力系數(shù)PSD 曲線的頻率峰值比Ur=6 工況時更多,主頻能量值卻少了兩個數(shù)量級,致使DC 的橫流向振幅減小.當(dāng)Ur=8 時,雷諾數(shù)對PSD 曲線影響十分顯著,當(dāng)Re=100 時,如圖13(a)~圖13(c)所示,升力系數(shù)PSD 曲線光滑且中出現(xiàn)四階頻率峰值.當(dāng)Re=160 時,升力系數(shù)PSD 曲線由繁多的峰值構(gòu)成整體呈下降趨勢的曲線.研究發(fā)現(xiàn),上述的3 種不同的曲線形態(tài)是串列布置三圓柱體的PSD 曲線的最基本形態(tài).C 區(qū)間內(nèi),升力系數(shù)PSD曲線中出現(xiàn)多階頻率峰值,包含較高的能量值,導(dǎo)致下游圓柱產(chǎn)生強(qiáng)烈的動力學(xué)響應(yīng).當(dāng)Ur>13,PSD 曲線的波動性逐漸增強(qiáng),雜頻占據(jù)的能量值增大,圓柱體獲取的能量值減少,導(dǎo)致下游圓柱的橫流向振幅逐漸減小.與UC 和MC 不同,當(dāng)Re=80 時,在Ur>15后,升力系數(shù)PSD 曲線中出現(xiàn)三階頻率峰值,主頻能量值隨Ur增大而變大,使得DC 橫流向振幅不斷增大.隨Ur增大,阻力系數(shù)PSD 曲線的主頻能量值逐漸降低,導(dǎo)致DC 順流向振幅減小.當(dāng)r=2.0 時,升力系數(shù)PSD 曲線隨約化速度的變化規(guī)律與r=1.0工況基本一致,如圖13(d)和圖13(f)所示,僅在頻率峰值數(shù)量及主頻能量值有細(xì)微差距.由此可見,流體力系數(shù)PSD 曲線隨約化速度的變化規(guī)律:Ur較小時PSD 曲線光滑,隨Ur增大PSD 曲線波動性增強(qiáng),最終PSD 曲線逐漸光滑,與圓柱的運(yùn)動軌跡的變化規(guī)律類似.

圖13 不同工況下,下游圓柱升力系數(shù)功率譜密度(紅線)與阻力系數(shù)功率譜密度(黑線)Fig.13 PSD of lift coefficient(red)and drag coefficient(black)of downstream cylinder under different cases

3.4 結(jié)構(gòu)頻譜機(jī)理分析

結(jié)構(gòu)功率譜密度能反映單位頻帶內(nèi)激勵荷載功率隨頻率的分布,可以得到各個峰值頻率成分的幅值,從而衡量不同峰值頻率對結(jié)構(gòu)響應(yīng)的貢獻(xiàn),是結(jié)構(gòu)頻譜分析的重要手段,渦激振動中激勵荷載包括流體力、位移等[32].通過結(jié)構(gòu)功率譜密度分析,可以判斷激勵荷載信號響應(yīng)的強(qiáng)弱,從而對結(jié)構(gòu)動力學(xué)響應(yīng)進(jìn)行分析.圓柱的橫流向位移與升力系數(shù)功率譜密度曲線包圍的面積Ay和Acl分別代表位移與升力系數(shù)的平均功率值.本節(jié)以Re=100,r=1.0 工況為例,對結(jié)構(gòu)功率譜密度曲線和激勵荷載平均功率值與結(jié)構(gòu)振動的關(guān)系進(jìn)行分析.

由圖14(a)可知,圓柱順流向位移平均功率值A(chǔ)y隨約化速度的變化趨勢與其在y方向的振幅變化趨勢類似.當(dāng)Ur=5 時,上游圓柱位移PSD 曲線出現(xiàn)fsy與3fsy二階頻率峰值,如圖15(a)所示,其主頻能量值較Ur=3 時更大,增大了位移平均功率值,激勵UC 發(fā)生大振幅振動.另一方面,此時渦脫落頻率fsy=0.18 與圓柱的固有頻率(fn,y=0.20)接近,UC 發(fā)生橫流向鎖定現(xiàn)象,產(chǎn)生較大振幅.當(dāng)Ur=5 和8 時,UC 的PSD 曲線表現(xiàn)為光滑的多峰譜特征,其運(yùn)動軌跡表現(xiàn)為規(guī)則的“8”字形;但Ur=9 時,PSD 曲線的波動性增強(qiáng)并在主頻附近出現(xiàn)許多鋸齒形的窄帶峰值,削弱了主頻能量值,導(dǎo)致UC 運(yùn)動的不規(guī)則性及隨機(jī)性增強(qiáng),誘發(fā)產(chǎn)生小振幅順流向運(yùn)動,如圖15(a)所示.進(jìn)一步研究發(fā)現(xiàn),在某約化速度區(qū)間內(nèi),圓柱橫流向振幅Ymax/D與位移信號平均功率值A(chǔ)y的大小成正比.

圖14 Re=100,r=1.0 工況下圓柱位移與升力系數(shù)平均功率值隨約化速度的變化Fig.14 Variation of average power value of cylinder displacement and lift coefficient with reduced velocity at Re=100 and r=1.0

通過對比發(fā)現(xiàn),在區(qū)間A 和B 內(nèi),升力系數(shù)平均功率值隨Ur的變化規(guī)律與升力系數(shù)均方根值的變化規(guī)律類似,且圓柱升力系數(shù)均方根值隨升力平均功率值A(chǔ)cl的增大而增大.與Ur=3 工況相比,如圖15(b)所示,當(dāng)Ur=5 時,UC 的升力系數(shù)PSD 曲線的頻譜成分更豐富,出現(xiàn)fscl,3fscl與5fscl三階頻率峰值,UC 的升力系數(shù)平均功率值增大,進(jìn)而增強(qiáng)了UC 的振動響應(yīng).Ur=9 工況下PSD 曲線的窄帶峰值數(shù)量明顯增多,且分布的頻帶更寬.在區(qū)間C 內(nèi),升力系數(shù)均方根值變化的波動加強(qiáng),但維持了升力系數(shù)均方根值的變化趨勢.

圖15 Re=100 與r=1.0 時,上游圓柱功率譜密度Fig.15 PSD of upstream cylinder at Re=100 and r=1.0

需要注意的是,對不同約化速度區(qū)間的工況進(jìn)行頻譜分析時,必須優(yōu)先考慮振動頻率比(fs/fn,y)的影響.當(dāng)Ur=8 時,UC 升力系數(shù)平均功率值比Ur=3 工況更小,而此時UC 的Ymax/D達(dá)到0.271是Ur=3 工況下的近9 倍.這是由于Ur=8 時的主頻fscl=0.128 接近固有頻率fn,y=0.125,UC 發(fā)生橫流向共振并產(chǎn)生較大振幅.由此可見,對升力系數(shù)功率譜密度分析時,振動頻率比與激勵荷載平均功率值是影響結(jié)構(gòu)振動響應(yīng)強(qiáng)弱的兩個重要參數(shù).

4 結(jié)論

基于四步半隱式特征線分裂算子有限元方法,對串列布置三圓柱雙自由度結(jié)構(gòu)體系渦激振動問題進(jìn)行了數(shù)值模擬,分析了雷諾數(shù)、頻率比和約化速度的變化對結(jié)構(gòu)振動響應(yīng)及頻譜特性的影響,主要結(jié)論如下:

(1)UC 兩個方向的振幅均比較小,雷諾數(shù)、頻率比對振幅的影響較小.MC 的鎖定區(qū)間的范圍隨著雷諾數(shù)的增大而擴(kuò)大,其動力響應(yīng)受頻率比的影響很小.本文的臨界約化速度在6.0 附近,且受雷諾數(shù)的影響較小.DC 的順流向振幅受頻率比和雷諾數(shù)的影響更為顯著.

(2)UC 的CD-mean 與CL-rms 隨約化速度的變化出現(xiàn)先增大后降低直至平穩(wěn)的趨勢,受雷諾數(shù)和頻率比影響較小.由于MC 受到UC 尾流的影響,導(dǎo)致其CL-rms 與CD-mean 與UC 完全不同,約化速度較小時受雷諾數(shù)和頻率比的影響較大.DC 的CL-rms 受雷諾數(shù)和頻率比的影響十分顯著.中下游圓柱CD-mean的變化趨勢與UC 類似.

(3) 圓柱流體力系數(shù)PSD 曲線的波動性隨雷諾數(shù)和頻率比的增大而增強(qiáng),主峰能量值越大圓柱振動響應(yīng)越積極,PSD 曲線波動性越大,圓柱運(yùn)動軌跡越不規(guī)則.當(dāng)升、阻力系數(shù)PSD 曲線的主峰重疊時,圓柱沿著非對稱軌跡運(yùn)動.約化速度較小時PSD 曲線光滑,在大振幅區(qū)間內(nèi)PSD 曲線波動性增強(qiáng)且頻譜成分增多,最終PSD 曲線趨于光滑.

(4)激勵荷載平均功率值隨Ur的變化趨勢與對應(yīng)的結(jié)構(gòu)動力響應(yīng)的變化類似.在同一約化速度區(qū)間內(nèi),結(jié)構(gòu)振動響應(yīng)的強(qiáng)弱與位移的平均功率值成正比.對升力系數(shù)功率譜密度分析時,區(qū)間A 和B 內(nèi)升力系數(shù)均方根值與升力平均功率值成正比.在其他區(qū)間振動頻率比對結(jié)構(gòu)振動響應(yīng)的影響更大.

猜你喜歡
雷諾數(shù)升力振幅
高速列車車頂–升力翼組合體氣動特性
無人機(jī)升力測試裝置設(shè)計(jì)及誤差因素分析
基于自適應(yīng)偽譜法的升力式飛行器火星進(jìn)入段快速軌跡優(yōu)化
基于Transition SST模型的高雷諾數(shù)圓柱繞流數(shù)值研究
十大漲跌幅、換手、振幅、資金流向
十大漲跌幅、換手、振幅、資金流向
滬市十大振幅
升力式再入飛行器體襟翼姿態(tài)控制方法
失穩(wěn)初期的低雷諾數(shù)圓柱繞流POD-Galerkin 建模方法研究
基于轉(zhuǎn)捩模型的低雷諾數(shù)翼型優(yōu)化設(shè)計(jì)研究
文成县| 长乐市| 鄄城县| 庆阳市| 贵定县| 阿勒泰市| 调兵山市| 清徐县| 浑源县| 黄平县| 明光市| 溧阳市| 车险| 桂林市| 杭锦后旗| 漾濞| 开平市| 承德县| 龙山县| 苗栗县| 汶上县| 出国| 金山区| 舒城县| 龙州县| 广南县| 大邑县| 德钦县| 得荣县| 南宁市| 孟州市| 青冈县| 南通市| 图们市| 应城市| 宁阳县| 邢台市| 丰顺县| 辛集市| 忻城县| 绥阳县|