DIFFUSION PIECEWISE HOMOGENIZATION VIA FLUX DISCONTINUITY RATIOS

Nuclear Engineering and Technology.
2013.
Nov,
45(6):
707-720

- Received : September 26, 2013
- Accepted : October 02, 2013
- Published : November 25, 2013

Download

PDF

e-PUB

PubReader

PPT

Export by style

Share

Article

Metrics

Cited by

TagCloud

We analyze piecewise homogenization with flux-weighted cross sections and preservation of averaged currents at the boundary of the homogenized domain. Introduction of a set of flux discontinuity
ratios
(FDR) that preserve reference interface currents leads to preservation of averaged region reaction rates and fluxes. We consider the class of numerical discretizations with one degree of freedom per volume and per surface and prove that when the homogenization and computing meshes are equal there is a unique solution for the FDRs which exactly preserve interface currents. For diffusion submeshing we introduce a Jacobian-Free Newton-Krylov method and for all cases considered obtain an ‘exact’ numerical solution (eight digits for the interface currents). The homogenization is completed by extending the familiar full assembly homogenization via flux discontinuity factors to the sides of regions laying on the boundary of the piecewise homogenized domain. Finally, for the familiar nodal discretization we numerically find that the FDRs obtained with no submesh (nearly at no cost) can be effectively used for whole-core diffusion calculations with submesh. This is not the case, however, for cell-centered finite differences.
where
x
= (
r
,
E
,
Ω
) is a point in phase space
X
=
D
×
ε
×
S
_{2}
of entering boundary Γ_= {
x
,
r
∈
∂D
,
E
∈
ε
,
Ω
∈
S
_{2}
,
Ω
⋅
n
< 0},
D
= geometric domain,
ε
= energy domain,
S
_{2}
= unit sphere and
∂D
= boundary of
D
with outward normal
n
. Also,
B
=
Ω
⋅▽+Σ–
H
is the transport operator with
H
representing scattering, λ is a positive real number,
P
is the fission production operator and
ψ_{in}
is the angular flux entering the sub domain. We note that Eq. (1) is a source problem with fixed value λ and a surface source given by the incoming flux
ψ_{in}
. Furthermore, if we take λ and
ψ_{in}
equal to the values λ
_{core}
and (
ψ_{in}
)
_{core}
, that would be obtained from the hypothetical solution of the whole-core transport equation, then the solution of this local source problem is the exact angular flux in the subdomain that would be obtained from the whole-core transport equation.
(8)
Although the need for the exact λ
_{core}
and (
ψ_{in}
)
_{core}
invalidates the use of local problem (1) to compute the
ψ_{core}
in the sub domain, this problem provides the basis for a homogenization paradigm. By supplying a λ for the core eigenvalue and an approximate model for the boundary condition
ψ_{in}
, the physicist can create a Reference Homogenization Problem (RHP) and, by requesting that the solution of a homogenized version of the RHP with a low-order operator preserves the reaction rates provided by the transport solution, he can implement a method to determine the homogenization parameters to be used with the low-order operator in the final whole-core calculation. The intuitive idea is that the RHP gives a good approximation for the true transport flux in the reactor core and that small ‘errors’ in λ
_{core}
and (
ψ_{in}
)
_{core}
result in small errors for the homogenized parameters, i.e., the homogenization procedure is continuous on the data. In practice the sub domain is an assembly or a colorset, the calculations are carried out in two-dimensions and an infinite lattice approximation is often used to model the incoming flux. This defines an RHP which is customarily complemented with an approximate critical buckling term in order to ensure λ = 1, as for a critical core. Although the infinite-lattice approximation is often far from the actual boundary condition of the sub domain in the reactor core, there are theoretical arguments that justify this choice for uniform cores
(8)
and this model is traditionally adopted.
In any case, regardless of the approximations adopted for λ and
ψ_{in}
, classical homogenization theory is based on constructing homogenized parameters from the solution of the RHP for use in low-order, whole-core calculations. The former are obtained by constraining the solution of the homogenized RHP to preserve reference reaction rates. When diffusion is the low-order operator, the solution Φ
^{h}
of the homogenized RHP
where
h
denotes homogenized data and results,
x^{h}
= (
r
,
E
),
D^{h}
is the diffusion coefficient and
Q^{h}
= [
H^{h}
+(1/λ)
P^{h}
]Φ
^{h}
is the source term, must preserve the few-group reaction rates over a set of homogeneous regions
where
is the volume-averaged diffusion flux in region
i
and group
g
,
V_{i}
is the region volume and
τ^{g}_{i}
is the reference transport reaction rate. The boundary condition in the diffusion equation is an approximation of the boundary condition of the RHP transport problem which preserves the averaged normal transport current over each side on the boundary,
where
k
indicates a side of a homogenized region lying on the boundary of the domain.
Typically the sub domain
D
is rectangular, the regions are defined from an
N
×
M
Cartesian Homogenization Mesh (HM) and the multigroup diffusion problem is solved by a numerical spatial discretization.
The homogenized cross sections for all other reactions are proportional to the homogenized total cross sections via coefficients that are ratios of transport reaction rates: for a reaction
r
and for multigroup transfer we have
respectively. Hence, the source term
Q^{h}
depends only on the set of all homogenized total cross sections
RHP (Left) and Homogenized RHP (Right) Obtainedwith a 4×4 HM.
One can then formalize the problem of homogenization as that of computing a set
of homogenized parameters solution of nonlinear equations (3), where the
are the region-averages of the flux
solution of (2). The vector
comprises all the total cross sections
and possibly other parameters. Let
be a hypothetical solution of nonlinear equations (2,3) and replace in Eq. (2) the
Q^{h}
with the exact source term
The equation so obtained is unrealistic because the source cannot be known if the multigroup solution for
is not known, but together with the constraints in (3) these equations are equivalent to the initial ones. However, with multigroup coupling replaced by a known source, the initial nonlinear system of equations splits now into a set of independent nonlinear equations for energy group. Hence, from now on we shall consider Eqs. (2,3) for a fixed energy group
g
as a nonlinear source problem defining the
for the
NM
regions in the domain. In the following the group index
g
should be omitted except if otherwise necessary.
We are now in a position to show that the constraints in (3) are not independent and that there is in fact one degeneracy per group. By integration of Eq. (2) over the domain we see that the solution
of this equation satisfies the global balance
and this regardless of the
and therefore of the constraints in (3). To write this equation we have used the fact that the sources and boundary conditions of equations (2) preserve transport averages over regions and sides, respectively, to replace
by the transport values
Also, in the equation
A_{k}
is the area of side
k
and the
are calculated in the outward normal direction. Finally, comparison of Eq. (4) with the equivalent equation satisfied by the reference transport solution gives
and, since this relation is satisfied regardless of
, we can safely conclude that the number of linearly independent constraints in (3) is
NM
– 1.
In this work we have analyzed two homogenization models based on diffusion equations (2) and constraints (3). Both models share the following two properties: a) the diffusion coefficients are obtained from the transport data or from independent data, and b) both models enforce the continuity of the normal component of the current at regions interfaces.
But, prior to the discussion of these models, we shall analyze how the solution
of (2,3) constraints the averaged values of the interface currents
, a point which will help clarify the differences between the two models. We start by observing that, in a similar manner as we just did with global balance to derive Eq. (4), we can use local region balance together with constraints (3). The result is an equivalent constraint for the total current exiting each homogenized region,
where the sum in
k
is over the sides of region
i
, the
are the normal total currents leaving through the cell sides and the
are reference transport values. Hence, each reactionrate constraint translates into a linear constraint for the currents exiting the region. The number
N_{f}
of degrees of freedom for the internal
is then obtained by subtracting the number
NM
– 1 of linearly independent reaction-rate constraints from the total number
of interface currents at internal region sides. The result,
N_{f}
= (
N
– 1)(
M
– 1), shows that with
N
= 1 and/or
M
= 1 all interface currents are constrained to their transport values,
. However, for
N
,
M
> 1 the interface currents are not necessarily constrained to their transport values, even though the conservative solution
remains among the set of all admissible solutions.
We now return to the analysis of the two homogeneous models.
where right and left indicates the limit values at both sides of interface
k
. The homogenization vector is
and the number of unknowns is
NM
per group. Then, since the number of independent reaction-rate constraints is
NM
> 1, the system of equations for the
is undetermined and, to select a single solution, one adds a supplementary constraint, which is often that of preserving the domain-averaged transport flux,
The solution
is obtained from the minimization of a functional of the form
where
with the
solution of (2).
are directly defined by flux weighting,
and one seeks to introduce new homogenization parameters in order to preserve the averaged interface currents
which, given that the averaged region sources
are preserved, results in region reactions rates and regionaveraged fluxes conservation,
τ^{h}
=
τ
and
We argue next that the hypothetical solution of this problem does not necessarily preserve flux continuity at region interfaces. Take for example the class of numerical approximations with one degree of freedom at the interfaces, which without loss of generality we take to be the interface averaged current, and consider the case of an isolated region. The input values of the numerical model are the averaged interface currents
, that serve here as boundary conditions, and the averaged source
; therefore, the output side-averaged fluxes are a function of the input parameters, the region cross section and the region dimensions,
Since this applies also to the neighboring regions which, except for the common interface current, may have different parameters, it is clear that this class of numerical discretizations cannot support a solution which preserves interface currents and region reaction rates without allowing the fluxes to be discontinuous at region interfaces.
Hence, in order to support these types of solutions one introduces a set of Flux Discontinuity Ratios (FDR)
to enforce the flux interface condition
where right denotes the direction of increasing coordinate values. In the present case the homogenization vector
has
N_{int}
unknowns per group and the solution can be sought by minimization of the functional
where the
are the averaged interface currents obtained from the solution
of (2) with the flux interface conditions in (5).
We note that in this case the number of unknowns equals the number of constraints on the interface currents so the problem, although nonlinear, is well determined.
where + (–) denotes the right (left) sides of the region,
is the region averaged flux,
J
± are the surface averaged interface currents and
α_{k}
= 2
D^{h}
/ Δ
_{k}
with Δ
_{k}
= width of the cell in the
k
direction.
On the other hand, a finite approximation is used for the transverse averaged flux in the NEM. In the case of a parabolic approximation, the coefficients of the expansion can be calculated in terms of the transverse currents and the region averaged flux, and this leads to an explicit formula similar to the previous one:
However, this cannot be done with the usual quartic expansion, which brings in the source term and requires a parabolic buckling approximation based on the preservation of neighboring region transverse leakages.
^{1}
The result is a multigroup system of equations for the interface averaged fluxes
, coupled via the source term, and the inversion of this system gives the
’s in terms of the multigroup region-averaged fluxes of the region and the multigroup interface currents for the region and its neighbor regions.
(12)
Finally, with only one-degree of freedom per interface, this leads to a constructive proof of the existence of an exact solution of problem (2-3,5) with the ratios provided by the use in (5) of the single-region values for the
which are directly obtained from the one-region explicit calculations:
However, in the case of submeshing this construction is not possible because the shape of the current is not known at the interfaces.
^{2}
We have then applied a Jacobian-Free Newton Krylov (JFNK) method to the minimization of functional (6). The Krylov stage was based on GMRES and, in order to compute the finite-differences directional derivatives, we adopted double precision throughout our FORTRAN90 computer code NEM4. Among other problems, we have extensively analyzed the FDR homogenization of a colorset with reflective boundary conditions adapted from reference (13) and depicted in
Fig. 2
.
Figure 3
shows the reference two-group fluxes and the relative errors in the pin averaged fluxes obtained with the 4×4 homogenization with unit FDRs; i.e., with the familiar flux continuity interface condition. The strong flux gradient in the thermal group and the large errors in the thermal pin fluxes witness to the difficulty of the homogenization. Similar errors were also observed with the pin-by-pin 14×14 homogenization.
We have computed the FDRs at each interface of the four homogenized configurations shown in
Fig. 2
. We have run calculations with and without submeshing. In all cases all interface currents converged to the sought precision of 10
^{–8}
. For these calculations we required an absolute precision of 10
^{–10}
for the convergence of the functional and computed the diffusion eigenvalue and fission source with a relative precision of 10
^{–12}
and 10
^{–11}
, respectively; the relative precision in the region reaction rates varied from 10
^{–8}
for the 4×4 homogenization to 10
^{–6}
for the pin-by-pin one. The loss of precision is mainly due to inconsistencies in the number of significant digits considered for the transport input data.
Left: Model PWR Colorset with four DifferentAssemblies. A: Quarter of MOX Assembly with 2.359(Black), 2.993 (Dark Grey) and 4.252 (Light Grey) % ofFissile Pu. B, C and D: Quarters of UOX Assemblies with1.053, 3.400 and 1.349 % U235 Enrichment. The RHP is aCritical Leakage Lattice Problem. Right: Homogenized RHPMeshes that have been Analyzed.
Reference Averaged Fluxes (Top) and Relative Errors in the Pin Averaged Fluxes (Bottom) as Obtained from the 4×4 Homogenization Problem with Unit FDFs.
All the NEM calculations presented hereafter have been obtained with the quartic expansion of the nodal diffusion solver NEM4.
Table 1
reflects the performance of the JFNK minimization procedure and of a preconditioned JFNK. The sub meshing calculation for the 2×2 HM were run for all the six compatible CMs but only one intermediary mesh will be included in the tables. The number of JFNK steps depends only on the number of unknowns
n
of the problem, but the cost of the calculation increases with the size of the multigroup diffusion problems that are solved at each GMRES iteration and can be as high as 1.73 hours in a desktop computer for the pin-by-pin homogenization. In order to reduce the computing time we have implemented a crude left preconditioning using as preconditioner the inverse of the Jacobian matrix computed at the beginning of the Newton search. This reduces the previous 1.73 hours to 62 sec.
Next, we investigate how the FDRs vary with submeshing.
Table 2
gives the values for the FDRs for the 2×2 HM
Performance of Homogenization Procedure using NEM. p-JFNK = Preconditioned JFNK Iterations. HM = Homogenization Mesh, n = Number of Unknowns, CM = Diffusion Mesh, #GMRES = Total Number of GMRES Iterations,t (s) = Time in Seconds, t (%) = Percentage of t Saved.
FDRs Obtained with the FD and NEM Discretizations for the 2×2 HM. The Interface Notation S, N, W and E is the One in Fig. 2
for the FD and NEM discretizations. The results in this table show that the FDRs vary appreciably with the FD discretization but remain approximately constant with the NEM one. The reason is that the NEM discretization converges very fast. The lower FD discretizations are so bad that they need distorted FDR values in order to preserve interface currents, but for refined meshes the FD discretization becomes better and its FDRs are close to the NEM ones.
The implication is that one might use the FDRs obtained with a coarser CM rather than those computed with the mesh used for the final core calculation. To test this hypothesis we have calculated the error introduced in reaction-rate conservation. Our results show that calculations for a submeshing using the FDRs of the case with no submeshing introduces acceptable errors in the final reaction rates. For example for the pin-by-pin submesh of the 2×2 HM, the relative errors in the NEM fluxes did not exceed 0.15% and 0.4% for the fast and the thermal groups, respectively. This allows using the FDRs from the case with no submesh, computed nearly at no cost, for any other sub meshing. However, this does not apply to the FD discretization because in this case the errors in the fluxes are not acceptable.
^{1} As an aside we mention that this leakage approximation is equivalent to introducing a spatial discretization of the diffusion equation with extra ‘boundary’ terms, i.e., the transverse leakages for the adjacent regions.
^{2} This is also the case with no submeshing but with a spatial discretization with more than one degree of freedom per interface.
where
k
denotes the side of a region on the boundary between two separately homogenized sub domains, so that
The FDFs
f
are then calculated separately in their own RHPs. Because the latter are critical problems, the FDFs must be defined so they do not depend on normalization. In this work we have considered two boundary homogenization models. In the Generalized Equivalence Theory (GET) the FDFs are computed to preserve the ratio of the reference to the diffusion boundary fluxes,
and are basically intended to preserve the limits of the asymptotic flux at the boundary between the two homogenized sub domains. On the other hand, preservation of the local albedo is used to determine the FDFs for the Black-Box (BB) model. With the albedo defined as the ratio of the partial outgoing (+) to the partial incoming (–) currents, and after use of the diffusion boundary condition, 4
J
^{±}
_{k}
=
f_{k}
Φ
^{h}_{k}
± 2
J_{k}
, one obtains
We note that for full-assembly homogenization with a perfectly reflected RHP one has Φ
^{h}_{k}
= Φ
_{k}
and, therefore, the FDFs are independent of the diffusion mesh.
F
. We have proved by simple construction that, when the computing mesh is equal to the homogenization mesh, an exact solution
F
= 0 exists and that the FDRs can be obtained from the simple separated region problems. This has been demonstrated for the finite differences and for the nodal expansion method but we think that the result could also be proved for other numerical discretizations such as the nodal analytical method. We have applied a Jacobian-Free Newton Krylov technique to the minimization of the functional and for the case of submeshing we have obtained a numerical zero of the functional for all cases considered. Thus, contrarily to the equivalence technique for piecewise homogenization, we think that we can assume with confidence that the exact solution for the FDR homogenization problem exists (we recognize that this conclusion contradicts the expectations that two of the authors made in a previous conference paper.
(3)
) Moreover, we have completed our FDR homogenization technique to include the calculation of FDRs at interface between subdomains separately homogenized in a consistent way with the classical techniques of fullassembly homogenization with flux discontinuity factors. We have also shown that for coarse homogenization meshes with quartic nodal discretization one can use without significant loss of precision the FDRs and FDFs obtained with no submeshing. This means that the coefficients can be recovered from explicit, simple calculations performed independently for each homogenized region. Therefore, there is no need for costly JFNK nonlinear iterations and no need also for a version of the diffusion program in double precision. This also applies to pin-by-pin homogenization with finite differences.
The use of FDRs may have an important impact in storage requirements for core burnup calculations. For these calculations a library, containing homogenization data parameterized in terms of physical local variables (burnup, fuel and coolant temperatures, etc.), has to be stored prior to the core calculations. For each homogenization one has to store
cross sections, 2
NM
–
N
–
M
FDRs and 2(
N
+
M
) FDFs for a total of
homogenized parameters per group. Here
is the average number of isotopic data per homogenized region and per group, which depends on the average number of isotopes per region
N_{isot}
and on the average size of cross section data per isotope and per group
N_{dat}
. This results in the following proportion of FDR and FDF data:
which as
N
and
M
increase goes to
For precise calculations one can decide to store as many as 100 to 150 isotopes for which absorption, fission, scattering and transfer microscopic cross sections are required and therefore
p
is very small. However, for less precise calculations most of the isotopes in the regions may be represented by a single set of macroscopic cross sections while only a few isotopes are represented individually by their microscopic cross sections and the storage of FDRs and FDFs might significantly increase the size of the parameterized library. For full assembly homogenization this amounts to
and for coarse isotopic representations, i.e.,
small, the contribution of the FDRs and FDFs may be significant. This is a motivation for coarse piecewise homogenization, whereas a pin-by-pin homogenization will demand a very large library. It is important to note that other homogenization paradigms can be used to minimize FDR storage, albeit to a prize in the precision of the homogenization. For example, Herrero
(5)
has proposed a fitting based on the introduction of a linear approximation for the FDFs in terms of local dynamic variables whose coefficients are determined from a few simplified pin environments. These local variables comprise the ratios of the reference interface flux and current to the diffusion interface flux and a local leakage correction term. A problem though with this approach is that it requires nonlinear iterations for the FDFs in the solution of the multigroup diffusion core calculations, and this considerably increases computing time. We think that the computation of the FDRs' environments could be advantageously be done with FDR piecewise homogenization which assures a better homogenization approximation.
x
a generic direction (either
x
,
y
, or
z
) of unit vector
e
_{x}
,
x
∈ (
x
_{–}
,
x
_{+}
) within the region, Δ
_{x}
=
x
_{+}
–
x
_{–}
the region width in direction
x
,
A
_{x}
the area of the side perpendicular to
x
and
V
= Δ
_{x}
A
_{x}
the region volume. The diffusion equation is considered separately for each homogeneous region,
together with boundary and interface conditions. The latter may account for flux discontinuity. At an interface orthogonal to direction
x
we must have at any position
r
:
where
n
denotes the neighboring node,
α
is the flux discontinuity coefficient and Φ(
r
) and
J
_{x}
(
r
) =
e
_{x}
·
J
(
r
) are the limit values for the flux and the normal net current from the interior of the node. In practical applications these interface conditions are applied to the averaged transfer currents and fluxes defined in (14) and (16).
Two Neighboring Nodes showing Current Continuity, J _{x} = (J _{x})_{n} , and Flux Discontinuity, f Φ = (f Φ)_{n} , at the Interface.
We observe that the interface condition for the flux depends only on the ratio
between the discontinuity coefficients at the interface and, therefore, the solution of the homogenized diffusion equation, depend only on these ratios at each interface. Note that the ratios can be computed from the relation of the averaged interface fluxes resulting from the solution of Eqs. (7,8) with boundary condition (12) and that this fact leads to an iterative computation of the ratios.
Regarding boundary conditions, we consider a general condition of the form
where
J
_{x}
is the net current leaving the domain, Φ is the flux at the boundary and
J
_{0}
and
γ
are parameters. For
γ
= ∞ we have the zero flux boundary condition Φ = 0, for
γ
= (1/2) (1 –
β
)/(1 +
β
) and
J
_{0}
= 0 we get the albedo boundary condition
for the partial currents leaving (
out
) and entering (
in
) the domain,
J^{out}
= Φ/4 +
J
_{x}
/2 and
J^{in}
= Φ/4 –
J
_{x}
/2. Finally, for
γ
= 0 we have a net current condition
which for
J
_{0}
= 0 gives the reflection condition
J
_{x}
= 0. Boundary conditions (12), with
J
_{0}
equal to the reference transport boundary current, are the one typically applied for the homogeneization problem that it is used for the determination of the flux discontinuity ratios at interfaces between homogenized regions.
and the node surface currents
J
_{x}
_{±}
=
J
_{x}
(
x
_{±}
),
where
is the averaged transverse nodal current in direction
x
. In this equation
x
_{⊥}
(
x
) denotes the transverse node area orthogonal to direction
x
. By using Fick's law we can write
where
is the averaged transverse nodal flux. We note that we must have
The solutions Φ
_{x}
(
x
) and
J
_{x}
(
x
) for consecutive nodes along direction
x
are connected by the transverse averaging of the interface conditions in (8). Writing the nodal equations in a response matrix form facilitates the incorporation of the interface conditions. The response matrix formulation requires two relations per direction giving the partial currents exiting the node via the faces normal to the direction in terms of the incoming partial currents and the flux in the node. We write these two equations at
x
_{±}
:
which, together with the supplementary equations
completely specify the partial currents in terms of the interface flux and current. The sought response matrix equations will also require expressing the interface conditions in terms of partial currents. By using the precedent equations with interface conditions (8) we obtain
where
n
denotes the neighbor node at
x
± and
f
and
f_{n}
must be understood as
f
_{±}
and (
f
_{∓}
)
_{n}
, respectively. These equations allow writing the incoming current in terms of the exiting current from the neighbor node and either the outgoing current from the node or the entering neighbor current. We choose to minimize the coupling with the neighbor node:
In order to preserve relation (17) we select to have
It is then convenient to use a truncated expansion in Legendre polynomials with
where
This yields the expansion
This expansion can now be used to compute Φ
_{x}
_{±}
and
J
_{x±}
= –
D
Φ'
_{x}
(
x
±) in equations (19) and (18) that, after elimination of the flux expansion coefficients Φ
_{x,n}
for
n
> 0, will provide the response matrix formulation for the partial exiting currents which, together with nodal balance (13), will close the system of equations. We write then
where we have noted that
∂_{x} f
(
z
) = (2/Δ
_{x}
)
f '
(
z
) and used the properties of the Legendre polynomials
P_{n}
(±1) = (±)
^{n}
and
P'_{n}
(±1) = (±)
^{n+1}
n
(
n
+ 1)/2. Also, in these equations
We shall exclusively consider the case
N
= 4, for which elimination of the four flux components requires the two relations in (19) plus two additional equations. The latter are obtained from projection of the nodal transverse equation.
^{3}
The nodal transverse diffusion equation is a onedimensional diffusion equation obtained by transverse integration of the original diffusion equation (13). Noticing that ▽⋅ = (
∂_{x}
e
_{x}
+▽
_{⊥}
)⋅ and by averaging (13) over the transverse area we obtain
where Φ
_{x}
(
x
) and
J
_{x}
(
x
) are the averaged transverse nodal flux and current in (16) and (15),
Q
_{x}
(
x
) is the averaged transverse nodal source, similarly defined, and
is the averaged transverse leakage.
To derive two supplementary equations for Φ
_{x,3}
and Φ
_{x,4}
we use a projective technique where the nodal expansion (21) is introduced in the transverse nodal equation and the resulting equation is then multiplied by
P_{m}
(
z
) for
m
= 1, 2 and integrated over
x
.
^{4}
The replacement of expansion (21) in (15) results in
where it has been assumed that the source
Q_{x}
has also a Legendre expansion. Then, projecting and using the orthogonality of the Legendre polynomials leads to
where
Also, in this equation we have noticed that for
n
= 2, 4 and
m
= 1, 2 we have (
P_{m}
,
P'_{n}'
) =
γ_{m}δ_{m}
_{+2,n}
with
γ
_{1}
= 5 and
γ
_{2}
= 7. Note that the last equation allows to explicit Φ
_{x,2}
and Φ
_{x,4}
in terms of Φ
_{x,1}
and Φ
_{x,3}
, respectively:
where
and
with, according to the leakage model in (33),
L_{x,m}
= (2
_{m}
+ 1)(
P_{m}
,
L_{x}
).
n
= 1, 4. By using (26) we can recast (22) as
Then, Eqs. (19) are used to determine Φ
_{x,1}
and Φ
_{x,2}
,
and the resulting expressions are then used back into the second equation in (28) to yield:
where
c
_{x,1}
= (1 + 6
α
_{1}
σ
)/(1 +
α
_{1}
σ
),
c
_{x,2}
= (3 + 10
α
_{2}
σ
)/(1 +
α
_{2}
σ
) and
Next, (29) is fed into Eqs. (18) to yield a system of two equations of the form
where
and where we have noted
a_{x}
= (2
D
/Δ
_{x}
)(
c
_{x,1}
+
c
_{x,2}
),
b_{x}
= (2
D
/Δ
_{x}
)(
c
_{x,1}
–
c
_{x,2}
).
Finally, the coupling between a node and its neighbors is mediated by interface conditions (20) which in vector notation read
with
where
r
_{±}
denote the ratios of the discontinuity coefficient of the node divided by that of the neighbor node at the node interfaces
x
_{±}
.
By replacing this last relation in (30) we obtain a tridiagonal system of equations with two-dimensional unknowns,
which can be solved by a forward and backward iteration. In the last equation:
l_{x,k} = -B_{x}M_, u_{x,k} = -B_{x}M_{+}, d_{x,k} = A_{x} - B_{x}M.
and introduce the forward elimination
where
α_{k}
and
β_{k}
are 2×2 matrices, leading to the recursive formulas
α_{k}
= –
D_{k}u_{k}
and
β_{k}
=
D_{k}
(
q_{k}
–
l_{k}β_{k}
_{–1}
), where
D_{k}
= (
l_{k}α_{k}
_{–1}
+
d_{k}
)
^{–1}
with initial values,
D
_{1}
=
d
^{1}
_{–1}
,
α
_{1}
= –
D
_{1}
u_{k}
and
β_{k}
=
D
_{1q1}
.
where for example
Note that the average node value is
We now introduce a quadratic expansion of the form
and determine the coefficients
L_{xy,n}
for
n
= 1, 2 by imposing the condition that the node averaged values of this expansion in the neighboring nodes
k
±1 (noted hereafter with ±) along direction
x
equal the corresponding average leakages
This gives the conditions
where, by adopting the change of variable
which implies the change
z
→ ±
z'
, we have
with α
_{±}
= 2(Δ
_{x}
)
_{±}
/Δ
_{x}
. This results in (
γ
_{x,1}
)
_{±}
= ±(2 + α
_{±}
)/2 and (
γ
_{x,2}
)
_{±}
=(α
^{2}
_{±}
+ 3α
_{±}
+ 2)/2, which for a constant mesh size gives (
γ
_{x,1}
)
_{±}
= ±2 and (
γ
_{x,2}
)
_{±}
= 6.
The final result for the total transverse leakage is
where
with
Another particular case of the general albedo condition in (11) is the vacuum boundary condition
J^{in}_{x}
= 0 resulting from
β
= 0.
To complete the equations it remains to compute the interpolate for the transverse leakage. For the case of reflection or, in general, for a geometrical motion, there is no problem because the assembly has a ‘neighbor’ via the motion, which can then be used for the interpolation. This is not the case, however, for the general albedo boundary condition and in the numerical implementation we have followed the rules given by Kord Smith for the vacuum case.
(14)
where
d
is the reduced coefficient in (23). Also, in Eq. (34) we have used a first-order, finite-differences approximation, between the center of the cell
and the interface
x
_{+}
, for the derivative of the flux and we have assumed a linear variation so that
With the help of (34), the interface conditions in (8) lead to an expression for the interface averaged flux Φ and current
J
_{x}
, which, in turn, yields the inner interface flux and the current at the interface:
where
r
is the ratio of the discontinuity coefficient of the cell to that of the neighboring cell in (9).
J
_{x}
be the side averaged flux and current. Replacing the finite-difference approximation for
J
_{x}
in Eq. (34) in the general expression (10) for the boundary condition yields
and
For
γ
= 0 we have
and
J
_{x}
=
J
_{0}
so that
J
_{0}
stands for the total boundary current. For
γ
= ∞ we have the zero-flux boundary condition Φ = 0 with
Finally for the albedo boundary condition with
J
_{0}
= 0 and
γ
= (1/2) (1 –
β_{ab}
)/(1 +
β_{ab}
) we have
and
J
_{x}
=
γ
Φ.
where
denotes the cell-averaged flux in cell
a
, the sum is over the neighboring cells and
In these equations the sums in
b
∈
int
and
b
∈
bd
are over the neighboring cells and over the faces of the cell on the boundary of the domain, respectively, the lower index
ab
represent the value of a quantity at interface
ab
,
A
_{(ab)}
stand for the area of interface
ab
and
J
_{0ab}
is the averaged current at interface
ab
on the boundary of the domain.
The equation for cell
a
contains a diagonal term and 2
n_{d}
off diagonal terms, where 2
n_{d}
is the number of neighbors and
n_{d}
is the dimension of the geometry (1, 2 or 3). Note that the off-diagonal terms are negative while all the contribution to the diagonal terms are positive. We also have the relations
Finally, note that
Hence, system of equations (39) is neither symmetric nor has diagonal dominance, except if
r_{ad}
= 1 ∀
a
,
b
. We note that the solution of diffusion equations (39) depends only on the ratios and not on the values of the discontinuity coefficients and that, once the cell-averaged fluxes have
Two Neighboring Homogenized Cells a and b . The Dots Indicate the Intersection with the Interface of the RefinedComputing Mesh.
been calculated, cell interior interface fluxes and interface currents can be recovered from Eqs. (35,37) and (36,38) using only the ratios of the discontinuity coefficients.
Assume cells
a
and
b
are homogenized cells and that
r_{ab}
is the ratio of discontinuity coefficients along their common side
ab
. When a submesh is used to decompose further the cells for the final diffusion calculation one has
where the sum in
r
is for all neighbor calculation regions along the common boundary
ab
,
A
_{(ab)}
(
r
) is the area on the interface common to regions of index
r
with
A
_{(ab)}
= Σ
_{r}
A
_{(ab)}
(
r
), and
are the respective region-averaged fluxes.
^{3} It should be noted that if one uses a quadratic expansion (N = 2), then there is no need of adding two new equations and therefore there is no need to introduce or discuss the transverse nodal equation, or to introduce a model for the transverse leakage.
^{4} Note that the moment n = 0 gives Eq. (13) and does not provide a relation for the flux expansion coefficients.

