Difference between revisions of "Unstable 3D pin-cell burnup problem"

From Serpent Wiki
Jump to: navigation, search
(Created page with "In certain burnup problems such as long, axially symmetric 3D assemblies or fuel rods or large symmetric core geometries traditional Monte Carlo burnup schemes may run into in...")
 
m
 
(12 intermediate revisions by one other user not shown)
Line 1: Line 1:
−
In certain burnup problems such as long, axially symmetric 3D assemblies or fuel rods or large symmetric core geometries traditional Monte Carlo burnup schemes may run into instabilities<ref name="DufekNSE09"/><ref name="DufekANE13a"/>.
+
In certain burnup problems such as long, axially symmetric 3D assemblies or fuel rods or large symmetric core geometries traditional Monte Carlo burnup schemes may run into instabilities<ref name="DufekNSE09"/><ref name="DufekANE13a"/>. This page describes the simulation of one such case first using the traditional burnup schemes (to showcase the instability) and then using a stable burnup scheme.
  
 +
== Problem description ==
 +
 +
<div class="toccolours mw-collapsible mw-collapsed" style="width:40em;">
 +
'''Base input for unstable 3D pin-cell problem'''
 +
<div class="mw-collapsible-content">
 +
<nowiki>/*****************
 +
** Run options **
 +
*****************/
 +
 +
% --- Neutron population to be used
 +
 +
set pop 10000 500 200
 +
 +
% --- 200 W/cm linear power
 +
 +
set power 60000
 +
 +
/*************************
 +
** Geometry definition **
 +
*************************/
 +
 +
% --- Fuel Pin definitions:
 +
 +
pin p1
 +
fuel  0.47
 +
void  0.48
 +
clad  0.54
 +
cool
 +
 +
% --- Lattice
 +
 +
lat l1  1  0.0 0.0 1 1 1.5
 +
p1
 +
 +
% --- Surrounding surfaces:
 +
 +
% Boundary of geometry:
 +
 +
surf 3 cuboid -0.75 0.75 -0.75 0.75 -160 160
 +
 +
% Lower boundary of fuel
 +
 +
surf 4 pz -150
 +
 +
% Upper boundary of fuel
 +
 +
surf 5 pz  150
 +
 +
% --- Cell definitions:
 +
 +
% Active fuel pin
 +
 +
cell  3  0  fill l1  -3  4 -5
 +
 +
% Coolant below active fuel (bottom reflector)
 +
 +
cell  4  0  cool    -3 -4
 +
 +
% Coolant above active fuel (top reflector)
 +
 +
cell  5  0  cool    -3  5
 +
 +
% outside world
 +
 +
cell 99  0  outside  3    % Outside world
 +
 +
% --- Reflective boundary conditions in XY, black in Z:
 +
 +
set bc 3 3 1
 +
 +
/**************************
 +
** Material definitions **
 +
**************************/
 +
 +
% --- Fuel material (4.85 % enrichment):
 +
 +
mat fuel    -10.283 vol 208.19 rgb 200 200 125
 +
92235.09c      0.016166667
 +
92238.09c      0.317166667
 +
8016.09c      0.666666667
 +
 +
% --- Cladding (Zr-4)
 +
 +
mat clad -6.56000E+00 rgb 180 180 180
 +
8016.06c    -1.19276E-03
 +
8017.06c    -4.82878E-07
 +
24050.06c    -4.16117E-05
 +
24052.06c    -8.34483E-04
 +
24053.06c    -9.64457E-05
 +
24054.06c    -2.44600E-05
 +
26054.06c    -1.12572E-04
 +
26056.06c    -1.83252E-03
 +
26057.06c    -4.30778E-05
 +
26058.06c    -5.83334E-06
 +
40090.06c    -4.97862E-01
 +
40091.06c    -1.09780E-01
 +
40092.06c    -1.69646E-01
 +
40094.06c    -1.75665E-01
 +
40096.06c    -2.89038E-02
 +
50112.06c    -1.27604E-04
 +
50114.06c    -8.83732E-05
 +
50115.06c    -4.59255E-05
 +
50116.06c    -1.98105E-03
 +
50117.06c    -1.05543E-03
 +
