双层非线性耦合反应扩散系统中复杂Turing斑图*

李新政1) 白占国1)† 李燕1) 赵昆1) 贺亚峰2)

1)(河北科技大学理学院,石家庄 050018)

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

1 引言

斑图是自然界中广泛存在的一种典型的非线性自组织现象,是在时间或空间以及时空上具有某种规律性的非均匀结构.在化学反应系统、非线性光学系统以及生物系统中均可观测到种类丰富的斑图[1-7],虽然各个系统斑图形成的具体物理机制各不同,但斑图的形成与演化规律却具有一定的共性,Turing在1952年首先利用化学反应扩散方程组对斑图的形成进行了理论解释[8].近年来,介质阻挡放电系统作为一个新兴的斑图系统,以其本身发光可视、形成斑图所需时间短、构成斑图的单元(放电丝)可被单独测量、控制参数调节方便等优点备受关注[9-14].气体放电系统属于一种远离平衡的非线性系统,能在高级分岔点附近研究复杂斑图的形成和演化,在该系统中已从实验获得了种类丰富的复杂斑图(例如:蜂窝、白眼、超四边、超六边、点线斑图和螺旋波等),经研究发现复杂斑图往往由两套及两套以上子结构嵌套而成,归因于壁电荷斑图与放电丝斑图相互作用的结果[15,16].构成复杂斑图的简单子斑图结构是按照一定的波长关系组合的,只有当系统内出现的失稳模具有两个或多个波长时,才会产生具有复杂的空间结构和对称性的超晶格斑图,因此怎样获得多个不同波长的失稳模以及不同波长模间的相互作用是研究的热点[17-19].人们开始利用双层耦合模型对实验中观察到的超点阵斑图进行理论模拟,主要源于Barrio等[18]提出超点阵斑图的形成是因为两个模间的相互作用,并获得了许多有意义的研究成果.例如Yang等[6]用耦合两个Brusselator扩散反应模型得到了黑眼和白眼超点阵斑图,Berenstein等[20]用双层耦合的Lengel-Epstein模型加四边形或六边形挡板数值模拟了在化学反应扩散系统观察到的多种超点阵和振荡斑图等.前人在该方面的研究主要集中于模与模的线性相互作用,对非线性相互作用研究较少.本文采用双层耦合Brusselator模型,研究模之间的非线性相互作用,获得了种类丰富的斑图,通过分析Turing模对不同斑图形成的影响,从理论上分析了复杂斑图的形成机理及其动力学行为.

2 模型

利用反应扩散方程研究Turing分岔斑图的主要有Lengyel-Epstein模型、Brusselator模型和Schnackenberg模型等,其中Lengyel-Epstein模型和Brusselator模型应用最广[6,20-28].Brusselator模型是以双变量形式来描述化学反应过程中化学元素变化的一类反应扩散方程组.文中所采用的是双层Brusselator四变量反应扩散模型,每层中包括一个活化子和一个禁阻子,在无量纲的情况下该模型形式如下:其中i和 j表示两子系统1(u1,v1)和2(u2,v2),u和v分别为活化子和禁阻子的浓度;Du和Dv分别为u和v变量的扩散系数,α,β为两子系统非线性耦合系数(且i/=j),动力学行为由控制参数a和b所决定.为便于研究与对比,本文中固定a=3和b=9,两个子系统具有相同的均匀定态解:(u0,v0)=(a,b/a).通过对该均匀定态解做线性稳定性分析可以得到斑图动力学中最常见的两类初级分岔现象:当控制参数b满足b＞bH=(1+a2)时,系统经√历Hopf分 岔;当控制参数b满足:b＞bT=时√,系统经历Turing分岔并产生特征波数为k2=a的Turing模.选择适当参数使系统处于Turing空间,从式中可以看出Turing模的波数大小反比于两个变量的扩散系数的大小.图1是两系统的色散关系:具有较大扩散系数的失稳模(长波模)具有较小的波数,而具有较小扩散系数的稳定模(短波模)具有较大的波数.通过合理选择参数使两个子系统分别具有长波模和短波模,进而讨论2模作用引起的复杂斑图.

