大渦模擬中模型系數(shù)對(duì)方柱繞流的影響

張童偉，聶 欣，陶雪峰

(杭州電子科技大學(xué)機(jī)械工程學(xué)院，浙江 杭州 310018)

大渦模擬中模型系數(shù)對(duì)方柱繞流的影響

張童偉，聶 欣，陶雪峰

(杭州電子科技大學(xué)機(jī)械工程學(xué)院，浙江 杭州 310018)

采用大渦模擬中3種不同Smagorinsky系數(shù)的標(biāo)準(zhǔn)Smagorinsky-Lilly模型和動(dòng)態(tài)模型，對(duì)雷諾數(shù)為22 000的三維方柱繞流進(jìn)行數(shù)值研究.對(duì)計(jì)算結(jié)果中的特征變量進(jìn)行了驗(yàn)證與分析；并對(duì)不同尺度的渦旋及其相互關(guān)系進(jìn)行研究；同時(shí)對(duì)流場(chǎng)不同位置處和不同Smagorinsky系數(shù)對(duì)計(jì)算結(jié)果的影響及其準(zhǔn)確性進(jìn)行了對(duì)比分析.結(jié)果表明：不同Smagorinsky系數(shù)對(duì)流場(chǎng)中的大尺度渦旋以及斯特勞哈爾數(shù)的計(jì)算結(jié)果影響較??；小尺度渦旋的含能量隨Smagorinsky值的減小而減小，且較小Smagorinsky系數(shù)的模型展現(xiàn)出的小尺度渦旋更多；不同模型的在方柱尾流區(qū)不同位置計(jì)算準(zhǔn)確度存在差異.整體來看，動(dòng)態(tài)模型相比于標(biāo)準(zhǔn)Smagorinsky-Lilly模型更為合適.

大渦模擬；方柱繞流；Smagorinsky系數(shù)；動(dòng)態(tài)模型；湍流能譜

0 引 言

方柱繞流廣泛存在于環(huán)境、海洋、水利、動(dòng)力機(jī)械及建筑等諸多工程領(lǐng)域，其繞流渦旋脫落產(chǎn)生的振動(dòng)對(duì)機(jī)械零部件、建筑和橋梁的穩(wěn)定具有重要影響，渦旋升阻力在海洋方面也得到廣泛應(yīng)用，另外方柱后方的回流現(xiàn)象對(duì)泥沙以及污染物的輸運(yùn)起到重要作用.與此同時(shí)，由于方柱幾何模型簡(jiǎn)單，且涉及到剪切層的轉(zhuǎn)捩，流動(dòng)的分離和附著，周期性的渦脫落等復(fù)雜的流動(dòng)現(xiàn)象，因此被廣泛應(yīng)用于實(shí)驗(yàn)和數(shù)值模擬的研究[1-5].其中大渦模擬作為一種效率高且計(jì)算精準(zhǔn)的模型被應(yīng)用于繞流問題的研究.

大渦模擬所用的亞格子模型中Smagorinsky系數(shù)一直是研究者們關(guān)注的對(duì)象.Smagorinsky系數(shù)Cs是對(duì)亞格子渦粘系數(shù)進(jìn)行描述時(shí)引出的一個(gè)參數(shù)，它與過濾尺寸的乘積稱為混合長(zhǎng)度.對(duì)于Cs系數(shù)的確定方法，文獻(xiàn)[6]根據(jù)經(jīng)典的局部各向同性湍流理論，在Kolmogorov能譜慣性子區(qū)范圍內(nèi)，推導(dǎo)得出的Cs取值范圍在0.17～0.21之間；文獻(xiàn)[7]通過對(duì)一些學(xué)者的研究結(jié)果進(jìn)行總結(jié)，得出Cs取值范圍在0.19～0.24之間；而根據(jù)文獻(xiàn)[8]對(duì)渠道湍流的大渦研究表明，較大的Cs值產(chǎn)生的耗散也大，特別是在近壁區(qū)表現(xiàn)得更明顯，因此對(duì)于內(nèi)流問題，其建議Cs值在0.1左右較合適；文獻(xiàn)[9]采用3種不同的Cs系數(shù)來模擬渠道湍流問題，與直接數(shù)值模擬DNS的結(jié)果比較表明，Cs=0.065時(shí)計(jì)算結(jié)果更準(zhǔn)確.近年來的研究表明，Cs不是一個(gè)常數(shù)，即使對(duì)于特定的流動(dòng)，它也是一個(gè)與物理模型和數(shù)值計(jì)算緊密相關(guān)的值，需要在使用中不斷調(diào)整才能得到合適的結(jié)果.而對(duì)于繞流問題，文獻(xiàn)[10]在進(jìn)行繞流大渦模擬時(shí)采用的Cs值為0.07；文獻(xiàn)[11]取Cs為0.12進(jìn)行計(jì)算.而關(guān)于Cs系數(shù)對(duì)三維方柱繞流模擬結(jié)果的影響，還沒有較為系統(tǒng)的研究.