1. INTRODUCTION

Traditionally, reactor core calculations are based in a simplified few-macrogroup 3D diffusion calculation with fully homogenized assemblies and a homogenized reflector. The homogenization problem consists of calculating in transport an assembly in an infinite lattice, often with critical leakage, and deriving homogenized parameters that reproduce assembly reaction rates. Most often these parameters are flux-weighted cross sections and Flux Discontinuity Factors (FDF) on the assembly surfaces, and the latter are obtained using a physical model that predicts how the assembly couples with neighboring assemblies in the core. Lately, the increase in computer power has awakened an interest in more-detailed core calculations based on perquarter or pin-by-pin assembly homogenization
(1)
, which until now were only used for limited applications. Calculation of FDFs for pin-by-pin diffusion calculations has been done either directly by using the reference transport currents on the boundary of the pins
(2)
,
(3)
,
(4)
or by zero-leakage calculations with corrections that account for the impact of the cell environment and led to an iterative calculation in the final 3D core diffusion calculation.
(5)
In this work we adopt the classical homogenization paradigm where homogenized parameters are sought so the reaction rates obtained from the diffusion calculation of the homogenization problem reproduce the transport reaction rates for each homogenized region of the assembly. To do so we extend the standard full-assembly homogenization technique based on Flux Discontinuity Factors for use in piecewise homogenization but with two differences. The first is that we recognize that the solution of the diffusion equation with flux discontinuity interface conditions does not depend on the two opposite FDFs across each interface but only on the ratios of these values or, more precisely on the ratios of the diffusion fluxes at both sides of the interfaces.
(4)
Hence, we adopt as homogenization parameters the Flux Discontinuity Ratios (FDR) at the interfaces. The second difference is that, instead of directly aiming to preserve region reaction rates, we opt for preserving the averaged normal currents at the interfaces. Preservation of reaction rates follows from the use of flux-weighted cross sections, and this also results in the preservation of region averaged fluxes. We can then formalize the homogenization problem as the minimization of a functional that depends on the sum of the squares of the differences between the reference and the diffusion normal interface currents, where the latter are obtained from the solution of the diffusion equation with the homogenization parameters.
We realize that the functional to be minimized does not include any physical model and that, therefore, the interface ratios are the solution of a purely mathematical problem. Although this is a nonlinear problem, we prove that there is a unique solution which can be easily calculated by considering a diffusion problem separately for every homogenized region, with boundary currents equal to the reference ones. The proof, which is done by construction and does not require any mathematics, is carried out for two typical discretizations of the diffusion equation, finite differences and the nodal transverse method, but it could be extended to other techniques, such as finite elements and analytical nodal methods. However, this solution holds only if the homogenized regions are not submeshed for the diffusion calculation, which would be always the case for pin-by-pin calculation. In the case of submeshing we have used a Jacobian-Free Newton Krylov method
(6)
(JFNK) to minimize the functional and obtain the ratios. For all the cases we have considered we have found that this minimization problem had always a unique exact solution and we have confidence that this is always the case.
Classical equivalence theory is reviewed in Section 2 for use in piecewise homogenization. Is here that we introduce the flux discontinuity ratios and the nonlinear functional that helps determine them. The object of the following section is the calculation of the FDRs and the illustration of these calculations via the exhaustive analysis of a colorset homogenization problem. We focus in particular on the existence of a solution for general subdomain homogenization and prove that an exact solution exists for the cell-centered finite differences as well as for the nodal element spatial discretizations for the diffusion equation with no submeshing. For cases with submeshing we also have found that the JFNK minimization converges to the numerical zero of the functional. Finally, we analyze the behavior of the solution with submeshing and show that for nodal calculations the ratios do not change much with submeshing and, therefore, there is no need to use the time-consuming JFNK search; however, this is not the case for finite differences. In Section 4 we complete our FDR-based piecewise homogenization model by generalizing the classical FDF homogenization method to determine the FDFs for the sides of the regions lying on the boundary of the piecewise homogenized domain. This is done for surface FDFs determined either with the well-known Kord Smith's generalized equivalence theory
(7)
or with the black-box model.
(8)
Conclusions are given in the last section where, in particular, we analyze storage requirements for the FDRs. The equations for the nodal and for the finite-difference discretizations of the diffusion equation with FDR interface conditions, as they have been implemented for the calculations in this work, are detailed in the appendices.
2. CLASSICAL EQUIVALENCE THEORY

