基于浸入邊界法和編程的圓柱繞流的數值模擬

牛朝

摘要：浸入邊界法是在非結構網格下對N-S方程進行求解。本文使用浸入邊界法和有限差分法，運用C++編程實現對靜止圓柱繞流的模擬，在不同雷諾數下求解不可壓縮N-S方程和驗證投影浸入邊界法的可靠性，并分析圓柱周圍速度場、壓力場和渦結構的分布情況。同時，得到了最優的網格劃分范圍。

Abstract： In immersed boundary method （IBM）， N-S equations are solved with unstructured grid. In this article， we used IBM and the finite difference method， used C++ codes to simulate the flow over a static cylinder， to solve incompressible N-S equation and verify the reliability of projection of IB for the different Reynolds number. The distribution of velocity field， pressure field and vortex structure around circular cylinder were calculated and analyzed. At the same time， the best division of Cartesian grid is obtained.

關鍵詞：浸入邊界法;有限差分法;圓柱繞流

Key words： immersed boundary method;finite difference method;flow around a circular cylinder

中圖分類號：O174.63 文獻標識碼：A 文章編號：1006-4311（2018）07-0155-03

0 引言

隨著科技的發展和時代的進步，越來越多的水工建設和科學研究都涉及到流固耦合作用研究。流固耦合現象包含著劇烈的旋渦運動和壓力變化，會增加結構體的不穩定負載，產生噪音和振動，嚴重的情況甚至會造成破壞和災難，所以研究流體-結構耦合問題[1-3]變得非常重要。圓柱繞流是流固耦合研究中經典的繞流問題，在流場中結構體的研究中，圓柱繞流幾乎包含了鈍體繞流的所有復雜現象。在流體力學中，這類問題的研究已經成為基礎研究和工程應用研究中最有意義的課題之一。

本文結合浸入邊界法的思想，基于有限差分離散控制方程，主要研究以下三個方面的內容，首先，通過在不同網格尺寸劃分情況下，比較流場的變化情況，驗證浸入邊界投影法的網格獨立性;其次，計算分析圓柱繞流特性。最后對圓柱繞流的動態特性研究。

1 數值計算模型

1.1 浸入邊界法的基本控制方程

對黏性不可壓縮的流體，在笛卡爾坐標系下包含浸入結構的矩形區域中，流動的無量綱控制方程為：

其中：u是速度矢量，t是時間，p是壓強;f表示在形狀規則的固體邊界上的速度場分布產生的外部體力，雷諾數定義為Re=UL/v，U、L分別為無量綱特征速度和特征長度尺度，其中v表示流體的運動粘度。

1.2 數值計算方法

本文是基于浸入邊界法的模式，采用有限差分法求解流場，采用均勻交錯網格，最后通過編程實現圓柱繞流流場的求解。使用曲率系數對流方案[4]處理對流項，本文中采用Quick法來做計算;時間差分項以顯式Adams-Bashforth 法來離散;壓力項的處理采用投影法（projection method），導出含有混合邊界的壓力場Poisson方程式進行離散求解。

1.3 計算模型及邊界條件設置

物理模型為一矩形求解區域，使用無量綱計算，計算區域長L=20D，寬W=10D，圓柱直徑D=1.0，建立直角坐標系x-y，圓柱圓心坐標（x，y）=（3D，5D），假設入口邊界流體x方向速度符合拋物線分布，最大速度為1.0，y方向上速度為v=0，u=u（y）=-2/50·y2+4·y，v=0;流場上下邊壁為無滑移邊界，u=v=0，整個流場為均勻流場，時間步長Δt=0.001。

2 計算結果分析

2.1 浸入邊界法網格獨立性驗證

由圖1不同雷諾數下的速度場變化圖可看出：Re≤50時，流場還處于對稱尾流區，旋渦還沒有脫落，隨著雷諾數的增大，旋渦逐漸從圓柱體兩側周期性的脫落，形成規則性排列的卡門渦街。當雷諾數Re>300后，遠離圓柱的區域旋渦的排列變得不再是規則的周期排列。

表1給出了圓柱繞流在三種不同網格劃分下，低雷諾數Re=25，40，50時的計算結果的旋渦比較，可以看出在相同雷諾數下隨著網格密度的增加，旋渦的長度L、旋渦內的最大速度U和兩個對稱旋渦的渦核間距d不斷增加，在網格為800×400時所得結果與實驗結果[5]最為接近。

2.2 圓柱繞流的特性

圖2為網格劃分為200×100，Re=1時的速度等值線圖和流線圖，由圖可知圓柱兩側的流線前后上下對稱，流體繞圓柱流動沒有形成旋渦，仍處于層流狀態，與圓柱繞流流場性質相吻合。

endprint

以下結果為網格劃分800×400。如圖3所示為不同雷諾數下的圓柱附近的流線圖，Re>10后繞圓柱分離后的流體在圓柱后形成一對固定不動的對稱旋渦并且不斷被拉長，隨著Re的增大，旋渦逐漸失去對稱性，對稱旋渦消失旋渦開始從兩側不斷地輪流脫落，由Re>100時，可以看出圓柱體的下游的流場不再具有規則性，流場變得較復雜，由此得出的計算結果與圓柱繞流的流場特性相吻合。

表2為Re=25時，旋渦分離的物理參數和文獻的對比，與所得結果吻合得很好，說明了該方法求解圓柱繞流流場的有效性和正確性。

3 結論

通過利用浸入邊界法對單圓柱繞流進行數值模擬計算，得出如下結論：①使用浸入邊界法模擬圓柱繞流，不同網格劃分下所得結果都符合流場的一般性質，只是精確度不同，網格劃分越密，精度越高，在格劃分間距d=0.025D時能很好的滿足流場求解的需求。②使用浸入邊界法編程來實現圓柱繞流流場的數值模擬，所得結果與圓柱繞流的流場性質相吻合，該方法能有效、正確的模擬圓柱繞流流場。

參考文獻：

[1]施耀華.橋塔繞流與流固耦合的數值模擬[D].北京交通大學，2008，06，01.

[2]周光坰，嚴宗毅，許世雄，等.流體力學[M].北京：高等教育出版社，2003.

[3]邢景棠，周盛，崔爾杰.流固耦合力學概述[J].力學進展， 1997，27（1）：19-38.

[4]Goharian N， Jain A， Sun Q. Comparative analysis of sparse matrix algorithms for information retrieval[J]. Computer， 2003， 2： 0-4.

[5]M. Coutanceau and R. Bouard， J. Fluid Mech. 79（2）， 231 （1977）.

[6]P. M. Gresho， R. Chan， C. Upson， and R. Lee， Int. J. Numer. Methods Fluids 4， 619 （1984）.

[7]E. M. SAIKI ANDS. B IRINGEN，Numerical Simulation of a Cylinder in Uniform Flow： Application of a Virtual Boundary Method，JOURNAL OF COMPUTATIONAL PHYSICS ， 1996， 123：450-465.

[8]Taneda S. Experimental investigation of the wakes behind cylinders and plates at low reynolds Numbers[J]. Journal of the Physical Society of Japan， 1956， 11： 302.

[9]Takami H， Keller H B. Steady Two-dimensional Viscous Flow of an in-compressible fluid past a circular cylinder[J]. Physics of Fluids， 1969， 12（12）： II-51-II-56.

[10]姜知亨.浸入邊界法及在大渦模擬中的應用[D].清華大學.endprint