低雷諾數(shù)下的釘螺繞流特征數(shù)值模擬

，，，　

(1.長江水利委員會(huì) 綜合管理中心，武漢　 430010；2.河海大學(xué) 水利水電學(xué)院，南京　210098；3.長江科學(xué)院 院長辦公室，武漢　430010)

1　研究背景

釘螺是日本血吸蟲病傳播的唯一中間宿主，通常隨水流遷移擴(kuò)散。釘螺的在動(dòng)水、靜水中的起動(dòng)、沉降和漂移規(guī)律是設(shè)計(jì)水利阻螺措施的主要依據(jù)。張威等[1-2]對(duì)釘螺的動(dòng)水沉降規(guī)律進(jìn)行了觀測及研究，長江科學(xué)院等單位通過實(shí)驗(yàn)觀察建立了釘螺的靜水沉降公式[3-6]?；谏鲜鲅芯?，長江科學(xué)院、湖北省血吸蟲病防治研究所、武漢大學(xué)等單位提出了沉螺池等有效的水利阻螺設(shè)施。通過多年的科學(xué)防控，我國釘螺分布得到有效控制，為進(jìn)一步提高水利血防工程的阻螺效率，更有效地控制釘螺擴(kuò)散，需對(duì)釘螺在水流中的運(yùn)動(dòng)規(guī)律做更精細(xì)的研究。

圖1　肋殼釘螺

釘螺的外形對(duì)釘螺在水中的運(yùn)動(dòng)規(guī)律有較大的影響。如圖1所示釘螺的外形近似圓錐體，代表釘螺形體的特征值主要是螺長h和螺徑d。通過螺長h與螺徑d的比值φ將釘螺劃分為幼螺(φ<2.0)、中螺(2.0≤φ≤2.5)、大螺(φ>2.5)3級(jí)，不同級(jí)別的釘螺力學(xué)性質(zhì)也有所差別[5]。根據(jù)螺殼表面縱紋的不同,釘螺分為光殼和肋殼2類，光殼釘螺與肋殼釘螺對(duì)水流阻力的影響也各有不同，分布在湖沼和水網(wǎng)區(qū)域釘螺一般為肋殼釘螺[6]。據(jù)觀測，釘螺一般在岸邊水流較緩的區(qū)域活動(dòng)，釘螺在靜水沉降時(shí)一般能保持平穩(wěn)，螺長方向多與水流方向垂直，但在動(dòng)水中漂移的姿態(tài)則復(fù)雜得多。

本文使用開源流體力學(xué)計(jì)算軟件OpenFOAM對(duì)較低雷諾數(shù)肋殼釘螺在螺徑與水流平行(簡稱“橫向水流”)和螺徑與水流垂直(簡稱“縱向水流”)2個(gè)典型條件下的繞流開展數(shù)值模擬研究，對(duì)釘螺繞流中的尾流結(jié)構(gòu)、卡門渦街的形成等現(xiàn)象進(jìn)行分析，對(duì)比2個(gè)方向上的差異，為進(jìn)一步掌握釘螺擴(kuò)散規(guī)律、改進(jìn)水利阻螺措施提供理論依據(jù)。

2　數(shù)學(xué)模型、網(wǎng)格劃分與數(shù)值方法

鈍體繞流的數(shù)值模擬一直以來都是流體力學(xué)領(lǐng)域的研究熱點(diǎn)，國內(nèi)外眾多學(xué)者在數(shù)值模擬方面已有深入細(xì)致的研究，但對(duì)釘螺外形這樣復(fù)雜邊界的鈍體繞流現(xiàn)象的研究卻十分少見。