本文采用不同Cs系數(shù)下的標(biāo)準(zhǔn)Smagorinsky-Lilly模型和動(dòng)態(tài)模型對(duì)三維方柱繞流進(jìn)行大渦模擬，計(jì)算結(jié)果與文獻(xiàn)[3-4]的結(jié)果進(jìn)行了對(duì)比.研究了方柱下游流場(chǎng)的時(shí)均速度、湍流特性及3D流動(dòng)狀態(tài)，對(duì)比分析了不同Cs系數(shù)下Smagorinsky-Lilly模型的差異，同時(shí)驗(yàn)證了動(dòng)態(tài)模型在湍流中計(jì)算的準(zhǔn)確性，對(duì)方柱繞流的流動(dòng)特性進(jìn)行分析，為湍流問題研究模型的選擇和實(shí)際工程應(yīng)用提供有價(jià)值的參考.

1 控制方程與模型

1.1 控制方程

運(yùn)用大渦模擬(Large Eddy Simulation，LES)方法，對(duì)連續(xù)性和不可壓縮流動(dòng)的Navier-Stokes方程進(jìn)行濾波，得到如下的控制方程：

(1)

(2)

(3)

1.2 亞格子模型

要實(shí)現(xiàn)大渦數(shù)值模擬，必須構(gòu)造亞格子應(yīng)力的封閉模式.其中以Boussinesq假設(shè)為基礎(chǔ)的計(jì)算公式如下：

(4)

(5)

1.2.1 標(biāo)準(zhǔn)Smagorinsky-Lilly模型

標(biāo)準(zhǔn)Smagorinsky-Lilly模型最早是由J. S. Smagorinsky于1963年基于普朗特混合長(zhǎng)度模式的原理而提出的[12]，將渦粘系數(shù)寫成以下形式：

(6)

1.2.2 動(dòng)態(tài)Smagorinsky-Lilly模型

動(dòng)態(tài)Smagorinsky-Lilly模型是D. K. Lilly于1992年基于Germano等式提出的動(dòng)態(tài)確定Smagorinsky系數(shù)的方法，簡(jiǎn)稱動(dòng)態(tài)模型DCs[13-14].將Germano等式用Smagorinsky模式代入得到：

(7)

采用最小誤差法來解決超定問題，對(duì)結(jié)果進(jìn)行系統(tǒng)平均得：

(8)

2 數(shù)值試驗(yàn)

2.1 算例描述

圖1 計(jì)算域幾何尺寸簡(jiǎn)圖

本文算例以高雷諾數(shù)下的方柱繞流為研究對(duì)象，采用不同的Smagorinsky系數(shù)及動(dòng)態(tài)模型對(duì)其進(jìn)行大渦模擬，計(jì)算結(jié)果與相關(guān)文獻(xiàn)的結(jié)果進(jìn)行對(duì)比驗(yàn)證.計(jì)算區(qū)域?yàn)?0.5D×14D×4D(D為方柱直徑)，其中方柱上游4.5D，下游15D，如圖1所示.入口速度U=0.535 m/s，均勻流動(dòng)，出口條件為自由流動(dòng)，四周表面均設(shè)為對(duì)稱邊界條件.工作介質(zhì)為水(密度ρ=998.2 kg/m3，動(dòng)力粘度μ=1.005×10-3Pa·s)，方柱的繞流雷諾數(shù)Re=ρUD/μ=22 000.

