Chemical Kinetics in Objectica

Stephan Leibbrandt

({{Library,  , ←,  , Application,  , Concepts of Ojectica}, { , ↘,  , ↙,  ,  , Atoms}, { ,  , Objectica,  ,  ,  , Molecules}, { ,  , ↓,  ,  ,  , Reaction}, { ,  , Mathematica,  ,  ,  , Reactor}, { ,  ,  ,  ,  ,  , Dry}, { ,  ,  ,  ,  ,  , Wet}, { ,  ,  ,  ,  ,  , Mass Balance}})

Concepts of Ojectica

Standard Package

Get["Class`Class`"]

Class, Object and Access

Class[Person]:={
name=0,
Virtual.kind[]:="someone",
Person[name_]:=This.name=name,
info[]:=kind[]<>" called "<>name
}

per=New.Person["Miller"]

per.info[]

Inheritance and Polymorphism

Class[Employee,Person]:={
Private.salary=0,
Employee[name_,sal_]:=(Super[name];salary=sal),
Virtual.Override.kind[]:="worker with salary "<>ToString[salary]
};

emp=New.Employee["Smith",10000]

emp.info[]

Some Details

Clone

per2=CloneObject[per]

Save for External Use

Export[ToFileName[$TemporaryPrefix,"per"],Person,LicenseID->"L1234-5678"]

Export[ToFileName[$TemporaryPrefix,"per"],per]

Thread

{per,per2}.name

Variable Generation

Class[Var]:={
get[]:=$.x
};
Table[New.Var[],{5}].get[]

Straightforward Operator Usage

Class[Field]:={conts=0};
obj=New.Field[]
obj.conts++
obj.conts

Class[Field]:={conts={1,2,3}};
obj=New.Field[]
obj.conts[[2]]=27
obj.conts

Documentation

Class[Var]:={
generated=DateList[],
"Generate object symbol.".get[]:=$.x
};
Var

New.Var[]

({{Atoms,  }, { ,  }})

Definitions

