双层耦合非对称反应扩散系统中的超点阵斑图*

刘富成 刘雅慧 周志向 郭雪 董梦菲

(河北大学物理科学与技术学院,保定 071002)

通过线性耦合Brusselator模型和Lengyel-Epstein模型,数值研究了双层耦合非对称反应扩散系统中图灵模之间的相互作用以及斑图的形成机理.模拟结果表明,合适的波数比以及相同的对称性是两个图灵模之间达到空间共振的必要条件,而耦合强度则直接影响了图灵斑图的振幅大小.为了保证对称性相同,两个图灵模的本征值高度要位于一定的范围内.只有失稳模为长波模时,才能对另一个图灵模产生调制作用,并形成多尺度时空斑图.随着波数比的增加,短波模子系统依次经历黑眼斑图、白眼斑图以及时序振荡六边形斑图的转变.研究表明失稳图灵模与处于短波不稳定区域的高阶谐波模之间的共振是产生时序振荡六边形的主要原因.

1 引 言

自组织斑图是非线性动力学系统的一个典型特征,普遍存在于自然界以及众多非线性系统中,例如物理、化学以及生物系统等[1-6].在过去的半个世纪内,人们对斑图形成机理的认识有了突飞猛进的发展,其中比较著名的就是反应扩散系统的图灵机制[7],一个均匀稳定态经历图灵失稳后会产生不随时间变化的空间周期性分布的结构,例如条纹斑图和六边形斑图[8,9].

近年来,多层耦合系统中不同模式之间的相互作用越来越受到人们的关注[10-18],其中不同图灵模之间的相互作用可以产生各种复杂的多尺度时空斑图.2002年,Yang等[11]通过线性耦合双层反应扩散系统研究了图灵模之间的相互作用,结果表明图灵模之间的波数比、子斑图对称性以及本征值高度是空间共振的三个必要条件.该研究小组还发现仅仅是双层间的耦合作用也可以导致一种新的叉形分叉,进而产生两个图灵模式.这两个图灵模相互作用可以形成对称、非对称以及反相图灵斑图[12].在此基础上,多个研究小组对图灵模相互作用机理进行了进一步研究.例如Liu等[13]在更广的本征值高度范围内获得了超点阵斑图,并且在波数比接近4∶1时观察到一种新类型的白眼斑图.Catlla等[14]发现当两个图灵模的波数比例为时,系统可以产生稳定的四边形斑图.白占国等[15-18]详细研究了各种耦合条件对斑图的影响,发现了多种超点阵四边形斑图,包括大小点、白眼、点线和复杂超四边等.

以前的研究主要集中在相同反应扩散模型[19-27]之间的耦合,例如两个Brusselator模型[19,20]、两个Lengyel-Epstein模型[21-23]或者是两个Oregonator模型[24]等.在大多数非线性系统中,同时存在着多种不同的非线性物理过程,这些过程各自具有不同的斑图形成机理.例如在介质阻挡放电系统中,放电层与电介质层相互耦合,可以产生种类繁多的时空斑图,这些时空斑图通常由两套或者多套子结构相互嵌套而成[25-27],因此有必要对不同物理过程之间的耦合进行详细研究.然而到目前为止,不同反应模型之间的耦合以及斑图选择机制还尚未见报道,人们对不同斑图形成系统之间的耦合机制还不甚清楚.

本文通过线性耦合两个不同的反应扩散模型数值研究了两个图灵模之间的相互作用,并分析了两个图灵模的波数比、相对强度以及耦合强度对斑图选择和形成的影响.

2 模型分析

为了研究不同系统形成的图灵模之间的相互作用,本文采用线性耦合两个不同的反应扩散模型,构建了一个双层耦合非对称反应扩散系统模型,在无量纲的情况下,该模型具有以下形式:

该模型包含两个子反应扩散系统 (u1,v1)和(u2,v2),α为两个子系统活化子之间的耦合强度;方程f(u,v)和 g (u,v) 为每个子反应扩散系统的局部动力学方程,不同的反应系统有着不同的具体表达式.本文选取Brusselator模型和Lengyel-Epstein模型来充当两个子系统.在无量纲的情况下,Brusselator模型的局部动力学方程具有下列形式:

其中系数a和b是系统的控制参数.

Lengyel-Epstein模型在无量纲情况下的形式为

其中系数c和d是系统的控制参数.对于子系统1和2,其均匀定态解分别为 (u10,v10)=(a,b/a) 和本文中选取c=5a,则双层耦合系统的均匀定态解可以表示为

