Sample Problems > Usage > Optimization > criticality

{ CRITICALITY_SIZE.PDE

This problem demonstrates the use of FlexPDE in the solution of optimization problems.

FlexPDE implements the Nelder-Mead "amoeba" algorithm to minimize or maximize an objective function.

This is not the method of greatest speed, but it is very flexible, and allows FlexPDE to perform

optimization searches in a wide range of problem environments.

A simple model of nuclear criticality can be made using a Fick's-Law diffusion equation for single-velocity neutrons.

In this model, the criticality relation can be stated as

Div(D*grad(N)) -sigmar*N + (1/k)*nper*sigmaf*N = 0,

where

N = neutron density

beta = Fick's Law proportionality factor

D = beta/sigmat

k = 1/lambda = "criticality eigenvalue"

mix = fractional content of fissionable material in the source region

nper = number of neutrons produced at each absorption

sigmat = transport cross-section = 1/(transport mean free path)

sigmar = removal cross-section

sigmaf = fission cross section

In this probem, we use the "mix" parameter discovered in CRITICALITY.PDE and solve for the critical size of the assembly.

This should return us 5.0, the size used in CRITICALITY.PDE.

}

title 'Nuclear Criticality'

variables

N { neutron Density }

select

modes=1 { calculate only the smallest eigenvalue }

cell_limit=2000

definitions

source { name the material parameters, values will be declared by region }

sigmat { transport cross-section }

sigmar { removal cross-section }

sigmaf { fission cross-section }

beta = 1/3 { Fick's Law proportionality factor }

nper = 2 { number of neutrons produced per fission }

convergence = (lambda-1)^2 { optimization parameter: a function with a smooth minimum at lambda=1 }

Hw { the outer halfwidth of the structure }

box = 3*Hw/5 { the halfwith of the reactive inclusion }

mix = 0.086052 { as reported by CRITICALITY.PDE }

! Here is the optimization request:

! Modify the halfwidth "Hw" until the eigenvalue "lambda" is 1.0

! "Hw" starts at 10, and the initial range of "Hw" samples is 1.0

! Iteration continues until "convergence" is less than OPTERRLIM (default 2e-6) or "Hw" is bracketed to less than OPTERRLIM*Hw.

minimize(convergence) vs Hw (10,1)

equations { The neutron density equation: }

N : div(beta/sigmat*grad(N)) + source - sigmar*N + lambda*mix*nper*sigmaf*N = 0

boundaries

region 1 { the bounding region is absorbing }

source=0 sigmar=1 sigmat=2 sigmaf=0

start(-Hw,-Hw)

natural(N)= 0

line to (Hw,-Hw)

line to (Hw,Hw) to (-Hw,Hw)

line to close

region 2 { this region has fission }

source=0 sigmar=0.1 sigmat=2 sigmaf=1

start(-box, -box)

line to (box, -box) to (box, box) to (-box, box) to close

monitors

contour(N) as "Neutron Density"

history(lambda,Hw,convergence)

plots

grid(x,y)

contour(N) as 'Neutron Density' report(Hw) report(lambda)

surface(N) as 'Neutron Density'

vector(-beta/sigmat*grad(N)) as 'Neutron Flux'

contour(sigmaf*nper*N) as "Fission Source"

contour(sigmar*N) as "Neutron Absorption"

history (lambda)

history(Hw)

history(convergence)

history (lambda) vs Hw

end

criticality_size

criticality_size