50118.06c    -3.35688E-03
 +
50119.06c    -1.20069E-03
 +
50120.06c    -4.59220E-03
 +
50122.06c    -6.63497E-04
 +
50124.06c    -8.43355E-04
 +
 +
% --- Coolant:
 +
 +
mat cool    -0.75 moder lwtr 1001 rgb 50 50 255
 +
1001.06c    0.666666667
 +
8016.06c    0.333333333
 +
 +
% --- Thermal scattering data for light water:
 +
 +
therm lwtr lwj3.11t</nowiki>
 +
</div>
 +
</div>
 +
 +
The base input for the problem is given above. The input describes a 300 cm long fuel rod in infinite lattice geometry. Axially the fuel rod is reflected from top and bottom with 10 cm water layers after which a black boundary condition is applied. The radial geometry is shown here:
 +
 +
[[File:Unstable3Dpin_geom1.png|frameless|200px|xy-plot of the pin-cell geometry.]]
 +
 +
Let's say that we want to calculate the axial power distribution and flux distribution using 100 axial bins over the active fuel length at burnups between 0 MWd/kgU and 20 Mwd/kgU. In order to capture the axial burnup distribution we will divide the fuel material axially into a rather small number of division zones (10) using the [[Input syntax manual#div|div card]].
 +
 +
== Unstable solution using explicit Euler's scheme ==
 +
 +