Present day reactor calculations are based on the use of homogenized cross sections in three-dimensional, few-group diffusion calculations. The homogenization is obtained from detailed calculations of a series of local transport problems over representative sub domains of the core. For any one of these sub domains one can write a local transport source problem,
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Nint= 2NM– (N+M)

Lager Image

Lager Image

- 2.1 Equivalence

For this method
(9)
,
(10)
the second interface condition is classical flux continuity,
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- 2.2 Flux Discontinuity Ratios

In this case the
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

3. CALCULATION OF FDRS

The solution of problem (2-3,5) depends on the computing mesh (CM) used to solve the diffusion equation. The CM decomposes the domain into homogeneous regions and can be equal to the HM or can be a submesh of the HM. The first case would be typically that of pin-by-pin homogenization, but submeshing can be used to increase the precision of the final core calculation with coarser piecewise homogenizations requiring less storage.
In our calculations we have considered two typical spatial approximations: the mesh-centered Finite-Differences (FD) and the Nodal Expansion Method (NEM). These two-methods belong to the class of methods considered in Section (2.2) and one can directly calculate the FDRs from the averaged boundary fluxes obtained from the individual region solutions. Although the proof has been formally argued in that section, here we give a detailed account for both spatial discretizations. In the case of the FD method the finite-differences expression for the surface currents provides the formula
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Performance of Homogenization Procedure using NEM. p-JFNK = Preconditioned JFNK Iterations. HM = Homogenization Mesh,n= Number of Unknowns, CM = Diffusion Mesh, #GMRES = Total Number of GMRES Iterations,t(s) = Time in Seconds,t(%) = Percentage oftSaved.