本文運(yùn)用有限體積法，利用高斯公式對(duì)二維不可壓縮流體的控制方程進(jìn)行離散，將控制方程轉(zhuǎn)換為代數(shù)方程來求解，其中，對(duì)流擴(kuò)散項(xiàng)、壓力梯度項(xiàng)采用二階精度的中心格式離散，非恒定項(xiàng)使用一階精度的隱式離散，時(shí)間步長根據(jù)雷諾數(shù)的不同選取適當(dāng)?shù)臅r(shí)間間隔。壓力與速度耦合采用PISO算法，從而在數(shù)值計(jì)算中能夠比較準(zhǔn)確地滿足邊界條件提高解的精度。離散得到的代數(shù)方程采用Gauss-Seidel迭代法求解。

2.1　控制方程

在笛卡爾坐標(biāo)系下二維不可壓縮黏性流動(dòng)控制方程的無量綱形式為：

(1)

(2)

(3)

式中：u,v分別是x,y方向上的流速：Re為雷諾數(shù);p為流體的靜壓強(qiáng)。

2.2　計(jì)算區(qū)域及模擬對(duì)象

本文選取較有代表性的肋殼中螺進(jìn)行數(shù)值模擬研究。本文模擬的中螺特征值φ=2.27，螺高為h，螺徑為d，其等面積圓的直徑為D，計(jì)算區(qū)域?yàn)闊o界流動(dòng)狀態(tài)。為了使計(jì)算邊界對(duì)數(shù)值模擬的影響最小，本文嘗試使用足夠大的計(jì)算區(qū)域來消除邊界對(duì)計(jì)算的影響，當(dāng)釘螺距來流入口邊界超過20D時(shí)，上游及上下邊界對(duì)流動(dòng)的影響幾乎可以忽略，因此本文的計(jì)算區(qū)域高達(dá)40D×40D。釘螺位于計(jì)算區(qū)域的中部，橫向水流條件下螺徑與x軸平行，縱向水流條件下螺徑與x軸垂直。

各邊界條件如下所述。①來流面(inlet)：U為均勻來流的流速，來流方向沿x軸正方向。②出流面(outlet)：指定為壓力出口。③固墻(釘螺)：采用無滑動(dòng)條件。④上下邊界(wall)：采用無滑動(dòng)條件，速度與來流面一致。

圖2　釘螺附近網(wǎng)格劃分

2.3　網(wǎng)格劃分

如何擬合螺體的復(fù)雜邊界是本研究要解決的重要問題之一。此研究使用非結(jié)構(gòu)化的三角形網(wǎng)格貼合釘螺的復(fù)雜邊界，并對(duì)釘螺表面附近的網(wǎng)格進(jìn)行加密處理，以保證物理量變化較劇烈區(qū)域的計(jì)算精確性，如圖2所示。作比對(duì)的二維圓柱體盡管具有規(guī)則的幾何邊界，在網(wǎng)格劃分時(shí)有曲線網(wǎng)格、結(jié)構(gòu)化網(wǎng)格等多種選擇，但為保證計(jì)算結(jié)果更具可比行，在計(jì)算二維圓柱繞流時(shí)也使用相同算法劃分為三角形網(wǎng)格。

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

顆粒雷諾數(shù)Red、阻力系數(shù)CD、斯特羅哈數(shù)St等無量綱數(shù)是描述低雷諾數(shù)下鈍體繞流規(guī)律的常見參數(shù)，其定義分別為：

(4)

(5)

(6)

式中：U為均勻來流的流速；D為特征長度，本文取圓的直徑或釘螺等面圓直徑；ν為運(yùn)動(dòng)黏滯系數(shù)；FD為鈍體受到的阻力；ρ為流體密度；S為迎流面面積；f為漩渦脫落頻率，本文通過阻力系數(shù)CD的周期變化來得到漩渦脫落的周期。

3.1　算法有效性驗(yàn)證

