基于格子Boltzmann 方法的釘扎螺旋波反饋控制*

賴瑤瑤 陳鑫夢 柴振華 施保昌?

1) (華中科技大學數學與統計學院,武漢 430074)

2) (華中科技大學,工程建模與科學計算湖北重點實驗室,武漢 430074)

3) (華中科技大學,數學與應用學科交叉創新研究院,武漢 430074)

螺旋波是心室跳動過速和纖維性顫動的根源,釘扎螺旋波相對于自由螺旋波來說更難消除.本文采用格子Boltzmann 方法求解,以FitzHugh-Nagumo 模型為對象,研究了使用反饋控制法消除釘扎螺旋波.數值結果表明,無論釘扎螺旋波釘在圓形障礙物還是矩形障礙物上,反饋控制法對其都具有很好的控制作用.此外,通過數值模擬系統研究了可激性系數、反饋控制信號幅度、記錄反饋信號時間和障礙物的大小對釘扎螺旋波的控制情況.研究表明,釘扎螺旋波消除有三種情況.首先,反饋控制信號幅度和可激性系數與釘扎螺旋波消除所需的時間有關,反饋控制信號幅度越大或可激性系數越小,釘扎螺旋波消除越快.其次,障礙物大小和可激性系數影響著能成功消除釘扎螺旋波下記錄反饋信號時間與加入反饋控制時間之間對應的時間間隔.最后,在保持加入反饋控制時間不變的情況下,記錄反饋信號時間影響著能成功消除釘扎螺旋波所需的最小反饋控制信號幅度.

1 引言

在許多非線性反應擴散系統的斑圖中,螺旋波動力系統一直受到廣泛關注.螺旋波廣泛存在于自然界中,例如Belousov-Zhabotinsky (BZ)反應[1]、催化表面過程[2]和心臟系統[3]等.其中,BZ 反應是最早發現螺旋波現象的自然反應系統,對BZ 反應中螺旋波的研究最早開始于20 世紀70 年代[4].在心臟系統中,螺旋電波與心動過速和危及生命的纖顫有關,它的存在會危及人體生命[5].然而,如果心臟中存在缺陷區域,螺旋波可能會被釘住,其核心將停留在缺陷區域,使它們停留的時間更長[6,7].當自由螺旋波固定在缺陷上時,就產生了所謂的釘扎螺旋波,從而導致一類生理性心律失常.因此,尋求有效的釘扎螺旋波控制方法是十分重要的.對于釘扎螺旋波,許多控制方法被研究人員提出.例如Fu 等[8]通過外部局部刺激產生靶波,然后靶波與釘住的螺旋波相互作用來消除釘扎螺旋波.Ponboonjaroenchai 等[9]提出了在給定位置進行周期性刺激產生波列控制釘扎螺旋波,研究了波列在大小和方向不同的不可激矩形障礙物上螺旋波的終止問題.Chen 等[10]研究了同步復極化能有效消除螺旋波和時空混沌.Yuan 等[11]通過邊界附近的周期性擾動方法成功地解除和消除釘扎螺旋波,并且該方法對于大障礙物同樣適用.

反饋法是一種控制自由螺旋波的經典方法,它具有持續時間短和反饋信號弱的優點.Yuan 等[12]利用全局反饋法消除了螺旋波和時空混沌.2017年,Hou 等[13]利用局部反饋法對螺旋波進行了控制,并且研究了反饋控制區域和位置等參數的影響.但是,現有的研究都是利用反饋法對自由螺旋波進行控制,還沒有利用此方法研究釘扎螺旋波的消除.基于反饋法的優勢,本文使用反饋法對釘扎螺旋波進行控制.此外,還從反饋控制信號幅度、可激性系數、障礙物的形狀及大小等方面研究了反饋控制法對消除釘扎螺旋波的影響.