Lager Image

FDRs Obtained with the FD and NEM Discretizations for the 2×2 HM. The Interface Notation S, N, W and E is the One inFig. 2

Lager Image

4. FLUX DISCONTINUITY FACTORS

The final core calculation requires homogenized cross sections and FDRs at all interfaces. However, since the FDR homogenization process only provides FDRs at internal interfaces, all the interfaces between sub domains separately homogenized have to be provided with FDRs. These FDRs are not obtained from homogenization and we have to introduce a new boundary homogenization model which is provided by the classical Flux Discontinuity Factor (FDF) theory earlier introduced for whole-assembly homogenization.
(7)
The discontinuity between two interface fluxes is now written as
Lager Image

Lager Image

Lager Image

Lager Image

5. CONCLUSIONS

A strong motivation for piecewise homogenization is its potential for improvement in pin power reconstruction. In this work we have considered piecewise homogenization based on a set of flux discontinuity ratios that preserve interface currents between homogenized regions and therefore region reaction rates. The FDRs are the result of the minimization of a nonlinear functional
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- APPENDIX

- Appendix A: The diffusion Equation with Flux Discontinuity Conditions

We consider a domain composed of Cartesian homogeneous regions. We note
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Appendix B: Diffusion Nodal Equations

Here we adopt the approach in Ref. (11). Integration of the diffusion equation in a node gives the basic nodal relation between the node averaged flux
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Nodal Expansion for the Transverse Flux