2.2 網(wǎng)格與計(jì)算過程

根據(jù)大渦模擬的物理原則，模型第一層網(wǎng)格的y+應(yīng)小于1.因此，在粘性底層(y+<5)之內(nèi)，設(shè)置了3～5層網(wǎng)格，根據(jù)最終計(jì)算的收斂結(jié)果，壁面處第一層網(wǎng)格y+在0.5～0.9之間.流域截面的網(wǎng)格形態(tài)如圖2所示，所有模型的網(wǎng)格總數(shù)均為320萬(wàn).采用FLUENT14.0作為求解器，有限體積法離散控制方程，對(duì)于壓力與速度的耦合采用SIMPLE算法，時(shí)間項(xiàng)采用二階隱式差分，壓力項(xiàng)采用二階迎風(fēng)格式進(jìn)行離散，對(duì)流和擴(kuò)散項(xiàng)采用二階中心差分格式.根據(jù)網(wǎng)格尺度和速度尺度設(shè)定時(shí)間步長(zhǎng)，保證每一次迭代都在一個(gè)網(wǎng)格范圍內(nèi)[16]，此處Δt設(shè)置為1×10-3s，每個(gè)子步迭代8次后達(dá)到所設(shè)定的殘差值(1×10-3).本算例使用曙光5 000計(jì)算機(jī)(單節(jié)點(diǎn)32核)并行計(jì)算.先進(jìn)行穩(wěn)態(tài)分析，獲得較合適的初場(chǎng)后進(jìn)行瞬態(tài)計(jì)算，在取樣前計(jì)算10個(gè)的流動(dòng)周期(以D/U作為一個(gè)流動(dòng)周期)，取樣后再計(jì)算10個(gè)以上的流動(dòng)周期(約為30個(gè)渦旋脫落周期)，以獲得穩(wěn)定的統(tǒng)計(jì)平均結(jié)果.其中標(biāo)準(zhǔn)Smagorinsky-Lilly模型的計(jì)算時(shí)間為140 h，動(dòng)態(tài)模型的計(jì)算時(shí)間為180 h，由于此模型需要進(jìn)行二次過濾且在計(jì)算過程中不斷改變Cs的值，因此需要消耗更多的計(jì)算資源，約為標(biāo)準(zhǔn)模型的1.3倍.

圖2 計(jì)算域網(wǎng)格

3 結(jié)果與分析

不同Cs系數(shù)和模型下，方柱繞流特征變量的模擬結(jié)果與相關(guān)文獻(xiàn)對(duì)比如表1所示，包括時(shí)均且展現(xiàn)平均的阻力系數(shù)Cd，回流區(qū)長(zhǎng)度L/D，斯特勞哈爾數(shù)St=fD/U，f為通過對(duì)升力系數(shù)隨時(shí)間的變化量進(jìn)行快速傅里葉變換求得的渦旋頻率.

從表1可以看出，所有模型的Cd值和實(shí)驗(yàn)結(jié)果相比都偏大，與前人的研究結(jié)果相似，而DCs模型和Cs=0.04時(shí)的壓力系數(shù)值與實(shí)驗(yàn)值更接近，其余模型的Cd值隨著Cs系數(shù)的增加而變大.對(duì)于回流區(qū)長(zhǎng)度的預(yù)測(cè)，不同的研究者得出的結(jié)果有較大差異[1]，本文的計(jì)算結(jié)果與實(shí)驗(yàn)值更接近，其中Cs=0.04時(shí)的計(jì)算結(jié)果和實(shí)驗(yàn)值對(duì)比差異最小，這與阻力系數(shù)的值有一定的對(duì)應(yīng)關(guān)系.由于亞格子模型只是對(duì)小渦進(jìn)行建模計(jì)算，而對(duì)渦旋頻率起主導(dǎo)作用的大渦在不同Cs系數(shù)下的求解方法相同，因此不同模型的繞流St數(shù)計(jì)算結(jié)果基本相同.

表1 特征變量計(jì)算結(jié)果的對(duì)比

3.1 不同模型時(shí)均流動(dòng)特性對(duì)比與分析