確定以二維圓柱繞流驗(yàn)證計(jì)算方法的準(zhǔn)確性及有效性的理由有3點(diǎn):①二維圓柱繞流計(jì)算是驗(yàn)證外部流求解的事實(shí)標(biāo)準(zhǔn)之一;②釘螺繞流的水流結(jié)構(gòu)實(shí)驗(yàn)觀測及數(shù)值模擬資料都較少，而二維圓柱繞流有大量的試驗(yàn)觀測和模擬計(jì)算數(shù)據(jù)作比對(duì);③目前在研究釘螺運(yùn)動(dòng)研究時(shí)多將其概化圓形或球形。

二維低速定常繞流的流型與顆粒雷諾數(shù)Red有關(guān)。Red>4時(shí)流體脫離圓柱表面，在下游形成一對(duì)固定不動(dòng)的對(duì)稱回流區(qū)域；Red>47時(shí)圓柱后緣上下兩側(cè)有渦周期性地輪流脫落，形成規(guī)則排列的卡門渦街；Red>180后尾流呈現(xiàn)出三維形態(tài)[7]。

本文首先計(jì)算Red=40，80時(shí)二維圓柱繞流與已有數(shù)值模擬數(shù)據(jù)比對(duì)以驗(yàn)證算法:Red=40時(shí)，隨著計(jì)算進(jìn)行阻力系數(shù)CD、回流區(qū)長度Lw逐漸趨于穩(wěn)定；Red=80時(shí)，阻力系數(shù)、升力系數(shù)出現(xiàn)周期性震蕩，同時(shí)下游回流區(qū)不再存在，取而代之的是渦周期性脫落而形成的卡門渦街，斯特羅哈數(shù)St為描述卡門渦街的脫落頻率的無量綱數(shù)。阻力系數(shù)、回流區(qū)長度和斯特羅哈數(shù)等與Mittal等[7-10]的數(shù)值模擬結(jié)果基本一致，驗(yàn)證了算法的有效性。

3.2　橫向水流條件模擬

本文首先模擬了釘螺處于橫向水流時(shí)，在不同雷諾數(shù)(Red=1，2，3，4，10，20，30，40，50)下的流場，并對(duì)尾流的形態(tài)及變化規(guī)律作了初步分析。計(jì)算表明，釘螺繞流流動(dòng)同二維圓柱繞流一樣，隨著雷諾數(shù)的增加，展現(xiàn)出2個(gè)重要的流態(tài)變化：一是雷諾數(shù)達(dá)到一定值后流體脫離螺體在下游形成穩(wěn)定的附著渦;二是雷諾數(shù)進(jìn)一步增大后釘螺下游附著渦瓦解形成卡門渦街。

圖3(a)給出了Red=1時(shí)，釘螺周圍流場的跡線圖，釘螺表面流動(dòng)無分離現(xiàn)象，流體也未脫離螺面；圖3(b)為Red=10時(shí)，流體經(jīng)過釘螺表面后開始脫落并形成固定不動(dòng)的非對(duì)稱旋渦，渦內(nèi)流體自成封閉回路而形成“死水區(qū)”，尾流區(qū)長度Lw為釘螺末端到“死水區(qū)”末端的長度。

圖4　Red=40時(shí)不同時(shí)刻渦量云圖

Red≥30時(shí)，附著渦瓦解，流體通過釘螺后下游不再是定常流，而是周期性上下脫落形成卡門渦街。圖4選取了Red=40時(shí)一個(gè)脫落周期[0，T]內(nèi)3個(gè)具有代表性時(shí)刻的渦量云圖(T為一個(gè)脫落周期時(shí)長)。渦脫落的周期圖中清晰地表現(xiàn)了流場運(yùn)動(dòng)的周期特性，在t=0和T時(shí)刻，渦量圖實(shí)際上是相同的，0時(shí)刻上部渦即將脫落，而t=T/2時(shí)刻的渦量圖是t=0和T時(shí)刻的鏡像，下部渦即將脫落。實(shí)際上，t=T/4與3T/4時(shí)刻也互為鏡像。

