三維有限長圓柱繞流數值模擬

王曉聰，桂洪斌，劉洋

哈爾濱工業(yè)大學（威海）船舶與海洋工程學院，山東威海264209

0 引 言

繞流現象廣泛存在于自然界中，無論是陸上結構物（如橋梁上的纜索），還是海洋中各類柱狀結構（如石油平臺的樁腿），都會受到繞流的影響。圓柱繞流一般會產生漩渦脫落并誘發(fā)流體的脈動載荷與結構振動［1］，從而發(fā)生渦激振動現象。當振動嚴重時，會對結構物本身造成疲勞等結構破壞［2］，所以對經典的圓柱繞流現象進行研究具有重要意義。截至目前，盡管學術界對于經典的圓柱繞流流體力學問題進行了諸多研究，但對于有限長圓柱（Finite Cylinder，FC）的研究卻并不多見。鑒于圓柱上端端面的下洗作用，有限長圓柱的繞流現象與無限長圓柱（Infinite Cylinder，IFC）的繞流現象相比區(qū)別較大，其三維特性顯著［3］，故非常有必要對其開展研究。

首先，在實驗方面，國內外學者主要從圓柱長徑比AR、自由端斯托羅哈爾數St、自由端阻力系數的變化規(guī)律等方面進行了研究。Park等［4］利用粒子圖像測速（Particle Imagine Velocimetry，PIV）技術對雷諾數Re=7.5×103的有限長圓柱進行了實驗分析，發(fā)現其自由端成對出現梢渦，經分析后發(fā)現梢渦是由于端面的下洗流動和圓柱的繞流流動相互作用產生。Benitz等［5］在雷諾數 Re=2.9×103和傅汝德數Fr=0.65的條件下對有限長圓柱同時進行了一系列長徑比的實驗和數值模擬研究。研究發(fā)現，相比于無限長圓柱，有限長圓柱繞流的阻力系數更小，且在AR=2時阻力系數減小得最為明顯。此外，當AR＜3時，由于自由端效應，下洗的流動會將卡門渦街完全抑制，如圖1所示。

其次，在數值模擬方面，Fr?hlich 等［6］使用LESOCC2代碼，通過大渦模擬（Large Eddy Simula?tion，LES）方法求解不可壓縮的N-S方程，并計算了AR=2.5的有限長圓柱。模擬中，分別使用了SM（Smagorinsky model）和 DSM（Dynamic Smago?rinsky model）亞網格尺度（Subgrid Scale，SGS）模型。結果發(fā)現，DSM對于模擬當前長徑比下的有限長圓柱缺陷明顯，它會導致計算中渦粘性偏大，致使下壁面邊界條件從邊界層向對稱面轉變，而SM模型則能更好地完成仿真模擬。Palau-Salvador等［7］基于浸入邊界法，使用大渦模擬對2種有限長圓柱（AR=2.5，5）在 2種雷諾數（Re=4.3×104，2.2×104）工況下進行了實驗和仿真分析。結果發(fā)現，較短圓柱的交替漩渦脫落僅發(fā)生在接近壁面處，且不穩(wěn)定。

本文基于前人的研究結果，擬分別對具有相同幾何尺寸且在相應計算域下為有限長和無限長的圓柱，選取前人研究較多的Re=3.9×103工況進行模擬，以驗證本文所提方法的可行性。然后，分別對有限長圓柱流場的時均特性和瞬時特性進行分析，以模擬出有限長圓柱特有的“馬蹄”渦、梢渦及“蘑菇”渦。最后，將得到的數據和流場特征與無限長圓柱的模擬結果進行對比，以了解有限長圓柱流場的三維特性和自由端面以及固定壁面對流場的影響，并期望從流場的渦泄角度來解釋有限長圓柱阻力系數損失的原因，為半空間結構聲振問題研究提供參考。

1 數學模型

1.1 控制方程

本文將采用大渦模擬的湍流模型進行數值模擬，該方法的基本原理是將漩渦在空間過濾，其中大漩渦采用直接模擬方法，而小漩渦則采用SGS模型來模擬其對大渦的影響。當考慮粘性、不可壓縮性的水作為流體介質時，N-S方程可寫為

式中：ρ為流體密度；p為壓力；ν為流體的運動粘性系數；ui，uj為速度分量；xi，xj為位移分量；t為時間。從湍流瞬態(tài)運動方程中將尺度比濾波函數尺度小的渦濾掉，分解出描寫大渦流場的運動方程。采用濾波函數處理下的N-S方程及連續(xù)性方程，即大渦模擬的控制方程為：