圖3為不同Cs系數(shù)和模型下，方柱后方中心線處時(shí)均速度u和湍動(dòng)能K的分布情況.由于在相同網(wǎng)格精度的情況下，亞格子模型的過濾尺度Δ相同，由Cs與Δ的乘積所構(gòu)成的混合長(zhǎng)度值LS不同，使得亞格子的渦粘系數(shù)值隨Cs值的增大而增大，而較大的渦粘系數(shù)會(huì)使得流動(dòng)耗散偏大.因此，在方柱附近位置流動(dòng)為層流向湍流過渡的區(qū)域，應(yīng)變率張量大而湍流耗散較小[13]，較小的Cs值更適合對(duì)回流區(qū)長(zhǎng)度及阻力系數(shù)等參數(shù)的預(yù)測(cè)，同時(shí)DCs模型在此種流動(dòng)狀態(tài)處會(huì)自動(dòng)減小模型系數(shù)來修正流動(dòng)行為[17].從圖3中可以看出，Cs=0.04時(shí)方柱附近的時(shí)均速度與實(shí)驗(yàn)值最接近，DCs模型的計(jì)算結(jié)果也較為精準(zhǔn).由于同時(shí)，從回流區(qū)最低速度的大小和位置可以看出，不同模型最低速度都位于方柱后方0.85D的距離處，而DCs和Cs=0.04時(shí)的最小值分別為-0.098 m/s和-0.009 m/s，與實(shí)驗(yàn)值偏差更小.而在遠(yuǎn)離方柱的位置，由于渦旋的脫落和破碎，湍流耗散較大，此時(shí)較大的Cs值更適合此種流動(dòng)狀態(tài).因此對(duì)比不同模型在方柱后方較遠(yuǎn)處的速度恢復(fù)值可以發(fā)現(xiàn)，Cs=0.18時(shí)計(jì)算結(jié)果與實(shí)驗(yàn)值更為接近.

圖3 中心線上時(shí)均速度和湍動(dòng)能的分布

從圖3明顯可以看出，在X/D=1.5的位置存在一個(gè)峰值，而此位置為兩側(cè)渦旋交替產(chǎn)生的區(qū)域，因此速度脈動(dòng)大.同時(shí)，在方柱背面的回流區(qū)，由于靠近壁面處的流體以沿壁面的縱向流動(dòng)為主，在層流底層中由壁面效應(yīng)而猝發(fā)的卷吸渦使得此處速度脈動(dòng)值較大[18]，因此該處的湍動(dòng)能值也有一個(gè)明顯的峰值，而Cs=0.18時(shí)得到的湍流粘度較大，在此處峰值也更為明顯，這與文獻(xiàn)[2]和[11]在其他問題的研究中提出的觀點(diǎn)類似.綜合分析表1中的統(tǒng)計(jì)值可以得出，方柱后方較大的速度脈動(dòng)和較低速度恢復(fù)值使得其回流區(qū)長(zhǎng)度偏小而阻力系數(shù)值偏大，從DCs和Cs取0.04和0.12中也可以得到相對(duì)應(yīng)的結(jié)論，此結(jié)論與文獻(xiàn)[2]提出的觀點(diǎn)一致.

總之，對(duì)于較高雷諾數(shù)下方柱繞流的大渦模擬，Cs=0.04時(shí)對(duì)方柱后方回流區(qū)速度場(chǎng)及阻力系數(shù)的預(yù)測(cè)更準(zhǔn)確，而較遠(yuǎn)處的速度恢復(fù)值在Cs=0.18時(shí)與實(shí)驗(yàn)值更吻合，因此，對(duì)于局部位置或流動(dòng)狀態(tài)較為單一的問題(如渠道流動(dòng)，管內(nèi)流動(dòng)等)進(jìn)行研究，合理地選擇Cs系數(shù)對(duì)計(jì)算結(jié)果的精確性至關(guān)重要.而對(duì)于方柱繞流這種不同位置流動(dòng)狀態(tài)差異較大的問題，標(biāo)準(zhǔn)Smagorinsky-Lilly模型很難做到對(duì)整個(gè)流場(chǎng)的各個(gè)位置進(jìn)行精準(zhǔn)預(yù)測(cè)，此時(shí)動(dòng)態(tài)模型更合適.

