3.5. Transient Heat Transfer after the Start of Boiling

After the start of boiling, the coolant temperatures are calculated in the coolant routines, rather than being calculated simultaneously with fuel, cladding, and structure temperatures. Coupling between the boiling calculations and the non-coolant heat-transfer calculation takes place in two parts for each time step. First, the boiling routines us extrapolated cladding and structure temperatures to calculate the heat fluxes to the coolant for the boiling calculation. Then the heat fluxes actually used in the coolant routines are passed to the heat-transfer routines to be used as boundary conditions at the cladding, structure, and reflector surfaces. The net results of this procedure are that energy is conserved, and a fully implicit boiling calculation can be made without requiring a direct simultaneous solution of all of the fuel-pin temperatures in the boiling model. The coupling through extrapolated cladding and structure temperatures and heat fluxes at the cladding and structure surfaces imposes numerical stability limitations on the heat-transfer time-step sizes. Currently, fuel-pin temperatures are calculated at the end of the fuel-pin heat-transfer time step, whereas structure and reflector temperatures are calculated at every coolant time step. The coolant time step can be no longer than the heat-transfer step, and the coolant step is often much shorter. For typical fuel pins, the stability limit for the heat-transfer time step is of the order of .02 s. With a thin structure, the stability limit for structure temperature calculations could be less than .02 s, although for typical duct wall thicknesses (.12 in. or .003 m) the stability limit would be closer to one second. Should timing studies indicate that the structure and reflector temperature calculations account for a significant fraction of the total computing time, then the code will be modified so as to do these calculations less often than once every coolant time step.

3.5.1. Fuel and Cladding Temperatures in the Core and Axial Blanket

The equations used for fuel and cladding temperatures after the switch to the boiling module are the same as those used in the non-voiding module, except that in the boiling module the fuel-pin heat-transfer calculations stop at the cladding outer surface rather than carrying through to the structure outer node. The finite difference equations for radial nodes 1-NE are the same as Eq. (3.3-17) to Eq. (3.3-37). For node $\text{NE}'$ , the cladding outer node, the heat flux to the coolant at the cladding surface must be accounted for. Also, in a boiling region, a film of liquid sodium can be left on the cladding. Since the film is in intimate contact with the cladding, the heat capacity of the film is added to the heat capacity of the cladding outer node, rather than being accounted for in the boiling calculation. Thus, the finite difference equation for node $\text{NE}'$ becomes

$$\begin{multline} \left\lbrack \frac{m_{\text{e}}c_{\text{e}}}{4} + 2\pi r\left( \text{NE}' \right)\rho_{\text{c}}c_{\text{c}}w_{\text{fe}} \Delta z \left( j \right) \right\rbrack \cdot \left\lbrack \frac{T_{2}\left( \text{NE}' \right) - T_{1}\left( \text{NE}' \right)}{ \Delta t } \right\rbrack \\ = - \frac{2\pi{\overline{r}}_{\text{NE}'} \Delta z \left( j \right)\ {\overline{k}}_{\text{NE,NE}'}}{\Delta r_{\text{NE,NE}'}} \left\{ \theta_{21} \left\lbrack T_{1}\ \left( \text{NE}' \right) - T_{1} \left( \text{NE} \right) \right\rbrack \right. \\ \left. + \theta_{2}\ \left\lbrack T_{2}\ \left( \text{NE}' \right) - T_{2} \left( \text{NE} \right) \right\rbrack\ \right\} - 2\pi r\left( \text{NE}' \right) \Delta z \left( j \right) \frac{E_{\text{ec}}\left( j \right)}{ \Delta t } \\ + Q\left( \text{NE}' \right) \end{multline}$$

where $w_{\text{fe}}$ is the thickness of liquid sodium film left on the cladding after voiding occurs, and $E_{\text{ec}}$ is the integrated heat flux from cladding to coolant.

The value of $E_{\text{ec}}$ is computed in the coolant routines as