The basic nodal approximation consists of introducing an expansion for the transverse nodal fluxes,
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Response Matrix Equations

We are now ready to eliminate all the fluxes for
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Solution by Forward Elimination and Backward Substitution

We write (31) as
Lager Image

Lager Image

- Model for Transverse Leakage

The leakage term (25) can be separate into two different contributions
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Lx,n=Lxy,n+Lxz,n

Lager Image

- Nodal Boundary Conditions

For the nodal solution we have considered geometrical motions and albedo boundary conditions. Geometrical motion (rotation, translation and specular reflection) can be written as interface conditions between the nodes put into contact by the geometrical motion. Rotation and translation boundary condition can put two different assemblies in contact and one should use the general formulation for the interface condition. For the particular case of specular reflection (for which the assemblies in contact are identical and, therefore, the discontinuity ratio is 1) the condition is particularly simple,
Lager Image

- Appendix C: Finite-differences discretization of the diffusion equation

Here we shall use cell instead of node. The basic equations are the cell-integrated balance equations in (13). But, to compute the currents we apply a finite-difference approximation to Fick's law in (15):
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Boundary conditions

Consider one side of a cell on the boundary of the domain and let Φ and
Lager Image

Lager Image

Lager Image

Lager Image

Lager Image

- Final equations

By replacing expressions (36) and (38) in (13) we obtain a system of equations for the cell-averaged fluxes:
Lager Image