在本文的数值模拟中,采用欧拉向前差分的方法进行积分,数值模拟在一个含有N×N(128×128)个格子的二维平面上进行,时间积分步长为Δt=0.01个时间单位,初始条件为均匀定态上加一个很小的随机扰动,微扰幅度Δu=0.03u0,Δv=0.01ν0,边界条件选用周期性边界条件,为确保其稳定性所有结果的积分时间均超过10000个时间单位.

图1 两耦合系统色散关系

3 模拟结果与分析

斑图的选择完全取决于两系统Turing模之间的相互作用,在一定条件下两个子系统Turing模会产生空间共振,若长波模和短波模满足共振条件,在短波模所处的子系统中会形成复杂的斑图.由于两子系统间的相互作用较为复杂,为简化分别通过改变系统反应参数对Turing模间的相互作用进行研究.

3.1 波数比值的影响

斑图的形成首先是因为系统内禁阻子和活化子浓度在时间及空间上分布不均匀造成的,扩散系数的大小对于斑图的形成及类型起着决定性作用.本文通过改变短波模与失稳模的波数比,观察波数比的改变对系统1斑图的形成及类型选择的影响.

图2中是两子系统波数比为1:1时,两子系统出现相同的简单六边形与条纹斑图.波数相等时两子系统具有相同的Turing模,且地位相同,此时系统处于初级Turing分岔点附近,两波模无耦合发生,不会产生新的模,此时系统对失稳模的波长选择是单一的或一条较窄的波段,同时斑图的自组织行为受空间共振规律的约束,因此系统仅有单一空间尺度的简单斑图形成,从图形相应的二维傅里叶频谱中可以看出此时系统仅有一种波数的波矢,根据镜像对称图2(a3)中的两个波矢可以写成q1和-q1,图2(b3)中的六个波矢可以约化为三个大小相等基矢q1,q2,q3(|q1|=|q2|=|q3|).

图2 波数比1:1时,(a)Du1=51.7,Dv1=155.5,Du2=51.7,Dv2=155.5,α=0.03,β=0.03;(b)Du1=27.75,Dv1=73.75,Du2=27.75,Dv2=73.75,α=0.045,β=0.045(其中 1 为子系统1图,2为子系统2图,3为傅里叶频谱图)

两子系统波数比大于1时,两子系统Turing模波长不同,地位不同,短波模受失稳模调制,系统1中斑图受到系统2中斑图的影响,而系统2中的斑图不受系统1中斑图的影响.通过调整参数发现:短波模受失稳模调制不仅可以出现系统2中简单形状斑图(如六边、条纹、四边形等斑图);失稳模与短波膜发生共振偶合时,系统有新的不同波长失稳模产生,若系统处于高级Turing分岔点附近时还可以出现丰富的超点阵斑图.在调整波数比的过程中始终保持子系统2中呈现稳定的简单六边形斑图,分别考察了波数比为整数2:1,3:1,4:1时子系统1中产生复杂斑图情况.

