Simulating Radiation Chemistry with Mathematica

Chemical Reaction Equations Into ODEs

Karen Colins

Atomic Energy of Canada Limited

Prologue

In water-cooled nuclear reactors, the phenomena of radiolysis, chemical reactions, and transport take place simultaneously in irradiated, aqueous environments.

Computer models of these phenomena can be used to calculate the evolving concentrations of radiolysis products and radioactive species.  Developed over many decades, such models are the subject of ongoing correction and exhibit increasing sophistication and complexity.

Primary variables of interest

Primary variables of interest are the concentrations of hydrogen gas and iodine-containing compounds.

Both substances are potentially hazardous and must be accounted for in safety analyses and systems design.

Atomic Energy of Canada Limited (AECL)

AECL develops and produces (CANDU^®) heavy water reactors, in which radiation chemistry is important and must be taken into account.

Chalk River Laboratories

Computer program for radiation chemistry

Grid[({{Style["LIRIC",Bold,Blue,Larger], Style["Radiation-chemical kinetics program, >450 reactions",Smaller]}, {Style["↑",Blue,Larger], " "}, {Style["FACSIMILE",Bold,Blue,Larger], Style["Automatic differential equation solver \nand problem description language",Smaller]}, {Style["↑",Blue,Larger], " "}, {Style["FORTRAN",Bold,Blue,Larger], Style["High-level programming language for problems \nthat can be expressed algebraically",Smaller]}}),Alignment->{{Center,Left},Center},Spacings->2]//Style[#,FontFamily->"Arial"]&

How LIRIC works

The programmer first takes a set of >450 coupled chemical reaction equations

{(Subscript[H, 2])N(Subscript[O, 3]) (Underoverscript[⇌, (1.5×(10^10) (((s^(-1))(mol·(dm^(-3))))^(-1))), (2.4×(10^5) (s^(-1)))]) HN(Subscript[O, 3]^(-)) + (SuperPlus[H]), …}

and translates them into a seqence of FACSIMILE language statements

% 2.4E5 % 1.5E10: H2NO3 = H+ + HNO3-;
⋯

to create the LIRIC program library.

When LIRIC executes, the FACSIMILE engine converts all of the translated chemical reaction equations into a set of Ordinary Differential Equations (ODEs), with time as the independent variable, and solves them numerically.

Physical systems simulated in LIRIC

Organizational structure of LIRIC

Grid[({{Style["SUBMODELS",Bold,Blue], Style["ADDITIONAL SETS OF EQUATIONS",Bold,Blue]}, {Style["Water radiolysis",Smaller], Style["Hydrogen leakage",Smaller]}, {Style["Inorganic iodine chemistry",Smaller], Style["Acid addition",Smaller]}, {Style["Radiolytic iodine reactions",Smaller], Style["Transport",Smaller]}, {Style["Generic metal ion radiolysis",Smaller], Style["Acid-base equilibria:",Smaller]}, {Style["Iron ion radiolysis",Smaller], Style["  ∘ Water",Smaller]}, {Style["Copper ion radiolysis",Smaller], Style["  ∘ Water radiolysis products",Smaller]}, {Style["Organic impurity radiolysis",Smaller], Style["  ∘ Carbonate",Smaller]}, {Style["Organic iodide production and destruction",Smaller], Style["  ∘ Nitric Acid and products",Smaller]}, {Style["Nitrate radiolysis",Smaller], Style["  ∘ Ammonia",Smaller]}, {Style["Ammonia radiolysis",Smaller], Style["  ∘ Hydrazine",Smaller]}, {Style["Hydrazine radiolysis",Smaller], Style["  ∘ Phosphate",Smaller]}, {Style["Nitric acid production by air radiolysis",Smaller], Style[" ",Smaller]}}),Alignment->Left,Dividers->{{False,True},{False,True,{False}}},Spacings->2]//Style[#,FontFamily->"Arial"]&

Submodels and additional sets of equations are many and complex!

Nitrate radiolysis submodel of LIRIC

LIRIC theory

Problem: not explicit

The basis of LIRIC is an engineered collection of submodels etcetera, not a comprehensive theory

Theory is implicit to some extent in FACSIMILE, in terms of the symbolic equations for generating Ordinary Differential Equations (ODE) from chemical reaction equations

Theory is implicit to a large extent within LIRIC itself, in terms of submodels and sets of additional equations and physical systems.

Question:  How to synthesize an explicit theory for LIRIC?

Answer:  By partitioning LIRIC phenomenologically, in as general terms as possible.

Phenomenological expression of theory

The concentration of a chemical species depends on time t and physical system σ.

ω≡ω(t,σ)=time and system dependent concentration of chemical species Ω

The temporal derivative of the concentration of a chemical species can be expressed in terms of its additive, phenomenological components.  (Phenomena are assumed to affect chemical species concentrations independently.)

(∂_t) ω=(Subscript[((∂_t) ω), ho])+(Subscript[((∂_t) ω), ra])+(Subscript[((∂_t) ω), he])+(Subscript[((∂_t) ω), pp])+(Subscript[((∂_t) ω), gr])

