Stability of reaction fronts in random walk simulations

(1)

simulations NoémiNagy

∗

Feren Izsák

†

January19, 2011 Abstra t

Amodelofpropagatingrea tionfrontsisgivenforsimpleauto atalyti

rea tionsandthestabilityofthepropagatingrea tionfrontsarestudied

inseveralnumeri alexperiments. The orrespondingrandomwalk

simu-lations-extendingofare entalgorithm-makepossiblethesimultaneous

treatment of movingparti les. A systemati omparison withthe

stan-dard deterministi simulations highlight the advantages of the present

sto hasti approa h. The main favor of the random walk simulation is

that the initial perturbation hasnostrong ee tonthe stability of the

frontunlikeindeterministi ases.

1 Introdu tion

Frontpropagation in hemi aland bio hemi alsystems hasbeeninvestigated

inseveralstudiesand the orrespondingmathemati almodelsprovideus

hal-lengingproblems [16℄. In aseof ertain hemi alrea tionsthepropagatingof

rea tionfronts anexhibitinstabilitieswhi his alledshortlyngering. Typi al

examplesareprovidedby ombustion[1℄andauto atalyti rea tions[2℄,[3℄but

inawider ontext,growthofba teriaensembles[4℄,ormotionofbiolms[5℄are

frequentlyshowsimilarbehavior. Thenatureof thisinstabilitiesaredes ribed

independen eofsomephysi alparameterssu hastheRayleighnumber[6℄,the

Lewisnumber[6℄,[7℄, theZeldovi hnumber[8℄andtheDamköhlernumber[9℄,

[10℄.

Abovetothedetailed experimentalstudies,severalmathemati almodelsof

thesephenomenahavebeenelaborated. The hiefquantityto omputehereis

the on entration of the reagentswhi h indi ate dire tly the geometry of the

rea tionfront. However,in aseofrea tionsin the ombustiontheoryalso the

∗

DepartmentofAppliedAnalysisandComputationalMathemati s,EötvösUniversity,

H-1117Budapest,Pázmánysétány1/C,Hungary

†

DepartmentofAppliedAnalysisandComputationalMathemati s,EötvösUniversity,

H-1117Budapest,Pázmánysétány1/C,Hungary,

&

DepartmentofAppliedMathemati s, Uni-versityofTwente,P.O.Box217,7500AEEns hede,TheNetherlands(izsakf s.elte.hu),

(2)

fun tions. If the rea tiono urs in porous media, hydrodynami instabilities

havebeenreportedaswell(seee.g. [7℄and[9℄),wherealsothepressurehastobe

determined. Inthemathemati al ontextrea tionfronts,thestabilityofwhi h

isinvestigated,areidentiedwithtravellingwavesolutionsofthe orresponding

systemof PDE's. The existen e and qualitative properties of these solutions

arestudied in detailforsimpleauto atalyti rea tionse.g. in [12℄and further

generalizedin [13℄. The stabilityof thetravellingwavesisusually determined

usingalinearstabilityanalysis. Themainquestionhereisto al ulatea riti al

value orresponding to the above physi al parameters at whi h the rea tion

front be omes unstable. The on rete omputations lead to time onsuming

numeri alproblems involvingthedis retelevelthesolutionoflargeeigenvalue

problems[10℄,where ertainparameters(des ribingtheshapeandthevelo ity

ofthe front)havebeenapproximated previously. Consequently, the resultsin

deterministi simulationsmaybepollutedheavilywith omputationalerrors.

In real life situations the ngering of the rea tion front [14℄ is evoked by

some inhomogeneities, e.g. in the temperature, in the material stru ture, in

the velo ity eld of the reagentsor in an external (ele tri ) eld. Therefore,

it anbeusefultoin orporate theseu tuations intothemodelandexe utea

orrespondingMonte Carlosimulation. At thesametime, weshould keepthe

omputational osts of the method relatively low and also, a lear

orrespon-den ehastobeprovidedbetweenthematerialparametersandtheonesusedin

thesimulations.

In this work we are fo using to ubi auto atalyti rea tions in a

two- omponentsystem,wherethestabilityofaplanarrea tionfrontdependsonthe

ratioofthediusion oe ientsofthetwospe ies. Thisphenomenonhasbeen

observedin real experiments and the ngering of the rea tionfront has been

reportedinvarious hemi alsystems,see[2℄,[3℄and[15℄. Basedonasto hasti

approa h we providea Monte Carlosimulation formodelingthe instabilityof

thefrontandwedeterminethe riti alratioofthediusion oe ientsatwhi h

therea tionfrontbe omesunstable. Themovingrea tionfrontsaresimulated

in a three dimensional setup and their stability is studied by in orporating a

randomnoiseinthediusionpro ess. Wewillusehereatwodimensional

gen-eralizationofare entee tivealgorithm[17℄,andwepointouttheadvantages

ofthesto hasti simulationin thestabilityanalysisoftherea tionfronts.

2 Stability of the rea tion front

2.1 The rea tion-diusion equations

Weinvestigatea ubi auto atalyti rea tionwhi hisgivenwiththes heme

(3)

tions

a(t, x, y)

and

b(t, x, y)

satisfythesystemofpartialdierentialequations

∂a

∂t

= D

A

∇

2

a

− kab

2

,

∂b

∂t

= D

B

∇

2

b

+ kab

2

,

(1)

for

(t, x, y) ∈ (0, T ) × (0, X) × (0, Y )

,wherethe

D

A

and

D

B

denotethediusion oe ientsoftherea tantsand

k

referstotherea tionrate oe ient.

A ordingtoareal experimentintheiodate-arsenous a idsystem(see [2℄)

weusethefollowingparameterset:

D

A

= D

B

= 2 · 10

−

5

m

2

s

−

1

_,

_k

_{= 9 · 10}

3

M

−

3

s

−

1

_.

Corresponding to the me hanism of the rea tion, we rewrite (1) regarding to

the on entration

ˆb(t, x, y)

(whi h denotesthat ofthe free iodide ions),su h

that(1)transformsinto

∂a

∂t

= 2 · 10

−

5

_∇

2

_a

_{− 9 · 10}

3

_aˆ

_b

2

_,

∂ˆ

b

∂t

=

2 · 10

−

5

ν

∇

2

_{ˆb +}

9 · 10

3

ν

aˆ

b

2

,

(2)

for the details on the me hanism we refer to [2℄. Introdu ing the quantities

˜

a

=

1 ν

· a

and

˜b = ˆb

wearriveattheequation

∂˜

a

∂t

= 2 · 10

−

5

_∇

2

_˜

_a

_{− 9 · 10}

3

_{· ˜}

_a˜

_b

2

_,

∂˜

b

∂t

= D · 2 · 10

−

5

_∇

2

˜b + 9 · 10

3

_{· ˜}

_a˜

_b

2

_,

(3) where

D

=

D

B

D

A

=

1 ν

istheratioofthediusion oe ients. Basedonthesetup ofrealexperimentsweapplied theinitial onditions

˜

a(0, x, y) =

(

0

for

x

≤ x

0

0 _.048

ν

Mfor

x > x

0

,

˜b(0, x, y) =

(

0.048

M

,

for

x

≤ x

0

0,

for

x > x

0

,

(4)

wherethematerialsareinitiallyseparatedat

x

= x

0

. Attheleftandtheright endofthetesttubethe on entrationsof

A

and

B

,respe tively,arekeptxed

duringtherea tion:

˜

a(t, X, y) =

0.048 ν

M and

˜b(t, 0, y) = 0.048

M

.

(5) Itisassumedalsothattheyvanishattheotherend:

a(t, 0, y) = b(t, X, y) = 0.

(6)

Atthe rest of the boundary homogeneous Neumannboundary onditionsare

(4)

proa h

In the literature, the stability of the rea tion fronts are investigated usually

usingthefollowingsteps.

•

The rea tion front is identied with a one-dimensional travelling wave solutionofthePDE(3)withslightlymodiedboundary onditions: The

travellingwave(andthe orrespondingODE)isdenedonthewholereal

axis, su h that instead of (5) and (6) the limits of the two spe ies are

givenat innityandminusinnity,respe tively. Forthedetails werefer

to[12℄.

•

Thetravellingwaveis omputedasasolutionofanODE system. Using realparametersthissystemisusuallystiandaspe ialmethodisrequired

for an a urate approximation of the solution on a ne grid. Also, the

speedoftherea tionfronthastobedetermined. Theanalyti predi tions

forthisquantityexhibitalargedeviationfromthepropervalues,see[12℄.

Re entlytheseboundshavebeensharpenedin[11℄.

•

ANeumanntypestabilityanalysisgivesthestabilityoftheabove travel-ling wave. Here we arelooking for thoseperiodi ( omplex exponential)

perturbationsofthetravellingwavesolutionwhi hsolvetheoriginal

equa-tions. Therealpartoftheexponentsareobtainedbysolvinganeigenvalue

problemwhi h ontainstheapproximatedtravelingfrontsolutionandthe

propagationspeed.

A ordingtotheprevioussteptheeigenvalueproblemshouldbesolvedona

negrid,whi hresultsinalargelinearsystem. Aspe ialmethodisrequiredto

solvethese eigenvalueproblems. Forthe algorithmi details in amoregeneral

ontextwereferto[10℄.

2.3 Deterministi simulation

Tobypassthe di ulties aboveweapplied adire t, spatially two-dimensional

simulation of the dynami s of the rea tion front. Based on (3)-(6) we have

simulatedthepropagationoftherea tionfrontrstinadeterministi way. The

rea tionspa ewasdis retizedinagridof

63×L

y

points,wheretheunitdistan e orrespondsto

0.0625

mand

L

y

,thenumberofgrid olumns maydependon thenumberoftimesteps(seeFig. 1). Inea hgridpointweapproximatedthe

reagents' on entrations.

Theinitial onditions orrespond to(4): in ea hgrid ellat theleft handside

of the grid onsisting of

63 × 5

points we have initiated the simulation with

b(0, x, y) = 0.048

M,andontherighthandside onsistingof

63 × L

y

− 5

points wedene

a(0, x, y) = 0.048D

M.Asasmallperturbation,wehaveadditionally

xed

b(t, x, y) = 0.048

MofBinthemiddle

17

unitgridsofthe

6

thgrid olumn.

Followingamethodoflineste hnique,anexpli itEulermethodwasusedforthe

(5)

Investigating the stability of the rea tionfront it has been reported in many

arti les (see e.g. [2℄, [3℄, [18℄) that there exist a riti al value

D

c

su h that for

D < D

c

therea tionfrontbe omesunstable,while for

D > D

c

itremains stable. We onrmthis observation: ifwe set upthedeterministi simulation

with

D

= 0.1

and

D

= 0.2

theinitial perturbation wasgrowing,a ngering

ofthe rea tionfront evolved, whileat theparametervalue

D

= 0.3

theshape

of therea tionfront didnot exhibit any majorvariation; moreover,in along

runtheperturbationhasbeendisappeared. Therefore,usingtheparameterset

aboveand perform moreruns we obtainthe empiri alvalue

D

c

∼

= 0.25

. Note thatthemagnitudeoftheu tuationswasdierentforthedierentratioofthe

diusion oe ients. Theresultofthesimulationsareshownin Figures1-2.

0

10

20

30

40

50

0

15

30

45

60

0

30

60

90

120

150

0

15

30

45

60

Figure1: Contourlinesof the on entrationof

B

in deterministi simulations

yieldingunstablerea tionfronts. Left: ontourlinesofthe on entrationsat

b

=

0.025 · D

usingthediusion oe ientratio

D

= 0.1

afterthetimesteps

t

max

=

15000, 45000

and

75000

,respe tively. Right: ontourlinesofthe on entrations at

b

= 0.025 · D

using the diusion oe ient ratio

D

= 0.2

after time steps

t

max

= 15000, 45000

and

75000

,respe tively.

As an important output of the simulations we have also determined the

speed of the rea tion fronts. This quantity arises not only in the theoreti al

investigations (see the explanation in 2.2 above), but also in the pra ti e: in

the ombustiontheorytheestimationofthisparameterisof entralimportan e

(6)

0

50

100

150

200

250

0

15

30

45

60

0 ₁₅

30 ₄₅

60

0

30

60

0

0.03

0.06

Figure 2: Left: ontourlines of thedeterministi simulation yielding astable

rea tionfrontat

b

= 0.025 · D

D

= 0.3

after

timesteps

t

max

= 15000, 45000

and

75000

. Right: athree dimensionalviewof the on entration prole of

B

orresponding to the unstable ase in Figure 1

leftafter75000timesteps.

3 Stability of rea tion fronts using a sto hasti

approa h: modeling

We have applied also a sto hasti model a ording to the rea tion-diusion

system in (3)-(6). The benet of this approa h is that we an in orporate

theu tuations whi h are present during the rea tionand areresponsiblefor

the ngering in the experimental observations. Also, the random deviations

duringthe rea tionsareusually mu hsmallerthanthe initialdeviation in the

usualsimulations. Inthetheory,theappli ationoftheNeumanntypestability

analysispresumesthepresen eofaperturbationwithanarbitrarywavelength,

Table1: Dynami s of themotionof the rea tionfrontin the numeri al

simu-lationsfor somevaluesof thediusion oe ientratio

D

. At ertain number

oftimestepsthefrontpositionisgivenin theunits ofthenumeri algrid. The

orrespondingfrontvelo ityisgivenin m/s.

D 0.05 0.1 0.15 0.2 0.25 0.3

velo ity

7.04 · 10

(7)

Motivated by the aboveaspe ts, we des ribe herea sto hasti model and

developane ientsimulationoftheunderlying pro essusingageneralization

oftheMonteCarlosimulation,whi hintrodu edre entlyin [17℄.

The model is adis retesto hasti one, i.e. we onsider anite numberof

parti les, whi h are denoted again with

a

and

b

. These are distributed in the

subdomains of rea tion spa e

(0, X) × (0, Y )

, whi h is dis retized again into

uniform ells.

Therea tionoftheparti leso ursinadeterministi waya ordingto(3):

a(t + δ

t

, x, y) = a(t, x, y) − δ

t

k

D

· a(t, x, y)b

2

(t, x, y)

b(t + δ

t

, x, y) = b(t, x, y) + δ

t

k

D

· a(t, x, y)b

2

(t, x, y),

(7)

where

δ

t

denotesthetimestep,

k

D

isthedis reterea tionrate oe ient,and

x

and

y

denotethe oordinatesofthemidpointsofthe ells.

Thediusionof theparti lestakespla e a ordingto asymmetri random

walk,wheretheparti les anjumpintotheneighboring ell(intofourdire tions

denoted withtheve tors

(±δ

x

,

0)

and

(0, ±δ

y

)

)orresideat theiroriginal site. Aslongasthediusionisisotropi ,one anassumethat

δ

x

= δ

y

. Formally,for allsites

(x, y)

wehave

P

p

((t + δ

t

, x

+ δ

x

, y)|(t, x, y)) = P

p

((t + δ

t

, x

− δ

x

, y)|(t, x, y))

= P

p

((t + δ

t

, x, y

+ δ

y

)|(t, x, y)) = P

p

((t + δ

t

, x, y

− δ

y

)|(t, x, y)),

(8)

where

P

p

(t, x, y) > 0

denotestheprobabilitythattheparti le

p

residesin

(x, y)

at time

t

. It is lear that the number of the parti les jumping between the

neighboring ellshasabinomialdistribution. Usingthisobservation,themotion

oftheparti les anbe simulatedsimultaneously, andpush moreparti lesinto

thesamedire tion.

3.1 A Monte Carlo simulation

Thepro esswasinitiatedbypla ing

1.5·10

5

parti lesoftype

B

inea hgridpoint

inthelefthandsideof thegrid onsistingof

63 × 5

points. Into theremaining

grid points

1 _.5·10

5 ν

parti les of type

A

are pla ed. This number of parti les providesagoodbalan e: we anavoidhighlo al on entrationdeviationsarising

from the sto hasti approa h and the simulations an be exe uted within a

reasonabletime. For a perfe t omparison with the deterministi simulation

a small perturbation has been applied: we haveadditionally pla ed

1.5 · 10

5

parti lesof

B

in the middle

17

ellsofthe6thgrid olumn. A ording to the

initial onditionsin (4)anumberof

1 _.5·10

5

0 _.048

parti les orrespondsto1M.Then the in rease

k

D

1 ν

(1.5 · 10

5

)

3

in the number of parti les yields an in rease of

0 .048

1 _.5·10

5

k

_D

1 ν

(1.5 · 10

5

)

3

M,andtherefore,

0.048 1.5 · 10

5

k

D

1 ν

(1.5 · 10

5

)

3

M

= 9 · 10

3

·

1 ν

0.048

3

M

(8)

k

D

= 9 · 10

3

0.048

2

(1.5 · 10

5

₎

2

≈ 9 · 10

3

· 10

−

13

_.

(9)

Usingunittimestepsagainwegivedire tlythetimestepping orrespondingto

(7),(8)and(9):

a(t + δ

t

, x, y) = a(t, x, y) + δ

t

∇

2 d

a(t, x, y) − 9 · 10

−

10

· b

2

(t, x, y)a(t, x, y) ,

b(t + δ

t

, x, y) = b(t, x, y) + δ

t

∇

2 d

b(t, x, y) + 9 · 10

−

10

· b

2

(t, x, y)a(t, x, y) .

(10) Here

∇

2 d

yieldsthedis retediusionoperatorgivenin(8)while

a

and

b

refer to thenumberof parti les, whi h are given in the ellsat ea h timestep, i.e.

x

∈ {1, 2, . . . , L

y

}, y ∈ {1, 2, . . . , 63}

and

t

= 1, 2, . . . , t

max

. Thedis retediusionoperator

∇

2 d

isgivenasaspe ialtypeofrandomwalk basedontheglobalrandomwalkalgorithmin[17℄. Inthispro ess,inea hnode

apartoftheparti lesjumps totheneighboringnodes,andtherestremainsin

theprimarynode. Thisproportionisgivenwith

r

A

, r

B

∈ [0, 1]

forthespe ies

A

and

B

,respe tively: in expe tation

⌈(1 − r

A

) · N ⌉

ofparti les

A

stayinanode, ifwehad there

N

parti lesoftype

A

originally. Similarrelationholds for

r

B

. Inthedis retelevelthediusion oe ients

D

A

and

D

B

aregivenby

D

A

=

r

A

δ

2 x

4δ

t

D

B

=

r

B

δ

2 x

4δ

t

,

(11)

see[17℄. Consequently,the hoi e

r

A

=

2 · 10

−

5

· 4δ

t

0.0625

2

and

r

B

=

D

· 2 · 10

−

5

· 4δ

t

0.0625

2

,

ontherighthandsideof (10)givesthe orresponden ebetweentheparameters

in(3)and (10). If thenumberofthemovingparti lesis determined,weapply

aMonte Carlosimulationto s atter the parti lesinto the neighboring nodes.

Using the

(i, j)

th entry of the random matrix, we rst al ulate the number

oftheparti les,whi h movehorizontally(therest ispushedverti ally)andwe

determinethenumberofthose,whi haremovingtotheright(therestispushed

to the left). The simulations havebeenexe uted using MATLAB. As behind

therea tionfront we have

D

· 150000

parti les in ea h node, the ontourline

representing the rea tion front has been omputed at

0.5 · D · 150000

. This

orrespondstothelevelappliedinthedeterministi simulations.

Inthese ondseriesofnumeri alexperimentswedidnotperturbtheinitial

onditions,thepro esswasinitiatedbypla ing

1.5 · 10

5

parti lesoftype

B

into

the rst ve grid olumns. This orresponds to the initially planar rea tion

frontsin the real experiments. This ase providesalso abetter experimental

studyoftheinstability: we antra ktheevolutionoftherea tionfrontdriven

alsoby the random u tuations. Theresults of both series of simulationsfor

(9)

0

50

100

150

200

250

0

12

24

36

48

60

0

50

100

150

200

250

0

15

30

45

60

Figure3: Contourlinesofthe on entrationof

B

in thesto hasti simulations

D

= 0.3

. The ontourlinesofthe

on entra-tionsaredrawnat

b

= 0.5 · D · 150000

aftertimesteps

t

max

= 15000, 45000

and

75000

, respe tively. Left: small initial perturbation has beenapplied. Right: thesimulationwithoutanyinitialperturbation.

4 Results and dis ussion

The stability in the sto hasti ase is studied out based on the variation of

the rea tionfront. For a true omparison between the deterministi and the

sto hasti approa h,wehave hosenthisquantitytomeasuretheinstabilityalso

inthedeterministi ase. Wehaveappliedthesametypeofinitialperturbations:

the on entration/numberof parti les at themiddle 17 gridpointsin the 6th

grid olumn have been in reasedto the level in the rst5 grid olumns. The

variationofthepositionof ontourlineatthelevel

b

= 0.5·D ·0.048

after15000,

45000,75000and90000timestepsinthedeterministi aseare omparedwith

the ones at

b

= 0.5 · D · 150000

in the sto hasti ase. The summary of the

simulationsfor several valuesof thediusion oe ient ratio

D

anbefound

in Tables 2 and 3, respe tively. In ase of the parameter values

D

= 0.15

and

D

= 0.2

thedeterministi simulationsexhibitamoreintensivein reaseof

the standardvariation of the rea tionfront ompared to the sto hasti ones.

An interpretation of this observation is that the presen e of inhomogeneities

stabilizetherea tionfront.

Initially, in theabsen e of anyperturbation, the standardvariation of the

(10)

0

30

60

90

120

150

0

15

30

45

60

0

30

60

90

120

150

0

15

30

45

0

60

B

D

= 0.2

b

= 0.5 · D · 150000

aftertimesteps

t

max

= 15000, 45000

and

75000

grow initially by random u tuations [14℄. This anbeperfe tly modelled in

the framework of the sto hasti simulations. In ase of stability this growth

is stopped or anturn into ade rease, while in aseof an unstable frontthis

growthgoeson. Thevariationoftherea tionfrontisshownafter15000,45000,

75000and90000 timesteps at severalvalues ofthe diusion rateratio

D

. In

ase of the deterministi simulations, it is ne essary to in lude the arti ial

initialperturbationgivenabove(anyway,theplanarfrontisonlymoving). The

resultsofthe omputationsforthis aseareshowninTable3. Forthesto hasti

simulation,the orrespondingresults anbe foundin Table2. Similarresults

without any initial perturbations are presented in Table 4. One an observe

thattheinitialperturbationshavenosigni antinuen eonthegrowthofthe

variationoftherea tionfront: wefore aststabilityat thesamevaluesof

D

in

bothsto hasti ases.

Toobtainanevenmoredetailed omparisonbetweenthestabilityprovided

bythe two dierentapproa hes,the ee t of randomperturbations are

inves-tigated at a xed rate of the diusion oe ients. Sin e this behavior is of

interestat theinterfa e ofthe parametersets yieldingstable and unstable

re-a tionfronts, we(arbitrarily) hoose su h a parameter:

D

= 0.17

. With this

(11)

0

10

20

30

40

50

0

15

30

45

60

0

10

20

30

40

50

0

15

30

45

60

B

D

= 0.1

b

= 0.5 · D · 150000

aftertimesteps

t

max

= 15000, 45000

and

75000

standardvariationsinthedeterministi andsto hasti simulations. Inall ases

wehave hangedthe6th olumnofthegridasfollows:

•

P0: noinitialperturbationhasbeenapplied.

•

P1: in the middle7 ellsthe numberofparti les havebeen in reasedto

1.5 · 10

5

in thesto hasti andto0.048inthedeterministi ase.

•

P2: inthemiddle17 ellsthenumberofparti leshavebeenin reasedto

1.5 · 10

5

•

P3: in the ellsfrom 10to 20andin theones from43 to53 thenumber ofparti leshavebeenin reasedto

1.5 · 10

5

in thesto hasti andto0.048

inthedeterministi ase.

•

P4: inthemiddle45 ellsthenumberofparti leshavebeenin reasedto

1.5 · 10

5

TheresultsarepresentedinTables6and5forthe orrespondingsto hasti

anddeterministi simulations,respe tively.

Theresultsinthetablesshowthatthestabilityin aseofthedeterministi

(12)

Table2: Standardvariationoftherea tionfront(the

x

oordinateofthe ontour

lineatthelevel

0.5 · D · 150000

)after15000,45000,75000and90000timesteps

inthesto hasti simulationatseveralvaluesof

D

withtheperturbationofthe

entral17gridsinthe6thgrid olumn. Thevalueshavebeenaveragedin four

independentruns. Theunit inthetable orrespondsto

6.25 · 10

−

5

m.

Diusion oe ientrate(

D

)

Timesteps 0.05 0.1 0.15 0.2 0.25 0.3

15000 531 888 1055 899 739 518

45000 1067 2812 1569 1082 834 563

75000 2114 4178 1678 1124 876 620

90000 2718 4618 1723 1105 913 648

x

lineat thelevel

0.5 · D · 0.048

)after 15000,45000,75000and90000time steps

inthedeterministi simulationatseveral valuesof

D

withtheperturbationof

the entral17gridsin the6thgrid olumn.

D

)

Timesteps 0.05 0.1 0.15 0.2 0.25 0.3

15000 488 676 869 719 663 192

45000 685 2282 2569 1842 1394 659

75000 1038 4304 4116 2950 1832 999

90000 1271 5121 4666 3312 2112 1025

initial perturbation. This an be a real barrier if the instability has to be

predi ted: dependingontheinitialperturbation(whi hhastobeappliedinany

ase)therateof theinstability anbeverydierent. Also, the orresponding

theory predi ts the exponential growth of the perturbations in the rea tion

front. Inrealexperiments,this anholdonlyforashorttime,theperturbations

grow afterwards only slowly. At the same time, the dieren e between the

variations in Table 6 arises only from the dierent standard variation of the

initialperturbations.

In the framework of the dis rete sto hasti model dis ussed herewe have

presentedMonteCarlosimulationswhi hfor ertainparametervaluesstabilize

therea tionfrontin ontrasttothedeterministi simulationsandprovidemore

realisti results: theinstabilitydoesnotdependontheinitialperturbation.

A knowledgements

(13)

x

lineatthelevel

0.5·D·150000

)after15000,45000,75000and90000timestepsin

thesto hasti simulationatseveralvaluesof

D

withoutanyinitialperturbation.

D

)

Timesteps 0.05 0.1 0.15 0.2 0.25 0.3

15000 14 309 655 562 452 451

45000 299 1754 1231 709 521 475

75000 837 2655 1289 760 629 484

90000 1240 2954 1316 792 643 523

Table5: Standardvariationofthepositionof therea tionfront(the

x

oordi-nateofthe ontourlineatthelevel

0.5·D·150000

)after15000,45000,75000and

90000timestepsinthedeterministi simulationsforsomeinitialperturbations

at

D

= 0.17

. Perturbation Timesteps P0 P1 P2 P3 P4 15000 0 797 810 960 655 45000 0 2371 2370 1912 1431 75000 0 3806 3667 1967 1512 90000 0 4403 3984 1943 1528 Referen es

[1℄ S.B. Margolis, H.G. Kaper, G.K. Leaf, B.J. Matkowsky, Bifur ation of

pulsatingandspinningrea tionfrontsin ondensedtwo-phase ombustion.

Combust.S i.andTe hnol.43(1985),127-165.

[2℄ D. Horváth, K. Showalter, Instabilities in propagating rea tion-diusion

frontsoftheiodate-arsenousa idrea tion,J.Chem.Phys.103(1994)

2471-2478.

[3℄ D. Horváth, Á. Tóth, Diusion-driven front instabilities in the

hlorite-tetrathionaterea tion,J.Chem.Phys.108(1998)1447-1451.

[4℄ J. Müller, W. van Saarloos, Morphologi al instability and dynami s of

fronts in ba terialgrowthmodels withnonlinear diusion, Phys.Rev. E:

Stat.Nonlin.SoftMatter.Phys.65,(2002)No.061111.

[5℄ J.Do kery,I.Klapper,Fingerformationinbiolmlayers,SIAMJ.ofAppl.

Math. 62(2001)853-869.

[6℄ M.Rappaz,M.Bellet,M.Deville,Numeri alModelinginMaterialsS ien e

(14)

Table6: Standardvariationofthepositionof therea tionfront(the

x

oordi-nateofthe ontourline atthelevel

0.5 · D · 150000

)after15000,45000,75000

and90000timestepsinthesto hasti simulationsforsomeinitialperturbations

at

D

= 0.17

. Theresultsyieldsanaverageof fourindependentruns.

Perturbation Timesteps P0 P1 P2 P3 P4 15000 648 967 1004 1113 905 45000 975 1267 1240 1448 1187 75000 1030 1339 1373 1531 1291 90000 1043 1363 1414 1558 1322

[7℄ S. Kalliadasis, J. Yang, A. De Wit, Fingering instabilities of exothermi

rea tion-diusion frontsin porous media, Physi s of Fluids16 (5),

1395-1409,2004.

[8℄ K.Allali,A.Du rot,A.Taik,V.Volpert,Conve tiveInstabilityofRea tion

FrontsinPorousMedia,Math. Mod. Nat.Phenom. 2(2007)20-39.

[9℄ A.DeWit, Fingeringof hemi alfrontsinporousmedia,Phys.Rev.Lett.

87(2001),No.054502.

[10℄ J. Yang, A. D'Onofrio, S. Kalliadasis, A. De Wit, Rayleigh - Taylor

in-stability of rea tion-diusion a idity fronts, J. Chem. Phys. 117 (2002),

9395-9408.

[11℄ , Qi, Y.: Travelling fronts of rea tion diusion systems modeling

auto- atalysis,Dis reteContin.Dyn.Syst.Ser.A,7thAIMSConferen e,suppl.,

(2009),622629.

[12℄ J. Billingham, D.J. Needham, The Development of Travelling Waves in

Quadrati and Cubi Auto atalysiswithUnequal DiusionRates. I.

Per-manent Form Travelling Waves,Phil. Trans. R. So . Lond.A 334 (1991)

1-24.

[13℄ D.J. Needham, J.H. Merkin, The ee ts of geometri al spreadingin two

andthreedimensionsontheformationoftravellingwavefrontsinasimple,

isothermal, hemi alsystem,Nonlinearity5(1992),413-452.

[14℄ M.Bö kmann,Fingeringofanas endingfrontin aHele-ShawCell,

avail-ableathttp://www.uni-magdeburg.de/abp/resear h/ d / d .htm

[15℄ A.Hanna,A.Shaul,K.Showalter,DetailedStudiesofPropagatingFronts

in the Iodate Oxidationof ArsenousA id, J.Am.Chem. So .104(1982)

(15)

(2000),161230.

[17℄ C.Vamos,N.Su iu,H.Veree ken,Generalizedrandomwalkalgorithmfor

thenumeri almodelingof omplexdiusionpro ess,J.Comput.Phys.186

(2003)527-544.

[18℄ J.H. Merkin, I.Z. Kiss, Dispersion urves in the diusional instability of

auto atalyti rea tionfronts,Phys.Rev.E 72,(2005)No.026219.

[19℄ J.Fort,D.Campos,J.R.González,J.Velayos,Boundsforthepropagation