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

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

1) (華中科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,武漢 430074)

2) (華中科技大學(xué),工程建模與科學(xué)計(jì)算湖北重點(diǎn)實(shí)驗(yàn)室,武漢 430074)

3) (華中科技大學(xué),數(shù)學(xué)與應(yīng)用學(xué)科交叉創(chuàng)新研究院,武漢 430074)

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

1 引言

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

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

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

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

2 模 型

2.1 反饋法的格子Boltzmann 建模

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

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

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

其中,kfb為反饋控制的幅度;Ifeedback為反饋控制項(xiàng),其具體形式為

這里,D為控制的區(qū)域,t0為記錄反饋信號(hào)的時(shí)間,t1為加入反饋控制的時(shí)間,t2為反饋控制結(jié)束的時(shí)間,即反饋控制持續(xù)的時(shí)間為t1—t2.

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

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

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

其中,ωi為權(quán)系數(shù),滿足=1 .

源項(xiàng)為

宏觀量u(x,t) 為

對(duì)于方程(2),該方程沒有擴(kuò)散項(xiàng),則通過有限差分法進(jìn)行求解:

2.2 邊界條件

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

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

邊界的處理在數(shù)值模擬中至關(guān)重要,對(duì)采用無流邊界處理的螺旋波進(jìn)行研究,無流邊界格式:

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

其中,xb為邊界格點(diǎn),xf為系統(tǒng)內(nèi)部格點(diǎn),為碰撞后的分布函數(shù),i-為i的相反方向.

2.3 模型驗(yàn)證

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

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

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

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

其中,λ 表示1200×1200網(wǎng)格下的波長數(shù)據(jù),λ′表示1600×1600 網(wǎng)格下的波長數(shù)據(jù).

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

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

3 模擬結(jié)果分析

3.1 全局反饋控制

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

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

圖3 在點(diǎn)(4,20)處u 和v 的值(a) u 與時(shí)間步t 的關(guān)系;(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在系統(tǒng)內(nèi)順時(shí)針旋轉(zhuǎn),如圖4(a)所示.t=20001時(shí)間步時(shí)加入反饋控制,如圖4(b)所示.由于系統(tǒng)內(nèi)部分點(diǎn)處于易損期,會(huì)使系統(tǒng)內(nèi)激發(fā)形成第二個(gè)釘扎螺旋波S2,且S2的波寬會(huì)經(jīng)歷逐漸增大后減小的過程,如圖4(c)—(f)所示.隨后,釘扎螺旋波S1順時(shí)針旋轉(zhuǎn),而釘扎螺旋波S2逆時(shí)針旋轉(zhuǎn).當(dāng)兩個(gè)波陣面相遇時(shí),根據(jù)激發(fā)態(tài)波的特性,激發(fā)態(tài)會(huì)碰撞消失,沒有進(jìn)行湮滅的殘余部分形成激發(fā)態(tài)S3,如圖4(g)—(j)所示.此時(shí),系統(tǒng)內(nèi)沒有點(diǎn)缺陷行為,不會(huì)形成新的螺旋波.當(dāng)S3被逐出系統(tǒng)時(shí),系統(tǒng)恢復(fù)靜息態(tài),如圖4(k)—(l)所示.

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

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

其中,u(x,y) 為系統(tǒng)內(nèi)全局的快變量u的值.當(dāng)uall≤10-4時(shí),認(rèn)為系統(tǒng)恢復(fù)了靜息狀態(tài).

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

圖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.

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

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

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

3.2 參數(shù)對(duì)釘扎螺旋波控制的影響

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

3.2.1 反饋控制的幅度kfb

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

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

3.2.2 可激性系數(shù)ε

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

圖8 (a)釘扎螺旋波的周期隨可激性系數(shù) ε 的變化;(b)記錄反饋信號(hào)時(shí)間 t0 隨可激性系數(shù) ε 的變化;(c) kfbmin 隨可激性系數(shù)ε的變化;(d)釘扎螺旋波的控制時(shí)間t 隨可激性系數(shù) ε 的變化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 ε .

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

圖9 不同 ε下,能夠成功控制釘扎螺旋波的記錄反饋信號(hào)時(shí)間 t0 與加入反饋控制時(shí)間 t1 之間時(shí)間間隔對(duì)應(yīng)的范圍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 記錄反饋信號(hào)時(shí)間t0及圓形障礙物半徑

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

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

圖11 不同R下,能夠成功控制釘扎螺旋波的記錄反饋信號(hào)時(shí)間 t0 與加入反饋控制時(shí)間 t1 之間時(shí)間間隔對(duì)應(yīng)的范圍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 結(jié)論

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

附錄A

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

由(A2)式可以得到:

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

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

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

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