不可壓縮繞流的并行格子Boltzmann模擬

張建影, 閆廣武

(1. 長春工業(yè)大學 基礎科學學院, 長春 130012； 2. 吉林大學 數(shù)學學院, 長春 130012)

0 引 言

格子Boltzmann方法(LBM)是一種模擬不可壓縮黏性流的工具, 已得到廣泛關注[1-2]. 與傳統(tǒng)計算流體動力學(CFD)的方法不同, LBM從細觀尺度上的Maxwell方程出發(fā), 通過特殊的離散方式給出粒子束的演化方程, 進一步得到宏觀物理量. 因此, 最初的對流擴散方程, 如Navier-Stokes方程, 被一個格子Boltzmann方程取代. 由于標準格子Boltzmann方法具有如下優(yōu)點： 1) 由狀態(tài)方程給出壓力, 不需要求解壓力方程[3]； 2) 公式簡單易于實現(xiàn)編程； 3) 算法容易實現(xiàn)并行[4]； 4) 復雜邊界處理較方便[5]. 因此, LBM在流體力學和非線性偏微分方程等領域應用廣泛[6]. 但當流動區(qū)域較大或需要處理流動細節(jié)的計算問題時, 計算網格數(shù)量的要求已超出普通計算機的能力. 在這種情況下, 需要使用并行模式的格子Boltzmann 算法[7-8], 利用LBM與其他算法的結合實現(xiàn)并行計算[9-12].

本文提出一個簡單的并行模式, 通過在流動區(qū)域內劃分若干子塊以及與周圍子塊之間的重疊邊界, 每個子塊的計算是并行的[13], 利用MPI(message passing interface)并行平臺實現(xiàn)全場的并行計算, 進而達到增加計算網格數(shù)量的目的.

1 格子Boltzmann 模型

在標準的二維網格空間中, 定義fα(x,t)為位置x、 在時刻t沿方向α的粒子分布函數(shù). 流場中的密度和動量定義為

(5)

平衡態(tài)分布函數(shù)為

(6)

(7)

其中τ為松弛因子.

引入Knudsen數(shù)ε, 其定義為平均自由程l與特征長度L之比, 即

(8)

假設ε與格子尺度Δx同量級, 并令Δt=ε, 則Δx=c·Δt=c·ε, 利用Chapman-Enskog展開及時間多尺度展開, 可得不同時間尺度下的系列偏微分方程[14]：

將式(9)+ε×式(10)并求和得

在不可壓縮假設和低速度假設下, 可得不可壓縮流動的Navier-Stokes方程:

其中

(15)

2 管道內繞流的并行格子Boltzmann算法

圖1 矩形計算區(qū)域子塊劃分Fig.1 Sub-block partition of rectangular computing area

將矩形區(qū)域用D2Q9網格鋪滿, 記水平方向為x軸m個格點, 豎直方向為y軸n個格點. 為了采用MPI并行技術, 將計算區(qū)域分解成如圖1所示的若干子塊, 并要求一個子塊與相鄰子塊之間有重合的兩列網格, 第一個子塊和最后一個子塊只有一個相鄰的子塊, 子塊的數(shù)量根據全場y方向網格數(shù)決定. 在計算中, 這些子塊同時計算, 除第一個和最后一個子塊外, 每塊數(shù)據都需要調用5部分, 分別為左右兩列以及該塊的內部區(qū)域[13].

計算圖1中第一個子塊B1時, 需要調用的數(shù)據為： 列M1N1,M2N2,M5N5,M6N6和塊M3N3N4M4. 時刻(t+Δt)完成數(shù)據M2N2N5M5的計算, 列M1N1和M5N5作為該子塊的邊界使用. 考慮到子塊之間的對接, 采用雙邊界M5N5,M6N6, 即在計算第一個子塊時M6N6是邊界, 但在計算第二個子塊B2時M6N6為內部, 此時邊界為M5N5. 第二個子塊B2的計算需要調用列M5N5,M6N6,M9N9,M10N10和塊M7N7N8M8. 列M5N5,M10N10為邊界.

2.1 方柱繞流

流場為600×100網格, 運算分為3個子塊, 每塊為200×100. 用區(qū)域[92,108]×[42,58]表示邊長16個網格的方柱, 中心位于(100,50)處, 上下邊界采用固體邊界. 圖2為Re=100, 時間步分別為t=12 000,18 000,24 000,30 000時的渦線； 圖3為Re=300, 時間步分別為t=12 000,18 000,24 000,30 000時的渦線.

圖2 Re=100時方柱繞流渦量等值線Fig.2 Vorticity isolines of flow around a square cylinder at Re=100