对该均匀定态解附加一个均匀的线性微扰,去掉高阶项得到其线性本征方程:

其中A和D分别为雅克比矩阵和扩散系数矩阵,

解此本征方程可得到本征值λ,本征值λ与微扰模之间的关系被称之为系统的色散关系曲线.通过系统的色散关系可以分析失稳模的性质以及系统的分岔类型.

图1给出了两个耦合系统处于图灵空间时的色散关系.图1中蓝色实线为本征值的实部,代表模式的线性增长率,红色虚线为本征值的虚部.对于每一个子系统,当活化子u的扩散系数远小于禁阻子v的扩散系数时,系统经历图灵分岔并激发出图灵模.该图灵模的波数大小为其中Δ是与控制参数有关的常数.Du相比于 Dv越小,系统偏离分岔点越远,图灵模线性增长率越高,即图灵模的本征值高度越高.因此通过调节模型中各变量的扩散系数,可以使两个子系统均经历图灵分岔,从而各自激发相应的图灵模式.由于图灵模的波数与扩散系数成反比,所以具有较小扩散系数的子系统的图灵模波数大(短波模),称为 k2,本征值高度为 h2,而具有较大扩散系数的子系统的图灵模波数较小(长波模),称为 k1,相应的本征值高度为 h1.图1(a)和图1(b)给出了两个图灵模波数比分别为1∶2和1∶3时的色散关系图.研究发现,只有长波模 k1对短波模 k2有影响,而短波模 k2对长波模 k1没有影响.因此长波模 k1一般设置为失稳模.

图1 不同参数下耦合系统的色散关系 (a) Du1=12.6,Dv1=27.9,Du2=22,Dv2=420,α=0.1;(b) Du1=5.3,Dv1=20,Du2=22,Dv2=500,α=0.1Fig.1.Dispersion relations of coupled systems under different parameters:(a) Du1=12.6,Dv1=27.9,Du2=22,Dv2=420,α=0.1;(b) Du1=5.3,Dv1=20,Du2=22,Dv2=500,α=0.1.

3 模拟结果与讨论

采用欧拉向前差分的方法进行积分,数值模拟在一个含有 1 28×128 个空间格点的二维平面上进行,时间积分步长取为 Δ t=0.01 个时间单位,空间积分步长取为 Δ h=1.0 个空间单位,扩散项在数值模拟中采取5点差分格式,边界条件选用周期性边界条件,初始条件为均匀定态加上一个很小的随机扰动.如果没有特殊说明,一般设置子系统2具有较大的扩散系数,即激发长波模,而观察子系统1中形成的各种斑图.

接下来分别从波数比、相对强度以及耦合强度三个角度来讨论两个图灵模之间的耦合作用以及相应的斑图结构.

3.1 两个图灵模波数比的影响

固定两个图灵模的本征值高度(h1,h2)=(0.20,-0.55)以及耦合强度 α=0.1,通过改变两个图灵模的波数比值 K1:k2,研究不同波数比对斑图形成的影响.如前文所述,图灵模的波数大小取决于 Du和Dv的乘积,而其高度取决于 Du和Dv的相对大小,因此两个图灵模的波数 k1和k2的大小可以通过其各自的 DuDv来控制.波数大小确定后,保持 DuDv的乘积不变,改变二者的相对大小,调节其本征值高度h.

图2给出了波数比分别为1∶2,1∶3和1∶4三种比值下子系统1中形成的超点阵斑图及其相应的傅里叶频谱图.图2(a)是波数比为1∶2时得到的黑眼斑图,其中每个中心点的浓度最低,且周围环绕着一圈高浓度环,整体呈现六边形阵列排布.从其傅里叶频谱图可以看出,该黑眼斑图包含三套不同波数的模式,其波数值分别为 k1,k2和q.每一套波矢都组成一个完美的六边形结构,其中 k1和k2为两个子系统激发的图灵模 (K1:k2=1:2),而q是新生模式.这三种模式之间满足空间三波共振关系:K1+k2=q.图2(b)是波数比为1∶3时得到的白眼六边形斑图,与黑眼斑图不同,白眼斑图的中心为高浓度斑点,周围依次围绕一圈低浓度环和一圈高浓度环,然后被6个更高浓度的亮点包围,整体也呈现六边形阵列结构.分析其傅里叶频谱图可知,该白眼斑图同黑眼斑图类似,也包含上述三套满足空间共振关系的波矢,唯一的区别是新生模式q具有两套六边形结构,即12个波矢.通过几何关系可知,此两套六边形之间的夹角为θ=21.8°.图2(c)是波数比为1∶4时得到的白眼斑图.与图2(b)相比,每一个六边形晶胞中心的亮斑点浓度下降,同时,晶胞外围的六个亮点连接在一起,形成了一个蜂窝状六边形结构.它包含的波矢成分也较为复杂,除了基本的图灵模 k1和k2外,还产生了各种高阶模 2 k1,3 k1以及这些模式共振而生成的新生模式q.需要特殊说明的是,图灵模式 k2的强度较弱,而高阶模 3 k1的强度较强,这意味着当波数比较大时,子系统1中的斑图主要是失稳模 k1及其激发的高阶模相互作用而形成.