Where

Grid[{{HoldForm[(Subscript[((∂_t) ω), ho])],"= homogeneous chemical reaction term"},{HoldForm[(Subscript[((∂_t) ω), ra])],"= radiolysis term"},{HoldForm[(Subscript[((∂_t) ω), he])],"= heterogeneous chemical reaction term"},{HoldForm[(Subscript[((∂_t) ω), pp])],"= partial pressure term"},{HoldForm[(Subscript[((∂_t) ω), gr])],"= gravitation term"}},Alignment->Left]//Style[#,FontFamily->"Arial"]&

Each of these terms can be given explicit mathematical expression.

Mathematical expressions of LIRIC theory

(Subscript[((∂_t) ω), ho])=-(∑_j^( ))(((Subscript[α, j])-(Subscript[β, j]))((Subscript[K, fj])(∏_i^( ))(Subscript[ω, ji]^(Subscript[α, ji]))-(Subscript[K, bj])(∏_i^( ))(Subscript[ω, ji]^(Subscript[β, ji]))))

(Subscript[((∂_t) ω), ra])=(Subscript[y, p]) ξ (Subscript[ρ, r])

(Subscript[((∂_t) ω), he])=-(∑_h^( ))(((Subscript[α, h])-(Subscript[β, h]))((Subscript[K, ADh])(∏_i^( ))((((Subscript[s, i]^(Subscript[δ, i])) (c^(δ-(Subscript[δ, i])))(Subscript[ω, i])))^(Subscript[α, hi]))-(Subscript[K, DSh])(∏_i^( ))((((Subscript[s, i]^(Subscript[δ, i])) (c^(δ-(Subscript[δ, i])))(Subscript[ω, i])))^(Subscript[β, hi]))))

(Subscript[((∂_t) ω), pp])=(((((-1)/(Subscript[V, aq]))))^(1-η))(((Subscript[k, aq])(Subscript[k, g])(Subscript[A, int]))/(Subscript[V, g]^η))(([Ω(aq)]-[Ω(g)])/(H (Subscript[k, aq])+(Subscript[k, g])))

(Subscript[((∂_t) ω), gr])=-(α-β)((Subscript[K, DO])(∏_i)(Subscript[c, gr]^(ε-(Subscript[ε, i]))) (Subscript[ω, i]^(Subscript[α, i])))

Upon examination of these expressions, chemists demand, "Show us the chemistry!"

This is a job for Mathematica

Here is a small fraction of the homogeneous chemical reactions (reduced iodine set).