本文采用SM模型，定義亞網格應力為

式中：τkk為各向性部分的SGS應力；δij為Kronecker符號；為 SGS湍流粘度，其中Cs為Smagorinsky常數，在外流條件下通常取分別為沿x，y，z軸方向的網格尺寸，Δ為濾波尺寸，為可解尺度變形率的張量表達式。

1.2 幾何模型與網格劃分

為對比有限長圓柱繞流的流動特性，本文選取具有相同幾何尺寸的有限長及無限長圓柱模型進行計算和分析，如圖2所示。圖2（a）中，無限長圓柱計算域的展向高度等于圓柱高度；圖2（b）中，有限長圓柱計算域的展向高度大于圓柱高度。

上述幾何模型采用笛卡爾坐標系作為計算域參考坐標，圖中x，y，z分別表示順流向、橫流向和展向，坐標原點位于圓柱地面中心。其中，有限長圓柱計算域尺寸為10D×20D×6D，無限長圓柱計算域尺寸為10D×20D×πD，D為特征尺寸，即圓柱直徑）。此圓柱展向高度為πD，圓柱中心距離上游入口為5D，距離下游出口為15D，模型阻塞比為5.23%。由于圓柱自由端面與流體域上表面必須留出足夠的空間，以便精確捕捉自由端面附近的流動情況，考慮到自由端上部流場的發(fā)展，本文中自由端面距離流場域上表面約為3D。

圖3所示為圓柱計算域網格尺寸劃分，其中無限長圓柱劃分俯視圖與有限長圓柱的一致，采用O形網格對圓柱周邊網格區(qū)進行局部加密。為有效監(jiān)測有限長圓柱自由端的流動，該區(qū)域網格劃分也進行了局部加密，以精確捕捉流動特性。

1.3 邊界條件和時間步選擇

邊界條件設置為速度入口、壓力出口、壁面、對稱面，分別對應幾何體的入口面、出口面、側面、上下壁面。

對于n維問題的顯式計算，時間步應滿足Courant-Friedirchs-Lewy（CFL）條件，其表達式為

式中：C為庫朗數（Courant number），Cmax=1；ui為i方向流速；Δi在此處為最小計算網格尺寸；Δt為計算時間步。

對于本文的計算，取n=1，因此只需滿足u(Δt/Δx)≤1即可。此時，u即為順流向速度U0。設置無因次時間Δt*=ΔtD/U0=5.13×10-4，則計算時間步Δt選為2.5×10-3s，C=0.975，故滿足條件。

1.4 模型驗證

為有效驗證計算的準確性，在同等初始條件下，分別對有限長和無限長圓柱計算域進行計算。在Fluent計算流體力學軟件下，設置壓力速度耦合求解的SIMPLE格式，差分格式則采用中心差分。以遠點流速和圓柱直徑為特征參數的雷諾數Re=3.9×103。將計算結果與前人所做實驗和數值模擬結果進行了對比，如表1所示。表中，Sim1指文獻［8］得到的數值模擬結果，EXP1，EXP2和EXP3分別指文獻［9］、［10］和文獻［5］的實驗結果。

表1 數值計算結果與文獻實驗結果的對比Table 1 Comparisonof numericalsimulation and experimental results provided by the literatures

表中列出了時均阻力系數d、升力系數均方根值Clrms、漩渦脫落頻率相關參數斯托羅哈爾數St。斯托羅哈爾數St計算公式為

式中，f為漩渦脫落主頻率。

由表1可知，在同等SGS模型下，無限長圓柱算例IFC1和有限長圓柱算例FC1與實驗結果非常接近，盡管無限長圓柱算例IFC3與有限長圓柱算例FC3更加接近，但是相比計算資源消耗的增加，無限長圓柱算例IFC1和有限長圓柱算例FC1是更合理的方案。

為驗證文獻［6］中關于SGS應力模型的結論，本文針對有限和無限長圓柱，在SM及DSM模型中均計算了相關參數。經對比發(fā)現，對于當前工況，SM模型確實能更精確地模擬亞臨界區(qū)圓柱繞流。在后續(xù)計算中，無論是有限長還是無限長圓柱繞流的計算，本文均選用SM亞網格尺度應力模型，且SM模型常數Cs=0.2［11］。同時，本文將計算結果與前人的實驗結果進行了對比，發(fā)現在IFC1和FC1工況下，得到的結果與實驗結果非常接近，且與作為衡量準確性最強的物理量d最接近［6］，故本文后續(xù)計算均取IFC1和FC1的工況進行運算。

