相變傳熱與流體流動數(shù)值分析第114講課件

1、相變傳熱與流體流動數(shù)值分析(第11-14講)格子Boltzmann方法Lattice Boltzmann Method主要內(nèi)容 14.1 啟發(fā)式邊界條件 14.2 動力學(xué)邊界條件 14.3 外推邊界條件 14.4 曲面邊界條件概述 14.5 LBM模擬的基本步驟本章重點；邊界條件處理 上章回憶；理論推導(dǎo)邊界條件處理如何進行LBM模擬單相模型(D2Q9模型)多相和多組分模型（SC模型）14. LBM邊界處理每個時步之后，內(nèi)部流場節(jié)點上的分布函數(shù)均已經(jīng)獲得，但是邊界節(jié)點上的部分分布函數(shù)是未知的。需要根據(jù)已知的宏觀邊界條件確定出邊界節(jié)點上相應(yīng)的分布函數(shù)的取值。根據(jù)邊界處理格式的特性，對邊界處理的方式

2、進行分類，主要分為啟發(fā)式格式、動力學(xué)格式、外推格式以及其他復(fù)雜邊界處理格式。原理： 根據(jù)邊界上諸如周期性、對稱性、充分發(fā)展等宏觀物理學(xué)特性，通過微觀粒子的的運動規(guī)則直接確定邊界節(jié)點上的未知分布函數(shù)。優(yōu)點： 不需要較復(fù)雜的數(shù)學(xué)推導(dǎo)和公式求解類型： 周期性邊界、對稱邊界、充分發(fā)展以及用于固壁邊界的反彈格式、鏡面反彈格式、反彈與鏡面反彈混合格式等14.1 啟發(fā)式格式如果流場在空間呈現(xiàn)周期性變化或在某個方向無窮大，常常將周期性單元取出作為模擬區(qū)域，并在相應(yīng)邊界上采用周期性邊界。其邊界處理格式是指，當流體粒子從一側(cè)邊界離開流場時，在下一個時步就會從流場的另側(cè)邊界進入流場。周期性邊界能嚴格保證整個系統(tǒng)的質(zhì)

3、量和動量守恒。14.1.1 周期性邊界處理格式周期格式可以表示為：14.1.1 周期性邊界處理格式14.1.2 對稱邊界處理格式由于對稱性問題，為了節(jié)省計算資源，可以取物理模型的一半作為模擬的區(qū)域并在對稱軸上采取對稱邊界處理當流體在通道內(nèi)流動達到充分發(fā)展后，密度與速度等物理量在主流方向上不再發(fā)生變化，即他們的空間導(dǎo)數(shù)為0其處理格式可以表示為：這是格子Boltzmann方法中處理充分發(fā)展邊界最簡單和最常用的方法14.1.3 充分發(fā)展邊界處理格式處理靜止無滑移壁面的一種常用格式，包括標準反彈格式、半步長反彈格式及修正反彈格式標準反彈格式：對于靜止固體邊界，常用的處理方法就是對邊界上的粒子進行彈回處

4、理，即假設(shè)粒子與壁面碰撞后沿粒子原方向的反方向逆轉(zhuǎn)14.1.4 反彈格式標準反彈格式： 操作簡單，能夠嚴格保證系統(tǒng)的質(zhì)量和動量守恒； 一階精度，降低了格子玻爾茲曼方法的精度； LBM格子波爾茲曼方法在內(nèi)節(jié)點上具有二階精度 -降低方法精度 修正反彈格式或者半步長反彈格式-二階精度與標準反彈格式的最大區(qū)別：在碰撞過程中，邊界節(jié)點與流體區(qū)域的節(jié)點一樣，所有的離散速度方向的粒子均參與碰撞。得到的分布函數(shù)在進入下一次的碰撞和遷移的過程，這一格式稱為修正反彈格式。修正反彈格式： 14.1.4 反彈格式 半步長反彈格式與標準反彈格式的反彈原理相同。不同的是：固體邊界位于兩個格點中間。 半步長反彈格式： 物理

