，，，　

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

1　研究背景

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

圖1　肋殼釘螺

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

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

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

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

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

2.1　控制方程

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

(1)

(2)

(3)

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

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

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

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

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

2.3　網(wǎng)格劃分

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

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

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

(4)

(5)

(6)

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

3.1　算法有效性驗證

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

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

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

3.2　橫向水流條件模擬

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

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

圖4　Red=40時不同時刻渦量云圖

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

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

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

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

(3)Red≥30時,尾流出現(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為一個脫落周期內(nèi)的算術(shù)平均值

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

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

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

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

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

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

3.3　縱向水流條件模擬

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

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

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

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

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

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

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

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

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

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

4　結(jié)論與展望

4.1　結(jié)　論

通過對低雷諾數(shù)下的橫向和縱向2個流向下釘螺繞流的直接數(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ù)雜有其獨特的水動力學(xué)特性，在進(jìn)行釘螺的水力學(xué)研究時不能簡單概化為簡單的幾何體。

4.2　展　望

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

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

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

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