本文計(jì)算釘螺的雷諾數(shù)Red、阻力系數(shù)CD、斯特羅哈數(shù)St等無量綱數(shù)時(shí)，特征長度D均取釘螺等面圓直徑。從表1可知：

(1)Red≤3時(shí)，流體未從釘螺表面脫離,Lw/D=0，此雷諾數(shù)范圍的繞流稱為斯托克斯區(qū)。

(2)4≤Red≤20時(shí),流體開始從釘螺表面脫落形成尾流，且隨著Red的增加“回流區(qū)”長度Lw/D明顯增大，此雷諾數(shù)范圍的繞流稱為對(duì)稱尾流區(qū)。

(3)Red≥30時(shí),尾流出現(xiàn)非穩(wěn)定態(tài)，回流區(qū)演變?yōu)榭ㄩT渦街，描述旋渦脫落頻率的斯特羅哈數(shù)St隨著Red增加基本保持穩(wěn)定，此雷諾數(shù)范圍的繞流稱為卡門渦街區(qū)。

表1　橫向水流不同雷諾數(shù)下釘螺阻力系數(shù)、尾流區(qū)長度和斯特羅哈數(shù)

注：Red≥30后的阻力系數(shù)CD為一個(gè)脫落周期內(nèi)的算術(shù)平均值

在螺徑與水流方向平行時(shí)，與等容圓柱相比，二維釘螺繞流隨Red的增加呈現(xiàn)相同的3種尾流形態(tài)，不同之處在于：

(1)二維圓柱繞流分離點(diǎn)臨界雷諾數(shù)Rec約為4[9]，釘螺則在3～4之間。

(2)二維圓柱繞流由對(duì)稱尾流區(qū)過渡到卡門渦街區(qū)的臨界雷諾數(shù)Rec約為46[9]，釘螺則在20～30之間。

(3)相同雷諾數(shù)下釘螺下游回流區(qū)長度、斯特羅哈數(shù)均小于二維圓柱繞流。

圖5　 Red=40時(shí)CD-t關(guān)系

當(dāng)Red≤20時(shí)，水流通過釘螺后均為恒定流，阻力系數(shù)亦為恒定值；其中，Red<3是水流通過釘螺后幾乎無流動(dòng)分離，阻力以摩擦力為主；4≤Red≤20時(shí)水流通過釘螺后有流動(dòng)分離，釘螺后方回流形成一對(duì)駐渦，阻力由摩擦阻力和壓差阻力2部分組成；當(dāng)Red≥30后，流動(dòng)不再保持恒定，水流在釘螺后方交替脫落形成卡門渦街，此時(shí)阻力主要由壓差阻力產(chǎn)生，阻力隨旋渦脫落的規(guī)律周期性變化。圖5截取了Red=40時(shí)CD隨時(shí)間變化關(guān)系圖，明顯看出CD隨時(shí)間變化的波形圖是由2個(gè)波疊加而成，該現(xiàn)象在其他鈍體繞流現(xiàn)象(如：圓柱、方柱繞流)中均未出現(xiàn)。

3.3　縱向水流條件模擬

基于相同的算法和模型，本文模擬了釘螺處于縱向水流條件時(shí)，在不同雷諾數(shù)(Red=20，40，80，120，140，160)下的流場。計(jì)算表明，螺徑平行水流方向與螺徑垂直水流方向時(shí)一樣，隨著雷諾數(shù)的增加，同樣展現(xiàn)出3種重要的流態(tài)變化。圖6(a)給出了Red=40時(shí)釘螺周圍流場的跡線圖，釘螺表面流動(dòng)的流體未脫離螺面；圖6(b)為Red=80時(shí)，流體經(jīng)過釘螺表面后開始脫落并形成固定不動(dòng)的非對(duì)稱旋渦，渦內(nèi)流體自成封閉回路而形成“死水區(qū)”。

圖6　縱向水流條件下釘螺繞流跡線圖