為進一步驗證算例的正確性，分別列出有限長和無限長圓柱沿中心線(y=0，z=πD/2)上的時均順流向速度ux及圓周壓力系數Cp分布，并與文獻［12］的計算結果進行對比，如圖4所示。其中，壓力系數Cp的計算公式為

式中，p，p0分別為監(jiān)測點壓力和無窮遠處壓力。

由圖4可知，本文計算結果與文獻［12］結果誤差在3%以內，這進一步說明了本文計算工況選取合理。此外，通過對比無限長圓柱的計算結果還可發(fā)現，有限長和無限長圓柱的順流向速度在時變?yōu)檎?，這說明有限長圓柱的回流區(qū)更小，自由端面引起的下洗對回流區(qū)的影響導致順流向阻力系數Cd損失，這與表1的計算結果一致。

2 三維流動特性分析

2.1 時均流場特性分析

圖5所示為有限長和無限長圓柱的升、阻力系數穩(wěn)態(tài)時歷曲線。由圖可知，相比于無限長圓柱，有限長圓柱的時均阻力系數要小23.6%，這說明自由端的存在明顯改變了圓柱受力，進一步說明了下洗作用對圓柱繞流流場的干擾。此外，相比于無限長圓柱，有限長圓柱的升力系數Cl的峰值更小，同時也表現得更無規(guī)律性。鑒于2種圓柱的實際幾何尺寸一致，說明在有限長圓柱計算域下，下洗作用至少局部地擾亂了原有亞臨界區(qū)表現出的卡門渦街。

為研究圓柱存在的自由端面對流體下洗作用的影響，將對圓柱展向居中處水平截圓，以其前緣點為初始監(jiān)測點（θ=0°），每隔10°沿主流方向至180°均勻分布圓周監(jiān)測點，分別計算各監(jiān)測點處的壓力系數時均曲線。為與無限長圓柱進行對比，在同等監(jiān)測位置（即無限長圓柱的展向中間）處設置了同等數量的監(jiān)測點，檢測結果如圖6所示。

由圖6可知，有限長和無限長圓柱在中間位置處分離角相同，不同的是對于壓力系數Cp，有限長圓柱與無限長圓柱的差值達到了35%。由此可見，對于較小長徑比的有限長圓柱，端面下洗流可對流動造成巨大干擾，從而對流場帶來極大的改變。此外，由于下洗的作用，圓柱尾流近壁面處流動受到了一定的抑制。

為進一步研究自由端面和固定壁面對圓柱近壁區(qū)尾流的影響，在y=0，x=D的垂向直線上分別為有限長和無限長圓柱均勻設置了60余個監(jiān)測點，用于監(jiān)測順流向速度ux，無因次化處理后的結果如圖7所示（圖中，h為圓柱展向高度）。由圖可知，對于無限長圓柱，在x=D處，其順流向速度ux均為負值。結合對順流向速度ux沿x方向的分布曲線，說明圓柱的回流區(qū)大于D，與前文對壓力系數的分析一致。同時，在z≤0.8h的區(qū)域內，無論有限長圓柱還是無限長圓柱，ux同時為負，但當z≥0.8h時，有限長圓柱速度迅速躍升，從自由端面直到距離自由端面約1.2D的位置，速度變?yōu)闊o窮遠處流速。從本文進一步分析z/h=0和z/h=π時的有限長及無限長圓柱的相對差值可以發(fā)現，自由端對順流向速度ux的影響要大于固定壁面的影響。

2.2 瞬時流場特性分析

由于有限長圓柱的實際流動時，下壁面，即圓柱固定壁面并非光滑壁面，固定壁面的非滑移性質導致產生了“馬蹄”渦。為能同時得到梢渦和“馬蹄”形漩渦，本文將原有有限長圓柱計算中的上、下面對稱條件改為下壁面非滑移，并重新計算得到了有限長圓柱繞流的流場。

為能清晰得到圓柱后漩渦脫落的流動，本文采用Q準則（Q-criterion）定義渦量如下：

式中：Sij為應變張量；Ωij為計算渦量。

在三維笛卡爾坐標系下，簡化公式為