3.2 方柱繞流渦旋特性分析

圖4 Q準(zhǔn)則數(shù)分布

在湍流動(dòng)能輸運(yùn)過程中，大渦占湍動(dòng)能的絕大部分，因此對(duì)數(shù)坐標(biāo)下湍流能譜的含能量主要集中在小波數(shù)附近，而小渦占湍動(dòng)能耗散的絕大部分，在含能區(qū)與耗散區(qū)之間湍動(dòng)能通過慣性子區(qū)將能量逐級(jí)傳遞下去，此時(shí)的湍流能譜函數(shù)基本滿足波數(shù)k的-5/3冪次律[21].特征波數(shù)與渦旋尺度及含能量的關(guān)系如表2所示，從表2可以看出，大尺度渦旋(小k)在方柱附近位置的含能量最高，隨著脫體渦旋的不斷發(fā)展和破碎，遠(yuǎn)離方柱的位置，大渦的含能量降低，而小渦整體的含能量增加.由于較小的Cs系數(shù)得到的亞格子湍流耗散更小，同時(shí)對(duì)小尺度渦旋的描述更細(xì)致(數(shù)量更多)，單個(gè)渦旋的含能量也會(huì)降低，方柱中心線上不同位置的湍流能譜分布情況如圖5所示，可以看出小尺度渦的含能量隨著Cs系數(shù)的減小而減小.而DCs模型的小尺度渦旋含能量在方柱附近位置(X/D=0.5～2.0處)與Cs=0.04相近，在遠(yuǎn)離方柱位置(X/D=5.0～7.0處)與Cs=0.18相近，這也是DCs模型整體計(jì)算結(jié)果較為準(zhǔn)確的一個(gè)重要因素.

表2 特征波數(shù)與渦旋尺度、含能量的關(guān)系

圖5 中心線上湍流能譜的分布

4 結(jié)束語(yǔ)

本文通過對(duì)大渦模擬中3種不同模型系數(shù)Cs的Smagorinsky-Lilly及動(dòng)態(tài)模型對(duì)方柱繞流計(jì)算結(jié)果的影響進(jìn)行分析得出，不同Smagorinsky系數(shù)Cs對(duì)流場(chǎng)中的大尺度渦旋以及斯特勞哈爾數(shù)St的計(jì)算結(jié)果影響較小.在流場(chǎng)不同位置出，不同模型對(duì)小尺度渦旋及整體計(jì)算結(jié)果的準(zhǔn)確性存在差異.對(duì)于方柱繞流這種渦旋尺度與強(qiáng)度變化較大且雷諾數(shù)較高的繞流問題，動(dòng)態(tài)模型相比于標(biāo)準(zhǔn)Smagorinsky-Lilly模型更為合適.因此，在今后類似繞流問題的大渦模擬計(jì)算中，本文的研究為模型以及模型參數(shù)值的選擇提供了參考.

[1] RODI W, FERZIGER J H, BREUER M, et al. Status of large eddy simulation: results of a workshop[J]. Journal of Fluids Engineering,1997,119(2):248-262.

[2] SOHANKER A, DAVIDSON L, NORBERG L. Large eddy simulation of flow past a square cylinder: comparison of different subgrid scale models[J]. Journal of Fluids Engineering,1999,122(1):39-47.

[3] LYN DA, EINAV S, RODI W, et al. A laser Doppler velocimetry study of ensemble averaged characteristics of the turbulent near wake of a square cylinder[J]. Journal of Fluid Mechanics,1995,304:285-319.

[4] TRIAS FX, GOROBETS A, OLIVA A. Turbulent flow around a square cylinder at Reynolds number 22,000: A DNS study[J]. Computers & Fluids,2015,123:87-98.

[5] ELKHOURY M. Assessment of turbulence models for the simulation of turbulent flows past bluff bodies[J]. Journal of Wind Engineering & Industrial Aerodynamics,2016,154:10-20.

[6] LILLY D K. Representation of small scale turbulence in numerical simulation experiments[J]. Proceedings of IBM Scientific Computing Symposium on Environmental Sciences,1967:195-210.