5、意義： 在 格子節(jié)點處碰撞，速度指向壁面的粒子經(jīng)過 時間到達固體壁面，然后與壁面碰撞后反彈，速度逆轉(zhuǎn)，再經(jīng)過 時間沿原路返回到 格子節(jié)點處。反彈邊界的優(yōu)點在于操作簡單，對于處理復(fù)雜的不規(guī)則邊界（如多孔介質(zhì)）具有很大的優(yōu)勢;14.1.4 反彈格式對于光滑壁面，如果壁面對于流體沒有摩擦作用，即壁面與流體之間沒有動量交換，則通常采用鏡面反射格式來實現(xiàn)自由滑移邊界。流體節(jié)點(i-1,2)入射粒子 遷移到壁面節(jié)點(i,1)后，如果沿原路返回即 方向，即為反彈處理；如果以壁面法線方向 反射，則為鏡面反射處理。 14.1.5 鏡面反射格式壁面法向速度為零；切向速度分量不為零!適合既不能用簡單的反彈格式處理，

6、也不能用鏡面反射格式的情況，如微通道中的氣體流動；將反彈和鏡面反射相綜合起來，準確地實現(xiàn)真實的氣體與固體之間的相互作用。因此需要定義一個彈回系數(shù) （ ），用來表述粒子在和壁面作用時沿原路彈回所占的比例。14.1.6 反彈與鏡面反射混合格式 當 時，表示純反彈處理當 時，表示純鏡面反射處理當 時，表示理想的漫反射 越小，壁面滑移速度越大；對任意 ，壁面法向速度為零；14.1.6 反彈與鏡面反射混合格式 動力學(xué)格式主要利用邊界上宏觀物理量的定義，直接求解邊界點上未知分布函數(shù)的方程組以獲得邊界節(jié)點上待定的分布函數(shù)，動力學(xué)格式主要包括。 Nobel格式 非平衡反彈格式 反滑移格式 質(zhì)量修正格式14.2

7、 動力學(xué)格式 14.2.1 Nobel格式通常情況下，速度或壓力其中一個已知，以速度邊界條件為例：未知量：3個方程1995年，Nobel等人使用內(nèi)能作為補充條件方程封閉，聯(lián)立求解1995年，Inamuro等提出；邊界節(jié)點上的未知分布函數(shù)可以用一個新的平衡態(tài)分布函數(shù)來替代為了克服壁面滑移，需要在這個平衡態(tài)分布函數(shù)中加入一個反滑移速度進行修正，該反滑移速度可以抵消掉采用反彈格式處理無滑移邊界時所產(chǎn)生的滑移速度，從而使壁面上的流體速度與避免速度相等。設(shè)未知分布函數(shù)形式如下：其中 ， ， ，都是未知參數(shù)， 表示壁面的速度。14.2.2 反滑移格式通用性差1997年，Zou與He給出了另外一種補充條件：

8、假設(shè)分布函數(shù)的非平衡部分在垂直于邊界的方向上仍然滿足反彈格式速度邊界/壓力邊界14.2.3 非平衡態(tài)反彈格式流場內(nèi)部對于速度邊界條件：已知量：未知量：四個未知量，需要四個方程The macroscopic density formula is one equation:14.2.3 非平衡態(tài)反彈格式The macroscopic velocity formula gives two equations:x-direction:y-direction:Components of ea are all unit vectorsAssuming ux = 014.2.3 非平衡態(tài)反彈格式假設(shè)分布函數(shù)

9、的非平衡部分在垂直于邊界的方向上仍然滿足反彈格式14.2.3 非平衡態(tài)反彈格式第四個方程，方程封閉，聯(lián)立求解Two equations have the directional density unknowns f4, f7 and f8 in common, so rewrite them with those variables on the left hand side:14.2.3 非平衡態(tài)反彈格式Equating them gives: 14.2.3 非平衡態(tài)反彈格式Solving for r:Solving the bounceback equation for f4: Solvi