式中，u，v，w分別為3個方向的速度分量。

圖8所示為以Q準則為基準的渦量等值線分布圖。圖8（a）為2個時刻自由端面的梢渦渦量等值線圖。由圖可看出不同時刻的渦量發(fā)展，梢渦是成對出現，另外頂面處漩渦分離角約為70°；圖8（b）為圓柱展向中間位置處的平面渦量等值線圖。由圖可看出，此時的漩渦已不再是標準的卡門渦街漩渦分布，上、下漩渦也不再對稱?？梢娮杂啥说拇_已對中間處的漩渦發(fā)展產生影響，這也印證了本文之前對升、阻力系數分析時得到的結論；圖8（c）為固定圓柱壁面處的渦量等值線圖，由圖可看出，在底面處因同時考慮了底面的非滑移屬性，底面漩渦流動非常復雜，但在圓柱前端仍可清晰地看到“馬蹄”形漩渦的渦線。值得注意的是，底面處圓柱后因受到的下洗作用影響較小，故保留了一定的、類似于無限長圓柱的漩渦脫落屬性。這也間接說明，對于AR=π的圓柱，下洗并未將整個圓柱的流場完全擾亂，自圓柱底面至其展向高度的一半處，卡門渦街仍是流動的主要成分。

圖9所示為有限長圓柱自由端面的x方向速度分量云圖及其自由端面流線圖的局部放大圖。由圖可知，自由端面自前緣開始產生漩渦，漩渦逐漸發(fā)展，并覆蓋自由端面，形成“蘑菇”形漩渦，漩渦發(fā)展區(qū)域約占自由端表面長度的2/3。這說明，從圓柱前緣發(fā)生分離的流體在2/3處出現回旋后，其余流體繼續(xù)向前流動，未發(fā)生回貼；在自由端的最后向下，流體繼續(xù)分為2個部分：一部分向圓柱壁面流去，形成回旋，即梢渦，正是這種梢渦的存在與卡門渦街一起構成了有限長圓柱的主要流動特性；另一部分繼續(xù)向后流去，但其流動軌跡已不再平直，而是向著偏下的方向流動。

圖10所示為圓柱壁面極限流線圖。圖10（a）中，極限流線交匯成沿展向的某一直線，此直線即為流動繞圓柱壁面的分離點構成的線。由圖10（b）所示頂面自由端極限流線圖可明顯看到2個流線源點，該源點即為“蘑菇”形漩渦形成的源點，正是在此分離的流體分別沿中線向兩側回旋，形成漩渦，這與文獻［2］的研究結果一致。

為更清晰地了解圓柱三維渦結構，分別得到有限長和無限長圓柱基于Q準則的渦量等值線圖，如圖11所示。從圖中可清晰地看到，相比于無限長圓柱，有限長圓柱的漩渦結構更復雜，但卻相對較弱，其繞流被限制在相對較小的區(qū)域內，但有限長圓柱的漩渦在趨近下游區(qū)域逐漸變得光滑。從圖 11（a）和圖11（b）中可以看到，由于存在下洗作用，有限長圓柱的回流區(qū)更短，且渦截面形狀明顯呈“三角形”，這與前文流線圖中圓柱后的尾流流場中流線斜向下發(fā)展的現象一致；而無限長圓柱則不存在這種現象，圓柱后的流場比較均勻，截面呈標準的“矩形”。正是由于上述原因，造成了有限長圓柱阻力系數Cd的損失。除尾流渦結構外，從圖11（c）和圖11（d）可以看出，在圓柱柱體與地面交界處，可以清晰地觀察到“馬蹄”形漩渦，這正是由底面的非滑移屬性造成的。

3 結 論

本文基于大渦模擬數值方法，分別對Re=3.9×103工況下的有限長和無限長圓柱進行了數值模擬。通過與文獻研究結果的對比，證明了本文所提方法的正確性。在時域范圍內，本文分別從有限長圓柱的流場速度、阻力系數、升力系數和瞬態(tài)流場云圖等方面與無限長圓柱進行了對比分析，并得到如下結論：

1）相對于無限長圓柱，有限長圓柱的流動更為復雜，其阻力系數要小23.6%；而對于圓周壓力系數的監(jiān)測結果表明，有限長圓柱的壓力系數小35%，其主要原因是自由端面的下洗流動改變了原有卡門渦街的流動模式，造成其尾流流動減弱。同時，有限長圓柱的升力系數也比無限長圓柱的更小，這說明有限長圓柱對圓柱受力起到了明顯的抑制作用。在垂直方向上對有限長和無限長圓柱順流向速度的分析，得出自由端面對有限長圓柱的影響大于固定壁面。

