Hi all, attempting to make a Method of Characteristics solver in Matlab. I'm particularly hoping that there are some computational nuclear engineering guys about who might have a bit of experience programming a simple version of one of these.
I'm trying to create the solver to find the scalar neutron flux and criticality in a one-group, homogeneous, multiplying with slab subject to vacuum boundary conditions but, infuriatingly, while I can attain a neutron flux profile which is correct (verifying with other codes) my criticality is off by a factor of 10! I was hoping that someone may be able to provide a bit of insight into what I'm doing incorrectly - it may be a long shot but it's a problem I've been suffering with for far too long. If anyone would care to have a look at the code as it stands then do let me know!
For the method of characteristics, the fundamental equations that I'm iterating around in my solver are quite standard. I won't include them here as I would rather just share my code if you are familiar with the general problem.
The Attempt at a Solution
For my problem, I'm examining a geometry of length 100cm with a Σtr of 0.2758/cm, νΣf of 0.0366/cm and Σa of 0.0253/cm. As I say, I can get the correct normalised neutron flux profile across the slab but sadly k is approximately a factor of 10 smaller than the keff value of about 1.3876 that I am expecting. It is interesting to note that for perfectly reflective boundary conditions I attain a k∞ which is equal to νΣf/Σtr. As I say, feel free to ask for the code if you think you might be of any help - if you can identify the problem I'll owe you my dissertation!