Lager Image

Lager Image

- Mba=rabMab,
- rabrba= 1.

Lager Image

Lager Image

- Matrix coefficients and surface-integrated diffusion currents

In our implementation the matrix coefficients are computed using the expression.
Lager Image

Lager Image

Lager Image

Lautard J.J.
,
Loubière S.
,
Fedon-Magnaud C.
1990
‘Three Dimensional Pin by Pin Core Diffusion Calculation’
Proc. IAEA Specialist Meeting on Advanced Calc. Methods for Power Reactors
Cadarache, France
Sept. 10-14

Nichita E.
2010
‘Generalized Equivalence Theory for Checkerboard Configurations in Natural-Uranium CANDU Lattices’
Proc. Int. Mgt. PHYSOR 2010 -Advances in Reactor Physics to Power the Nuclear Renaissance
Pittsburgh, Pennsylvania, USA
May 9-14

Sanchez R.
,
Zmijarevic I.
2011
‘Diffusion piecewise homogenization via discontinuity factors’
Int. Conf. on Mathematics and Computational Methods Applied to Nuclear Science and Engineering
Rio de Janeiro, Brazil
May 8-12

Dante G.
,
Sanchez R.
,
Zmijarevic I.
2012
‘On pin-by-pin discontinuity factors’
Trans. Am. Nucl. Soc.
San Diego, CA, USA
Nov. 11-15
1121 -
1123