2）基于Q準則（Q=6）對有限長的流場進行了分析，在展向（z=πD，z=1/2πD，z=0）3個面分別進行了渦量等值線分析，發(fā)現自由端對圓柱后緣尾流影響的范圍在AR=3的工況下超過了圓柱高度的一半。在圓柱底面，由于存在固定壁面，圓柱底部前緣會出現有限長圓柱特有的“馬蹄”形漩渦，該漩渦對于圓柱后緣底部流動同樣具有一定影響，導致底部流場同樣非常復雜，但卡門渦街相對于頂面受到的影響相對較小。對圓柱中面(y=0)進行順流向速度云圖及流線分析可以發(fā)現，頂渦，即“蘑菇”形漩渦的發(fā)展區(qū)域約占自由端面長度的2/3，觀察自由端流線分布可清晰得到2個流線源點。此外，繞過自由端的流體會分為2個部分：一部分回貼向壁面，形成梢渦；另一部分沿斜向下的方向流去。后面的渦量等值面則間接證明了這一點。

3）同樣基于Q準則，分別對有限長和無限長圓柱作出渦量等值面，可以看到有限長圓柱的流動比無限長圓柱的弱，回流區(qū)較短，下洗流動的存在致使其界面流場并非像無限長圓柱那樣呈“矩形”，而是呈“三角形”形狀。

參考文獻：

［1］王露，李天勻，朱翔，等.浸沒邊界格子Boltzmann方法的改進及轉動圓柱繞流模擬［J］.中國艦船研究，2016，11（6）：97-103.WANG L，LI T Y，ZHU X，et al.Improved immersed boundary lattice Boltzmann method and simulation of flow over rotating cylinder［J］.Chinese Journal of Ship Research，2016，11（6）：97-103（in Chinese）.

［2］NAUDASCHER E，ROCKWELLD.Flow-induced vi?brations：an engineering guide：IAHR hydraulic struc?tures design manuals 7［M］.Rotterdam：Balkema Pub?lishers，1994.

［3］SUMNER D.Flow abovethefreeend of a sur?face-mounted finite-height circular cylinder：a review［J］.Journal of Fluids and Structures，2013，43：41-63.

［4］PARK C W，LEE S J.Effects of free-end corner shape on flow structure around a finite cylinder［J］.Journal of Fluids and Structures，2004，19（2）：141-158.

［5］BENITZ M A，CARLSON D W，SEYED-AGHAZA?DEH B，et al.CFD simulations and experimental mea?surements of flow past free-surface piercing，finite length cylinders with varying aspect ratios［J］.Comput?ers&Fluids，2016，136：247-259.

［6］FR?HLICH J，RODI W.LES of the flow around a cir?cular cylinder of finite height［J］.International Journal of Heat and Fluid Flow，2004，25（3）：537-548.

［7］PALAU-SALVADOR G，STOESSER T，FR?HLICH J，et al.Large eddy simulations and experiments of flow around finite-height cylinders［J］.Flow，Turbu?lence and Combustion，2010，84（2）：239-275.

［8］WORNOM S，OUVRARD H，SALVETTI M V，et al.Variational multiscale large-eddy simulations of the flow past a circular cylinder：Reynolds number effects［J］.Computers&Fluids，2011，47（1）：44-50.

［9］KAWAMURA T，HIWADA M，HIBINO T，et al.Flow around a finite circular cylinder on a flat plate：cylin?der height greater than turbulent boundary layer thick?ness［J］. Bulletin ofJSME，1984， 27（232） ：2142-2151.

［10］OKAMOTO S，SUNABASHIRI Y.Vortex shedding from a circular cylinder of finite length placed on a ground plane［J］.Journal of Fluids Engineering，1992，114（4）：512-521.

［11］WELLER H G，TABOR G，JASAK H，et al.A tenso?rial approach to computational continuum mechanics using object-oriented techniques［J］.Computers in Physics，1998，12（6）：620-631.

［12］ZHANG H，YANG J M，XIAO L F，et al.Large-ed?dy simulation of the flow past both finite and infinite circular cylinders at Re=3900［J］.Journal of Hydrody?namics（Ser.B），2015，27（2）：195-203.