10、ng for f7:14.2.3 非平衡態(tài)反彈格式Solving for f8: 14.2.3 非平衡態(tài)反彈格式/ Zou and He velocity BCs on north side. for( i=0; ini; i+) fi = ftempnj-1i; rho0 = ( fi0 + fi1 + fi3 + 2.*( fi2 + fi5 + fi6) / ( 1. + uy0); ru = rho0*uy0; fi4 = fi2 - (2./3.)*ru; fi7 = fi5 - (1./6.)*ru + (1./2.)*( fi1-fi3); fi8 = fi6 - (1./6.)

11、*ru + (1./2.)*( fi3-fi1); 14.2.3 非平衡態(tài)反彈格式壓力（密度）邊界Dirichlet boundary conditions constrain the pressure/density at the boundariesSolution is closely related to that for velocity boundaries A density r0 is specified and velocity is computedSpecifying density is equivalent to specifying pressure since t

12、here is an equation of state (EOS) relating them directlyFor single component D2Q9 model, the relationship is simply P = RTr with RT = 1/3. More complex EOS applies to single component multiphase modelsWe assume that velocity tangent to the boundary is zero and solve for the component of velocity no

13、rmal to the boundary. 14.2.3 非平衡態(tài)反彈格式Assume that velocity tangent to the boundary is zero and solve for the component of velocity normal to the boundary Need to solve for v, f4, f7 and f8 壓力（密度）邊界14.2.3 非平衡態(tài)反彈格式壓力（密度）邊界14.2.3 非平衡態(tài)反彈格式壓力（密度）邊界14.2.3 非平衡態(tài)反彈格式/ Zou and He pressure boundary on north sid

14、e. for( i=0; ini; i+) fi = ftempnj-1i; uy0 = -1. + ( fi0 + fi1 + fi3 + 2.*( fi2 + fi5 + fi6) / rho0; ru = rho0*uy0; fi4 = fi2 - (2./3.)*ru; fi7 = fi5 - (1./6.)*ru + (1./2.)*( fi1-fi3); fi8 = fi6 - (1./6.)*ru + (1./2.)*( fi3-fi1); 壓力（密度）邊界14.2.3 非平衡態(tài)反彈格式對于充分發(fā)展流動而言，如果邊界條件不能使流體在整個流動區(qū)域內(nèi)滿足質(zhì)量守恒要求，則數(shù)值結(jié)果的精度

15、就會降低，而且數(shù)值計算的收斂速度以及穩(wěn)定性都會受到影響；如果在處理充分發(fā)展邊界條件時，采取一定的方法，使得流體在整個流動區(qū)域內(nèi)的總體質(zhì)量守恒要求得到滿足，那么數(shù)值結(jié)果的精度以及計算過程的收斂狀況就會有較大程度的提高；何雅玲等基于此思路，提出了格子Bolzmann方法處理充分發(fā)展流動的質(zhì)量修正格式。所謂質(zhì)量修正，就是通過計算進口和出口的質(zhì)量流率，確定一個質(zhì)量修正系數(shù)，以此系數(shù)修正出口速度，并用修正后的速度更新邊界節(jié)點上的未知粒子分布函數(shù)。14.2.4 質(zhì)量修正格式具體做法如下：假設(shè)出口截面上各點速度x方向分量的相對變化率為常數(shù)，即然后，考慮整體質(zhì)量守恒，也即進出口的質(zhì)量流率相等，得質(zhì)量修正系數(shù)為

16、：利用質(zhì)量修正系數(shù)即可得到出口速度分布，利用該速度計算出口邊界的粒子分布函數(shù)。14.2.4 質(zhì)量修正格式啟發(fā)式格式和動力學(xué)格式的局限性依賴使用的格子Boltzmann模型， 通用性差；需要對邊界處的一些物理性質(zhì)（如速度、密度、壓力等）作假設(shè)，因而很難推廣到更一般的邊界處理上。對一些非常規(guī)邊界（包含梯度信息的邊界）存在誤差借鑒傳統(tǒng)計算流體力學(xué)方法中的邊界處理方法來構(gòu)造格子Bolzman方法邊界處理格式；外推格式非平衡外推格式14.3 外推格式1996年Chen等人提出；假定在物理邊界外有一層虛擬邊界，物理邊界和虛擬邊界格點上的分布函數(shù)值都按照標準的演化步驟進行演化，但虛擬邊界格點上指向流場的分布

17、函數(shù)值使用外推方法確定。14.3.1 外推格式流體層物理邊界虛擬邊界在碰撞過程中在物理邊界上使用給定的邊界條件（速度或壓力）來計算平衡態(tài)分布函數(shù)。二階精度優(yōu)點普適性好，可以根據(jù)這種方法容易地設(shè)計出一般邊界條件的格式；計算量也比較小，容易實現(xiàn)。在確定未知分布函數(shù)時，不需要額外的假設(shè)條件，這是外推格式的技巧所在；具有二階精度；缺點 數(shù)值穩(wěn)定性較差（在外推格式中，因為有負的系數(shù)出現(xiàn)，所以數(shù)值誤差有可能隨時間增加，從而導(dǎo)致數(shù)值不穩(wěn)定的發(fā)生）14.3.1 外推格式受外推方法和非平衡態(tài)反彈方法的啟發(fā)，吸收兩種方法的優(yōu)點，2002年Guo等人提出了一種新的邊界處理方法，即非平衡態(tài)外推方法。其基本思想是，將邊

18、界節(jié)點上的分布函數(shù)分解為平衡態(tài)和非平衡態(tài)兩部分，并根據(jù)具體的邊界條件定義新的平衡態(tài)分布來近似平衡態(tài)部分，而非平衡態(tài)部分則使用非平衡態(tài)外推（一階精度）確定，所得的分布函數(shù)的整體逼近精度為二階精度。14.3.2 非平衡外推格式流場內(nèi)部節(jié)點： G 邊界節(jié)點： O14.3.2 非平衡外推格式流體內(nèi)部物理邊界流場外部對于非平衡態(tài)部分，用鄰近的流體格子節(jié)點G出的非平衡態(tài)部分近似得到，即14.3.2 非平衡外推格式對于速度邊界：u已知，密度未知空間精度是二階的，時間精度也是二階的，這與格子Boltzmann方法的整體精度是一致的。由于低階外推格式的數(shù)值穩(wěn)定性一般要優(yōu)于高階外推格式的數(shù)值穩(wěn)定性，所以非平衡態(tài)外

19、推格式的數(shù)值穩(wěn)定性要好于根據(jù)線性外推得到的分布函數(shù)外推格式的數(shù)值穩(wěn)定性。同時，與外推格式一樣，非平衡態(tài)外推格式也具有適用范圍廣、計算簡單，容易實現(xiàn)的優(yōu)點。14.3.2 非平衡外推格式當所要處理的物理區(qū)域并不具有規(guī)則的幾何形狀時適體網(wǎng)格非結(jié)構(gòu)化網(wǎng)格在結(jié)構(gòu)化的直角正交網(wǎng)格，并在適當?shù)奈恢貌捎秒A梯逼近或者差值處理，以保證滿足物理邊界上的條件。 階梯逼近可與諸如反彈格式等前面提到的多種邊界處理格式配合使用，實施簡單，但是計算精度較低，常用的處理復(fù)雜邊界處理格式為Filippova與Hanel格式、Bouzidi格式、Lallemand與Luo格式、Guo格式等14.4 復(fù)雜邊界處理格式1998年， F

20、ilippova與Hanel提出構(gòu)建虛擬平衡態(tài)分布函數(shù)，并用其進行線性插值的方法來計算由壁面彈回流體的分布函數(shù)值。14.4.1 Filippova與Hanel格式定義變量q，用來表示流體節(jié)點中離壁面最近的點到實際物理邊界的相對距離，即式中， 為實際物理邊界所在的位置?？紤]到算法的穩(wěn)定性， Filippova與Hanel提出用以下關(guān)系式確定 和 當 時，當 時，14.4.1 Filippova與Hanel格式2001年，Bouzidi等提出了一種結(jié)合反彈格式和空間插值的邊界處理方法，用以處理任意曲率的復(fù)雜幾何邊界。該格式的基本思想是利用在彈回方向上用空間插值來獲得由邊界進入流體的分布函數(shù)值?？紤]