({{"Chemical Reaction Equation", "Forward Rate Constant", "Backward Rate Constant"}, {" (Subscript[I, 2]) + (SuperMinus[OH]) ⇌ (Subscript[I, 2])(SuperMinus[OH]) ", K3, Kneg3}, {" (Subscript[I, 2])(SuperMinus[OH]) ⇌ HOI + (SuperMinus[I]) ", 1.0×(10^7), 6.7×(10^9)}, {" (Subscript[I, 2]) + (SuperMinus[I]) ⇌ (Subscript[I, 3]^(-)) ", K6, Kneg6}, {" HOI ⇌ (SuperMinus[IO]) + (SuperPlus[H]) ", K7, Kneg7}, {" (Subscript[I, 2]) + (Subscript[•O, 2]^(-)) -> (Subscript[•I, 2]^(-)) + (Subscript[O, 2]) ", 1.6×(10^12)(E^(-(1662./TAB))), 0}, {" HOI + (Subscript[•O, 2]^(-)) -> (SuperMinus[OH]) + •I + (Subscript[O, 2]) ", 1.0×(10^6)EA, 0}, {" (Subscript[•I, 2]^(-)) + (Subscript[•O, 2]^(-)) -> (Subscript[O, 2]) + (SuperMinus[I]) + (SuperMinus[I]) ", 3.×(10^9)EA, 0}, {" HOI + •OH -> •IO + (Subscript[H, 2])O ", 7.×(10^9)EA, 0}, {" •IO + •OH -> (Subscript[HIO, 2]) ", KTEST EA, 0}, {" (SuperMinus[I]) + •OH ⇌ (SuperMinus[•HOI]) ", 1.1×(10^10)EA, 1.13×(10^6)EA}, {" (Subscript[I, 2]) + •OH -> HOI + •I ", 1.1×(10^10)EA, 0}, {" (Subscript[I, 2])(SuperMinus[OH]) + (Subscript[H, 2])(Subscript[O, 2]) ⇌ (Subscript[HIO, 2]) + (SuperMinus[I]) + (Subscript[H, 2])O ", 2.3×(10^6)(E^(31000`2tf)), 1.×(10^7)(E^(92000`2tf))}, {" (Subscript[HIO, 2]) + (SuperMinus[OH]) -> (SuperMinus[I]) + (Subscript[O, 2]) + (Subscript[H, 2])O ", (1/(1.×(10^(-10))+5.×(10^(-10)) (E^(-54000`2 tf)))), 0}, {" •I + (SuperMinus[I]) ⇌ \!(Subscript[•I, 2]^(-)) ", 1.2×(10^10)EA, 5.42×(10^8)(E^(-(2562./TAB)))}, {" (SuperMinus[•HOI]) -> •I + (SuperMinus[OH]) ", 1.2×(10^8)EA, 0}, {" (Subscript[•I, 2]^(-)) + (Subscript[•I, 2]^(-)) -> (Subscript[I, 3]^(-)) + (SuperMinus[I]) ", 4.5×(10^9)EA, 0}, {" (Subscript[•I, 2]^(-)) + •I -> (Subscript[I, 3]^(-)) ", 4.5×(10^9)EA, 0}, {" •I + •I -> (Subscript[I, 2]) ", 1.×(10^10)EA, 0}})

This table forms the sole input from which ODEs will be derived.

Extract chemical species formulae

Ω is a list of strings containing 2-D chemical formulae.

Ω=Union[Flatten[hoReactionFormulae]]

Parse chemical reaction equations

Style[Column[{"(1) Determine arrow position ordinals in chemical reaction equations","(2) Determine the set of chemical reaction equations that contain the given species Ω","(3) Determine the total number of chemical reaction equations that contain Ω","(4) Separate chemical reaction equations into reactants and products","(5) Determine formula occurence frequencies as reactants and products, (Subscript[α, j]) and (Subscript[β, j])"}],FontFamily->"Arial",Larger]

Function to transform 2-d formula strings "Ω" into 1-d concentration symbols ω and expressions ω[t]

makeConcentrationExpr[Ωs_String]:=
Module[{i=Ωs},i// (* Function module to make concentration symbols and expressions. *)
    StringTrim[#]& //(* Remove leading and trailing spaces. *) 
    StringCases[#,WordCharacter..|"•"|"-"|"+"]&  //(* Get sequences of characters. *)
Module[{j=#},j// (* This submodule gets rid of unwanted box-expr heads. *)
                    Not/@StringMatchQ[#,LetterCharacter..~~"Box"]& //(* Make True-or-False filter from box-expr heads. *)
                    Pick[j,#]&  (* Apply filter to get only sequences that are required. *)
]&//
     ReplaceAll[#,{"-"->"minus","+"->"plus"}]& //(* Convert + and - characters to make modified sequences. *)
    Apply[StringJoin,#]&  //(* Make string names of concentration symbols. *)
    StringReplace[#,s:RegularExpression["^[CDEIKNO]$"]:>s<>"1"]& //(* Append 1 to special single-character formulae. *)
    ToExpression //(* Make concentration symbols. *) 
    Apply[#,_[t]]&  (* Make concentration symbols variable functions of time t. *) 
]

(Programming style is Kristen Carlson's Postfix / Pure Function or "Functional-Procedural Fusion".)

Results of applying the function makeConcentrationExpr

Grid[Table[{Ω[[i]],ω[Ω[[i]]]=makeConcentrationExpr[Ω[[i]]]},{i,1,Length[Ω]}],Alignment->Left]

Differential equation generation

Here is a function to generate the differential equation for the concentration of the chemical species having chemical formula f.

genODE[f_]:=(∂_t) ω[f]==-(∑_(j=1)^(jL[f]))((α[f,j]-β[f,j])(Kf[f,j](∏_(i=1)^(iL[f,j]))((ω[f,j,i])^(α[f,j,i]))-Kb[f,j](∏_(i=1)^(iL[f,j]))((ω[f,j,i])^(β[f,j,i]))))

Here is a test of the function applied to a single 2-d chemical formula.

genODE[" (Subscript[I, 2])(SuperMinus[OH]) "]

Showing the chemistry

Using Mathematica, the equation

(Subscript[((∂_t) ω), ho])=-(∑_j^( ))((Subscript[((Subscript[α, j])-β), j])((Subscript[K, fj])(∏_i^( ))(Subscript[ω, ji]^(Subscript[α, ji]))-(Subscript[K, bj])(∏_i^( ))(Subscript[ω, ji]^(Subscript[β, ji]))))

applied to the reduced iodine set of 18 homogeneous chemical reactions, is shown to yield the following ODEs.

diffEqn=genODE[#]&/@Ω;
Grid[Transpose[{Ω,diffEqn}],Frame->All,Alignment->Left]

Conclusion

An arbitrary set of homogeneous chemical reaction equations have been automatically transformed into ODEs, by means of functions defined using Mathematica.

Once the "chemistry has been shown" for the remaining 4 equations in the phenomenologically based theory, it will be a simple matter to apply the defined functions to the complete set of >450 chemical reaction equations represented in LIRIC and generate the entire corresponding system of ODEs.  With the addition of initial conditions and embedded correlations, a Mathematica replacement for LIRIC is within easy reach.

Epilogue

Is a Mathematica radiation-chemistry program in the making?  Yes!

Benefits:

capability for direct use of 2-D chemical reaction equations without the need for manual translation, resulting in fewer errors

availability as a validation / testing tool for the existing FACSIMILE version of LIRIC

simpler configuration management

capability to provide a superior user interface