Herrero J.J.
,
Garcia-Herranz N.
,
Cuervo D.
,
Ahnert C.
2012
‘Neighborhood-corrected interface discontinuity factors for multi-group pin-by-pin diffusion calculations for LWR’
Annals of Nucl. Eng.
46
106 -
115
** DOI : 10.1016/j.anucene.2012.03.022**

KNOLL D.A.
,
KEYES D.E.
2004
‘Jacobian-free Newton-Krylov Method: a survey of approaches and applications’
J. Comp. Physics
193
357 -
397

Smith K.S.
1986
‘Assembly Homogenization Techniques for Light Water Reactor Analysis’
Progress in Nucl. Energy
17
(3)
303 -
335
** DOI : 10.1016/0149-1970(86)90035-1**

Sanchez R.
2009
“Assembly homogenization techniques for core calculations”
Progress in Nucl. Energy.
51
14 -
31
** DOI : 10.1016/j.pnucene.2008.01.009**

Hébert A.
1993
‘A Consistent Technique for the Pin-by-Pin Homogenization of a Pressurized Water Reactor Assembly’
Nucl. Sci. Eng.
113
227 -
238

Hébert A.
,
Mathonnière G.
1993
‘Development of a Third-Generation Superhomogenization Method for the Homogenization of a Pressurized Water Reactor Assembly’
Nucl. Sci. Eng.
115
129 -