螺旋波的控制方程主要是反應擴散系統,學者們通常使用有限差分法研究螺旋波動力學.格子Boltzmann 方法(LBM)是一種介觀數值方法,廣泛應用于求解線性及非線性偏微分方程[14,15].它具有算法簡單有效、計算并行、復雜邊界容易處理等優點[16].本文使用格子Boltzmann 方法對Fitz-Hugh-Nagumo (FHN)模型下的釘扎螺旋波進行模擬,結合反饋法對釘扎螺旋波進行控制.

本文的結構如下: 第2 節對基于FHN 的反饋法模型進行格子Boltzmann 建模;第3 節分析模擬結果,首先研究圓形障礙和矩形障礙下釘扎螺旋波的控制情況,然后以圓形障礙物為例分別對各參數進行研究;第4 節對全文總結.

2 模 型

2.1 反饋法的格子Boltzmann 建模

FHN 模型是一種簡單的雙變量形式.它是對Hodgkin-Huxley 模型的簡化修正[17].近年來,FHN方程已成為反應擴散系統中一種廣泛應用的模型,通常用于模擬可激介質中的傳播,FHN 模型也被用于模擬神經纖維和心臟組織以及化學反應中[18,19].此外,它還可用于生物學、化學、傳熱傳質等方面的研究[20].FHN 模型的混沌性和同步性在激光技術、醫學和保密通信等領域具有潛在的應用價值[21].因此,研究FHN 模型的動力學特性在實際應用中具有重要意義.FHN 模型的形式如下[11]:

其中,a和b為無量綱常數,本文固定a=0.03 和b=2.0;u為激活變量,v為抑制變量;ε為可激性系數;α 為擴散系數.

為了使用反饋法控制釘扎螺旋波,在(1)式的右側加入反饋控制項[13]:

其中,kfb為反饋控制的幅度;Ifeedback為反饋控制項,其具體形式為

這里,D為控制的區域,t0為記錄反饋信號的時間,t1為加入反饋控制的時間,t2為反饋控制結束的時間,即反饋控制持續的時間為t1—t2.

針對方程(3),可建立LBM 的演化方程[15]:

其中,hi(x,t) 表示在位置x及時間t處具有離散速度ci的粒子分布函數,τ 表示松弛時間,表示分布函數hi(x,t) 的平衡態分布函數,Gi(x,t) 為源項分布函數.

通過直接泰勒展開[16],格子Boltzmann 方法可將演化方程恢復成方程(3),具體的恢復過程可以查看附錄A.

其中,ωi為權系數,滿足=1 .

源項為

宏觀量u(x,t) 為

對于方程(2),該方程沒有擴散項,則通過有限差分法進行求解:

2.2 邊界條件

以二維螺旋波為例,研究反饋法對釘扎螺旋波的控制.采用D2Q9 速度模型,速度及權重ωi大小如下:

其中,粒子速度c=?x/?t,?x為空間步長.

邊界的處理在數值模擬中至關重要,對采用無流邊界處理的螺旋波進行研究,無流邊界格式:

在使用格子Boltzmann 方法模擬螺旋波中,采用反彈格式對該宏觀邊界條件進行處理.反彈格式是一種啟發式格子Boltzmann 方法的邊界處理格式[22],包括標準反彈格式、半反彈格式和修正反彈格式.反彈格式通常被用來處理靜止無滑移邊界,其原理為粒子與邊界碰撞后沿原路返回,具有過程簡單、容易編程實現的優點.本文采用標準反彈格式進行邊界處理,邊界上的分布函數可表示為

其中,xb為邊界格點,xf為系統內部格點,為碰撞后的分布函數,i-為i的相反方向.

2.3 模型驗證