图2 不同波数比下的超点阵斑图及其傅里叶频谱图 (a) 1∶2下的黑眼斑图,Du1=13.5,Dv1=27.5,Du2=22,Dv2=400 ;(b) 1∶3下的白眼斑图,Du1=6,Dv1=12.3,Du2=22,Dv2=400 ;(c) 1∶4下的白眼斑图,Du1=3.4,Dv1=6.96,Du2=21.9,Dv2=400.α=0.1Fig.2.Superlattice pattern and fourier spectrum under different wave number ratios:(a) Black eye pattern at 1∶2,Du1=13.5,Dv1=27.5,Du2=22,Dv2=400 ;(b) white eye pattern at 1∶3,Du1=6,Dv1=12.3,Du2=22,Dv2=400 ;(c) white eye pattern at 1∶4,Du1=3.4,Dv1=6.96,Du2=21.9,Dv2=400.α=0.1.

图3给出了两个图灵模波数比为1∶5时所产生的振荡超六边形斑图.图3(a)为相应的色散关系曲线图,由图3(a)可知,在 k1和k2之间存在着一个虚部不为零的区域,即短波不稳定性区域,在此区域内的各个模式均随时间周期性振荡,失稳模k1的高阶模 2 k1和3 k1均位于此区域内.在此情况下,失稳图灵模 k1与其高阶模 2 k1和3 k1相互作用时,就导致了一种正弦型时间振荡斑图的形成.图3(b)给出了三个不同位置处(分别标记为A,B,C)活化子 u1随时间变化的图.其中B点为振荡六边形的中心斑点位置,其格点坐标为(50,76);A点为两个中心亮斑的中心,其格点坐标为(53,59);C点为三个中心亮斑组成的三角形中心,其格点坐标为(65,61),如图3(c)所示.由图3(c)可知,三个位置处的振荡周期相同,均为 T=t6-t1=46 个时间单位,但是相位不同,A,B,C相位依次延迟.另外,这三点的振荡幅度也不相同,B点处振荡幅度远高于其他两点.图3(c)显示了一个振荡周期内该振荡六边形斑图的时空演化过程.在t1=24时刻,六边形亮斑结构浓度稍高于其它位置,从而形成了一种具有晕结构的超六边形斑图.在t2=30时刻,中心斑点浓度下降,低于晕的浓度,此时斑图呈现一种类似白眼斑图的结构.在 t4=47 时刻,中心斑点处浓度降至最低点,远远低于晕的浓度,此时斑图呈现为简单的蜂窝六边形斑图.随着时间演化,中心斑点处浓度开始上升,在 t6=64 时刻开始显现出中心亮斑,有意思的是单元中心由原来的一个斑点变为两个斑点,并且每个单元外包围着6个暗点.该时序振荡六边形斑图与对称耦合反应扩散模型获得的眨眼斑图[11]有着很大的区别.首先是形成机理不同,该振荡六边形是失稳模与其高阶模相互作用而成,而眨眼斑图则是失稳模与其次谐波模相互作用而成;其次是振荡模式不同,该振荡六边形呈现正弦式的振荡,而眨眼斑图则呈现脉冲式的振荡.

图3 波数比为1∶5时的时间振荡超六边形斑图,Du1=2.2,Dv1=4.5,Du2=21.9,Dv2=400,α=0.1 (a) 色散关系曲线;(b) 三个位置处u1的时间变化关系图;(c) 一个振荡周期内的斑图演化过程Fig.3.Oscillatory super-hexagon pattern with wave number ratio of 1∶5,Du1=2.2,Dv1=4.5,Du2=21.9,Dv2=400,α=0.1:(a) Dispersion curve;(b) time variation of u1 at three positions;(c) evolution of pattern in an oscillating period.

3.2 图灵模强度对斑图的影响

图灵模的本征值高度h代表了图灵模的线性增长率,h越大,表示增长率越快,意味着图灵模强度越强.在本节中固定两个图灵模波数比为1∶3以及耦合强度 α=0.1,通过改变两个图灵模的本征值高度h来研究图灵模强度对斑图形成产生的影响.