Lawrence R.D.
1986
‘Progress in nodal methods for the solution of the neutron diffusion and transport equations’
Progress in Nuclear Energy
17
(3)
271 -
301
** DOI : 10.1016/0149-1970(86)90034-X**

Dante G.
,
master thesis
2012
‘Piecewise assembly homogenization via flux discontinuity factors’
Politecnico di Torino
master thesis

Halsall M.J.
1990
‘The CACTUS Transport Method in WIMS’
Proc. IAEA Specialist Meeting on Advanced Calc. Methods for Power Reactors
Cadarache, France
Sept. 10-14

Smith K.S.
,
Nuclear Eng. Thesis
1979
‘An analytic nodal method for solving the two-group, multidimensional, static and transient neutron diffusion equations
Massachusetts Institute of Technology
Nuclear Eng. Thesis

Citing 'DIFFUSION PIECEWISE HOMOGENIZATION VIA FLUX DISCONTINUITY RATIOS
'

@article{ OJRHBJ_2013_v45n6_707}
,title={DIFFUSION PIECEWISE HOMOGENIZATION VIA FLUX DISCONTINUITY RATIOS}
,volume={6}
, url={http://dx.doi.org/10.5516/NET.02.2013.518}, DOI={10.5516/NET.02.2013.518}
, number= {6}
, journal={Nuclear Engineering and Technology}
, publisher={Korean Nuclear Society}
, author={SANCHEZ, RICHARD
and
DANTE, GIORGIO
and
ZMIJAREVIC, IGOR}
, year={2013}
, month={Nov}