由圖2可見, 當Re=100時, 方柱后面開始出現(xiàn)Karmann渦街, 流動變?yōu)椴欢ǔ? 受其影響, 方柱所受的升力出現(xiàn)周期振蕩, 方柱的阻力也出現(xiàn)振蕩. 由圖3可見, 當Re=300時, 方柱后面的渦旋區(qū)域變大, 流動很快變成不定常. 這些結果與已有的數(shù)值結果[15-16]相符.

圖3 Re=300時方柱繞流渦量等值線Fig.3 Vorticity isolines of flow around a square cylinder at Re=300

阻力系數(shù)與升力系數(shù)分別定義為

(16)

其中:FD為阻力;FL為升力;ρ∞為來流密度;U∞為來流速度;d為方柱特征寬度.

圖4為Re=100時阻力系數(shù)與升力系數(shù)隨時間步的變化曲線. 由圖4可見, 當Re=100時, 阻力系數(shù)約為1.31, 而升力系數(shù)約為0, 與文獻[15-16]中的阻力系數(shù)和升力系數(shù)相同, 當t>20 000時出現(xiàn)振蕩. 圖5為Re=300時阻力系數(shù)與升力系數(shù)隨時間步的變化曲線. 此時方柱后面Karmann渦街的區(qū)域變大, 由圖5可見， 當t>13 000時, 升力與阻力的振蕩振幅變大, 與文獻[15-16]的結果吻合.

圖4 Re=100時方柱繞流阻力系數(shù)和 升力系數(shù)與時間步的關系曲線Fig.4 Relationship curves between drag coefficient and lift coefficient of flow around a square cylinder and time step at Re=100

圖5 Re=300時方柱繞流阻力系數(shù)和 升力系數(shù)與時間步的關系曲線 Fig.5 Relationship curves between drag coefficient and lift coefficient of flow around a square cylinder and time step at Re=300

2.2 圓柱繞流

流場為600×100網格, 運算分為3個子塊, 每塊為200×100. 圓柱中心位于(100,50), 半徑為12.5, 上下邊界采用固體邊界. 圖6為Re=300, 時間步分別為t= 20 000,30 000,40 000,50 000時的渦線. 由圖6可見, Karmann渦街中渦旋較密集. 圖7為Re=300時圓柱繞流的阻力系數(shù)與升力系數(shù)隨時間步的變化曲線. 由圖7可見, 升力和阻力均出現(xiàn)振蕩, 與文獻[17]中的LBM結果相符.

2.3 NACA0012機翼繞流

流場為600×300網格, 運算分為3個子塊, 每塊為200×300. 機翼前緣位于(150,150), 翼弦長為75個網格, 攻角12°, 上下邊界采用Neumann邊界. 圖8為Re=1 000時圓柱繞流的阻力系數(shù)與升力系數(shù)隨時間步的變化曲線. 由圖8可見, 升力和阻力均出現(xiàn)振蕩, 與文獻[18-19]的結果相符. 圖9為Re=1 000, 時間步分別為t=6 000,8 000,10 000,40 000時的渦線. 由圖9可見, 當t>6 000時機翼后面出現(xiàn)渦旋, 當t>10 000時出現(xiàn)較穩(wěn)定的Karmann渦街.

圖6 Re=300時圓柱繞流渦量等值線Fig.6 Vorticity isolines of flow around a circular cylinder at Re=300

圖7 Re=300時圓柱繞流阻力系數(shù)和 升力系數(shù)與時間步的關系曲線Fig.7 Relationship curves between drag coefficient and lift coefficient of flow around a circular cylinder and time step at Re=300

圖8 Re=1 000時圓柱繞流阻力系數(shù)和 升力系數(shù)與時間步的關系曲線 Fig.8 Relationship curves between drag coefficient and lift coefficient of flow around a circular cylinder and time step at Re=1 000

圖9 Re=1 000, 攻角α=12°時渦量等值線Fig.9 Vorticity isolines at Re=1 000 and attack angle α=12°

綜上, 本文應用并行模式的格子Boltzmann算法模擬了方柱、 圓柱以及NACA0012機翼繞流問題, 并給出了不同雷諾數(shù)下不同時間的渦線圖, 由這些渦線圖可觀察到被繞物體后面形成的Karmann渦街, 由升力阻力系數(shù)隨時間的變化曲線可觀察到升力阻力均出現(xiàn)振蕩, 這些結果都與經典的數(shù)值結果相符， 顯示了本文方法的有效性.

[1] CHEN Shiyi, Doolen G D. Lattice Boltzmann Method for Fluid Flows [J]. Annual Review of Fluid Mechanics, 1998, 30: 329-364.