首先固定本征值高度 h1=0.2,讨论本征值高度 h2对图灵斑图的选择及形成产生的影响.图4给出了 h2逐步升高而获得的各种复杂斑图及其相应的傅里叶频谱图,其中图4(a)—(d)是 h2＜0 的情况,图4(e) 是 h2＞0 的情况.当 h2=-2.56 时系统1呈现简单的蜂窝六边形斑图(图4(a)).由于h2较低,该系统的本征模并不能被激发,因此该六边形斑图完全是由系统2的失稳模调制形成.随着h2的升高,系统1中的本征模开始被激发,如前文所述,两个图灵模相互耦合,形成了白眼斑图(图4(b)和图4(c)).继续升高 h2至—0.31,白眼斑图失稳,转变为一种新的超六边形结构,该斑图中每个晶胞均由三个相邻的暗点排列而成,这些暗点的空间波数为 k2,晶胞之间以波数 k1周期性排列(图4(d)).当h2=0.67时,系统1呈现单一的条纹斑图,也就是说此时系统1的斑图完全由其自身的本征模所决定,系统2的长波模对它的调制作用仅仅存在于背景中,这可以通过傅里叶频谱图看出(图4(e)).

另外h还表征了系统远离初级分岔点的程度,h越大,偏离程度越高.早期的研究表明[13],在初级分岔点附近,系统一般选择为六边形斑图,而当远离分岔点时,选择为条纹斑图.因此通过改变h的大小还可以控制斑图的对称性类型.图5给出了 h1分别为0.2和0.4两种情况下子系统1中形成的复杂斑图,这里固定两个图灵模的波数比为1∶3,h2=-0.5,耦合强度 α=0.1.当 h1=0.2 时,两个子系统中的斑图均为六边形斑图,对称性相同,因此可以满足空间共振关系,从而形成了白眼斑图(图5(a)).但是当 h1增大到0.4时,该失稳模呈现出条纹结构,与另一个图灵模的对称性不匹配,不满足空间共振关系,因此在子系统1中仅仅呈现简单的重叠斑图,即条纹结构和六边形的叠加(图5(b)).

图4 不同本征值 h 2 下的复杂斑图及其傅里叶频谱图 (a) 蜂窝状六边形斑图 h 2=-2.56,Du1=8.5,Dv1=12.5,Du2=22,Dv2=400 ;(b) 白眼斑图 h 2=-1.3,Du1=7,Dv1=12.3,Du2=22,Dv2=400 ;(c) 白眼斑图 h 2=-0.55,Du1=6,Dv1=12.3,Du2=22,Dv2=400 ;(d) 超六边形斑图 h 2=-0.31,Du1=5.9,Dv1=12.7,Du2=22,Dv2=400 ;(e) 条纹斑图 h 2=0.67,Du1=5.5,Dv1=16,Du2=22,Dv2=400.α=0.1.Fig.4.Complex patterns and Fourier spectrum under different eigenvalues h2 :(a) Honeycomb hexagon pattern h 2=-2.56,Du1=8.5,Dv1=12.5,Du2=22,Dv2=400 ;(b) white-eye pattern h 2=-1.3,Du1=7,Dv1=12.3,Du2=22,Dv2=400;(c) white-eye pattern h 2=-0.55,Du1=6,Dv1=12.3,Du2=22,Dv2=400 ;(d) super-hexagon pattern h2=-0.31,Du1=5.9,Dv1=12.7,Du2=22,Dv2=400 ;(e) stripe pattern h 2=0.67,Du1=5.5,Dv1=16,Du2=22,Dv2=400.α=0.1.

详细研究发现,在该双层耦合非对称反应扩散系统中,两个子系统要实现空间共振,其斑图必须均为六边形结构,即两个模式的空间对称性必须相同.对于失稳模 k1,当本征值高度处于0＜h1＜0.33时,子系统自发选择六边形斑图;当 h1＜0 时,系统是稳定的,无法激发图灵模;当 h1＞0.33 时,子系统选择条纹斑图.对于图灵模 k2,其本征值高度范围需在-2＜h2＜0 ;若 h2＞0 时转变为失稳模,此时子系统选择条纹斑图;而当其本征值高度过小时,尽管受到失稳模 k1的激励作用,但是依然无法激发出该本征模,此时子系统处于均匀定态.