21、如下圖一般的二維情況 ，則線性插值格式給定如下：當 時當 時14.4.2 Bouzidi格式14.4.2 Bouzidi格式在Bouzidi格式中，差值公式包含了在不同時間層次上的分布函數(shù)。2003年, Lallemand與Luo針對運動邊界對Bouzidi格式進行了改進，在引用差值公式時做了變化，把所有分布函數(shù)統(tǒng)一到同一個時間層次上，同時考慮了碰撞和遷移過程，提出了Lallemand與Luo格式。2002年，Guo等提出了一種結(jié)合非平衡外推格式和空間差值曲線邊界處理方法，其基本思想是，將邊界節(jié)點的分布函數(shù)分為平衡態(tài)和非平衡態(tài)兩部分，平衡態(tài)部分由一個虛擬的平衡態(tài)分布函數(shù)代替，非平衡態(tài)部分則由相

22、鄰流體節(jié)點的非平衡部分插值獲得，提出了Guo格式。14.4.3 Lallemand與Luo格式、Guo格式14.5 LBM模擬的基本步驟Define velocity vectors012345678D2Q9 e1e2e3e4e5e6e7e8%define esex(0)= 0; ey(0)= 0ex(1)= 1; ey(1)= 0ex(2)= 0; ey(2)= 1ex(3)=-1; ey(3)= 0ex(4)= 0; ey(4)=-1ex(5)= 1; ey(5)= 1ex(6)=-1; ey(6)= 1ex(7)=-1; ey(7)=-1ex(8)= 1; ey(8)=-114.5 LB