首先,對網格無關性進行驗證.物理邊界長度為Lx=Ly=40,選取參數ε=0.004,α=1.0 .分別選取800×800,1200×1200 和1600×1600 三種網格對模型形成自由螺旋波的過程進行模擬.該過程保持松弛時間τ 不變,選取時間步長為?t=1/1800.圖1(a)給出了在t=15時,800×800,1200×1200 和1600×1600 網格的計算情況.其中,藍色虛線表示800×800網格下u=0.5 的等值線圖,紅色虛線表示1200×1200網格下u=0.5 的等值線圖,藍色實線表示1600×1600網格下u=0.5的等值線圖.由圖可以看出,與網格800×800 相比,1200×1200和1600×1600 網格誤差較小,且在波頭位置幾乎重合.圖1(b)描繪了在y=15時,在網格800×800,1200×1200和1600×1600 下u的變化情況.藍色虛線、紅色實線和藍色實線分別表示800×800網格、1200×1200網格和1600×1600網格下u值情況.由圖1(b)可以看出1200×1200 和1600×1600 網格結果存在一定差異,但差異較小.

圖1 (a) t=15時,u=0.5 的等值線情況;(b) t=15時,y=15處u 值情況Fig.1.(a) Contour case for u=0.5at t=15 ;(b) value of u for y=15at t=15 .

為了進一步驗證網格的無關性,設定物理邊界長度Lx=Ly=40,其他參數ε=0.004,α=1.0 .計算了 1200×1200 網格和 1600×1600 網格下,不同半徑下釘扎螺旋波的波長情況.并根據波長情況計算出相對誤差,如表1 所列.相對誤差的計算公式如下:

表1 不同半徑下釘扎螺旋波的波長及相對誤差(1200×1200 網格和1600×1600 網格下)Table 1.Wavelengths and relative errors of pinned spiral waves at different radius (under 1200×1200 grids and 1600×1600 grids).

其中,λ 表示1200×1200網格下的波長數據,λ′表示1600×1600 網格下的波長數據.

從表1 的數據可知,在1200×1200 網格和1600×1600 網格下,半徑不同釘扎螺旋波的相對誤差均為 10-3量級.因此,為了提高計算效率并驗證格子Bolzmann 方法的準確性,選取1200×1200網格進行后續釘扎螺旋波的模擬.下面討論圓形障礙物下釘扎螺旋波的周期隨障礙物半徑的變化情況.設定物理邊界長度為Lx=Ly=40,選取參數ε=0.004,α=1.0,模擬結果如圖2 所示.可以看出,隨著圓形障礙物半徑的增加,釘扎螺旋波的周期也隨之增加.與文獻[11]中的結果進行對比,符合較好.

圖2 與文獻[11]對比,釘扎螺旋波的周期隨圓形障礙物半徑的變化Fig.2.Compared with Ref.[11],the period of the pinned spiral wave varies with the radius of the circular obstacle.

3 模擬結果分析

3.1 全局反饋控制

使用含反饋控制項的模型,對順時針旋轉的釘扎螺旋波進行控制.采用半徑為20 的圓形邊界情況進行討論,邊界條件采用無流邊界.選取參數ε=0.004,α=1.0.空間步長為?x=?y=1/30,選取時間步長為 ?t=1/1800 .圓形障礙物的中心位置為(20,20),障礙物半徑大小為R=2.0,控制的區域D選取為整個系統.

先使釘扎螺旋波在圓形邊界系統內順時針周期旋轉至少5 次.為了方便描述,以下時間t均為時間步數.圖3(a)描述了點(4,20)處在加入反饋控制前不同時間步t下u的變化值,圖3(b)描述了該點下u和v的關系.從圖中可以看出,釘扎螺旋波在半徑為2.0 的圓形障礙上進行周期運動,且通過u的最大值計算出此時釘扎螺旋波的周期為T=3020.根據釘扎螺旋波的周期,選定記錄反饋信號時間t0與加入反饋控制時間t1之間時間間隔為T/6(t1-t0=T/6).如果反饋信號注入易損期的某點,該點處反饋信號會被激發.如果反饋信號注入絕對不應期的某點,該點處反饋信號不會被激發[12].故選取記錄反饋信號時間為t0=19497,加入反饋控制時間為t1=20001,反饋控制結束時間為t2=20090 .