(3.5-2) $$E_{\text{ec}}\left( jc \right) = \int_{t}^{t + \Delta t}{ \frac{T_{\text{eex}}\left( jc \right) - {\overline{T}}_{c}\left( jc \right)}{R_{\text{ec}}\left( jc \right)}\ \text{dt}'}$$

where

$T_{\text{eex}}$ = extrapolated cladding temperature at a point ¼ of the way from the outer cladding surface to the inner cladding surface:

$$\begin{multline} T_{\text{eex}}\left( j,t' \right) = f_{1} \frac{\left\lbrack \text{T}\left( \text{NE},j,t_{1} \right) + T \left( \text{NEP},j,t_{1} \right) \right\rbrack}{2} \\ + f_{2} \frac{\left\lbrack \text{T}\left( \text{NE},j,t_{2} \right) + T \left( \text{NEP},j,t_{2} \right) \right\rbrack}{2} \end{multline}$$

(3.5-4) $$f_{1} = \frac{t_{2} - t'}{t_{2} - t_{1}}$$

(3.5-5) $$f_{1} = 1 - f_{1}$$

(3.5-6) $$R_{\text{ec}} = \frac{1}{h_{\text{c}}} + R_{\text{ehf}}$$

and

(3.5-7) $$R_{\text{ehf}} = \frac{\left\lbrack r \left( \text{NE}' \right) - r\left( \text{NE} \right) \right\rbrack}{k_{\text{e}}}\left\{ \frac{r \left( \text{NE}' \right) \left\lbrack 1 - \frac{\gamma_{\text{e}}}{4\left( 1 - \gamma_{\text{c}} - \gamma_{\text{s}} \right)} \right\rbrack}{r \left( \text{NE} \right) + 3r \left( \text{NE}' \right)} \right\}$$

If a node is partly voided, then $E_{\text{ec}}$ is averaged over the length of the node as well as integrated over time.

The $\rho_{\text{c}}c_{\text{c}}w_{\text{fe}}$ term is supplied by the boiling routines. It represents the heat capacity of any liquid film left on the cladding after voiding occurs. This term is zero unless voiding has occurred at this axial node. Since the film temperature tends to follow the cladding surface temperature more closely than it follows the vapor temperature, the film heat capacity is lumped with the cladding outer node, and film temperatures are not explicitly calculated in the boiling routines.

The finite difference equations are again put in a matrix form like Eq. (3.3-49), except that in the voiding case there are only $\text{NE}'$ elements. The definitions of $\alpha$, $\beta$, and $D$ are the same as in Eq. (3.3-50) to 3.3-51 for nodes 1-$\text{NE}$. For node $\text{NE}'$, $\alpha$ is still given by Eq.

(3.5-8) $$\beta_{\text{NE}'} = 0$$

Eq. (3.3-53) or Eq. (3.3-51), except that the $\rho_{\text{c}}c_{\text{c}}w_{\text{fe}}$ term is added to it. Also,

and

(3.5-9) $$D_{\text{NE}'} = \frac{\theta_{1}}{\theta_{2}} \beta_{\text{NE}} T_{1}\left( \text{NE} \right) + T_{1} \left( \text{NE}' \right)\left\lbrack \alpha_{\text{NE}'} - \frac{\theta_{1}}{\theta_{2}} \beta_{\text{NE}} \right\rbrack + \psi_{\text{NE}'} - r\left( \text{NE}' \right) E_{\text{ec}} \left( j \right)$$

The equations are solved in the same manner as in the non-voiding case.

3.5.2. Structure Temperatures

The basic equations used for the two structure radial nodes at an axial node are

(3.5-10) $$\left( \rho c \right)_{\text{sto}} d_{\text{sto}} \frac{\text{dT}_{\text{sto}}}{\text{dt}} = H_{\text{stio}} \left( T_{\text{sti}} - T_{\text{sto}} \right) + Q_{\text{st}} f_{\text{o}}$$

and

(3.5-11) $$\left\lbrack \left( \rho c \right)_{\text{sti}} d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack \frac{\text{dT}_{\text{sti}}}{\text{dt}} = H_{\text{stio}} \left( T_{\text{sto}} - T_{\text{sti}} \right) + Q_{\text{st}}\ f_{\text{i}} + \frac{\left( T_{\text{c}} - T_{\text{sti}} \right)}{R_{\text{sc}}}$$

where

(3.5-12) $$f_{\text{i}} = \frac{d_{\text{sti}}}{d_{\text{sti}} + d_{\text{sto}}}$$

(3.5-13) $$f_{\text{o}} = \frac{d_{\text{sto}}}{d_{\text{sti}} + d_{\text{sto}}}$$

and $w_{\text{fst}}$ = thickness of the liquid-sodium film left on the structure after voiding occurs.

The heat capacity of the film in the boiling region, as represented by the $\rho_{\text{c}}c_{\text{c}}w_{\text{fst}}$ term, is supplied by the boiling routines to be added to the inner structure node.

The structure temperature calculation is either a semi-implicit or a fully implicit calculation, depending on the time-step in relation to an inner structure node heat-transfer time constant, $\tau_{\text{sti}}$, calculated as

(3.5-14) $$\tau_{\text{sti}} = \frac{\left( \rho c \right)_{\text{sti}}d_{\text{sti}}^{2}}{2k_{\text{sti}}}$$

If $\Delta t$ is less than $\tau_{\text{sti}}$ then the semi-implicit calculation is used. Otherwise the fully implicit calculation is used.

3.5.2.1. Semi-Implicit Calculations

For the semi-implicit calculation, finite differencing of Eq. (3.5-10) and Eq. (3.5-11) gives

(3.5-15) $$\left( \rho c \right)_{\text{sto}} d_{\text{sto}} \frac{\left( T_{\text{sto}2} - T_{\text{sto}1} \right)}{\Delta t} = \frac{H_{\text{stio}}}{2} \left\lbrack T_{\text{sti}2} - T_{\text{sto}2} + T_{\text{sti}1} - T_{\text{sto}1} \right\rbrack + Q_{\text{st}} f_{\text{o}}$$

and

$$\begin{multline} \left\lbrack \left( \rho c \right)_{\text{sti}} d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack \frac{\left( T_{\text{sti}2} - T_{\text{sti}1} \right)}{\Delta t} \\ = \frac{H_{\text{stio}}}{2}\left\lbrack T_{\text{sto}2} - T_{\text{sti}2} + T_{\text{sto}1} - T_{\text{sti}1} \right\rbrack + Q_{\text{st}}\ f_{\text{i}} - \frac{E_{\text{sc}}}{\Delta t} \end{multline}$$

where

(3.5-17) $$E_{\text{sc}}\left( jc \right) = \int_{\text{t}}^{t + \Delta t}{ \frac{T_{\text{stex}}\left( jc \right)\ - {\overline{T}}_{\text{c}}\left( jc \right)}{R_{\text{sc}}\left( jc \right)}}\ \text{dt}'$$

The value of $E_{\text{sc}}$ is computed in the boiling routines. $T_{\text{stex}}$ is the extrapolated structure inner node temperature, and

(3.5-18) $$R_{\text{sc}} = \frac{1}{h_{\text{c}}} + \frac{d_{\text{sti}}}{2k_{\text{sti}}}$$

Eq. (3.5-15) and Eq. (3.5-16) are put in the form

(3.5-19) $$a_{11}T_{\text{sto}2} + a_{12}T_{\text{sti}2} = b_{1}$$

(3.5-20) $$a_{21}T_{\text{sto}2} + a_{22}T_{\text{sti}2} = b_{2}$$

with the solutions

(3.5-21) $$T_{\text{sto}2} = \frac{a_{22}b_{1} - a_{12}b_{2}}{a_{11}a_{22} - a_{12}a_{21}}$$

and

(3.5-22) $$T_{\text{sti}2} = \frac{a_{11}b_{2} - a_{21}b_{2}}{a_{11}a_{22} - a_{12}a_{21}}$$

where

(3.5-23) $$a_{11} = \left( \rho c \right)_{\text{sto}}\ d_{\text{sto}} + \frac{ \Delta t }{2} H_{\text{stio}}$$

(3.5-24) $$a_{12} = - \frac{ \Delta t }{2}\ H_{\text{stio}}$$

(3.5-25) $$a_{21} = a_{12}$$

(3.5-26) $$a_{22} = \left( \rho c \right)_{\text{sti}} d_{\text{sti}} + \frac{ \Delta t }{2} H_{\text{stio}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}}$$

(3.5-27) $$b_{1} = \left\lbrack \left( \rho c \right)_{\text{sto}}d_{\text{sto}} - \frac{ \Delta t }{2} H_{\text{stio}} \right\rbrack T_{\text{sto}1} + Q_{\text{st}}f_{0}\Delta t - a_{12} T_{\text{sti}1}$$

and

(3.5-28) $$b_{2} = \left\lbrack \left( \rho c \right)_{\text{sti}}d_{\text{stio}} - \frac{ \Delta t }{2}\ H_{\text{stio}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack T_{\text{sti}1} - a_{12}\ T_{\text{sto}1} + Q_{\text{st}}f_{\text{i}}\Delta t - \Delta t Q_{\text{sc}}$$

3.5.2.2. Fully Implicit Calculations

Since the inner structure node may represent only a small fraction of the total structure thickness, $\tau_{\text{sti}}$ can be small. If the time-step size is appreciably larger than $\tau_{\text{sti}}$ then the semi-implicit calculation can become numerically unstable. Therefore, a different algorithm is used for larger time-step sizes. This algorithm uses two steps. First, a fully implicit calculation is made, using a coolant temperature and thermal resistance to the coolant as structure surface boundary conditions, rather than using the integrated heat flux. In this first step, the heat flux from the coolant to the structure will, in general, not match the heat flux form structure to coolant used in the coolant calculations. Therefore, in the second step, the inner node and outer node structure temperatures are both adjusted by the same amount so that the integrated heat flux from structure to coolant is matched.

For the first step the finite difference equations used for the two structure node temperatures are

(3.5-29) $$\left( \rho c \right)_{\text{sto}} d_{\text{sto}} \frac{\left( T_{\text{sto}2} - T_{\text{sto}1} \right)}{ \Delta t } = H_{\text{stio}} \left( T_{\text{sti}2} - T_{\text{sto}2} \right) + Q_{\text{st}} f_{o}$$

and

(3.5-30) $$\left\lbrack \left( \rho c \right)_{\text{sti}} d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack \frac{\left( T_{\text{sti}2} - T_{\text{sti}1} \right)}{ \Delta t } = H_{\text{stio}}\left( T_{\text{sto}2} - T_{\text{sti}2} \right) + \frac{T_{\text{c}} - T_{\text{sti}2}}{R_{\text{sc}}} + Q_{\text{st}} f_{\text{i}}$$

where the values for $T_{\text{c}}$ and $R_{\text{sc}}$ are supplied by the coolant routines. The solutions for these two equations again have the same form as Eq. (3.5-21) and Eq. (3.5-22) except that in this case the coefficients are defined as

(3.5-31) $$a_{11} = \left( \rho c \right)_{\text{sto}} d_{\text{sto}} + \Delta t H_{\text{stio}}$$

(3.5-32) $$a_{12} = a_{21} = - \Delta t H_{\text{stio}}$$

and

(3.5-33) $$a_{22} = \left( \rho c \right)_{\text{sti}} d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} + \Delta t H_{\text{stio}} + \frac{ \Delta t }{R_{\text{sc}}}$$

The temperature difference, $\Delta T_{\text{st}}$, between the outer and inner nodes is then defined as

(3.5-34) $$\Delta T_{\text{st}} = T_{\text{sto}2} - T_{\text{sti}2}$$

In the second step, $\Delta T_{\text{st}}$ is preserved but the temperatures are adjusted so as to match the value of $E_{\text{sc}}$ supplied by the coolant routines:

(3.5-35) $$\left\lbrack \left( \rho c \right)_{\text{sti}}d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack\left( T_{\text{sci}2} - T_{\text{sti}1} \right) + \left( \rho c \right)_{\text{sto}}d_{\text{sto}}\left( T_{\text{sto}2} - T_{\text{sto}1} \right) = Q_{\text{st}} - E_{\text{sc}}$$

The solution to Eq. (3.5-34) and Eq. (3.5-35) is

$$\begin{multline} T_{\text{sto}2} = \left\{ Q_{\text{st}} \Delta t - E_{\text{sc}} + \left( \rho c \right)_{\text{sto}} d_{\text{sto}}T_{\text{sto}1} \right. \\ \left. + \left\lbrack \left( \rho c \right)_{\text{sti}}\ d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack \cdot \left( T_{\text{sti}1} + \Delta T_{\text{st}} \right) \right\} \left\lbrack \left( \rho c \right)_{\text{sto}}\ d_{\text{sto}} + \left( \rho c \right)_{\text{sti}}d_{\text{sti}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fst}} \right\rbrack \end{multline}$$

and

(3.5-37) $$T_{\text{sti}2} = T_{\text{sto}2} - \Delta T_{\text{st}}$$

Note that the second step can still cause numerical instabilities if the time-step size is too large or the total structure thickness is too small, but in the fully implicit scheme the stability limit is based on the total structure thickness, whereas in the semi-implicit scheme the stability is based mainly on the inner node thickness.

3.5.3. Reflector Temperatures

In the boiling module the treatment of reflector temperatures is almost identical to the structure temperature treatment. The main difference is that in the reflector, the outer node is in contact with the coolant, whereas in the structure, the inner code is in contact with the coolant. Also, in the reflector the density, specific heat, and thermal conductivity are the same for both nodes, whereas in the structure these properties can vary from inner node to outer node.

The basic equations are

(3.5-38) $$\left( \rho c \right)_{\text{r}} d_{\text{ri}} \frac{\text{dT}_{\text{ri}}}{\text{dt}} = H_{\text{rio}} \left( T_{\text{ro}} - T_{\text{ri}} \right)$$

and

(3.5-39) $$\left\lbrack \left( \rho c \right)_{\text{r}} d_{\text{ro}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack \frac{\text{dT}_{\text{ro}}}{\text{dt}} = H_{\text{rio}}\left( T_{\text{ri}} - T_{\text{ro}} \right) + \left( T_{\text{c}} - T_{\text{ro}} \right) H_{\text{erc}}$$

where

$w_{\text{fr}}$ = thickness of the liquid-sodium film left on the reflector after voiding occurs.

An outer reflector node heat-transfer time constant, $\tau_{\text{ro}}$, is calculated as

(3.5-40) $$\tau_{\text{ro}} = \frac{\left( \rho c \right)_{\text{r}} d_{\text{ro}}^{2}}{2k_{\text{r}}}$$

and a fully implicit calculation is used if the time-step size is greater than $\tau_{\text{ro}}$. Otherwise a semi-implicit calculation is used.

3.5.3.1. Semi-Implicit Calculations

Finite differencing of Eq. (3.5-38) and Eq. (3.5-39) gives

(3.5-41) $$\left( \rho c \right)_{\text{r}} d_{\text{ri}} \frac{\left( T_{\text{ri}2} - T_{\text{ri}1} \right)}{ \Delta t } = \frac{H_{\text{rio}}}{2} \left( T_{\text{ro}2} - T_{\text{ri}2} + T_{\text{ro}1} - T_{\text{ri}1} \right)$$

(3.5-42) $$\left\lbrack \left( \rho c \right)_{\text{r}} d_{\text{ro}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack \frac{\left( T_{\text{ri}2} - T_{\text{ri}1} \right)}{ \Delta t } = \frac{H_{\text{rio}}}{2}\left( T_{\text{ri}2} - T_{\text{ro}2} + \ T_{\text{ri}1} - T_{\text{ro}1} \right) - \frac{E_{\text{rc}}}{ \Delta t }$$

where

(3.5-43) $$E_{\text{rc}} = \int_{t}^{t + \Delta t}{ \left( T_{\text{rex}} - {\overline{T}}_{\text{c}} \right)} H_{\text{erc}} \text{dt}'$$

These equations are put in the form

(3.5-44) $$a_{11}T_{\text{ri}2} + \ a_{12}T_{\text{ro}2} = b_{1}$$

(3.5-45) $$a_{21}T_{\text{ri}2} + \ a_{22}T_{\text{ro}2} = b_{2}$$

with the solution

(3.5-46) $$T_{\text{ri}2} = \frac{a_{22}b_{1} - a_{12}b_{2}}{a_{11}a_{22} - a_{12}a_{21}}$$

and

(3.5-47) $$T_{\text{ro}2} = \frac{a_{11}b_{2} - a_{21}b_{1}}{a_{11}a_{22} - a_{12}a_{21}}$$

The coefficients are

(3.5-48) $$a_{11} = \left( \rho c \right)_{\text{r}} d_{\text{ri}} + \frac{ \Delta t }{2}\ H_{\text{rio}}$$

(3.5-49) $$a_{12} = - \frac{ \Delta t }{2} H_{\text{rio}}$$

(3.5-50) $$a_{21} = a_{12}$$

(3.5-51) $$a_{22} = \left( \rho c \right)_{\text{r}} d_{\text{ro}} + \frac{ \Delta t }{2} H_{\text{rio}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}}$$

(3.5-52) $$b_{1} = \left\lbrack \left( \rho c \right)_{\text{r}}d_{\text{ri}} - \frac{ \Delta t }{2} H_{\text{rio}} \right\rbrack T_{\text{ri}1} + \frac{ \Delta t }{2}H_{\text{rio}} T_{\text{ro}1}$$

and

(3.5-53) $$b_{2} = \left\lbrack \left( \rho c \right)_{\text{r}}d_{\text{ro}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}} - \frac{ \Delta t }{2}\ H_{\text{rio}} \right\rbrack T_{\text{ro}1} + \frac{ \Delta t }{2} H_{\text{rio}} T_{\text{ri}1} - E_{\text{rc}}$$

3.5.3.2. Fully Implicit Calculations

As in the structure temperature case, a two-step process is used. In the first step, the finite difference equations used are

(3.5-54) $$\left( \rho c \right)_{\text{r}} d_{\text{ri}} \frac{\left( T_{\text{ri}2} - T_{\text{ri}1} \right)}{ \Delta t } = H_{\text{rio}} \left( T_{\text{ro}2} - T_{\text{ri}2} \right)$$

and

(3.5-55) $$\left\lbrack \left( \rho c \right)_{\text{r}} d_{\text{ro}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack \frac{\left( T_{\text{ro}2} - \ T_{\text{ro}1} \right)}{ \Delta t } = H_{\text{rio}}\left( T_{\text{ri}2} - T_{\text{ro}2} \right) + \left( {\overline{T}}_{\text{c}} - T_{\text{ro}2} \right) H_{\text{erc}}$$

The solutions again have the same form as Eq. (3.5-46) and Eq. (3.5-47), with the coefficients given by

(3.5-56) $$a_{11} = \left( \rho c \right)_{\text{r}} d_{\text{ri}} + \Delta t H_{\text{rio}}$$

(3.5-57) $$a_{12} = - \Delta t H_{\text{rio}}$$

(3.5-58) $$a_{21} = a_{12}$$

(3.5-59) $$a_{22} = \left( \rho c \right)_{\text{r}} d_{\text{ro}} + \Delta t H_{\text{rio}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}}$$

(3.5-60) $$b_{1} = \left( \rho c \right)_{\text{r}}d_{\text{ri}} T_{\text{ri}1}$$

and

(3.5-61) $$b_{2} = \left\lbrack \left( \rho c \right)_{\text{r}}d_{\text{ro}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack{T}_{\text{rol}} + {\overline{T}}_{\text{c}} \Delta t \ H_{\text{erc}}$$

The temperature difference between nodes, $\Delta T_{\text{r}}$, is defined as

(3.5-62) $$\Delta T_{\text{r}} = T_{\text{ri}2} - T_{\text{ro}2}$$

In the second step, $\Delta T_{\text{r}}$ is preserved and $E_{\text{rc}}$ is matched. The energy conservation equation is

(3.5-63) $$\left( \rho c \right)_{\text{r}}d_{\text{ri}}\left( T_{\text{ri}2} - T_{\text{ri}1} \right) + \left\lbrack \left( \rho c \right)_{\text{r}}d_{\text{ro}} + \rho_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack\left( T_{\text{ro}2} - T_{\text{ro}1} \right) = - E_{\text{rc}}$$

The solution to Eq. (3.5-62) and Eq. (3.5-63) is

(3.5-64) $$T_{\text{ri}2} = \frac{\left\{ - E_{\text{rc}} + \left( \rho c \right)_{\text{r}}d_{\text{ri}}T_{\text{ri}1} + \left\lbrack \left( \rho c \right)_{r}d_{\text{ro}} + p_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack \cdot (T_{\text{rol}} + \Delta T_{\text{r}}) \right\}}{\left\lbrack \left( \rho c \right)_{\text{r}}\left( d_{\text{ri}} + d_{\text{ro}} \right) + p_{\text{c}}c_{\text{c}}w_{\text{fr}} \right\rbrack}$$

and

(3.5-65) $$T_{\text{ro}2} = T_{\text{ri}2} - \Delta T_{\text{r}}$$

3.5.4. Gas Plenum Region

The basic equations used for the cladding and gas temperatures in the gas plenum region are Eq. (3.5-68) and the following equation:

(3.5-66) $$\rho_{\text{e}}c_{\text{e}}A_{\text{ep}} \frac{\text{dT}_{\text{e}}\left( jp \right)}{\text{dt}} = 2\pi r_{\text{erp}}H_{\text{erc}}\left\lbrack {\overline{T}}_{\text{c}}\left( jc \right) - T_{\text{e}}\left( jp \right) \right\rbrack + 2\pi r_{\text{brp}}H_{\text{eg}} \left\lbrack T_{\text{g}} - T_{\text{e}}\left( jp \right) \right\rbrack$$

Since Eq. (3.5-68) links all of the cladding nodes in the gas plenum, a direct semi-implicit or implicit solution of Eq. (3.5-68) and Eq. (3.5-66) would require a simultaneous solution for the gas temperature and all of the cladding node temperatures. Instead, the cladding temperatures are calculated first, using the gas temperature at the beginning of the time step. Then the gas temperature is calculated using the newly computed cladding temperatures.

Finite differencing of Eq. (3.5-66) gives

(3.5-67) $$\rho_{\text{e}}c_{\text{e}}A_{\text{ep}} \frac{\left( T_{\text{e}2} - T_{\text{e}1} \right)}{ \Delta t } = - 2\pi r_{\text{erp}} \frac{E_{\text{ec}}}{ \Delta t } + \pi r_{\text{brp}}H_{\text{eg}}\left( 2T_{\text{g}1} - T_{e2} - T_{\text{e}1} \right)$$

where $E_{\text{ec}}$ is calculated in the coolant routines in the same manner as indicated in Eq. (3.5-2), except that in the gas plenum only one radial node is used in the cladding, and $R_{\text{ehf}}$ becomes

(3.5-68) $$R_{\text{ehf}} = \frac{r_{\text{erp}} - r_{\text{brp}}}{2k_{\text{e}}}$$

The solution of Eq. (3.5-67) for $T_{\text{e}2}$ gives

(3.5-69) $$T_{\text{e}2} = \frac{\left\{ \left( \rho_{\text{e}}c_{\text{e}}A_{\text{ep}} - \pi r_{\text{brp}}H_{\text{eg}} \Delta t \right)T_{\text{el}} - 2\pi r_{\text{erp}}E_{\text{cc}} + 2\pi r_{\text{brp}}H_{\text{eg}} \Delta t T_{\text{g}1} \right\}}{\left( \rho_{\text{e}}c_{\text{e}}A_{\text{ep}} + \pi r_{\text{brp}}H_{\text{eg}} \Delta t \right)}$$

In the second step, Eq. (3.5-70) is used with $\theta_{1} = \theta_{2} = 1/2$. The solution for $T_{\text{g}2}$ is

(3.5-70) $$T_{\text{g}2} = \frac{\left\lbrack \left( \rho c \right)_{\text{g}} A_{\text{g}} - \pi r_{\text{brp}} H_{\text{eg}}\Delta t \right\rbrack T_{\text{g}1} + \pi r_{\text{brp}} H_{\text{eg}}\Delta t \frac{s_{1}}{s_{2}}}{\left( \rho c \right)_{\text{g}} A_{\text{g}} + \pi r_{\text{brp}} H_{\text{eg}}\Delta t}$$

where

(3.5-71) $$s_{1} = \sum_{\text{jp}}{ \left\lbrack T_{\text{e}1} \left( jp \right) + T_{\text{e}2} \left( jp \right) \right\rbrack} \Delta z \left( jp \right)$$

and

(3.5-72) $$s_{2} = \sum_{\text{jp}}{ \Delta z \left( jp \right)}$$

3.5.5. Coolant Temperatures in Liquid Slugs

Before the onset of coolant voiding, coolant temperatures are calculated at all node boundaries, as indicated in Figure 3.2.3. After the start of boiling, liquid coolant temperatures are calculated at all node boundaries outside vapor regions, as well as at moving nodes near the bubble interfaces. Two different types of calculations are made. Eulerian temperature calculations are made for fixed coolant nodes in the inlet and outlet liquid slugs. Lagrangian temperature calculations are made for the moving interface nodes and for any fixed nodes in liquid slugs between bubbles. There is also an option to use Lagrangian temperature calculations for all nodes, both fixed and moving.

The Eulerian calculation is probably more accurate for the fixed nodes. The main disadvantage of this method is that a sudden jump in inlet temperature can lead to a sawtooth temperature pattern, with the temperature high at one node, low at the next, and high again at the third node. The Lagrangian calculation does not exhibit this behavior. This sawtooth behavior is not unstable: the perturbation at any node is no larger than the jump in the inlet temperature, and the perturbations tend to die out in later time steps. Also, the coolant inlet and reentry temperature calculations described in Section 3.3.6 tend to eliminate sudden jumps in inlet and reentry temperatures.

3.5.5.1. Eulerian Temperature Calculation

The basic equation used in this calculation is again Eq. (3.3-5). The heat fluxes $Q_{\text{ec}}$ and $Q_{\text{sc}}$ are calculated as

(3.5-73) $$Q_{\text{ec}} = \frac{\left( T_{\text{e}} - {\overline{T}}_{\text{c}} \right)}{R_{\text{ec}}} \frac{2\pi r\left( \text{NE}' \right)}{A_{\text{c}}}$$

and

(3.5-74) $$Q_{\text{sc}} = \frac{\left( T_{\text{si}} - {\overline{T}}_{\text{c}} \right)}{R_{\text{sc}}} \frac{S_{\text{st}}}{A_{\text{c}}}$$

where $R_{\text{ec}}$ and $R_{\text{sc}}$ are given by Eq. (3.5-6) and Eq. (3.5-18), $T_{\text{e}}$ is the average of $T(\text{NE})$ and $T(\text{NE}')$, and $T_{\text{si}}$ is the inner structure node temperature. In reflector zones, $T_{\text{e}}$ is replaced by the reflector outer node temperature; and in the gas plenum region, the one radial cladding node temperature is used. In the boiling module, the coolant temperatures are calculated before the cladding and structure temperatures are, so linear extrapolation in time is used to obtain values of $T_{\text{e}}$ and $T_{\text{si}}$ at the end of a time step.

A semi-implicit finite differencing of Eq. (3.3-5) gives

$$\begin{multline} \overline{p}\left( jc \right) {\overline{c}}_{\text{c}}\left( jc \right)A_{\text{c}}\left( jc \right)\frac{\left\lbrack T_{\text{c}2}\left( jc + 1 \right) + T_{\text{c}2}\left( jc \right) - T_{\text{c}1}\left( jc + 1 \right) + T_{\text{c}1}\left( jc \right) \right\rbrack}{2\Delta t} \\ + {\overline{c}}_{\text{c}}\left( jc \right)w_{1}\frac{\left\lbrack T_{\text{c}1}\left( jc + 1 \right) + T_{\text{c}1}\left( jc \right) \right\rbrack}{2\Delta z\left( jc \right)} + {\overline{c}}_{\text{c}}\left( jc \right)w_{2}\frac{\left\lbrack T_{\text{c}2}\left( jc + 1 \right) - T_{\text{c}2}\left( jc \right) \right\rbrack}{2\Delta z\left( jc \right)} \\ = Q_{\text{c}}\left( jc \right)A_{\text{c}}\left( jc \right) + \frac{k_{5}\left( jc \right) A_{\text{c}}\left( jc \right)}{4}\left\{ \frac{2T_{\text{e}2}\left( jc \right) - T_{\text{e}2}\left( jc \right) - T_{\text{c}2}\left( jc + 1 \right)}{R_{\text{ec}2}\left( jc \right)} \right. \\ \left. + \frac{2T_{\text{e}1}\left( jc \right) - T_{\text{e}1}\left( jc \right) - T_{\text{c}1}\left( jc + 1 \right)}{R_{\text{ec}1}\left( jc \right)} + \gamma_{2}\left( jc \right)\left\lbrack \frac{2T_{\text{si}2}\left( jc \right) - T_{\text{c}2}\left( jc \right) - T_{\text{c}2}\left( jc + 1 \right)}{R_{\text{sc}2}\left( jc \right)} \right. \right. \\ \left. \left. + \frac{2T_{\text{si}1}\left( jc \right) - T_{\text{c}1}\left( jc \right) - T_{\text{c}1}\left( jc + 1 \right)}{R_{\text{sc}1}\left( jc \right)} \right\rbrack \right\} \end{multline}$$

where

(3.5-76) $$\begin{split}k_{5}\left( jc \right) = \begin{cases} \frac{2\pi r\left( \text{NE}',jc \right)}{A_{\text{c}}\left( jc \right)} & \text{in the core and blankets} \\ \frac{s_{\text{er}}\left( \text{kz} \right)}{A_{\text{c}}\left( jc \right)} & \text{in a reflector region} \\ \frac{2\pi r_{\text{erp}}}{A_{\text{c}}\left( jc \right)} & \text{in the gas plenum region} \end{cases}\end{split}$$

and

(3.5-77) $$\gamma_{2}\left( jc \right) = \frac{S_{\text{st}}\left( jc \right)}{k_{5}\left( jc \right) A_{\text{c}}\left( jc \right)}$$

Solving for $T_{\text{c}2}\left( jc + 1 \right)$ gives

$$\begin{multline} T_{\text{c}2}\left( jc + 1 \right) = \left\{ T_{\text{c}1}\left( jc + 1 \right) \left\lbrack {\overline{\rho}}_{\text{c}}\left( jc \right) - \frac{ \Delta t w_{1}}{ \Delta z \left( jc \right)A_{\text{c}}\left( jc \right)} - \frac{k_{5}\left( jc \right) \Delta t h_{\text{br}}\left( jc \right)}{2{\overline{c}}_{\text{c}}\left( jc \right)} \right\rbrack \right. \\ \left. + T_{\text{c}2}\left( jc \right) \left\lbrack - \rho\left( jc \right) + \frac{ \Delta t w_{2}}{ \Delta z \left( jc \right)A_{\text{c}}\left( jc \right)} - \frac{k_{5}\left( jc \right) \Delta t }{2{\overline{c}}_{\text{c}}\left( jc \right)} h_{\text{b}2}\left( jc \right) \right\rbrack \right. \\ \left. + T_{\text{c}1}\left( jc \right) \left\lbrack \rho_{\text{c}}\left( jc \right) + \frac{ \Delta t w_{1}}{ \Delta z \left( jc \right)A_{\text{c}}\left( jc \right)} - \frac{k_{5}\left( jc \right) \Delta t }{2{\overline{c}}_{\text{c}}\left( jc \right)} h_{\text{b}1}\left( jc \right) \right\rbrack \right. \\ \left. + \frac{4\Delta t}{2{\overline{c}}_{\text{c}}\left( jc \right)}\ k_{5}\left( jc \right) \varphi_{1}\left( jc \right) + \frac{2\Delta t Q_{\text{c}}\left( jc \right)}{{\overline{c}}_{\text{c}}\left( jc \right)} \right\}\ \big/\ \left\{ {\overline{\rho}}_{\text{c}}\left( jc \right) \right. \\ \left. + \frac{ \Delta t w_{2}}{ \Delta z \left( jc \right)A_{\text{c}}\left( jc \right)} + \frac{k_{5}\left( jc \right) \Delta t h_{\text{b}2}\left( jc \right)}{2{\overline{c}}_{\text{c}}\left( jc \right)}\ \right\} \end{multline}$$

with

(3.5-79) $$h_{\text{b}1}\left( jc \right) = \frac{1}{R_{\text{ec}1}\left( jc \right) } + \frac{\gamma_{2}\left( jc \right)}{R_{\text{sc}1}\left( jc \right) }$$

(3.5-80) $$h_{\text{b}2}\left( jc \right) = \frac{1}{R_{\text{ec}2}\left( jc \right) } + \frac{\gamma_{2}\left( jc \right)}{R_{\text{sc}2}\left( jc \right) }$$

and

(3.5-81) $$\varphi_{1}\left( jc \right) = \frac{1}{2} \left\{ \frac{T_{\text{e}2}\left( jc \right)}{R_{\text{ec}2}\left( jc \right)} + \frac{T_{\text{e}1}\left( jc \right)}{R_{\text{ec}1}\left( jc \right)} + \gamma_{2}\left( jc \right) \left\lbrack \frac{T_{\text{si}2}\left( jc \right)}{R_{\text{sc}2}\left( jc \right)} + \frac{T_{\text{si}1}\left( jc \right)}{R_{\text{sc}1}\left( jc \right)} \right\rbrack \right\}$$

If the inlet flow is positive, then the coolant temperature at node 1 is determined by the inlet temperature. Eq. (3.5-78) is then used to march up the channel through the inlet liquid slug, with $T_{\text{c}2}\left( jc + 1 \right)$ being computed after $T_{\text{c}2}\left( jc \right)$. Similarly, if the flow in the upper liquid slug is downward, then the assembly outlet reentry temperature determines the coolant temperature at the last coolant node. Then an equation similar to Eq. (3.5-79) is used to march down through the upper liquid slug, with $T_{\text{c}2} (jc)$ being computed after $T_{\text{c}2}\left( jc + 1 \right)$.

The Eulerian calculations always go from node to node in the direction of flow. An inlet slug expelling downward and an outlet liquid slug going upward are special cases, since in these cases the calculation starts at a liquid vapor interface rather than an end of the subassembly. The interface liquid temperatures are first calculated using the Lagrangian treatment described below. Then Eq. (3.5-78) or the equivalent equation for downward flow is used to calculate the temperatures at the fixed nodes within the liquid slug. For the first fixed node near the interface, some of the terms in Eq. (3.5-78) are modified. The moving interface node is treated as node $jc$. The interface temperature is used for $T_{\text{c}2}\left( jc \right)$, and an interpolated value is used for $T_{\text{c}1}\left( jc \right)$. The distance from the fixed node to the interface at the end of the step is used for $\Delta z \left( jc \right)$. Interpolated interface cladding and structure temperatures are used in calculating $\phi_1$ for the interface node.

3.5.5.2. Lagrangian Calculations for Interface Temperatures

For every liquid-vapor interface a vapor temperature is calculated at or very near the interface. The liquid temperature right at the interface would be close to the vapor temperature, but there can be strong axial temperature gradients in the liquid near the interface. These strong axial gradients would only extend a short distance into the liquid. The heat flow through the interface into a small vapor bubble is accounted for, as described in Chapter 12; but since only one liquid temperature node is used near the interface, the axial temperature distribution neat the interface is not represented. A liquid temperature is calculated for each interface, but axial conduction is neglected in this calculation. Thus, the liquid interface temperature can be considered as either the interface temperature that would occur if there were no axial conduction or the temperature a short distance from the interface where axial conduction is negligible.

A Lagrangian formulation, moving with the liquid, is used for the interface temperature calculation. The basic equation used is

(3.5-82) $$\rho c \frac{\text{DT}_{\text{c}}}{\text{Dt}} = Q_{\text{c}} + Q_{\text{ec}} + Q_{\text{sc}}$$

where the Lagrangian total derivative is used. After finite differencing this equation gives

$$\begin{multline} \rho_{\text{ci}}c_{\text{ci}}\frac{T_{\text{li}2}\left( k,L \right) - T_{\text{li}1}\left( k,L \right)}{ \Delta t } \\ = Q_{\text{c}} \left( jc \right) + \frac{k_{5}\left( jc \right)}{2} \left\lbrack \frac{T_{\text{ei}1}\left( k,L \right) - T_{\text{li}1}\left( k,L \right)}{R_{\text{eci}1}\left( k,L \right)} \right. \\ \left. + \frac{T_{\text{ei}2}\left( k,L \right) - T_{\text{li}2}\left( k,L \right)}{R_{\text{sci}1}\left( k,L \right)} \right\rbrack \\ + \gamma_{2}\left( jc \right)\left\lbrack \frac{T_{\text{si}1}\left( k,L \right) - T_{\text{li}1}\left( k,L \right)}{R_{\text{sci}1}\left( k,L \right)} + \frac{T_{\text{si}2}\left( k,L \right) - T_{\text{li}2}\left( k,L \right)}{R_{\text{sci}2}\left( k,L \right)} \right\rbrack \end{multline}$$

or

(3.5-84) $$T_{\text{li}2} \left( k,L \right) = \frac{T_{\text{li}1}\left( k,L \right)\left\lbrack 1 - d_{1}\ h_{\text{bi}1}\left( k,L \right) \right\rbrack + d_{1} \left\lbrack \varphi_{1\text{i}} \left( k,L \right) + \frac{2 Q_{\text{c}}\left( jc \right)}{k_{5}\left( jc \right)} \right\rbrack}{ \left\lbrack 1 + d_{1}h_{\text{bi}2}\left( k,L \right) \right\rbrack}$$

where

(3.5-85) $$d_{1} = \frac{k_{5}\left( jc \right)\Delta t}{2\rho_{\text{ci}}\ c_{\text{ci}}}$$

(3.5-86) $$h_{\text{bi}1}\left( k,L \right) = \frac{1}{R_{\text{eci}1}\left( k,L \right)} + \frac{\gamma_{2}\left( jc \right)}{R_{\text{sci}1}\left( k,L \right)}$$

(3.5-87) $$h_{\text{bi}2}\left( k,L \right) = \frac{1}{R_{\text{eci}2}\left( k,L \right)} + \frac{\gamma_{2}\left( jc \right)}{R_{\text{sci}2}\left( k,L \right)}$$

$k$ = bubble number

$L$ = 1 for lower bubble interface, 2 for upper bubble interface

$jc$ = coolant node containing the interface

(3.5-88) $$\theta_{1\text{i}}\left( k,L \right) = \frac{T_{\text{ei}2}\left( k,L \right)}{R_{\text{eci}2}\left( k,L \right)} + \frac{T_{\text{ei}1}\left( k,L \right)}{R_{\text{eci}1}\left( k,L \right)} + \gamma_{2}\left( jc \right) \left\lbrack \frac{T_{\text{si}2}\left( k,L \right)}{R_{\text{sci}2}\left( k,L \right)} + \frac{T_{\text{si}1}\left( k,L \right)}{R_{\text{sci}1}\left( k,L \right)} \right\rbrack$$

$T_{\text{ei}2}$, $T_{\text{si}2}$ = cladding and structure interface temperatures at the end of the step, extrapolated in time and interpolated between fixed cladding nodes.

$T_{\text{ei}1}$, $T_{\text{si}1}$ = same at the beginning of the time step.

$R_{\text{eci}2}$, $R_{\text{sci}2}$ = values of $R_{\text{ec}}$ and $R_{\text{sc}}$ at the interface at the end of the time step.

$R_{\text{eci}1}$, $R_{\text{sci}1}$ = same at the beginning of the time step.

3.5.5.3. Lagrangian Calculation for Fixed Nodes

The Lagrangian temperature calculations for fixed coolant nodes are similar to those for interface nodes. The fluid particle that ends up at coolant node $jc$ at the end of a time step is considered. During the time step, the particle travelled a distance

(3.5-89) $$\Delta z' = \frac{\left( w_{1} + w_{2} \right) \Delta t }{2\rho_{\text{c}}\left( jc \right) A_{\text{c}}\left( jj \right)}$$

where

(3.5-90) $$\begin{split}jj = \begin{cases} jc & \text{if}\ w_{1} + w_{2}\ < \ 0 \\ jc - 1 & \text{otherwise} \\ \end{cases}\end{split}$$

At the beginning of the time step, the particle was at $z'$, given by

(3.5-91) $$z' = z_{\text{c}}\left( jc \right) - \Delta z'$$

The coolant temperature, ${T'}_{\text{c}1}$, at $z'$ at the beginning of the step is obtained by linear interpolation between the nodes on either side of $z'$. Also, the cladding and structure temperatures, ${T'}_{\text{e}1}$ and ${T'}_{\text{s}1}$, at $z'$ at the beginning of the step are obtained by linear interpolation. The cladding and structure temperatures, ${T'}_{\text{e}2}$ and ${T'}_{\text{s}2}$, at $z_{\text{c}}\left( jc \right)$ at the end of the time step are also obtained by linear interpolation between the cladding and structure nodes.

The result of finite differencing of Eq. (3.5-82) for the particle at node $jc$ is

(3.5-92) $$T_{\text{c}2}\left( jc \right) = \frac{\left\{ {T'}_{\text{c}1} \left\lbrack 1 - {d'}_{1} h_{\text{b}1}\left( jj \right) \right\rbrack + {d'}_{1} \left\lbrack {\theta'}_{1} + 2 \frac{ Q_{\text{c}}\left( jj \right)}{k_{5}\left( jj \right)} \right\rbrack \right\} }{\left\lbrack 1 + {d'}_{1}\ h_{\text{b}2}\left( jj \right) \right\rbrack}$$

where

(3.5-93) $${d'}_{1} = \frac{k_{5}\left( jj \right)\Delta t}{\rho_{\text{c}}\left( jc \right)\ {\overline{c}}_{\text{c}}\left( jj \right)}$$

and

(3.5-94) $$\phi'_{1} = \frac{T_{\text{e}2}'}{R_{\text{ec}2}(jj)} + \frac{T_{\text{e}1}'}{R_{\text{ec}1}(jj)} + \gamma_{2}\left( jj \right)\left\lbrack \frac{T_{\text{s}2}'}{R_{\text{sc}2}(jj)} + \frac{T_{\text{s}1}'}{R_{\text{sc}1}(jj)} \right\rbrack$$

Again, $h_{\text{b}1}$ and $h_{\text{b}2}$are given by Eq. (3.5-79) and Eq. (3.5-80).