23、M模擬的基本步驟Problem definitionLY=10LX=20tau = 1( relaxation time)g=0.00001(gravity acceleration)%set solid nodesis_solid_node=zeros(LY,LX)for i=1:LX is_solid_node(1,i)=1 is_solid_node(LY,i)=1endLXLY% if is_interior_solid_node(j,i)14.5 LBM模擬的基本步驟Initialize density and fs (assuming zero velocity)%define i

24、nitial density and fsrho=ones(LY,LX);f(:,:,1) = (4./9. ).*rho; f(:,:,2) = (1./9. ).*rho;f(:,:,3) = (1./9. ).*rho;f(:,:,4) = (1./9. ).*rho;f(:,:,5) = (1./9. ).*rho;f(:,:,6) = (1./36.).*rho;f(:,:,7) = (1./36.).*rho;f(:,:,8) = (1./36.).*rho;f(:,:,9) = (1./36.).*rho;14.5 LBM模擬的基本步驟/ Computing macroscopi

25、c density, rho, and velocity, u=(ux,uy). for( j=0; jLY; j+) for( i=0; iLX; i+) u_xji = 0.0; u_yji = 0.0; rhoji = 0.0; if( !is_solid_nodeji) for( a=0; a9; a+) rhoji += fjia; u_xji += exa*fjia; u_yji += eya*fjia; u_xji /= rhoji; u_yji /= rhoji; 14.5 LBM模擬的基本步驟14.5 LBM模擬的基本步驟Periodic BoundariesOn bound

26、ary, neighboring point is on opposite boundary ip = ( i0 )?( i-1):( LX-1); jp = ( j0 )?( j-1):( LY-1);/ Compute the equilibrium distribution function, feq. f1=3.; f2=9./2.; f3=3./2.; for( j=0; jLY; j+) for( i=0; iLX; i+) if( !is_solid_nodeji) rt0 = (4./9. )*rhoij; rt1 = (1./9. )*rhoij; rt2 = (1./36.

27、)*rhoij; ueqxij = uxij; /Can add forcing here; see MATLAB code ueqyij = uyij; uxsq = ueqxij * ueqxij; uysq = ueqyij * ueqyij; uxuy5 = ueqxij + ueqyij; uxuy6 = -ueqxij + ueqyij; uxuy7 = -ueqxij + -ueqyij; uxuy8 = ueqxij + -ueqyij; usq = uxsq + uysq; 14.5 LBM模擬的基本步驟 feqij0 = rt0*( 1. - f3*usq); feqij1

28、 = rt1*( 1. + f1*ueqxij + f2*uxsq - f3*usq); feqij2 = rt1*( 1. + f1*ueqyij + f2*uysq - f3*usq); feqij3 = rt1*( 1. - f1*ueqxij + f2*uxsq - f3*usq); feqij4 = rt1*( 1. - f1*ueqyij + f2*uysq - f3*usq); feqij5 = rt2*( 1. + f1*uxuy5 + f2*uxuy5*uxuy5 - f3*usq); feqij6 = rt2*( 1. + f1*uxuy6 + f2*uxuy6*uxuy6 - f3*usq); feqij7 = rt2*( 1. + f1*uxuy7 + f2*uxuy7*uxuy7 - f3*usq); feqij8 = rt2*( 1. + f1*uxuy8 + f2*uxuy8*uxuy8 - f3*usq); 14.5 LBM