Off[General::stop]
Class[Atom]:={
Atoms={},
Name="",m=0,ρ=0,(Subscript[c, p])=0,(Subscript[ΔH, f])=0,(Subscript[ΔH, v]):=0,
Private.Atom[name_String]:=(
Name=name;
Atoms={ToExpression[ElementData[Name,"Abbreviation"]<>"->1"]};
{m,ρ,(Subscript[c, p]),(Subscript[ΔH, f]),(Subscript[ΔH, v])}=Map[ElementData[Name,#]&,{"AtomicWeight","Density","SpecificHeat","FusionHeat","VaporizationHeat"}];
),
Static.Private.Create[name_String]:=Set@@{
ToExpression[ElementData[name,"Abbreviation"],InputForm,Unevaluated],
New[Automatic,ElementData[name,"Abbreviation"]].Atom[name]
},
Static[]:=Map[Create,ElementData[]]
};
On[General::stop]

Usage

N.Name

O.m

I.(Subscript[c, p])

Fe.ρ

Fe
I
Object[I]

Re.Atoms

({{Molecules,  }, { ,  }})

Definitions

Class[Molecule]:={
Atoms={},
Private.Molecule[]:=(
Atoms=Rule@@PadRight[#,2,1]&/@Split[ToExpression[StringJoin[#]]&/@Split[Characters[ToString[This]],LowerCaseQ[#2]||DigitQ[#1<>#2]&],NumberQ[#2]&]
),
Static.Get[c_Symbol]:=If[DirectlyObjectQ[c],c,New[c].Molecule[]],
m:=Plus@@Map[Last[#]*First[#].m&,Atoms]
};

Usage

{Molecule.Get[C2H4],Molecule.Get[C2H4],Objects[Molecule]}

C2H4.Atoms

Molecule.Get[H]
%.Atoms

Molecule.Get[N2O3].m

({{Reaction,  }, { ,  }})

Definitions

Class[Reaction]:={
R=8.314472,Compounds={},k=0,
Analyze[side_]:=Block[{p,t},List@@Rule@@@Map[Take[Reverse[Flatten[t[1,#]]],2]&,Flatten[p[side/.{Plus->p,Times->t}]]]],
Reaction[lhs_->rhs_,k_]:=With[{me=This},
FormatValues[me]=FormatValues[me]/.HoldPattern[FormatRules[Object,me,Automatic]]:>FormatRules[Object,me,HoldForm[lhs->rhs]];
This.k=k/.Import.R->R;
Compounds=MapAt[Molecule.Get,#,1]&/@Join[MapAt[-#&,#,2]&/@Analyze[lhs],Analyze[rhs]]
]
};

Coal

Values of constants are only for demonstration!

Coal[Dry]={
New.Reaction[C+O2->CO2,(10^4)*(E^(-(180000/(R*T))))],New.Reaction[CO+1/2O2->CO2,(10^9)*(E^(-(250000/(R*T))))],
New.Reaction[C+CO2->2CO,(10^8)*(E^(-(250000/(R*T))))],New.Reaction[C+1/2O2->CO,(10^4)*(E^(-(180000/(R*T))))]
}

Coal[Wet]=Join[Coal[Dry],{
New.Reaction[H2+1/2O2->H2O,(10^8)*(E^(-(250000/(R*T))))],New.Reaction[CO2+H2->CO+H2O,(10^9)*(E^(-(100000/(R*T))))],
New.Reaction[C+H2O->CO+H2,(10^7)*(E^(-(180000/(R*T))))]
}]

Brute Force: 1-2 Educts + 1-3 Products -> 420 Reactions

rep={a___,m_->c1_Integer,b___,m_->c2_Integer,c___}:>{a,b,c,m->c1+c2};
m=Molecule.Get/@{O,O2,O3,N,N2,N2O,N2O2,NO,N2O3,NO2,N2O4,N2O5}
e=Join[List/@Select[m,Apply[Plus,Map[Last,#.Atoms]]>1&],Union[Sort/@Flatten[Outer[List,m,m],1]]];
p=Join[e,Union[Sort/@Flatten[Outer[List,m,m,m],2]]];
pa=Map[Sort,Join@@@Map[#.Atoms&,p,{-1}]//.rep];
reacs=Apply[Plus,Sort[Flatten[Map[Function[e,Map[e->#&,Select[Extract[p,Position[pa,Sort[Apply[Join,e.Atoms]//.rep]]],FreeQ[#,Alternatives@@e]&]]],e],1]],{2}];
Length[reacs]
TableForm[Partition[reacs,4,4,{1,1},""]]

({{Reactor,  }, { ,  }})

Definitions

Class[Reactor]:={
Reactions={},
Reactor[{reactions__Symbol?InstanceOf::Reaction}]:=Reactions={reactions},
Concentrations[t_]:=C[#][t]&/@Union[First/@Apply[Join,Reactions.Compounds]],
Equations[Temp_,t_]:=Block[{T=Temp},Module[{compounds,none,equations},
compounds=Union[First/@Apply[Join,Reactions.Compounds]];
none=Thread[compounds->0];
equations=Map[Function[c,D[C[c][t],t]==Apply[Plus,Map[(c/.#.Compounds/.none)*#.k*Times@@Map[C[First[#]][t]&,Select[#.Compounds,Last[#]<1&]]&,Reactions]]],compounds];
equations
]],
Solve[Temp_,{t_Symbol,from_,to_},{conc__}]:=Block[{T=Temp},NDSolve[Join[Equations[T,t],Thread[Concentrations[0]=={conc}]],Concentrations[t],{t,from,to}]]
};

Usage

Dry

Combustion[Dry]=New.Reactor[Coal[Dry]]

Combustion[Dry].Equations[1400,t]//TableForm

Wet

Combustion[Wet]=New.Reactor[Coal[Wet]]

Combustion[Wet].Equations[1400,t]//TableForm

Plot Function

Needs["PlotLegends`"]

Draw[type_,conc_,moleQ_,opts___]:=Module[{flows},
flows=Combustion[type].Solve[1400,{t,0,20},conc];
Plot[Evaluate[If[moleQ,1,Map[#[[0,1]]&,Combustion[type].Concentrations[t]].m]*Combustion[type].Concentrations[t]/.First[flows]],{t,0,20},opts,PlotRange->All,PlotLegend->Combustion[type].Concentrations[t],LegendPosition->{1.1,-0.4},PlotStyle->Thread[{AbsoluteThickness[3],{Black,Red,Green,Blue,Orange,Magenta}}],ImageSize->800,PlotLabel->If[moleQ,"Mole","Mass"]]
]

({{Dry,  }, { ,  }})

Combustion[Dry].Concentrations[t]

Stoichiometric

Draw[Dry,{10,0,0,10},True]

Draw[Dry,{10,0,0,10},False]

More Coal

Draw[Dry,{10,0,0,5},True]

More Oxygen

Draw[Dry,{10,0,0,20},True]

({{Very Wet,  }, { ,  }})

Combustion[Wet].Concentrations[t]

Stoichiometric

Draw[Wet,{10,0,0,0,2,10},True]

Draw[Wet,{10,0,0,0,2,10},True,PlotRange->{{0,1},All}]

More Coal

Draw[Wet,{10,0,0,0,2,5},True]

More Oxygen

Draw[Wet,{10,0,0,0,2,20},True]

({{Mass Balance,  }, { ,  }})

flows=Combustion[Dry].Solve[1400,{t,0,20},{10,0,0,10}];
Plot[Evaluate[Map[Function[atom,Apply[Plus,Map[(atom/.#[[1,0,1]].Atoms/.atom->0)*Last[#]&,flows,{2}],{0,1}]],{C,O}]-{10,20}],{t,0,20},PlotStyle->Thread[{AbsoluteThickness[2],{Black,Blue}}]]

flows=Combustion[Wet].Solve[1400,{t,0,20},{10,0,0,0,2,10}];
Plot[Evaluate[Map[Function[atom,Apply[Plus,Map[(atom/.#[[1,0,1]].Atoms/.atom->0)*Last[#]&,flows,{2}],{0,1}]],{C,O,H}]-{10,22,4}],{t,0,20},PlotStyle->Thread[{AbsoluteThickness[2],{Black,Blue,Red}}]]

({{Objectica, }})

({{({{Industry, 1000, Ringstraße 100}, {Education, 350, 52078 Aachen}, { ,  , Germany}, { ,  , Tel +49/241/566 08 28}, { ,  , sales@symbols-and-numbers.com}})}, {→ Special Conference Offer:-20% ←}})

({{abstraction layer, stepwise refinement -> iterative model development}, {OO libraries, black box models -> nonexpert can build models}, {one class template -> many objects, easy to change -> refinement driven by simulation results}, {encapsulation of functionality, incremental development -> bottom up design}})