<content> dsol </content>

06fc258d-d2ad-4c5f-85dd-0ff4b7c5738e 92 Reg Curry 4 5 false false false 0 3 decimal

&[DATE] &[TIME] - &[FILENAME]

NUCLEAR FISSION YIELD CALCULATIONS

<content> dsol </content> <content> Utils </content>

M Cols C 0 : k 1 M cols range : 12 mat C 1 k el M k col : 21 line :

a b n pR a b a - n / 0 n range * + eval 11 line :

A D C B TriDiagSys n B length : v n 1 matrix : y v : 13 mat w A 1 el 1 -^: y 1 el B 1 el w * : 12 mat k 2 n range v k 1 - el C k 1 - el w * eval : w A k el D k el v k 1 - el * - (1 -^eval : y k el B k el D k el y k 1 - el * - (w * eval : 31 line for k n 1 - (1 range y k el y k el v k el y k 1 + el * - eval : for y 51 line :

x y pStem XY 03 matrix : S 02 matrix : 12 mat k 1 x length range XY XY x k el y k el o 13 mat stack : S S x k el 0 x k el y k el x k el 032 mat stack : 21 line for XY S 21 sys eval 31 line :

e# a# b# Equrep equrep(e# num2str, a# num2str, b# num2str) concat str2num 11 line :

x# a# Find n# 1 x# length range : f# x# n# el n# augment eval a# 1 findrows : (0 ≡ 0 f# 2 col eval if 21 line :

fs# x# feval fs# num2str (x# num2str) concat str2num 11 line :

fs# x# y# feval fs# num2str (x# num2str, y# num2str) concat str2num 11 line :

fs# x# y# z# feval fs# num2str (x# num2str, y# num2str, z# num2str) concat str2num 11 line :

<content> Options </content>

O# o# Options Ω# 01 matrix eval : n# 0 : 12 mat ω# o# sys2mat.1 num2str vectorize eval ω# ω# " strrep : n# n# 1 + : 12 mat g# ω# ≡ strsplit : length 1 > f# 1 O# rows range O# augment eval g# 1 el 2 findrows : (0 ≠ (O# f# 1 el 2 el g# 2 el : O# O# g# 1 el g# 2 el augment 11 O# cols range el stack eval : if Ω# Ω# g# 2 el stack eval : 21 line f# O# ω# 1 findrows : (0 ≠ (Ω# Ω# f# 2 el stack eval : Ω# Ω#'NaN stack eval : if if 21 line for Ω# rows 2 > o# cols 1 ≡ Ω# Ω# transpose if Ω# 1 el if 31 line :

<content> Example Options </content>

dsol solver Rkadapt steps 100 stepsmax 0 eps 10^{-7} 42 mat :

dsol dummy Options'NaN

dsol steps dummy eps 31 mat Options 100'NaN 10^{-7} 31 mat

dsol steps dummy eps 13 mat Options 100'NaN 10^{-7} 13 mat

dsol steps solver rkfixed ≡ eps 31 mat Options 100 rkfixed 10^{-7} 31 mat

dsol solver rkfixed steps 100 stepsmax 0 eps 10^{-7} 42 mat

dsol Jacobian J# ≡ Options J#

dsol solver rkfixed steps 100 stepsmax 0 eps 10^{-7} Jacobian J# 52 mat

Check inside a procedure

x f dsol dummy x ≡ Options dsol 21 line :

1 f solver rkfixed steps 100 stepsmax 0 eps 10^{-7} Jacobian J# dummy 162 mat

dsol solver rkfixed steps 100 stepsmax 0 eps 10^{-7} Jacobian J# 52 mat

lbvp.2(f(x),[a,b],M,n) solves the linear BVP y'' + p(x)y' + q(x)y = φ(x) in a < x < b with    f(x) = [p(x), q(x), φ(x)],  M = [α1 β1 c1; α2 β2 c2]  and     α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2using central differences. 1 f x M n lbvp.2 h x 2 el x 1 el - n / : X x 1 el 0 n range h * + : ν n 1 + : P 0 : Q 0 : Φ 0 : α1 β1 c1 α2 β2 c2 23 mat M : 17 mat r 1 ν range P r el Q r el Φ r el 13 mat X r el f eval : for A D C B 14 mat 2 - h 2^Q * + (2 h P * - 2 / (2 h P * + 2 / (h 2^Φ * (14 mat eval : A 1 el D 1 el C 1 el B 1 el 14 mat β1 0 ≡ 100 c1 α1 / 14 mat A 1 el 2 h * α1 β1 / * D 1 el * + (02 B 1 el 2 h * c1 β1 / * D 1 el * + (14 mat if eval : A ν el D ν el C ν el B ν el 14 mat β2 0 ≡ 100 c2 α2 / 14 mat A ν el 2 h * α2 β2 / * C ν el * - (20 B ν el 2 h * c2 β2 / * C ν el * - (14 mat if eval : Y A D C B TriDiagSys : Y' Y : j 2 n range : 13 mat Y' j el Y j 1 + el Y j 1 - el - 2 h * / : Y' 1 el β1 0 ≡ Y 2 el Y 1 el - h / c1 β1 / α1 β1 / Y 1 el * - if : Y' ν el β2 0 ≡ Y ν el Y n el - h / c2 β2 / α2 β2 / Y ν el * - if : 13 mat X Y Y' augment eval 81 line :

bvp.2(f(x),[a,b],M,n,ε) solves the non-linear BVP y'' = f(x,y,y') in a < x < b with M = [α1 β1 c1; α2 β2 c2] where     α1*y(a) + β1*y'(b) = c1   and  α2*y(a) + β2*y'(b) = c2using central differences. fs# x Y M n ε bvp.2 x# y# y'# f fs# x# y# y'# feval : h x 2 el x 1 el - n / : X x 1 el 0 n range h * + : ν n 1 + : P 0 : Q 0 : Φ 0 : α1 β1 c1 α2 β2 c2 23 mat M : 17 mat Y 0 ≡ ψ x 1 el α1 * β1 + α1 x 2 el α2 * β2 + α2 22 mat (1 -^c1 c1 21 mat * : Y ψ 2 el ψ 1 el X * + : 12 mat 0 if Δy Y max Y min - ν / (2^h 2^+ sqrt eval : Δy' Δy h 1 -^* : 12 mat x y y' fy x y Δy + y' f x y Δy - y' f - 2 Δy * / : x y y' fy' x y y' Δy' + f x y y' Δy' - f - 2 Δy' * / : A ν 1 matrix eval : D A : C A : B A : Y' A : j 2 n range : 16 mat Iter 1 dsol MaxIters Options str2num range Y' j el Y j 1 + el Y j 1 - el - 2 h * / : Y' 1 el β1 0 ≡ Y 2 el Y 1 el - h / c1 β1 / α1 β1 / Y 1 el * - if : Y' ν el β2 0 ≡ Y ν el Y n el - h / c2 β2 / α2 β2 / Y ν el * - if : 13 mat r 1 ν range P r el Q r el Φ r el 13 mat X r el Y r el Y' r el fy X r el Y r el Y' r el fy' X r el Y r el Y' r el f 13 mat eval : for A D C 13 mat 2 - h 2^P * - (2 h Q * - 2 / (2 h Q * + 2 / (13 mat eval : B j el 2 Y j el * - Y j 1 + el + Y j 1 - el + h 2^Φ j el * - eval : A 1 el D 1 el C 1 el B 1 el 14 mat β1 0 ≡ 100 Y 1 el c1 α1 / - 14 mat A 1 el 2 h * α1 β1 / * B 1 el * + (022 Y 2 el * 21 - h α1 β1 / * + (* Y 1 el * + 2 h * c1 β1 / * - h 2^Φ 1 el * - (14 mat if eval : A ν el D ν el C ν el B ν el 14 mat β2 0 ≡ 100 Y ν el c2 α2 / - 14 mat A ν el 2 h * α2 β2 / * C ν el * - (202 Y n el * 21 h α2 β2 / * + (* Y ν el * - 2 h * c2 β2 / * + h 2^Φ ν el * - (14 mat if eval : Yo A D C B - TriDiagSys : Y Y Yo + eval : 12 mat Yo norme ε Y norme * < break continue if 81 line for Iter dsol MaxIters Options str2num < X Y Y' augment eval MaxIters  rached. error if 91 line :

dsol(DE,X,tf,N) try to solve the Differential Equation system DE for the functions in X. Ev# Yv# TF# N# dsol m kg s K A cd 16 mat 11111116 mat : u# v# t# EQ_F# IG_F# EQ_F# RK_D# Clear ky# 0 : kn# 0 : kv# 0 : 14 mat ys# 0 : U# 0 : V# 0 : K# 0 : NK# 0 : IC# 03 matrix eval : 16 mat M# 0 : ES# Ev# num2str : ES'# ES#' strsplit rows : 13 mat BIC# 0 : bic# 0 : too# : 13 mat s# Yv# sys2mat.1 num2str vectorize eval a# b# 12 mat s# (findstr 1 el s#) findstr 1 el 12 mat : ts# s# a# 1 + b# a# - 1 - substr : Assuming all ≡ 12 mat ys# ky# ky# 1 + : el s# 1 a# 1 - substr : k# 1 ES'# range ES# s# s# ('(strrep : findstr 1 - ≠ M# ky# el k# : if for 41 line for s# Ev# sys2mat.1 num2str vectorize eval g# s# ≡ strsplit : g'# s# ≈ strsplit : 12 mat dsol g# 1 el 1 findrows 0 ≡ NK# NK# 1 + : s# (ts#) concat findstr 1 - ≡ a# b# 12 mat g# 1 el (findstr 1 el g# 1 el) findstr 1 el 12 mat : to# s# a# 1 + b# a# - 1 - substr : too# ≡ (too# to# ≠ (| BIC# bic# bic# 1 + : el to# : too# to# : 12 mat if g'# length 1 ≡ ico# s# 1 a# 1 - substr : S# ico#' strrep (ts#) concat strrep : f# ys# S# Find : (0 ≠ IC# IC# f# 1 el ico#' strsplit rows g# 2 el 13 mat stack : if K# NK# 1 el IC : K# NK# 2 el g# 1 el str2num g# 2 el str2num - : bvp 13 mat 41 line K# NK# 1 el IG : K# NK# 2 el g'# 1 el str2num g'# 2 el str2num - : 12 mat if 41 line g# length 1 ≡ K# NK# 1 el EQ : K# NK# 2 el s# str2num : 12 mat K# NK# 1 el EQ : K# NK# 2 el g# 1 el str2num g# 2 el str2num - : 12 mat if if 21 line dsol s# Options if 21 line for ic# 1000 IC# 1 col * IC# 2 col + IC# 3 col augment : ic# ic# 1 csort 2 col eval : ic# ic# str2num vectorize : 12 mat n# 1 ky# range kn# kn# 1 + : a# ys# n# el (ts#) concat : 12 mat K# K# a# u# kn# el Equrep : m# 1 : m# M# n# el ≤ m# m# 1 + : a# a# ('(strrep : m# M# n# el 1 - ≤ kn# kn# 1 + : kv# kv# 1 + : V# kv# el u# kn# el : 13 mat K# K# a# u# kn# el Equrep : 21 line if 21 line for kv# kv# 1 + : U# n# el v# n# el : V# kv# el v# n# el : 13 mat K# K# a# v# n# el Equrep : 51 line for EQ# K# EQ 1 findrows 2 col : v# EQ_F# EQ# eval : 12 mat bic# 1 > IC# K# IC 1 findrows 2 col : u# FI# IC# : 12 mat M# u1# u2# 21 mat FI# u1# diff vectorize u1# u2# 21 mat FI# u2# diff vectorize 0021 mat FI# - augment eval : dsol dsolver Options lbvp.2 ≡ u# v# FE# EQ# 1 el : dv# u1# u2# 21 mat v1# 11 mat FE# v1# diff : 12 mat L# u1# u2# 21 mat v1# 11 mat FE# u2# diff dv# / u1# u2# 21 mat v1# 11 mat FE# u1# diff dv# / 0021 mat 011 mat FE# dv# / - 13 mat : lf#(ts#): L# num2str concat str2num v1# lf# BIC# str2num vectorize M# N# lbvp.2 41 line dsol dsolver Options bvp.2 ≡ vg# dsol DefaultGuess Options str2num : vg# vg# vg# stack : f# K# IG 1 findrows : num2str 0 ≠ IG# f# 2 col : u# IG_F# IG# eval : 12 mat vg# vg# 01020 -^IG_F# al_nleqsolve 2 el N# 1 + 1 matrix + eval : 21 line vg# 0 : if v# BF# EQ# : ε# dsol eps Options str2num : go# dsol DefaultGuess Options str2num 11 mat : t# y# y'# bf# ts# : t# num2str concat str2num u# y# y'# 12 mat : go# BF# al_nleqsolve 1 el 31 line : bf# BIC# str2num vectorize vg# M# N# ε# bvp.2 81 line Not Implemented. error if 31 line vg# dsol DefaultGuess Options str2num ky# 1 matrix + eval : f# K# IG 1 findrows : num2str 0 ≠ IG# f# 2 col : v# IG_F# IG# eval : vg# vg# 01020 -^IG_F# al_nleqsolve : 13 mat if v# vg# : dsol UpdateGuess Options 1 ≡ t# u# v# RK_D# ts# : t# num2str concat str2num v# v# EQ_F# al_nleqsolve : V# eval 31 line : t# u# RK_D# ts# : t# num2str concat str2num v# vg# EQ_F# al_nleqsolve : V# eval 31 line : if dsol dsolver Options (ic#,str2num(to#),TF#,N#,RK_D#) concat str2num 51 line if 121 line :

DE# Y# N# dsol DE# Y# 0 N# dsol 11 line :

dsol dsolver rkfixed steps 100 nsolver al_nleqsolve stepsmax 0 eps 10^{-7} DefaultGuess 1 UpdateGuess 1 MaxIters 5082 mat :

nsol nsolver al_nleqsolve stepsmax 0 eps 10^{-7} 32 mat :

<content> NEUTRON TRANSPORT INPUT TO YIELD CALCULATIONS </content>

ALPHA AND LAMDA MODEL

DEFINE SOME UNITS

CHOOSE TAMPER MATERIAL; "DU" OR "WC"

Tamper DU :

M.Tamper 231.0 kg * :

CHOOSE CORE MATERIAL: "U" OR "Pu"

CORE Pu :

M.Core 6.3 kg * :

COMPRESSION RATIO: \[REGION[9aa76a8d31f14602b1b72effb33e911a]]\

CR 2.5 :

AMU 1.66 1027 -^* kg * :

1 eV * 1.602 1019 -^* J *

1 MeV * 1.602 1013 -^* J *

Avogadro's Number

N.A 6.022 1023^* 1 mol / *

Ton.TNT 4.2 109^* J * :

KT 1000 Ton.TNT * :

CORE DEFINITION

HEU 235U PERCENT: \[REGION[ee5d80608d224c02b46b1fbc905fb865]]\

HEU% 100 :

HEUMASS M.Core :

M.core HEUMASS HEU% 100 / * : 6.300 kg *

GRAM MOLECULAR WEIGHT OF 235U: \[REGION[030ac39937c445e5971bc1bfa63e3a99]]\

U235.gmw 235 g * :

MolesU235.kg kg U235.gmw / : mol 4.255 1 mol / *

NumberAtomsU235.kg MolesU235.kg N.A * : 2.563 1024^* 1 mol / *

DATA ELEMENT r (gm/cm3) lf (cm) lt (cm) n U235 1.871 101^* 1.689 101^* 3.596 2.637 Pu239 1.560 101^* 1.414 101^* 4.108 3.172 DU 1.895 101^* 0.000 4.342 0.000 WC 1.480 101^* 0.000 3.159 0.000 55 mat :

ρ DATA 22 el g cm 3^/ * : g cm 3^/ 18.710

λ.f DATA 23 el g cm 3^/ * : g cm 3^/ 16.890

λ.t DATA 24 el g cm 3^/ * : g cm 3^/ 3.596

ρ.DU DATA 42 el g cm 3^/ * : g cm 3^/ 18.950

λ.TampDU DATA 44 el g cm 3^/ * : g cm 3^/ 4.342

ρ.WC DATA 52 el g cm 3^/ * : g cm 3^/ 14.800

λ.TampWC DATA 54 el g cm 3^/ * : g cm 3^/ 3.159

ρ.pu DATA 32 el g cm 3^/ * : g cm 3^/ 15.600

λ.puf DATA 33 el g cm 3^/ * : g cm 3^/ 14.140

λ.put DATA 34 el g cm 3^/ * : g cm 3^/ 4.108

ν DATA 25 el : 2.637

ν.pu DATA 35 el : 3.172

NOW DEFINE λs IN TERMS OF NUCLEAR CROSS SECTION AND NUCLEI NUMBER DENSITY

A 235.04 g mol / * :

n ρ N.A * A / : cm 3 -^4.794 1022^*

σ.f 1.235 barn * :

λ.f 1 n σ.f * / : 0.169 m *

σ.el 4.566 barn * :

σ.t σ.f σ.el + :

σ.t barn 5.801

λ.t 1 n σ.t * / : cm 3.596

A.WC 195.85 g mol / * :

n.WC ρ.WC N.A * A.WC / : cm 3 -^4.551 1022^*

σ.elWC 6.857 barn * :

λ.TampWC 1 n.WC σ.elWC * / : cm 3.205

A.DU 238.05 g mol / * :

n.DU ρ.DU N.A * A.DU / : cm 3 -^4.79 1022^*

σ.elDU 4.804 barn * :

λ.TampDU 1 n.DU σ.elDU * / : cm 4.34

σ.put σ.puel σ.puf + :

A.pu 239.05 g mol / * :

σ.puel 4.394 barn * :

σ.puf 1.80 barn * :

ν CORE U ≡ ν ν.pu if :

ρ.core CORE U ≡ ρ ρ.pu if :

ρ.tamper Tamper WC ≡ ρ.WC ρ.DU if :

R.c 3 M.core * 4 π * ρ.core * / 3 nthroot : cm 4.585

R.T 3 M.Tamper * 4 π * ρ.tamper * / R.c 3^+ 3 nthroot : cm 14.433

ρ.core g cm 3^/ 15.600

ρ.tamper g cm 3^/ 18.950

M.core 6.300 kg *

M.Tamper 231.000 kg *

NOW LET THE RADII OF THE CORE AND TAMPER EXPAND BY Δr \[REGION[40f85b2fef844511aa77e4d06a791dd3]]\

maxk 50 :

k 1 maxk range :

Δr 0.01 cm * :

Rexp R.c k 1 - (Δr * + :

Rexp max cm 5.08

RTexp R.T k Δr * + :

RTexp max cm 14.93

ρexp 3 M.core * 4 π * / Rexp (3 -^* vectorize :

nexp CORE U ≡ ρexp N.A * A / ρexp N.A * A.pu / if :

σt CORE U ≡ σ.t σ.put if :

σf CORE U ≡ σ.f σ.puf if :

ρexp 1 el g cm 3^/ 15.600

ρTexp 3 M.Tamper * 4 π * RTexp 3^Rexp 3^- (* / vectorize :

nTexp Tamper WC ≡ ρTexp N.A * A.WC / ρTexp N.A * A.DU / if :

σ.tamper Tamper WC ≡ σ.elWC σ.elDU if :

σ.tamper barn 4.804

ρTexp 1 el g cm 3^/ 18.909

λf 1 σf nexp * / vectorize :

λt 1 σt nexp * / vectorize :

λT 1 σ.tamper nTexp * / vectorize :

Λ λT λt / vectorize :

λT max λT min / 1.097

Rexp Rexp CR 3 nthroot / :

RTexp RTexp CR 3 nthroot / :

λT λT CR / :

λf λf CR / :

λt λt CR / :

α f1 e 23 α * λT λf * / sqrt * Rexp RTexp - (*^vectorize :

α f2 3 ν α - 1 - (* λf λt * / sqrt Rexp * 3 ν α - 1 - (* λf λt * / sqrt Rexp * cot * 1 - Λ 3 α * λT λf * / sqrt Rexp * 1 - (* - RTexp 23 / λT * 3 α * λT λf * / sqrt RTexp * 1 - (* + / vectorize :

α f3 3 ν α - 1 - (* λf λt * / sqrt Rexp * 3 ν α - 1 - (* λf λt * / sqrt Rexp * cot * 1 - Λ 3 α * λT λf * / sqrt Rexp * 1 + (* + RTexp 23 / λT * 3 α * λT λf * / sqrt RTexp * 1 + (* - / vectorize :

α F α f1 α f2 * α f3 - vectorize :

Ans α F α solve 2 el vectorize :

α Ans :

r Rexp :

slope α maxk el α 1 el - r maxk el cm / r 1 el cm / - / : 0.719 -

αfit slope r r 1 el - cm / (* α 1 el + :

alpha α 1 el slope 21 mat :

alpha 0.797 0.719 - 21 mat

x r cm / :

λslope λf maxk el cm / λf 1 el cm / - r maxk el cm / r 1 el cm / - / :

λfit λslope r r 1 el - cm / (* λf 1 el cm / + :

α 0.797 0.791 0.785 0.779 0.773 0.768 0.762 0.756 0.751 91 mat

x 3.379 3.386 3.393 3.401 3.408 3.415 3.423 3.430 3.437 91 mat

αfit 0.797 0.791 0.786 0.781 0.775 0.770 0.765 0.759 0.754 91 mat

λf cm 5.65 5.69 5.73 5.77 5.80 5.84 5.88 5.92 5.96 91 mat

lambda λf 1 el m / λslope 21 mat :

lambda 0.057 5.577 21 mat

plot x α augment :

Fit x αfit augment :

Lambda x λf cm / augment :

λFit x λfit augment :

λf cm / max 7.668

α max 0.797

x max 3.74

λf cm / min 5.655

α min 0.537

x min 3.38

plot Fit 2 1 sys Lambda λFit 2 1 sys

Params

α alpha : 0.797 0.719 - 21 mat

λ lambda : 0.057 5.577 21 mat

M.core 6.30 kg * :

M.Tamper 231.00 kg * :

Δr 1.00 104 -^* m * :

R.c 4.59 cm * :

CR 2.50 :

E.f 180 MeV * :

M.total M.core M.Tamper + : 237.30 kg *

N.o 1 cm 3 -^* :

γ 13 / :

KT 1000 4.2 109^* J * * :

NeutronFissionEnergy 2 MeV * :

m.n 1.67493 1027 -^* kg *

Nvel 2 NeutronFissionEnergy * m.n / sqrt : 1.96 107^* m s / *

Initials

Eo 1010 -^KT * :

E Eo 11 mat :

τ λ 1 el m * Nvel / 11 mat : 2.89 109 -^* s * 11 mat

αf 1 el α 1 el :

ro R.c CR 3 nthroot / : cm 3.38

r ro 11 mat :

V.core 43 / π * r 1 el (3^* 11 mat : cm 3^162.03 11 mat

v 0 m s / * 11 mat :

Ode sol

t V.core 43 / π * t r 3^* :

t τ λ 2 el t r ro - (* λ 1 el m * + Nvel / :

t.end 0.25 μs * :

RK t r t E t α 3 1 sys t.end kg 1 : m 1 : s 1 : 1 3 mat 500 2 1 line dsol :

t r'' 4 π * t r (2^* γ * t E * t V.core M.total * / ≡

0 r ro ≡

r.RK RK 2 col :

T.RK RK 1 col :

E.RK RK 4 col :

0 r' 0 ≡

α.RK RK 5 col :

t E' N.o t V.core * E.f * t τ / e t α t τ / t *^* ≡

0 E Eo ≡

Plot T.RK μs / E.RK KT / augment T.RK μs / α.RK augment 21 sys :

dsolver Adams ≡

t α' α 2 el t r' cm / * ≡

0 α α 1 el ≡

Plot

Euler

npoints 4000 :

j 1 npoints range :

Δt t.end npoints 1 + / : μs 6.25 105 -^*

T j el 0 j Δt * + :

j 2 npoints range r j el r j 1 - el v j 1 - el Δt * + eval : αf j el α 2 el r j 1 - el r 1 el - cm / (* α 1 el + eval : τ j el λ 2 el r j 1 - el r 1 el - (* λ 1 el m * + Nvel / eval : v j el v j 1 - el 4 π * r j 1 - el (2^* γ * E j 1 - el * V.core j 1 - el M.total * / Δt * + eval : V.core j el 43 / π * r j 1 - el (3^* eval : E j el E j 1 - el N.o V.core j 1 - el * E.f * τ j 1 - el / (e αf j 1 - el τ j 1 - el / (T j 1 - el *^* Δt * + eval : 61 line for

Yield T μs / E KT / augment :

Ngrowth T μs / αf augment :

Plot Yield Ngrowth 21 sys :

Plot