圖3 在點(4,20)處u 和v 的值(a) u 與時間步t 的關系;(b) u 和v 的相圖Fig.3.Values of u and v at the point (4,20): (a) Relationship between u and the time step t;(b) phase diagram of u and v.

圖4 給出了kfb=0.5下,使用LBM 模擬圓形障礙下釘扎螺旋波的控制情況.首先,釘扎螺旋波S1在系統內順時針旋轉,如圖4(a)所示.t=20001時間步時加入反饋控制,如圖4(b)所示.由于系統內部分點處于易損期,會使系統內激發形成第二個釘扎螺旋波S2,且S2的波寬會經歷逐漸增大后減小的過程,如圖4(c)—(f)所示.隨后,釘扎螺旋波S1順時針旋轉,而釘扎螺旋波S2逆時針旋轉.當兩個波陣面相遇時,根據激發態波的特性,激發態會碰撞消失,沒有進行湮滅的殘余部分形成激發態S3,如圖4(g)—(j)所示.此時,系統內沒有點缺陷行為,不會形成新的螺旋波.當S3被逐出系統時,系統恢復靜息態,如圖4(k)—(l)所示.

圖4 利用反饋法控制圓形障礙物下釘扎螺旋波的情況Fig.4.Suppression of a pinned spiral wave under a circular obstacle by feedback control approach.

為了更加方便地描述反饋法對釘扎螺旋波的控制效果及控制時間,采用以下形式來描述系統的行為:

其中,u(x,y) 為系統內全局的快變量u的值.當uall≤10-4時,認為系統恢復了靜息狀態.

圖5 給出了上述圓形障礙物下釘扎螺旋波控制過程中uall的變化情況,從加入反饋控制時間t1開始記錄uall.從uall的變化過程可以看出,uall的第一個小高峰值為控制結束時間t2,uall的第二個高峰值為釘扎螺旋波S2的波寬達到最寬時的狀態.注入反饋信號的1100 個時間單位后,此時uall≤10-4,系統內的釘扎螺旋波被完全消除,小于一個激發態周期.

圖5 ε=0.004, R=20, kfb=0.5 下反饋法控制釘扎螺旋波過程中 uall 的變化值Fig.5.Change of uall of unpinning spiral waves with ε=0.004, R=20, kfb=0.5.

接下來,為了討論反饋控制法下障礙物形狀的無關性,考慮矩形障礙物下釘扎螺旋波的控制情況.采用半徑為20 的圓形邊界情況進行討論,邊界條件采用無流邊界.選取參數ε=0.004,α=1.0.空間步長為?x=?y=1/30 .選取時間步長為?t=1/1800.正方形障礙物的中心位置為(20,20),障礙物邊長為L=4.0 .由數值模擬可知,釘扎螺旋波在正方形障礙邊長為4.0 時的周期為T=3621,則選定t0與t1之間時間間隔約為T/6 .因此,選取t0=19397,t1=20001,t2=20090 .

圖6 給出了kfb=0.5時,使用LBM 模擬正方形障礙下釘扎螺旋波的控制情況,可以發現方形障礙下利用反饋法控制釘扎螺旋波的情況跟圓形障礙下的控制情況類似.釘扎螺旋波S1在系統內繞著正方形障礙物順時針旋轉,加入反饋控制后,系統內形成第二個釘扎螺旋波S1,見圖6(a),(b).釘扎螺旋波S2會被系統激發成可激狀態,其波寬逐漸增大后減小,如圖6(c)—(f)所示.隨后,釘扎螺旋波S2繞著正方形障礙物順時針旋轉,而釘扎螺旋波S1繞著正方形障礙物逆時針旋轉,直到兩者相遇形成激發態S3,如圖6(g)—(j)所示.隨著S3被逐出系統,系統恢復靜息態,見圖6(k)—(l).

