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

牛朝

摘要：浸入邊界法是在非結(jié)構(gòu)網(wǎng)格下對(duì)N-S方程進(jìn)行求解。本文使用浸入邊界法和有限差分法，運(yùn)用C++編程實(shí)現(xiàn)對(duì)靜止圓柱繞流的模擬，在不同雷諾數(shù)下求解不可壓縮N-S方程和驗(yàn)證投影浸入邊界法的可靠性，并分析圓柱周?chē)俣葓?chǎng)、壓力場(chǎng)和渦結(jié)構(gòu)的分布情況。同時(shí)，得到了最優(yōu)的網(wǎng)格劃分范圍。

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.

關(guān)鍵詞：浸入邊界法;有限差分法;圓柱繞流

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

中圖分類(lèi)號(hào)：O174.63 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1006-4311（2018）07-0155-03

0 引言

隨著科技的發(fā)展和時(shí)代的進(jìn)步，越來(lái)越多的水工建設(shè)和科學(xué)研究都涉及到流固耦合作用研究。流固耦合現(xiàn)象包含著劇烈的旋渦運(yùn)動(dòng)和壓力變化，會(huì)增加結(jié)構(gòu)體的不穩(wěn)定負(fù)載，產(chǎn)生噪音和振動(dòng)，嚴(yán)重的情況甚至?xí)斐善茐暮蜑?zāi)難，所以研究流體-結(jié)構(gòu)耦合問(wèn)題[1-3]變得非常重要。圓柱繞流是流固耦合研究中經(jīng)典的繞流問(wèn)題，在流場(chǎng)中結(jié)構(gòu)體的研究中，圓柱繞流幾乎包含了鈍體繞流的所有復(fù)雜現(xiàn)象。在流體力學(xué)中，這類(lèi)問(wèn)題的研究已經(jīng)成為基礎(chǔ)研究和工程應(yīng)用研究中最有意義的課題之一。

本文結(jié)合浸入邊界法的思想，基于有限差分離散控制方程，主要研究以下三個(gè)方面的內(nèi)容，首先，通過(guò)在不同網(wǎng)格尺寸劃分情況下，比較流場(chǎng)的變化情況，驗(yàn)證浸入邊界投影法的網(wǎng)格獨(dú)立性;其次，計(jì)算分析圓柱繞流特性。最后對(duì)圓柱繞流的動(dòng)態(tài)特性研究。

1 數(shù)值計(jì)算模型

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

對(duì)黏性不可壓縮的流體，在笛卡爾坐標(biāo)系下包含浸入結(jié)構(gòu)的矩形區(qū)域中，流動(dòng)的無(wú)量綱控制方程為：

其中：u是速度矢量，t是時(shí)間，p是壓強(qiáng);f表示在形狀規(guī)則的固體邊界上的速度場(chǎng)分布產(chǎn)生的外部體力，雷諾數(shù)定義為Re=UL/v，U、L分別為無(wú)量綱特征速度和特征長(zhǎng)度尺度，其中v表示流體的運(yùn)動(dòng)粘度。

1.2 數(shù)值計(jì)算方法

本文是基于浸入邊界法的模式，采用有限差分法求解流場(chǎng)，采用均勻交錯(cuò)網(wǎng)格，最后通過(guò)編程實(shí)現(xiàn)圓柱繞流流場(chǎng)的求解。使用曲率系數(shù)對(duì)流方案[4]處理對(duì)流項(xiàng)，本文中采用Quick法來(lái)做計(jì)算;時(shí)間差分項(xiàng)以顯式Adams-Bashforth 法來(lái)離散;壓力項(xiàng)的處理采用投影法（projection method），導(dǎo)出含有混合邊界的壓力場(chǎng)Poisson方程式進(jìn)行離散求解。

1.3 計(jì)算模型及邊界條件設(shè)置

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

2 計(jì)算結(jié)果分析

2.1 浸入邊界法網(wǎng)格獨(dú)立性驗(yàn)證

由圖1不同雷諾數(shù)下的速度場(chǎng)變化圖可看出：Re≤50時(shí)，流場(chǎng)還處于對(duì)稱(chēng)尾流區(qū)，旋渦還沒(méi)有脫落，隨著雷諾數(shù)的增大，旋渦逐漸從圓柱體兩側(cè)周期性的脫落，形成規(guī)則性排列的卡門(mén)渦街。當(dāng)雷諾數(shù)Re>300后，遠(yuǎn)離圓柱的區(qū)域旋渦的排列變得不再是規(guī)則的周期排列。

表1給出了圓柱繞流在三種不同網(wǎng)格劃分下，低雷諾數(shù)Re=25，40，50時(shí)的計(jì)算結(jié)果的旋渦比較，可以看出在相同雷諾數(shù)下隨著網(wǎng)格密度的增加，旋渦的長(zhǎng)度L、旋渦內(nèi)的最大速度U和兩個(gè)對(duì)稱(chēng)旋渦的渦核間距d不斷增加，在網(wǎng)格為800×400時(shí)所得結(jié)果與實(shí)驗(yàn)結(jié)果[5]最為接近。

2.2 圓柱繞流的特性

圖2為網(wǎng)格劃分為200×100，Re=1時(shí)的速度等值線圖和流線圖，由圖可知圓柱兩側(cè)的流線前后上下對(duì)稱(chēng)，流體繞圓柱流動(dòng)沒(méi)有形成旋渦，仍處于層流狀態(tài)，與圓柱繞流流場(chǎng)性質(zhì)相吻合。

endprint

以下結(jié)果為網(wǎng)格劃分800×400。如圖3所示為不同雷諾數(shù)下的圓柱附近的流線圖，Re>10后繞圓柱分離后的流體在圓柱后形成一對(duì)固定不動(dòng)的對(duì)稱(chēng)旋渦并且不斷被拉長(zhǎng)，隨著Re的增大，旋渦逐漸失去對(duì)稱(chēng)性，對(duì)稱(chēng)旋渦消失旋渦開(kāi)始從兩側(cè)不斷地輪流脫落，由Re>100時(shí)，可以看出圓柱體的下游的流場(chǎng)不再具有規(guī)則性，流場(chǎng)變得較復(fù)雜，由此得出的計(jì)算結(jié)果與圓柱繞流的流場(chǎng)特性相吻合。

表2為Re=25時(shí)，旋渦分離的物理參數(shù)和文獻(xiàn)的對(duì)比，與所得結(jié)果吻合得很好，說(shuō)明了該方法求解圓柱繞流流場(chǎng)的有效性和正確性。

3 結(jié)論

通過(guò)利用浸入邊界法對(duì)單圓柱繞流進(jìn)行數(shù)值模擬計(jì)算，得出如下結(jié)論：①使用浸入邊界法模擬圓柱繞流，不同網(wǎng)格劃分下所得結(jié)果都符合流場(chǎng)的一般性質(zhì)，只是精確度不同，網(wǎng)格劃分越密，精度越高，在格劃分間距d=0.025D時(shí)能很好的滿足流場(chǎng)求解的需求。②使用浸入邊界法編程來(lái)實(shí)現(xiàn)圓柱繞流流場(chǎng)的數(shù)值模擬，所得結(jié)果與圓柱繞流的流場(chǎng)性質(zhì)相吻合，該方法能有效、正確的模擬圓柱繞流流場(chǎng)。

參考文獻(xiàn)：

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

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

[3]邢景棠，周盛，崔爾杰.流固耦合力學(xué)概述[J].力學(xué)進(jìn)展， 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]姜知亨.浸入邊界法及在大渦模擬中的應(yīng)用[D].清華大學(xué).endprint