[2] Aidun C K, Clausen J R. Lattice-Boltzmann Method for Complex Flows [J]. Annual Review of Fluid Mechanics, 2010, 42: 439-472.

[3] HE Xiaoyi, CHEN Shiyi, Doolen G D. A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit [J]. J Comput Phys, 1998, 146(1): 282-300.

[4] Schmieschek S, Shamardin L, Frijters S, et al. LB3D: A Parallel Implementation of the Lattice-Boltzmann Method for Simulation of Interacting Amphiphilic Fluids [J]. Computer Physics Communications, 2017, 217: 149-161.

[5] Rosis A, de. Ground-Induced Lift Enhancement in a Tandem of Symmetric Flapping Wings: Lattice Boltzmann-Immersed Boundary Simulations [J]. Comput Struct, 2015, 153: 230-238.

[6] ZHANG Jianying, YAN Guangwu. Lattice Boltzmann Model for Complex Ginzburg-Landau Equation in Curvilinear Coordinates [J]. Computers and Mathematics with Applications, 2015, 70(12): 2904-2919.

[7] Clausen J R, Reasor D A, Jr, Aidun C K. Parallel Performance of a Lattice-Boltzmann/Finite Element Cellular Blood Flow Solver on the IBM Blue Gene/P Architecture [J]. Comput Phys Commun, 2010, 181(6): 1013-1020.

[8] GAO Jinfang, XING Huilin, Rudolph V, et al. Parallel Lattice Boltzmann Computing and Applications in Core Sample Feature Evaluation [J]. Transport in Porous Media, 2015, 107(1): 65-77.

[9] Campos J O, Oliveira R S, Santos R W, dos， et al. Lattice Boltzmann Method for Parallel Simulations of Cardiac Electrophysiology Using GPUs [J]. Journal of Computational and Applied Mathematics, 2016, 295: 70-82.

[10] Zudrop J, Masilamani K, Roller S, et al. A Robust Lattice Boltzmann Method for Parallel Simulations of Multicomponent Flows in Complex Geometries [J]. Computers and Fluids, 2017, 153: 20-33.

[11] Randles A, Kaxiras E. Parallel in Time Approximation of the Lattice Boltzmann Method for Laminar Flows [J]. Journal of Computational Physics, 2014, 270: 577-586.

[12] SU Yan, Ng T, Davidson J H. A Parallel Non-dimensional Lattice Boltzmann Method for Fluid Flow and Heat Transfer with Solid-Liquid Phase Change [J]. International Journal of Heat and Mass Transfer, 2017, 106: 503-517.

[13] 閆廣武. 二維直角彎道黏性流動的格子氣體仿真 [D]. 長春: 吉林大學, 1988. (YAN Guangwu. Lattice Gas Simulation of 2-D Viscous Flow in a Right-Angled Bend [D]. Changchun: Jilin University, 1988.)

[14] 閆廣武. 二維剪切流的格子Boltzmann模擬 [J]. 水動力學研究與進展, 2003, 18(5): 521-525. (YAN Guangwu. Lattice Boltzmann Simulation of the Two-Dimensional Shear Flows [J]. Journal of Hydrodynamics, 2003, 18(5): 521-525.)

[15] Breuer M, Bernsdorf J, Zeiser T, et al. Accurate Computations of the Laminar Flow Past a Square Cylinder Based on Two Different Methods: Lattice-Boltzmann and Finite-Volume [J]. International Journal of Heat and Fluid Flow, 2000, 21: 186-196.

[16] Jafari S, Salmanzadeh M, Rahnama M, et al. Investigation of Particle Dispersion and Deposition in a Channel with a Square Cylinder Obstruction Using the Lattice Boltzmann Method [J]. Journal of Aerosol Science, 2010, 41(2): 198-206.

[17] Shu C, Liu N, Chew Y T. A Novel Immersed Boundary Velocity Correction-Lattice Boltzmann Method and Its Application to Simulate Flow Past a Circular Cylinder [J]. Journal of Computational Physics, 2007, 226: 1607-1622.

[18] 邱全輝, 王健平. 不可壓縮機翼繞流的有限譜法計算 [J]. 計算力學學報, 2005, 22(5): 518-523. (QIU Quanhui, WANG Jianping. Finite Spectral Computation of Incompressible Flows around Airfoil [J]. Chinese Journal of Computational Mechanics, 2005, 22(5): 518-523.)

[19] Ashraf M A, Young J, Lai J C S. Reynolds Number, Thickness and Camber Effects on Flapping Airfoil Propulsion [J]. Journal of Fluids and Structures, 2017, 27(2): 145-160.