圖6 利用反饋法控制正方形障礙物下釘扎螺旋波的情況Fig.6.Suppression of a pinned spiral wave under a square obstacle by feedback control approach.

3.2 參數對釘扎螺旋波控制的影響

本節以圓形障礙物為例,考慮控制系統中各參數對釘扎螺旋波控制的影響.從(3)式可以看出,反饋控制的幅度kfb、可激性系數ε和記錄反饋信號時間t0對反饋控制有影響.此外,還將討論障礙物大小對釘扎螺旋波控制的影響.

3.2.1 反饋控制的幅度kfb

kfb的大小是反饋控制的關鍵因素,其取值一般介于0—1 之間.對于心臟系統來說,若kfb取值過大,則對系統損傷嚴重.取定擴散系數α=1.0,可激性系數ε=0.004,圓形障礙半徑R=2.0,在0—1 范圍內均勻增大kfb,考慮不同kfb下消除釘扎螺旋波的時間,如圖7 所示.從圖中可以看出釘扎螺旋波的控制時間隨著kfb的增大而減小.當kfb越大時,注入系統內反饋項的值越大,被激發成新的釘扎螺旋波所需的時間越少,因此影響釘扎螺旋波的控制時間.

圖7 不同 kfb 下消除釘扎螺旋波的時間步Fig.7.Time step of unpinning spiral wave with different .

3.2.2 可激性系數ε

可激性系數ε影響著釘扎螺旋波的動力學行為.接下來,考慮在不同的可激性系數ε下釘扎螺旋波的控制情況.控制擴散系數α=1.0,圓形障礙半徑R=2.0不變,討論能使不同的可激性系數ε下釘扎螺旋波得到控制的kfb最小值,這里用kfbmin來表示.加入反饋控制時間t1=20000,反饋控制結束時間t2=20090.圖8(a)記錄了不同ε下,釘扎螺旋波的周期T.根據不同ε下的周期T,選取記錄反饋信號時間t0,使其與t1之間的時間間隔為T/10,如圖8(b)所示.圖8(c)描繪了能成功控制釘扎螺旋波所需的kfbmin.可以看出kfbmin隨著可激性系數ε的增大呈現先減小后增大的現象.在保持kfb=0.5 的情況下,討論了不同可激性系數ε下的控制時間,如圖8(d)所示.反饋法下釘扎螺旋波的控制時間隨著可激性參數的增大而增大.

圖8 (a)釘扎螺旋波的周期隨可激性系數 ε 的變化;(b)記錄反饋信號時間 t0 隨可激性系數 ε 的變化;(c) kfbmin 隨可激性系數ε的變化;(d)釘扎螺旋波的控制時間t 隨可激性系數 ε 的變化Fig.8.(a) Variation of the period of the pinned spiral wave with ε ;(b) variation of the time of recording the feedback signal with ε ;(c) variation of kfbmin with ε ;(d) variation of the time required to eliminate the pinned spiral wave with ε .

此外,還討論了在不同的可激性系數ε下能夠成功控制釘扎螺旋波的記錄反饋信號時間t0與加入反饋控制時間t1之間時間間隔的范圍.在其余系統參數不變的情況下,釘扎螺旋波的周期跟障礙物的周長有關.并將該時間間隔與釘扎螺旋波的周期進行換算(例如,若t0與t1之間的時間間隔為60,某參數下釘扎螺旋波周期為600,則該時間間隔為周期的1/10).為了方便觀察記錄,將不同ε下的t1都選取為釘扎螺旋波的可激態剛好運動到圓形障礙物的正上方(即圓形障礙物的十二點鐘方向),并以T/360 為時間間隔進行討論,其中T為釘扎螺旋波的周期.從圖9 可以看出,隨著ε的增大,t0與t1之間時間間隔的范圍減小,這是由加入反饋信號時系統是否處于易損期決定的.并且在不同的ε下,t0與t1之間時間間隔最小均可為T/360,即當t0與t1之間的時間間隔為相應釘扎螺旋波周期的1/360時,也能成功消除釘扎螺旋波.