2種流向下的差異也十分明顯，從表2可以看出：

(1)縱向水流條件下，得到2個(gè)臨界雷諾數(shù)Rec分別在區(qū)間(40，80)和區(qū)間(120，140)范圍內(nèi)，遠(yuǎn)大于橫向水流條件下的臨界值。

(2)處于第2種形態(tài)(對(duì)稱尾流區(qū))時(shí)，隨著雷諾數(shù)的增加尾流區(qū)長度并無明顯增長。

(3)在第3種尾流形態(tài)(卡門渦街區(qū))，渦的脫落頻率也基本保持穩(wěn)定，但縱向水流條件下渦的脫落頻率遠(yuǎn)大于橫向水流條件下的頻率。

表2　縱向水流不同雷諾數(shù)下釘螺阻力系數(shù)、尾流區(qū)長度和斯特羅哈數(shù)

注：Red≥140后的阻力系數(shù)CD為一個(gè)脫落周期內(nèi)的算術(shù)平均值

不同流向下Red-CD關(guān)系對(duì)比如圖7所示。從圖7看出，2種流向下阻力系數(shù)總體上隨Red的增加而減小，但自流體從釘螺表面開始脫落后，阻力系數(shù)CD減小的速度明顯變緩，且卡門渦街區(qū)的CD變化規(guī)律同橫向水流條件下一樣，也是由2個(gè)波形疊加而成。Red=20時(shí)，縱向水流條件下繞流尚處于斯托克斯區(qū)，但橫向水流條件下繞流已處于卡門渦街區(qū)。盡管此時(shí)縱向水流條件下阻力系數(shù)>橫向水流條件的阻力系數(shù)，隨著Red的增加，縱向水流條件下的阻力系數(shù)迅速下降，Red>35.2后阻力系數(shù)反而低于橫向水流條件下的阻力系數(shù)。這一反常現(xiàn)象說明釘螺繞流在卡門渦街區(qū)比較特殊，阻力系數(shù)的算術(shù)平均值并不能正確反映水流對(duì)釘螺的阻力作用。

圖7　不同流向下Red-CD關(guān)系對(duì)比

4　結(jié)論與展望

4.1　結(jié)　論

通過對(duì)低雷諾數(shù)下的橫向和縱向2個(gè)流向下釘螺繞流的直接數(shù)值模擬(DNS)可知，同其他鈍體繞流現(xiàn)象一樣，釘螺繞流在不同的雷諾數(shù)下呈現(xiàn)出3種不同的尾流形態(tài)，不同在于：

(1)3種尾流形態(tài)轉(zhuǎn)變的臨界雷諾數(shù)Rec有差別，釘螺橫向水流<圓柱<釘螺縱向水流。

(2)釘螺在橫向水流條件下，下游回流區(qū)長度隨雷諾數(shù)的增加明顯，但在縱向水流條件幾乎不隨雷諾數(shù)變化。

(3)阻力系數(shù)CD隨雷諾數(shù)增加而減小。

綜上所述，釘螺繞流和其他鈍體繞流現(xiàn)象相比，有共性也有不同。釘螺外形復(fù)雜有其獨(dú)特的水動(dòng)力學(xué)特性，在進(jìn)行釘螺的水力學(xué)研究時(shí)不能簡單概化為簡單的幾何體。

4.2　展　望

今后還需開展以下方面的研究：

(1)提高釘螺數(shù)值模擬的精度，進(jìn)一步提高釘螺模擬效果。

(2)開展釘螺繞流試驗(yàn)觀測研究以更有效地驗(yàn)證數(shù)值模擬結(jié)果，進(jìn)一步豐富對(duì)釘螺隨水流擴(kuò)散規(guī)律的認(rèn)識(shí)。

(3)釘螺在卡門渦街區(qū)CD的變化規(guī)律較為特殊，需進(jìn)一步深入研究。