Dynamic behavior of non-isothermal 
chemical system
Authors' Data: Housam Binous & Nasri Zakia
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :
S. K. Scott and A. S. Tomlin, Period doubling and other complex bifurcations in non-isothermal chemical systems, Phil. Trans. R. Soc. Lond. A (1990) 332, 51-68.

Governing equations. α, β and ϕ are dimensionless concentrations of species A and B and dimensionless temperature. μ is a bifurcation parameter.

sol[μ_]:=NDSolve[{D[α[t],{t,1}]== μ Exp[δ ϕ[t]]-α[t] β[t]^2-(Subscript[κ, u])α[t],
 D[β[t],{t,1}]== α[t] β[t]^2+(Subscript[κ, u])α[t]-β[t],
D[ϕ[t],{t,1}]== β[t]-γ ϕ[t],α[0]== 1,β[0]== 0.5,ϕ[0]== 0},{α[t],β[t],ϕ[t]},{t,0,10000},MaxSteps-> 1000000]

Time series for μ = 0.498. Simple period one behavior is observed.

(Subscript[κ, u])=10^-3;γ=0.05;δ=0.025;

f=sol[0.498]

Plot[β[t]/.f,{t,500,800},PlotRange-> All]

Time series for μ = 0.5586. Chaotic behavior is obtained.

f=sol[0.5586]

Plot[β[t]/.f,{t,200,500},PlotRange-> All]

Time series for μ = 0.5586. Spiral attractor is observed

ParametricPlot3D[{α[t],β[t],ϕ[t]}/.f,{t,200,500},PlotRange->{{0,25},{0,25},{0,25}},ViewPoint->{0,5,1}]

Time series for μ = 0.558. Exponential decay is altered is this non-isothermal system.

(Subscript[κ, u])=10^-3;γ=0.05;δ=0.025;

f=sol[0.558]

Plot[ϕ[t]/.f,{t,200,500},PlotRange-> All]

Time series for μ = 0.695. Chaotic behavior is observed

(Subscript[κ, u])=5.5 10^-3;γ=0.5;δ=0.1;

f=sol[0.695]

Plot[β[t]/.f,{t,2215,2265},PlotRange-> All]

Time series for μ = 0.695. chaotic attractor is observed

ParametricPlot3D[{α[t],β[t],ϕ[t]}/.f,{t,200,500},PlotRange->All,ViewPoint->{0,5,1},PlotPoints-> 10000]

Computing the next maxima map for μ=0.695. We get a plot that is typical of such chaotic attractors.

k=0;i=0;

solt[0]=FindRoot[D[First[β[t]/.f],{t,1}]==0,{t,219,220}]
pt[k]=First[β[t]/.f]/.solt[0]
(D[First[β[t]/.f],{t,2}]/.t->219.069)

pt[0]

k=1;

ta=.

While[i<2501,{t0=220+i,i+=1,solt[i]=Quiet[FindRoot[D[First[β[t]/.f],{t,1}]==0,{t,t0,t0+1}]],If[(D[First[β[t]/.f],{t,2}]/.solt[i])<0 && t0<( t/.solt[i])< t0+1,{ta[k]=t/.solt[i];pt[k]=First[β[t]/.f]/.solt[i]}],If[pt[k]≠ pt[k-1],k+=1]}]

k

tbl1=Table[{ta[k],pt[k]},{k,1,204}]

tbl=Table[{pt[i],pt[i+1]},{i,1,202}]

ListPlot[tbl,PlotRange->{{0,10},{0,10}}]