圖9 不同 ε下,能夠成功控制釘扎螺旋波的記錄反饋信號時間 t0 與加入反饋控制時間 t1 之間時間間隔對應的范圍Fig.9.Corresponding range of the time interval between the recording feedback signal time t0 and the adding feedback control time t1 with different ε .

3.2.3 記錄反饋信號時間t0及圓形障礙物半徑

由前面的討論可知,在不同的可激性系數ε下,能夠成功控制釘扎螺旋波的記錄反饋信號時間t0與加入反饋控制時間t1之間都有對應的一個時間間隔.若保持ε及t1不變,在對應的時間間隔內討論t0,則可以得到不同t0下能成功控制釘扎螺旋波的情況.控制變量不變,保持t0持續變化,討論不同t0下能夠成功控制釘扎螺旋波所需的最小的kfb,如圖10 所示.kfbmin隨著t0的增大呈現先減小后增大的情況,這與文獻[13]反饋法控制自由螺旋波的情況一致.圓形障礙物的大小對使用反饋法控制螺旋波也有重要的影響.圖11 給出了不同半徑下能夠成功控制釘扎螺旋波的記錄反饋信號時間t0與加入反饋控制時間t1之間對應的時間間隔.隨著圓形障礙半徑R的增大,t0與t1之間時間間隔的范圍減小.這是由加入反饋信號時系統是否處于易損期決定的.并且R不同時,能成功控制釘扎螺旋波下t0與t1之間時間間隔最小均可為T/360 .

圖10 在不同的t0下,能 成功控制釘扎螺 旋波所需的kfbmin 的變化情況Fig.10.Change of kfbmin of unpinning spiral wave successfully with different t0 .

圖11 不同R下,能夠成功控制釘扎螺旋波的記錄反饋信號時間 t0 與加入反饋控制時間 t1 之間時間間隔對應的范圍Fig.11.Corresponding range of the time interval between the recording feedback signal time t0 and the adding feedback control time t1 with different R.

4 結論

本文利用反饋控制法研究了釘扎螺旋波的消除情況,并研究了反饋控制中參數的影響.模擬結果表明,反饋法能夠有效控制釘扎螺旋波.首先,分別在圓形障礙物和矩形障礙物下驗證了反饋法對于釘扎螺旋波消除的有效性.其次,研究了反饋控制信號幅度和可激性系數對控制釘扎螺旋波的影響,并研究記錄了能夠成功控制釘扎螺旋波的記錄反饋信號時間與加入反饋控制時間之間對應的時間間隔.可以發現釘扎螺旋波的控制時間隨著反饋控制信號幅度的增大而減小.而釘扎螺旋波的控制時間隨著可激性系數的增大而增大,這是由于可激性系數增大時,螺旋波的波寬隨之減小,從而影響控制時間.最后,研究了圓形障礙的半徑和記錄信號時間的影響.在不同的記錄信號時間下,能成功消除釘扎螺旋波的最小反饋控制信號幅度不同,呈現先減小后增大的趨勢.基于本文的討論可知,反饋法對于自由螺旋波和釘扎螺旋波都具有很好的控制作用.基于該方法具有反饋信號弱和反饋持續時間短的優點,對該方法進一步研究,可以大幅度提高控制螺旋波的效果.

附錄A

在本附錄中,采用直接泰勒展開(DTE)從現有的格子Boltzmann 模型中恢復帶有反饋控制項的FHN 模型.對方程(5)進行泰勒展開,可以得到:

由(A2)式可以得到:

將(A4)式代入(A3)式,有

為了恢復正確的帶有反饋控制項的FHN 模型,分布函數hi,和Gi滿足以下條件:

通過上述關系,對具有O(?t)和O(?t2) 階的(A2)式和(A5)式求零階矩,可得:

將(A10)式代入(A9)式即可恢復宏觀方程(3),且α=(τ-1/2).