图3(a)为波数比为2:1时产生蜂窝状六边形超点阵斑图,每个蜂窝单元边缘均由短线构成,中心区是相同的圆形点状斑图;从相应二维傅里叶频谱中可以看出蜂窝斑图有三种不同波数的波矢,每组相同波数的波矢均可构成一六边形,中心区为一组等波数的谐振波矢qi(|q1|=|q2|=|q3|),此外存在两种次谐振波矢Ki,其中|Ki|=2|qi|,=31/2|qi|满足三波共振关系:K1+q2=因此蜂窝斑图有三套子结构构成.图3(b)为波数比3:1时系统1中所呈现的超六边斑图,任何一个单元与周围均匀分布的六个相同单元都可组成正六边形的形状,形成六边形的每个单元与蜂窝超六边形的中心斑图相似,从二维傅里叶频谱中可以看出也是由两组Turing模qi和Ki相互作用形成的(其中|Ki|=3|qi|),共有三套次生模:第一套新生次谐振模有2组,其中3个ki模的大小为2|qi|(|Ki|-|qi|=2|qi|,ki//qi),另3个模的波数为31/2|qi|,且与谐振模满足共振关系:k1+q2=第二套为2组大小相等的新生次谐振模Qi和,分别满足三波共振关系:K1+q2=Q1和-K3-q2=第三套共有4组新生次谐振模Mi,Ni,Ri,其中Mi模的大小为4|qi|(|Ki|+|qi|=4|qi|,Mi//qi),Ri模与第一套新生次谐振模平行且大小为另两组模Ni和的大小相等且满足三波共振关系:K1-q3=N1和-K3+q1=.图3(c)为波数比为4:1时系统1中形成的黑眼斑图,每一个六边形晶胞中心为暗斑点,向外依次是一个亮环和暗环,形成类似于眼睛的结构,每三个单元中心处各有一个浓度不同的点状区域,每一个单元周围有六个这样的点状区域,且这六个点成正六边分布,二维傅里叶频谱显示由两组Turing模qi和Ki相互作用形成的(其中|Ki|=4|qi|,Ki//qi).Turing斑图随系统2中变量扩散系数的增大,超点阵单元间的距离逐渐增大,相同区域面积内的晶胞减少,即超点阵斑图的波数逐渐减小,其原因为系统2中失稳模波数随扩散系数的增加而减小,短波模受失稳模调制所引起,但失稳模与短波模共振耦合时新生次谐振模数量会逐渐增加,且新生模与固有模或新生次谐振模与谐振模均可形成三波共振关系.

3.2 耦合系数的影响

两子系统波数相同时无论耦合系数取0—1的任何值,两系统均出现相同的简单斑图.图4中分别是耦合系数为0.045与0.15时子系统1中所形成简单六边形斑图,可以看出两图形(图4(a)和(b))基本一致,说明此时耦合系数与系统斑图的类型无关,在波数比为1:1时系统始终处于初级Turing分岔点附近.两子系统波数比大于1时,两波模地位不同,短波模受失稳模调制,两波模随耦合系数的改变耦合共振形式也会相应发生变化,因此斑图形式表现出明显的变化.图4(c)和4(d)分别是耦合系数为0.05与0.3时子系统1中所形成复杂斑图与简单六边形斑图,系统1中斑图出现明显的变化说明耦合系数的大小也可以影响系统Turing分岔行为.

图3 (a)波数比2:1时,Du1=6.78,Dv1=15.75,Du2=25.3,Dv2=68.9,α=0.05,β=0.05;(b)波数比3:1时,Du1=6.48,Dv1=15.55,Du2=55.3,Dv2=147.5,α=0.035,β=0.035;(c)波数比4:1时,Du1=1.55,Dv1=3.505,Du2=22.5,Dv2=65.5,α=0.05,β=0.01

图4 耦合强度对斑图的影响 波数比1:1时,Du1=27.75,Dv1=73.75,Du2=27.75,Dv2=73.75:(a)α=0.045,β=0.045;(b)α=0.15,β=0.15;波数比2:1时,Du1=6.5,Dv1=14.15,Du2=22.75,Dv2=65.5:(c)α=0.05,β=0.05;(d)α=0.3,β=0.3

模拟过程中发现在耦合系数小于0.1时,随耦合系数强度的增加系统1中不仅呈现简单和复杂斑图,还可以呈现简单斑图向复杂斑图的变化(如图5(a)—(e)所示,其中1,2代表系统1,2中斑图).可以看出两系统无耦合时子系统2中为大点简单六边形斑图,系统1中呈现小点简单六边形的斑图(图5(a));耦合系数在0—0.0015时两子系统直接形成无耦合时简单六边形斑图(图5(b));耦合系数在0.0015—0.04时,系统2稳定于简单六边形斑图,系统1出现超六边形复杂图形(图5(c)),耦合系数在0.0015—0.02时系统1中超点阵斑图不稳定,随模拟时间延长逐渐演化成简单六边形斑图(图5(1b))或区域点状斑图(图5(2d));耦合系数在0.04—0.065间取值时两子系统同时呈现并稳定于简单四边形斑图(图5(e)),但随耦合系数增大,系统形成正四边形的图形效果会逐渐降低;耦合系数大于0.065时系统只出现短暂的点状斑图,而未能形成稳定的Turing斑图.其中耦合系数在0.025—0.04时,系统1中斑图出现简单斑图到复杂斑图以及简单斑图间的转化(如图6(a)—(e)所示),两子系统会首先同时呈现简单正四边形斑图(图6(a));之后随模拟时间增加系统2中相邻斑点间开始融合,中间背景的面积增大,这种变化在子系统1中表现尤为明显(图6(1b));随模拟时间的继续增加,在系统1增大的背景区中心位置呈现新的斑点,此时系统斑图的形状开始变形(图6(c));随着背景区新生斑点的增多,系统逐渐失去正四边形的形状(图6(d));经过一段时间从无序到有序的演化后系统2开始呈现简单六边形斑图,系统1中复杂斑图单元也向六边形发展(图6(e)所示);最终系统1稳定于超六边斑图,而系统2则稳定于简单六边形斑图.从图形随耦合强度的变化可以看出:在耦合系数很小时失稳模对短波模的调控作用不明显,两系统相互作用较弱,系统斑图未发生改变,随着耦合强度的增大,失稳模对短波模的调控作用逐渐增强,两系统相互作用增大,进而影响到系统所呈现斑图;非线性模拟过程中失稳模对短波模的影响较线性耦合时更为明显[20],且随系统耦合系数增大,在高级Turing分岔与初级Turing分岔相邻区域,系统斑图可以自发地转变.系统斑图出现简单斑图向复杂斑图以及简单斑图间的相互转化,在介质阻挡放电实验过程中通过逐步升高电压的方法已实现[29,30].分析其原因,驱动电压的升高使系统变量扩散浓度增大的同时,壁电荷斑图与放电丝斑图相互作用会随之增强,两子系统间的相互作用增强在模拟过程中可以类比成耦合系数的增大,这与模拟结果相符.

图5 系统斑图随耦合强度的变化 Du1=6.5,Dv1=15.85,Du2=55.3,Dv2=147.5;(a)α=0,β=0;(b)α=0.001,β=0.001;(c)α=0.0375,β=0.0375;(d)α=0.019,β=0.019;(e)α=0.04,β=0.04

图6 系统参数不变斑图随模拟时间自发转变过程 Du1=6.5,Dv1=15.85,Du2=55.3,Dv2=147.5,α=β=0.0375;(a)t=236;(b)t=302.3;(c)t=312.1;(d)t=316.8;(e)t=323

从系统斑图的转变过程可以看出:短波模与失稳模产生共振耦合时,有新的失稳模形成,系统处于初级Turing分岔点附近时,由于系统对失稳模的选择具有惟一性,耦合新产生的失稳模与原失稳模存在竞争,初始失稳模被取代时系统出现新的简单斑图;系统处于高级Turing分岔点附近时,新生失稳模与固有失稳模共存,由于失稳模具有两个不同的波长,新生模与固有模间满足三波共振从而子系统1中产生超六边形复杂斑图.参数不变情况下系统1由简单正四边形向超六边形转变,系统2由四边形向六边形转变,证明系统在稳定外界条件下可以自发由初级分岔向高级分岔转变,在转变时系统的暂态平衡被打破,系统1中的耦合新生失稳模会与固有模形成三波共振,从而子系统1中形成复杂斑图;系统2中由于固有模的对称性高于新生模,系统斑图向高对称性转变,因此固有失稳模会逐步取代新生模,从而斑图由四边形逐渐演化成简单六边形.

4 结论

通过利用双层四变量反应扩散方程,研究了两子系统在非线性耦合强度下Turing模间的相互作用以及复杂斑图的生成机理.模拟结果表明子系统间的波数比值以及耦合强度的大小可以影响两子系统波矢的共振状态,两子系统波矢发生共振耦合会有新的不稳定波矢产生,在高级分岔点附近由于不稳定模具有两个或两个以上不同的波长,且不同波矢间满足三波共振关系,从而引起短波模子系统产生复杂斑图.模拟过程中系统斑图呈现出简单正四边形向六边形以及超六边形点阵斑图转化,这一结果从理论验证了不同斑图间可自发相互转化的实验结论.本文模拟结果为深入研究反应扩散系统超点阵复杂斑图有重要的借鉴作用.

[1]Schenk C P,Or-Guil M,Bode M,Purwins H G 1997 Phys.Rev.Lett.783781

[2]Berenstein I,Dolnik M,Yang L,Zhabotinsky A M,Epstein I R 2004 Phys.Rev.E 70046219

在认真总结多年来对危险货物申报人员和集装箱装箱现场检查员的管理实践经验的基础上，结合党中央、国务院关于建设“诚信中国”的总体要求，建立了“两员”诚信管理制度，进一步加强对“两员”日常从业情况的监管，具体诚信管理制度由交通运输部海事局制定。

[3]Arbell H,Fineberg J 2002 Phys.Rev.E 65036224

[4]Sharpe J P,Ramazza P L,Sungar N,Saunders K 2006 Phys.Rev.Lett.96094101

[5]Bois J S,Jülicher F,Grill S W 2011 Phys.Rev.Lett.106028103

[6]Yang L F,Dolnik M,Zhabotinsky A M,Epstein I R 2002 Phys.Rev.Lett.88208303

[7]Berenstein I,Yang L F,Dolnik M,Zhabotinsky A M,Epstein I R 2003 Phys.Rev.Lett.91058302

[8]Turing A M 1952 Phil.Trans.R.Soc.London B 23737

[9]Dong L F,Fan W L,He Y F,Liu F C,Li S F,Gao R L,Wang L 2006 Phys.Rev.E 73066206

[10]Duan X X,Ouyang J T,Zhao X F,He F 2009 Phys.Rev.E 80016202[11]Stollenwerk L,Laven J G,Purwins H G 2007 Phys.Rev.Lett.98255001

[12]Shirafuji T,Kitagawa T,Wakai T,Tachibana K 2003 Appl.Phys.Lett.832309

[14]Liu C Z,Brown N,Meenan B J 2006 Appl.Surf.Sci.2522297

[15]Dong L F,Liu W L,Wang H F,He Y F,Fan W L,Gao R L 2007 Phys.Rev.E 76046210

[16]Nie Q Y,Ren C S,Wang D Z,Li S Z,Zhang J L 2007 Appl.Phys.Lett.90221504

[17]Cross M C,Hohenberg P C 1993 Rev.Mod.Phys.65851

[18]Barrio R A,Varea C,Aragón J L,Maini P K 1999 Bull.Math.Biol.61483

[19]Kytta K,Kaski K,Barrio R A 2007 Physica A 385105

[20]Berenstein I,Yang L F,Dolnik M,Zhabotinsky A M,Epstein I R 2005 J.Phys.Chem.A 1095382

[21]Bai Z G,Dong L F,Li Y H,Fan W L 2011 Acta Phys.Sin.60118201(in Chinese)[白占国,董丽芳,李永辉,范伟丽2011物理学报60118201]

[22]M´ıguez D G,Dolnik M,Epstein I R,Muuzuri A P 2011 Phys.Rev.E 84046210

[23]Rogers J L,Schatz M F,Brausch O,Pesch W 2000 Phys.Rev.Lett.854281

[24]Ni W M,Tang M X 2005 Trans.Amer.Math.Soc.3573953

[25]Mikhailova A S,Showalter K 2006 Physics Reports 42579

[26]Yuan X J,Shao X,Liao H M,Ouyang Q 2009 Chin.Phys.Lett.26024702

[27]Liu H Y,Yang C Y,Tang G N 2013 Acta Phys.Sin.62010505(in Chinese)[刘海英,杨翠云,唐国宁2013物理学报62010505]

[28]Wang W M,Liu H Y,Cai Y L,Li Z Q 2011 Chin.Phys.B 20074702[29]Dong L F,Li S F,Liu F,Liu F C,Liu S H,Fan W L 2006 Acta Phys.Sin.55362(in Chinese)[董丽芳,李树锋,刘峰,刘富成,刘书华,范伟丽2006物理学报55362]

[30]Dong L F,Yang Y J,Fan W L,Yue H,Wang S,Xiao H 2010 Acta Phys.Sin.591917(in Chinese)[董丽芳,杨玉杰,范伟丽,岳晗,王帅,肖红2010物理学报591917]