图5 不同本征值 h1 下的复杂斑图及其傅里叶频谱图 (a) 白眼斑图,h 1=0.2,Du1=6.1,Dv1=12.7,Du2=22.3,Dv2=403;(b) 条纹点状斑图,h 1=0.4,Du1=6.1,Dv1=12.6,Du2=20.3,Dv2=464,α=0.1Fig.5.Complex patterns and fourier spectrum under different eigenvalues h1 :(a) White-eye pattern,h 1=0.2,Du1=6.1,Dv1=12.7,Du2=22.3,Dv2=403 ;(b) stripe-spot pattern,h 1=0.4,Du1=6.1,Dv1=12.6,Du2=20.3,Dv2=464,α=0.1.

图6 不同耦合强度下的超六边形斑图 (a) 白眼斑图,α=0.01,Du1=5.8,Dv1=11.4,Du2=22,Dv2=367 ;(b) 白眼斑 图,α=0.1,Du1=6.3,Dv1=12.9,Du2=21.8,Dv2=395 ;(c) 新型超六边形斑图,α=0.2,Du1=6.1,Dv1=13,Du2=22,Dv2=432 ;(d) 新白眼斑图,α=0.3,Du1=6.3,Dv1=13.98,Du2=22,Dv2=460Fig.6.Super-hexagon patterns with different coupling strength:(a) White-eye pattern,α=0.01,Du1=5.8,Dv1=11.4,Du2=22,Dv2=367 ;(b) white-eye pattern,α=0.1,Du1=6.3,Dv1=12.9,Du2=21.8,Dv2=395 ;(c) new superhexagon pattern,α=0.2,Du1=6.1,Dv1=13,Du2=22,Dv2=432 ;(d) new white-eye pattern,α=0.3,Du1=6.3,Dv1=13.98,Du2=22,Dv2=460.

3.3 耦合强度的影响

耦合强度也是影响两个图灵模相互作用的一个重要因素.由于各个参量共同影响两个图灵模的大小和强度,本文中改变耦合强度的同时对扩散系数进行微调,以保证两个图灵模的大小 K1=0.2,k2=0.6和本征值高度 (h1,h2)=(0.2,-0.55) 始终保持恒定,从而研究耦合强度对斑图形成的影响.图6给出了不同耦合强度下获得的超六边形点阵斑图.模拟研究发现,耦合强度不仅可以影响斑图的振幅大小,还影响了其空间结构的具体形式.如果耦合系数足够小,失稳模对子系统1的影响可以忽略不计,此时子系统1呈现稳定的均匀状态.当耦合系数 α=0.01 时,子系统1开始受到失稳模k1的调制和激发,系统形成白眼斑图,但是此时的活化子浓度的变化范围非常小(图6(a)).图6(b)为耦合系数 α=0.1 时模拟得到的白眼斑图,从图6(b)中可以看出,活化子浓度的变化范围变大,图灵模k2的强度也开始增强.继续增大耦合强度α,子系统1中活化子浓度的变化范围持续增加,并且其空间形貌逐渐开始变化.图6(c)是模拟获得的一种新型超六边形斑图,此超点阵斑图虽然整体也呈现六边形阵列结构,但是每一个晶胞单元中心不再是一个圆斑点,而是形成了一个线段,这可能是六边形斑图和条纹斑图共同竞争的结果.图6(d)给出了耦合强度 α=0.3 时得到的新白眼斑图,与图6(b)不同,此白眼斑图在中心亮点周围还包围着一圈高浓度环,环外排列着一组六边形斑点.从傅里叶频谱图中可以看出,耦合强度很大,此斑图也是由两个图灵模 k1和k2直接耦合而成.

4 结 论

通过线性耦合两个不同的反应扩散模型(Brusselator模型和 Lengyel-Epstein模型)研究了不同系统间两个图灵模之间的相互作用以及斑图的形成和选择过程.模拟结果表明合适的波数比以及相同的对称性是不同反应扩散系统的图灵模之间相互共振的必要条件.只有长波模子系统对短波模子系统有调制和激发作用,因此超点阵斑图只能出现在短波模子系统中.在不同的波数比下,短波模子系统内可以形成包括黑眼斑图和白眼斑图在内的多种超点阵斑图.为了保证两个子系统斑图具有相同的对称性,两个图灵模的本征值高度必须处于 0＜h1＜0.33和-2＜h2＜0.当失稳模的高阶谐波处于短波不稳定区域时,失稳模与其高阶模相互共振可以形成一种正弦型时序振荡六边形斑图.耦合强度直接影响了图灵斑图的振幅大小.本文的研究结果有助于人们进一步认知非线性系统中各种模式之间的相互作用机理.