Nem izoterm kémiai reakció egyszerű modellje

Egyéb képletek

w = k(T)ϕ(η)
k(T) = (Subscript[k, 0]) exp(-(E/RT))
Φ(η) = (((1-η))^m)
Φ(η) = ((Subscript[η, 0]) + η)(1 - η)

cp ((δT)/(δτ)) = λΔT + Qw
((δη)/(δτ)) = w
((δη)/(δτ)) = DΔη + w

Θ = (E/(Subscript[RT, *]^2))(T - (SubStar[T]))
t = τ (Q/cp) (E/(Subscript[RT, *]^2)) (Subscript[k, 0]) exp(-(E/(SubStar[RT])))
β = ((SubStar[RT])/E)
γ = (cp/Q) ((Subscript[RT, *]^2)/E)

x = (1/r) (x^(,))
y = (1/r) (y^(,))
z = (1/r) (z^(,))
δ = (Q/λ) (E/(Subscript[RT, *]^2)) (r^2) (Subscript[k, 0]) exp(-(E/(SubStar[RT])))

((δθ)/(δt)) = (1/δ)Δθ + ϕ(η) exp (θ/(1 + βθ))
((δη)/(δt)) = γϕ(η) exp (θ/(1 + βθ))
((δη)/(δt)) = μΔη + γϕ(η) exp(θ/(1 + βθ))
μ = (1/δ) · (cpD/λ)

x = (cx^(,))
y = (cy^(,))
z = (cz^(,))
c = ((Q/λ)(E/(Subscript[RT, *]^2))(Subscript[k, 0])exp(-(E/(SubStar[RT]))))

θ = 0 or ((δθ)/(δν)) + σθ = 0
σ = ((αr)/λ)
((δη)/(δν)) = 0

θ = η = 0 when t = 0

((δθ)/(δt)) = (1/δ)Δθ + (e^θ)
θ(Subscript[|, t = 0]) = 0
(((δθ)/(δν)) + σθ)(Subscript[|, S]) = 0

(t^(,)) = γt

(((d^2)θ)/(dx^2)) + ((n + 1)/x) · ((dθ)/dx) + ((δe)^θ) = 0
(((dθ)/dx) + σθ)(Subscript[|, x = 1]) = 0

((δθ)/(δt)) = (1/δ)Δθ + (e^θ)
θ(Subscript[|, t = 0]) = - (Subscript[θ, 0])
θ(Subscript[|, x = 0]) = 0
(((δθ)/(δν)) + σ(θ + (Subscript[θ, 0]))) (Subscript[|, S=0]) = 0

Előkészületek

Első lépésben inicializáljuk a differenciálegyenletekben szereplő szükséges konstansokat. Ezek nagyságrendileg megegyeznek létező anyagok hasonló tulajdonságaival, annyiban módosítva, hogy az időskálán szemmel követhető legyen a folyamat:

E0=1000
R=8.314
k[T_]=k0 Exp[-E0/(R T[x,y,z,t])]
η0=0.2
c=918
ρ=1.429
λ=0.50
Q=10000000
k0=1

eq1=c*ρ*D[T[x,y,z,t],t]==λ*(D[T[x,y,z,t],x,x]+D[T[x,y,z,t],y,y]+D[T[x,y,z,t],z,z])+Q*D[η[x,y,z,t],t]
eq2=D[η[x,y,z,t],t]==k0 Exp[-E0/(R*T[x,y,z,t]) ]*(1-η[x,y,z,t])^2
boundary1={T[x,y,0,t]==5000*(2-(x^2+y^2)),T[x,y,100,t]==300,T[x,-1,z,t]==T[x,1,z,t]==300,T[-1,y,z,t]==T[1,y,z,t]==300,T[x,y,z,0]==1000}
boundary2={η[x,y,z,0]==1-(x^2+y^2)/Sqrt[2]}
eq3=c*ρ*D[T[x,y,z,t],t]==λ*(D[T[x,y,z,t],x,x]+D[T[x,y,z,t],y,y]+D[T[x,y,z,t],z,z])

Összehasonlítás a szimpla hőterjedéssel

sol1=NDSolve[{eq3,boundary1},T,{x,-1,1},{y,-1,1},{z,0,100},{t,0,10}]
sol2=NDSolve[{eq1,eq2,boundary1,boundary2},{T,η},{x,-1,1},{y,-1,1},{z,0,100},{t,0,10}]

Plot3D[{Evaluate[T[x,y,1,4]/.sol1],Evaluate[T[x,y,1,4]/.sol2]},{x,-0.7,0.7},{y,-0.7,0.7},ColorFunction->Hue]

Látható, hogy az anyag égését is figyelembe vevő modell eredménye az, hogy magasabb hőmérsékletet érünk el a hővezetéshez képest, valamint az anyag jelenléte miatt más a hőeloszlás. Nézzük meg ezzel párhuzamosan az anyag fogyását:

Plot3D[Evaluate[1-η[x,y,1,4]/.sol2],{x,-0.7,0.7},{y,-0.7,0.7}]

Látjuk, hogy az anyag itt már majdnem teljesen elfogyott.

Animált modell

A láthatóság kedvéért új peremfeltételeket adunk meg:

boundary3={T[x,y,0,t]==10000*(2-(x^2+y^2)),T[x,y,100,t]==300,T[x,-1,z,t]==T[x,1,z,t]==300,T[-1,y,z,t]==T[1,y,z,t]==300,T[x,y,z,0]==1000}

boundary4={η[x,y,z,0]==0}

sol3=NDSolve[{eq1,eq2,boundary3,boundary4},{T,η},{x,-1,1},{y,-1,1},{z,0,100},{t,0,10}]

Az alábbi animáción a hőmérsékletet ábrázoljuk síkmetszetben az idő függvényében. A sötétpiros szín jelzi a legforróbb állapotot. Jól megfigyelhető a hőmérséklet emelkedése, illetve a szélek és a hasáb közepe közti hőmérsékletkülönbség.

Animate[Plot3D[Evaluate[T[x,y,2,t]/.sol3],{x,-0.7,0.7},{y,-0.7,0.7},PlotRange->{2000,12000},PlotPoints->5,ColorFunction->Hue],{t,0,3,0.01}]

Az alábbi animáción az anyagmennyiséget ábrázoljuk síkmetszetben az idő függvényében. A piros szín jelzi a legkevesebb mennyiséget. Megfigyelhető, hogy kicsi a mennyiségbeli különbség a tartály különböző részein, szabad szemmel alig látható, inkább csak a sznekből következtethetünk rá.

Animate[Plot3D[Evaluate[1-η[x,y,2,t]/.sol3],{x,-0.7,0.7},{y,-0.7,0.7},PlotPoints->5,ColorFunction->Hue,PlotRange->{0,1}],{t,0.5,1.5,0.01}]

Az alábbi animáción a hőmérsékletet ábrázoljuk síkmetszetben a hőforrástól való távolság függvényében. Jól látható, hogy minél távolabb megyünk a hőforrástól, annál kisebb a hőmérséklet.

Animate[Plot3D[Evaluate[T[x,y,z,1.5]/.sol3],{x,-0.7,0.7},{y,-0.7,0.7},PlotPoints->5,ColorFunction->Hue,PlotRange->{6000,15000}],{z,0,4,0.1}]

Forrás

A.I. Vol’pert, S.I. Hudjaev: Analysis in classes of discontinuous functions and equations of mathematical physics