[[File:UnstableCEmeshAnimation.gif|thumb|150px|Yz-meshplot of the fission heat deposition and thermal flux as a function of burnup (explicit Euler's scheme). The instability of the burnup algorithm is clearly visible.]]
 +
 +
We'll first run the calculation using the explicit Euler's scheme for the discretization of the Bateman equations. This scheme is equal to using the constant values for the beginning-of-step (BOS) flux and cross sections to burn the materials through each of the burnup steps. Explicit Euler's scheme can be chosen with
 +
 +
set pcc ce
 +
 +
Where "ce" refers to "constant extrapolation".
 +
 +
The input used for this calculation is given below:
 +
 +
<div class="toccolours mw-collapsible mw-collapsed" style="width:40em;">
 +
'''Input for explicit Euler's scheme (constant extrapolation) burnup calculation'''
 +
<div class="mw-collapsible-content">
 +
<nowiki>/*****************
 +
** Run options **
 +
*****************/
 +
 +
% --- Neutron population to be used
 +
 +
set pop 10000 500 200
 +
 +
% --- 200 W/cm linear power
 +
 +
set power 60000
 +
 +
/******************
 +
** Burn options **
 +
******************/
 +
 +
% --- Divisor for the fuel material
 +
 +
div fuel subz 10 -150 150
 +
 +
% --- Use constant extrapolation burnup scheme
 +
%    (explicit Euler's method)
 +
 +
set pcc ce
 +
 +
% --- Depletion history
 +
 +
dep butot 0.01 0.25 0.5 1.0 2.0 3.0 5.0 7.5 10.0 12.5 15.0 17.5 20.0
 +
 +
% --- Fission yield and decay data libraries
 +
 +
set nfylib "sss_jeff311.nfy"
 +
set declib "sss_jeff311.dec"
 +
 +
% --- Depletion output for divided materials
 +
 +
set depout 1
 +
 +
% --- Nuclides included in depletion output
 +
 +
set inventory
 +
    531350    % Iodine 135
 +
    541350    % Xenon 135
 +
    621490    % Samarium 149
 +
    922350    % Uranium 235
 +
    942390    % Plutonium 239
 +
 +
/***************
 +
** Detectors **
 +
***************/
 +
 +
% --- Fission heat deposition detector
 +
 +
det power dr -8 void dz -150 150 10
 +
 +
% --- Neutron flux detector
 +
 +
det flux            dz -150 150 10
 +
 +
/*************************
 +
** Geometry definition **
 +
*************************/
 +
 +
% --- Fuel Pin definitions:
 +
 +
pin p1
 +
fuel  0.47
 +
void  0.48
 +
clad  0.54
 +
cool
 +
 +
% --- Lattice
 +
 +
lat l1  1  0.0 0.0 1 1 1.5
 +
p1
 +
 +
% --- Surrounding surfaces:
 +
 +
% Boundary of geometry:
 +
 +
surf 3 cuboid -0.75 0.75 -0.75 0.75 -160 160
 +
 +
% Lower boundary of fuel
 +
 +
surf 4 pz -150
 +
 +
% Upper boundary of fuel
 +
 +
surf 5 pz  150
 +
 +
% --- Cell definitions:
 +
 +
% Active fuel pin
 +
 +
cell  3  0  fill l1  -3  4 -5
 +
 +
% Coolant below active fuel (bottom reflector)
 +
 +
cell  4  0  cool    -3 -4
 +
 +
% Coolant above active fuel (top reflector)
 +
 +
cell  5  0  cool    -3  5
 +
 +
% outside world
 +
 +
cell 99  0  outside  3    % Outside world
 +
 +
% --- Reflective boundary conditions in XY, black in Z:
 +
 +
set bc 3 3 1
 +
 +
/**************************
 +
** Material definitions **
 +
**************************/
 +
 +
% --- Fuel material (4.85 % enrichment):
 +
 +
mat fuel    -10.283 burn 1 vol 208.19 rgb 220 220 115
 +
92235.09c      0.016166667
 +
92238.09c      0.317166667
 +
8016.09c      0.666666667
 +
 +
% --- Cladding (Pure zirconium)
 +
 +
mat clad -6.56000E+00 rgb 180 180 180
 +
40090.06c    -3.32125E+00
 +
40091.06c    -7.32348E-01
 +
40092.06c    -1.13172E+00
 +
40094.06c    -1.17187E+00
 +
40096.06c    -1.92818E-01
 +
 +
% --- Coolant:
 +
 +
mat cool    -0.75 moder lwtr 1001 rgb 50 50 255
 +
1001.06c      0.666666667
 +
8016.06c      0.333333333
 +
 +
% --- Thermal scattering data for light water:
 +
 +
therm lwtr lwj3.11t
 +
 +
mesh 2 300 900</nowiki>
 +
</div></div>
 +
 +
Running the calculation will produce one yz-meshplot for each of the burnup steps (animation shown on right) and already based on these meshplots it is easy to see that the solution is neither symmetric or stable.
 +
 +
In any case, we can plot the axial power distribution and (e.g.) the axial <sup>135</sup>Xe distribution for some of the steps:
 +
 +
[[File:UnstableCEpowers.png|frameless|400px|Axial power distribution of the 3D pin-cell, explicit Euler's method (constant extrapolation).]]
 +
[[File:UnstableCExenons.png|frameless|400px|Axial <sup>135</sup>Xe distribution of the 3D pin-cell, explicit Euler's method (constant extrapolation).]]
 +
 +
These plots show that the solution really is quite unstable. Even a small asymmetry in the fission power produces an asymmetry in the <sup>135</sup>Xe distribution, which in turn affects the power distribution. Setting the equilibrium xenon calculation on (set xenon 1) typically offers some relief, but similar asymmetries in other nuclide densities can destabilize the solution.
 +
 +
== Unstable solution using predictor-corrector scheme ==
 +
 +
== Stable solution using SIE burnup scheme ==
 +
 +
== Stable solution using a symmetry boundary ==
  
 
== References ==
 
== References ==
Line 8: Line 310:
 
<ref name="DufekANE13a">Dufek, J. ''et al.'' "''Numerical stability of the predictor–corrector method in Monte Carlo burnup calculations of critical reactors''", Ann. Nucl. Energy, [http://www.sciencedirect.com/science/article/pii/S0306454913000352 56 (2013) 34-38]</ref>
 
<ref name="DufekANE13a">Dufek, J. ''et al.'' "''Numerical stability of the predictor–corrector method in Monte Carlo burnup calculations of critical reactors''", Ann. Nucl. Energy, [http://www.sciencedirect.com/science/article/pii/S0306454913000352 56 (2013) 34-38]</ref>
 
</references>
 
</references>
 +
 +
[[Category:Example input files]]

Latest revision as of 14:36, 29 November 2017

In certain burnup problems such as long, axially symmetric 3D assemblies or fuel rods or large symmetric core geometries traditional Monte Carlo burnup schemes may run into instabilities[1][2]. This page describes the simulation of one such case first using the traditional burnup schemes (to showcase the instability) and then using a stable burnup scheme.

Problem description

Base input for unstable 3D pin-cell problem 

/*****************
 ** Run options **
 *****************/

% --- Neutron population to be used

set pop 10000 500 200 

% --- 200 W/cm linear power

set power 60000

/*************************
 ** Geometry definition **
 *************************/

% --- Fuel Pin definitions:

pin p1
fuel  0.47
void  0.48
clad  0.54
cool

% --- Lattice 

lat l1  1  0.0 0.0 1 1 1.5
p1

% --- Surrounding surfaces:

% Boundary of geometry:

surf 3 cuboid -0.75 0.75 -0.75 0.75 -160 160

% Lower boundary of fuel

surf 4 pz -150

% Upper boundary of fuel

surf 5 pz  150

% --- Cell definitions:

% Active fuel pin 

cell  3  0  fill l1  -3  4 -5

% Coolant below active fuel (bottom reflector)

cell  4  0  cool     -3 -4

% Coolant above active fuel (top reflector)

cell  5  0  cool     -3  5

% outside world 

cell 99  0  outside   3     % Outside world

% --- Reflective boundary conditions in XY, black in Z:

set bc 3 3 1

/**************************
 ** Material definitions **
 **************************/

% --- Fuel material (4.85 % enrichment):

mat fuel    -10.283 vol 208.19 rgb 200 200 125
92235.09c      0.016166667
92238.09c      0.317166667
8016.09c       0.666666667

% --- Cladding (Zr-4)

mat clad -6.56000E+00 rgb 180 180 180
 8016.06c     -1.19276E-03
 8017.06c     -4.82878E-07
24050.06c     -4.16117E-05
24052.06c     -8.34483E-04
24053.06c     -9.64457E-05
24054.06c     -2.44600E-05
26054.06c     -1.12572E-04
26056.06c     -1.83252E-03
26057.06c     -4.30778E-05
26058.06c     -5.83334E-06
40090.06c     -4.97862E-01
40091.06c     -1.09780E-01
40092.06c     -1.69646E-01
40094.06c     -1.75665E-01
40096.06c     -2.89038E-02
50112.06c     -1.27604E-04
50114.06c     -8.83732E-05
50115.06c     -4.59255E-05
50116.06c     -1.98105E-03
50117.06c     -1.05543E-03
50118.06c     -3.35688E-03
50119.06c     -1.20069E-03
50120.06c     -4.59220E-03
50122.06c     -6.63497E-04
50124.06c     -8.43355E-04

% --- Coolant:

mat cool     -0.75 moder lwtr 1001 rgb 50 50 255
 1001.06c    0.666666667
 8016.06c    0.333333333

% --- Thermal scattering data for light water:

therm lwtr lwj3.11t

The base input for the problem is given above. The input describes a 300 cm long fuel rod in infinite lattice geometry. Axially the fuel rod is reflected from top and bottom with 10 cm water layers after which a black boundary condition is applied. The radial geometry is shown here:

xy-plot of the pin-cell geometry.

Let's say that we want to calculate the axial power distribution and flux distribution using 100 axial bins over the active fuel length at burnups between 0 MWd/kgU and 20 Mwd/kgU. In order to capture the axial burnup distribution we will divide the fuel material axially into a rather small number of division zones (10) using the div card.

Unstable solution using explicit Euler's scheme

Yz-meshplot of the fission heat deposition and thermal flux as a function of burnup (explicit Euler's scheme). The instability of the burnup algorithm is clearly visible.

We'll first run the calculation using the explicit Euler's scheme for the discretization of the Bateman equations. This scheme is equal to using the constant values for the beginning-of-step (BOS) flux and cross sections to burn the materials through each of the burnup steps. Explicit Euler's scheme can be chosen with 

set pcc ce

Where "ce" refers to "constant extrapolation".

The input used for this calculation is given below:

Input for explicit Euler's scheme (constant extrapolation) burnup calculation 

/*****************
 ** Run options **
 *****************/

% --- Neutron population to be used

set pop 10000 500 200

% --- 200 W/cm linear power

set power 60000

/******************
 ** Burn options **
 ******************/

% --- Divisor for the fuel material

div fuel subz 10 -150 150

% --- Use constant extrapolation burnup scheme
%     (explicit Euler's method)

set pcc ce

% --- Depletion history

dep butot 0.01 0.25 0.5 1.0 2.0 3.0 5.0 7.5 10.0 12.5 15.0 17.5 20.0

% --- Fission yield and decay data libraries

set nfylib "sss_jeff311.nfy"
set declib "sss_jeff311.dec"

% --- Depletion output for divided materials

set depout 1

% --- Nuclides included in depletion output

set inventory
    531350    % Iodine 135
    541350    % Xenon 135
    621490    % Samarium 149
    922350    % Uranium 235
    942390    % Plutonium 239

/***************
 ** Detectors **
 ***************/

% --- Fission heat deposition detector

det power dr -8 void dz -150 150 10

% --- Neutron flux detector

det flux             dz -150 150 10

/*************************
 ** Geometry definition **
 *************************/

% --- Fuel Pin definitions:

pin p1
fuel  0.47
void  0.48
clad  0.54
cool

% --- Lattice

lat l1  1  0.0 0.0 1 1 1.5
p1

% --- Surrounding surfaces:

% Boundary of geometry:

surf 3 cuboid -0.75 0.75 -0.75 0.75 -160 160

% Lower boundary of fuel

surf 4 pz -150

% Upper boundary of fuel

surf 5 pz  150

% --- Cell definitions:

% Active fuel pin

cell  3  0  fill l1  -3  4 -5

% Coolant below active fuel (bottom reflector)

cell  4  0  cool     -3 -4

% Coolant above active fuel (top reflector)

cell  5  0  cool     -3  5

% outside world

cell 99  0  outside   3     % Outside world

% --- Reflective boundary conditions in XY, black in Z:

set bc 3 3 1

/**************************
 ** Material definitions **
 **************************/

% --- Fuel material (4.85 % enrichment):

mat fuel    -10.283 burn 1 vol 208.19 rgb 220 220 115
92235.09c      0.016166667
92238.09c      0.317166667
8016.09c       0.666666667

% --- Cladding (Pure zirconium)

mat clad -6.56000E+00 rgb 180 180 180
40090.06c     -3.32125E+00
40091.06c     -7.32348E-01
40092.06c     -1.13172E+00
40094.06c     -1.17187E+00
40096.06c     -1.92818E-01

% --- Coolant:

mat cool     -0.75 moder lwtr 1001 rgb 50 50 255
 1001.06c      0.666666667
 8016.06c      0.333333333

% --- Thermal scattering data for light water:

therm lwtr lwj3.11t

mesh 2 300 900

Running the calculation will produce one yz-meshplot for each of the burnup steps (animation shown on right) and already based on these meshplots it is easy to see that the solution is neither symmetric or stable.

In any case, we can plot the axial power distribution and (e.g.) the axial 135Xe distribution for some of the steps:

Axial power distribution of the 3D pin-cell, explicit Euler's method (constant extrapolation). Axial 135Xe distribution of the 3D pin-cell, explicit Euler's method (constant extrapolation).

These plots show that the solution really is quite unstable. Even a small asymmetry in the fission power produces an asymmetry in the 135Xe distribution, which in turn affects the power distribution. Setting the equilibrium xenon calculation on (set xenon 1) typically offers some relief, but similar asymmetries in other nuclide densities can destabilize the solution.

Unstable solution using predictor-corrector scheme

Stable solution using SIE burnup scheme

Stable solution using a symmetry boundary

References

  1. ^ Dufek, J. and Hoogenboom, E. "Numerical Stability of Existing Monte Carlo Burnup Codes in Cycle Calculations of Critical Reactors", Nucl. Sci. Eng., 162 (2009) 307-311
  2. ^ Dufek, J. et al. "Numerical stability of the predictor–corrector method in Monte Carlo burnup calculations of critical reactors", Ann. Nucl. Energy, 56 (2013) 34-38
Retrieved from "https://serpent.vtt.fi/mediawiki/index.php?title=Unstable_3D_pin-cell_burnup_problem&oldid=4005"