SAFPWR home page
Fuel rod thermal conduction model in safpwr
Differential form
pic_13 represents the fuel rod, made of the fuel pellet "
u", the gap "
g", the clad "
c" and the coolant of temperature
t.
r denotes the current radius; u(r) is the temperature distribution.
A linear power source f (W/m) is generated in u. The heat rate generated in clad and the thermal capacity c of the gap are neglected.
Let
q(r) be the linear heat rate flowing out of the circular region of radius
r.
(1) is the Fourier heat conduction equation where
k is the heat conductivity (W/m/K).
In terms (2) of the current section area
s, (1) becomes (3).
(4) is the heat balance equation, where
χ(s) is the power distribution in the pellet; uniform distribution would correspond to χ = 1/su.
c is the thermal capacity (J/K/kg).
k and c in the fuel and clad may be temperature-dependent.
Heat transfer coefficient hg (W/mē/K) of gap is imposed.
In the clad-to-water film "f" the heat transfer coefficient may depend on the local thermodynamic properties.
Combining (3) and (4) leads to the heat balance (5).
In steady-state, for uniform χ, (5) accepts a solution such that
k ∂u/∂s = const,
so that by integrating (5) we get (7).
(7) is the "conductivity integral", it allows, by integrating the fuel conductivity curve, to obtain the center temperature as a function of linear power and pellet border temperature.
For uniform χ and k, the u distribution is linear with s, parabolic with r.
Space discretization
The pellet is decomposed in 4 annular regions si of equal area su/4 and uniform thermal properties each.
In steady state, the distribution u(s) satisfies (6) which means that the distribution is strictly polygonal in the pellet.
The clad, owing its small thickness, can be modeled as a plate; hence, its temperature distribution is linear versus r (no heat source in clad).
Stationary solution
The steady-state solution is obtained by completing the balance equation (6) with the boundary conditions (zero temperature gradient at center, automatically satisfied with a linear distribution u(s), and u6=t, the current water temperature) and the u and q continuity conditions at each region interface.
This finally results in a tridiagonal symmetric linear system
(8:15) whose resolution is straightforward.
As the elements of
K are temperature-dependent, the resolution must, however, be iterated until convergence, controlled by the criterum (?)
The steady state solution is exact.
Time discretization
In transient condition, the
u solution is no longer strictly polygonal because the derivative
∂u/∂s in (5) is not longer uniform over each region.
However, if the curvature of
u(s) within each region is not too pronounced, it makes sense to keep calculating the
kij from the stationary relations
(13 to 15).
The linear system to be solved is (9).
It is time-discretized as (18), which means that the thermal properties in C and K are taken at bos "1" but, to insure stability, the heat leakage (K1 u2) is calculated from the eos u temperature distribution u2, by means of
qij = kij (u2i - u2j).
In order to overcome numerical precision problems arising when clad and water temperatures are very close, the problem (21) is solved in terms of the temperature
deviation u2 - t over the currently known water temperature
t (bos
t1).
This is formally obtained by setting
t to 0 in (10).
The linear heat rate qlu out of the fuel rod is then simply (17).
The stepping transition statements for u and t2 must be adapted accordingly as (22, 23).
In the course of Doppler feed-back iterations, u2 must be corrected to account for heat source f variation over the step, which may be significant in the case, for example, of a control rod ejection transient.
Instead of repeating the whole problem (21) at each iteration, we evaluate the u2 correction by means of the derivative (24), which imply that variations of coolant temperature t and the thermal properties in the course of the feedback iterations are neglected.
The accuracy of the solution scheme has been tested against high order, finite elements solution.