[7] ROGALLO R S, MOIN P. Numerical simulation of turbulent flows[J]. Annual Review of Fluid Mechanics,1984,16(1):99-137.

[8] DEARDORFF J W. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[J]. Journal of Fluid Mechanics,1970,41(2):453-480.

[9] UDDIN M, MALLIK M S I. Large eddy simulation of turbulent channel flow using smagorinsky model and effects of smagorinsky constants[J]. British Journal of Mathematics & Computer Science,2015,7(5):375-390.

[10] DOOLAN C J. Large eddy simulation of the near wake of a circular cylinder at sub-critical Reynolds number[J]. Engineering Applications of Computational Fluid Mechanics,2014,4(4):496-510.

[11] MURAKAMI S. Overview of turbulence models applied in CWE-1997[J]. Journal of Wind Engineering & Industrial Aerodynamics,1998,74:1-24.

[12] SMAGORINSKY J S. General circulation experiments with the primitive equations, the basic experiment[J]. Monthly Weather Review,1963,91(3):99-164.

[13] LILLY D K. A proposed modification of the Germano subgrid-scale closure method[J]. Physics of Fluids A Fluid Dynamics,1992,4(3):633-635.

[14] GERMANO M, PIOMELLI U, MOIN P, et al. A dynamic subgrid-scale eddy viscosity model[J]. Physics of Fluids A Fluid Dynamics,1991,3(3):1760-1765.

[15] 張兆順，崔桂香，許春曉.湍流大渦數(shù)值模擬的理論和應(yīng)用[M].北京:清華大學(xué)出版社,2008:98-99.

[16] 鄭力銘.ANSYS Fluent 15.0流體計(jì)算從入門到精通[M].電子工業(yè)出版社,2015：43-49.

[17] TU J Y, GUAN H Y, LIU C Q.計(jì)算流體力學(xué):從實(shí)踐中學(xué)習(xí)[M].沈陽(yáng):東北大學(xué)出版社,2014:279-281.

[18] 唐鵬，韓省思，葉桃紅，等.聯(lián)合RANS/LES方法數(shù)值模擬方柱繞流[J].中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào),2010,40(12):80-85.

[19] HUNT J C R. Eddies Stream, and Convergence Zones in Turbulent Flows[J]. Center for Turbulence Research CTR-S 88,1988:193-209.

[20] PIOMELLI U, CHASNOV J R. Large-eddy simulations: theory and applications[M]//Turbulence and transition modelling. Springer Netherlands,1996:269-336.

[21] KOLMOGOROV A N. The local structure of turbulence in an incompressible viscous fluid for very large Reynolds numbers[J]. Proceedings of the Royal Society of London,1941,1890(1890):9-13.

LargeEddySimulationofFlowaroundaSquareCylinder:EffectsofSmagorinskyConstants

ZHANG Tongwei, NIE Xin, TAO Xuefeng

(SchoolofMechanicalEngineering,HangzhouDianziUniversity,HangzhouZhejiang310018,China)

Large eddy simulation of flow around a three-dimensional square cylinder at Re=22000 is performed by using standard model of three Smagorinsky constants and dynamic model, the calculation results of characteristic variable are analysed and validated in this paper. Thus, the different scales vortex and their interrelation are analysed, also the effect of different values of Cs on the calculated result and accuracy is performed in different regions. The results demonstrate that the results of large scale vortex and St number are less affected by the value of Cs, but the energy of small scale vortex decrease with the decrease of the Cs and the model of smaller Cs exhibit more small scale vortex. The calculation accuracy of models is different in different positions of square cylinder. The dynamic model is more appropriate than standard Smagorinsky-Lilly model from overall results.

large eddy simulation; square cylinder; Smagorinsky constant; dynamic model; turbulent spectrum

10.13954/j.cnki.hdu.2017.06.011

2017-05-08

國(guó)家自然科學(xué)基金資助項(xiàng)目(11472095，50909031)

張童偉(1994-)，男，安徽馬鞍山人，碩士研究生，流體力學(xué).通信作者：聶欣副教授，E-mail:xin_nie2000@163.com.

O353

A

1001-9146(2017)06-0055-07