.. _section-12.6:
Vapor Pressure Gradient Model: Large Vapor Bubbles
--------------------------------------------------
As mentioned above, whenever a bubble grows to a length greater than a
user-specified value (usually 5-50 cm), the vapor bubble calculation is
switched from the uniform-vapor-pressure model to a
vapor-pressure-gradient model. Saturation conditions are assumed in this
model, so the vapor energy equation can be combined with the mass
continuity equation. The vapor continuity and momentum equations are
solved simultaneously for each node in a bubble using a fully implicit
finite differencing scheme.
The equations described in this section are only for the vapor. Any
liquid in a bubble region is assumed to be a film on the cladding and
structure, and this film is treated separately, as described at the
beginning of :numref:`section-12.1`. The vapor calculation is influenced by the
liquid film in only two ways: heat flow through the film provides the
vapor source for the vapor calculations; and if a two-phase friction
factor multiplier for the friction between the vapor and the film is
used, then the film thickness determines the value of this multiplier.
.. _section-12.6.1:
Continuity and Momentum Equations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The vapor continuity equation can be written as
(12.6-1)
.. _eq-12.6-1:
.. math:: \frac{\partial W}{\partial z} + \frac{\partial\left( \rho_{v}A_{c} \right)}{\partial t} = \frac{Q}{\lambda}A_{C}
while the vapor momentum equation is given by
(12.6-2)
.. _eq-12.6-2:
.. math:: \frac{\partial p}{\partial z} = - \frac{1}{A_{c}}\frac{\partial}{\partial z}\frac{W^{2}}{\rho_{v}A_{c}} - \frac{f_{v}W\left| W \right|}{2\rho_{v}D_{h}A_{c}^{2}} + F_{c} - \frac{1}{A_{c}}\frac{\partial W}{\partial t} - \frac{W\left| W \right|}{2\rho_{v}A_{c}^{2}}\frac{\partial k_{or}}{\partial z}
The symbols used in Eqs. :ref:`12.6-1 <eq-12.6-1>` and :ref:`12.6-2 <eq-12.6-2>` are defined as follows:
:math:`W =` vapor mass flow rate
:math:`z =` axial height
:math:`\rho_{v} =` sodium vapor density
:math:`A_{c} =` coolant flow area
:math:`t =` time
:math:`Q =` heat flow rate per unit volume
:math:`\lambda =` heat of vaporization
:math:`p =` pressure
:math:`f_{v} =` friction factor
:math:`D_{h} =` hydraulic diameter
:math:`F_{c} =` condensation momentum loss term
:math:`k_{or} =` orifice pressure drop term
Note that, because the flow area :math:`A_{c}` is included in the mass and
momentum equations as a function of both space and time, the model is
able to treat problems involving variable flow areas.
The heat flow rate per unit volume is given by
(12.6-3)
.. _eq-12.6-3:
.. math:: Q = \gamma\left\lbrack \frac{T_{e} - T}{R_{ec}} + \gamma_{2}\frac{T_{s} - T}{R_{sc}} \right\rbrack
where
:math:`\lambda =` ratio of cladding surface area to coolant volume
:math:`T_{e} =` cladding temperature
:math:`T_{s} =` structure temperature
:math:`T =` vapor temperature
:math:`R_{ec} =` thermal resistance between the cladding and the vapor
:math:`R_{sc} =` thermal resistance between the structure and the vapor
:math:`\gamma_{2} =` ratio of structure surface are to cladding surface area
The thermal resistances :math:`R_{ec}` and :math:`R_{sc}` account for the resistance
to heat flow of the cladding, liquid film, and vapor and are computed
according to the expressions discussed in :numref:`section-12.5.1`.
The expression for the heat flow rate :math:`Q` can be rewritten as the sum of
the heat flow rate through the cladding, :math:`Q_{e}`, and the heat flow rate
through the structure, :math:`Q_{s}`, where
(12.6-4)
.. _eq-12.6-4:
.. math:: Q_{e} = \gamma\frac{T_{e} - T}{R_{ec}}
and
(12.6-5)
.. _eq-12.6-5:
.. math:: Q_{s} = \gamma\gamma_{2}\frac{T_{s} - T}{R_{sc}}
These two quantities are used to compute the condensation momentum loss
term :math:`F_{c}`. This term accounts for momentum lost from the system
when flowing sodium vapor condenses on cladding or structure that is
colder than the vapor. Only a momentum loss term is included, since it
is assumed that the system does not gain momentum when vaporization
occurs from cladding or structure that is hotter than the sodium film.
The condensation momentum loss is simply the product of the rate of mass
condensing per unit volume, :math:`\frac{Q}{\lambda}`, times the average
velocity, :math:`\frac{W}{2\rho_{v}A_{c}}`. Splitting the condensation
momentum loss term into separate contributions from the cladding and
structure then gives.
(12.6-6)
.. _eq-12.6-6:
.. math:: F_{c} = F_{ce} + F_{cs}
where
(12.6-7a, b)
.. _eq-12.6-7a:
.. _eq-12.6-7b:
.. math::
F_{ce} = \Bigg\{\begin{matrix}
0,Q_{e} \geq 0, \\
\frac{Q_{e}W}{2\rho_{v}A_{c}\lambda},Q_{e} < 0 \\
\end{matrix}
and
(12.6-8a,b)
.. _eq-12.6-8a:
.. _eq-12.6-8b:
.. math::
F_{cs} = \Bigg\{ \begin{matrix}
0,Q_{s} \geq 0, \\
\frac{Q_{s}W}{2\rho_{v}A_{c}\lambda},\ Q_{s} < 0 \\
\end{matrix}
These definitions reflect the assumption that momentum is lost when
condensation of vapor occurs (:math:`Q_{e}`, :math:`Q_{s}` negative) but that no
momentum change occurs when liquid is vaporized (:math:`Q_{e}`, :math:`Q_{s}` positive).
The friction factor, :math:`f_{v}` is expressed as
(12.6-9)
.. _eq-12.6-9:
.. math:: f_{v} = A_{fr}\left( \frac{D_{h}\left| W \right|}{\mu A_{c}} \right)^{b_{fr}}f_{tp}
where :math:`\mu` is the vapor viscosity and :math:`A_{fr}` and :math:`b_{fr}` are input
constants. The quantity in parentheses, :math:`D_{h}|W|/\left( \mu A_{c} \right)`,
is just the Reynolds number of the sodium vapor. The two-phase friction factor
multiplier,
:math:`f_{tp}`, is based on a correlation by Wallis :ref:`[12-10] <section-12.references>`:
(12.6-10)
.. _eq-12.6-10:
.. math:: f_{tp} = 1 + 300 f_{2\phi}(w_{fe} + \gamma_{2}w_{fs})/\lbrack\left( 1 + \gamma_{2} \right)D_{h}\rbrack
where :math:`f_{2\phi}` is set to 1.0 when the two-phase multiplier is used and is set to 0
by the user if a smooth tube friction factor is described.
Two approximations are made to Eqs. :ref:`12.6-1 <eq-12.6-1>` and :ref:`12.6-2 <eq-12.6-2>`
before finite differencing. First, in the continuity equation, the term
(12.6-11)
.. _eq-12.6-11:
.. math:: \frac{\partial(\rho_{v}A_{c})}{\partial t} = \rho_{v}\frac{\partial A_{c}}{\partial t} + A_{c}\frac{\partial\rho_{v}}{\partial t}
is taken to be
(12.6-12)
.. _eq-12.6-12:
.. math:: \frac{\partial(\rho_{v}A_{c})}{\partial t} = A_{c}\frac{\partial\rho_{v}}{\partial t}
So that :math:`A_{c}` is assumed constant in time (but not space) over the
time step for purposes fo solving the mass and momentum equations for
the vapor pressures and mass flow rates. This assumption is also made
with the variables :math:`D_{h}`, :math:`\gamma`, and :math:`\gamma_{2}`. This simplification
is a reasonable one, since one of these four quantities changes rapidly
with time. Second, in the momentum equations, the momentum convection
term can be expanded as
(12.6-13)
.. _eq-12.6-13:
.. math:: \frac{\partial}{\partial z}\left( \frac{W^{2}}{\rho_{v}A_{c}} \right) = \frac{1}{A_{c}}\frac{\partial}{\partial z}\left( \frac{W^{2}}{\rho_{v}} \right) + \frac{W^{2}}{\rho_{v}}\frac{\partial}{\partial z}\left( \frac{1}{A_{c}} \right)
Any changes in area over a spatial node are incorporated into the
orifice term, and so
:math:`\frac{W^{2}}{\rho_{v}}\frac{\partial}{\partial z}(\frac{1}{A_{c}})` is
included in the :math:`\frac{\partial k_{or}}{\partial z}` term.
Therefore,
(12.6-14)
.. _eq-12.6-14:
.. math:: \frac{\partial}{\partial z}\left( \frac{W^{2}}{\rho_{v}A_{c}} \right) = \frac{1}{A_{c}}\frac{\partial}{\partial z}\left( \frac{W^{2}}{\rho_{v}} \right)
so :math:`A_{c}` appears piecewise constant in space. This is also true for
:math:`D_{h}`, :math:`\gamma`, and :math:`\gamma_{2}`.
Therefore, the continuity and momentum equations to be differenced are
(12.6-15)
.. _eq-12.6-15:
.. math:: \frac{\partial W}{\partial z} + A_{c}\frac{\partial\rho_{v}}{\partial t} = \frac{Q}{\lambda}A_{c}
and
(12.6-16)
.. _eq-12.6-16:
.. math:: \frac{\partial p}{\partial z} = - \frac{1}{A_{c}^{2}}\frac{\partial}{\partial z}\left( \frac{W^{2}}{\rho_{v}} \right) - A_{fr}\frac{W\left| W \right|^{1 + b_{fr}}}{2\rho_{v}D_{h}^{1 - b_{fr}}A_{c}^{2 + b_{fr}}} + F_{c} - \frac{1}{A_{c}}\frac{\partial W}{\partial t} - \frac{W\left| W \right|}{2\rho_{v}A_{c}^{2}}\frac{\partial k_{or}}{\partial z}
.. _section-12.6.2:
Finite Differencing
~~~~~~~~~~~~~~~~~~~
The mass and momentum equations must now be reduced to algebraic
equations that can be coded into SASSYS-1. This is accomplished in three
steps: (1) the equations are discretized using finite differences, (2)
the advanced time terms are linearized, and (3) saturation conditions
are imposed to limit the independent variables to the change in mass
flow rate and the change in pressure. The details of these three steps
are presented below. The result is the set of algebraic equations
:ref:`12.6-45 <eq-12.6-45>`, :ref:`12.6-79 <eq-12.6-79>`, :ref:`12.6-87 <eq-12.6-87>`
(discretized form of the continuity equation at an interior node, lower
interface, and upper interface of a bubble, respectively), :ref:`12.6-128 <eq-12.6-128>`,
:ref:`12.6-143 <eq-12.6-143>`, and :ref:`12.6-149 <eq-12.6-149>`
(discretized form of the momentum equation in the interior, the
lower interface, and the upper interface). These equations are then solved
according to the procedure discussed in :numref:`section-12.6.3`.
As indicated in :numref:`figure-12.6-1`, some of the variables in the mass and momentum
equations are calculated at node boundaries, whereas others are calculated at
node mid-points. In general, midpoint quantities are used as if they were
constant over the node. The quantities defined at node boundaries are
:math:`z`, :math:`T`, :math:`p`, :math:`W`, :math:`\rho_{v}`, :math:`\mu`, and the velocity :math:`V`.
At liquid-vapor interfaces, these quantities are defined at the interface, and
the vapor velocity is set to the value of the interface velocity.
Midpoint quantities include :math:`T_{e}`, :math:`T_{s}`, :math:`R_{ec}`, :math:`A_{c}`, :math:`D_{h}`, :math:`\gamma`, :math:`\gamma_{2}`, :math:`k_{or}`, :math:`\gamma`, :math:`Q`, :math:`F_{e}`,
and the spatially averaged coolant temperature :math:`\overline{T}`. With all
variables, a subscript 1 indicates that the variable is evaluated at the
beginning of the time step; a subscript 2 means the variable is taken to be at
the end of the time step.
.. _figure-12.6-1:
.. figure:: media/Figure12.6-1.png
:align: center
:figclass: align-center
SASSYS-1 Voiding Model Axial Coolant Mesh Variable Placement
.. _section-12.6.2.1:
Finite Differencing of the Continuity Equation for a Mesh Segment Contained Entirely in One Bubble
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
This section gives a detailed explanation of how the continuity equation
is expressed in finite difference form. The extra factors required to
model bubble-liquid interfaces correctly are not included here; they
will be discussed in :numref:`section-12.6.2.2`. Therefore, the equations to be
developed in this section apply as they stand only to axial mesh
segments which are entirely contained within the bubble.
.. _section-12.6.2.1.1:
Differencing of the Continuity Equation
'''''''''''''''''''''''''''''''''''''''
The terms making up the continuity equation are differenced as follows:
(12.6-17)
.. _eq-12.6-17:
.. math:: \left. \frac{\partial W}{\partial z} \right|_{J + 1/2} = \frac{W_{2}(J + 1) - W_{2}(J)}{z(J + 1) - z(J)}
and
(12.6-18)
.. _eq-12.6-18:
.. math:: A_{c}\left. \frac{\partial\rho_{v}}{\partial t} \right|_{J + 1/2} = A_{c}(J)\frac{\rho_{2}(J) - \rho_{1}(J) + \rho_{2}(J + 1) - \rho_{1}(J + 1)}{2\Delta t}
These differencing schemes are illustrated in :numref:`figure-12.6-2`. Also,
(12.6-19)
.. _eq-12.6-19:
.. math:: Q_{2}(J) = \gamma(J)\left\lbrack \frac{T_{e_2}(J) - {\overline{T}}_{2}(J)}{R_{{ec}_2}} + \gamma_{2}(J)\frac{T_{s_{2}}(J) - {\overline{T}}_{2}(J)}{R_{sc_{2}}(J)} \right\rbrack
.. math:: = Q_{e_{2}}(J) + Q_{s_{2}}(J)
Therefore, the differenced continuity equation becomes
(12.6-20)
.. _eq-12.6-20:
.. math:: \frac{W_{2}(J + 1) - W_{2}(J)}{z(J + 1) - z(J)} + A_{c}(J)\frac{\rho_{2}(J) - \rho_{1}(J) + \rho_{2}(J + 1) - \rho_{1}(J + 1)}{2\Delta t}
.. math:: = A_{c}(J)\frac{Q_{2}(J)}{\lambda_{2}(J)}
.. _figure-12.6-2:
.. figure:: media/Figure12.6-2.jpeg
:align: center
:figclass: align-center
Evaluation of Derivatives at an Interior Node in the SASSYS-1
Voiding Model
.. _section-12.6.2.1.2:
Linearization of Advanced Time Terms
''''''''''''''''''''''''''''''''''''
Equation :ref:`12.6-20 <eq-12.6-20>` can be linearized by expressing the advanced time terms
in the form
(12.6-21)
.. _eq-12.6-21:
.. math:: x_{2}(J) = x_{1}(J) + \Delta x(J)
where *x* is any quantity. This produces
(12.6-22)
.. _eq-12.6-22:
.. math:: \rho_{2}(J) = \rho_{1}(J) + \Delta\rho(J)
(12.6-23)
.. _eq-12.6-23:
.. math:: \lambda_{2}(J) = \lambda_{1}(J) + \Delta\lambda(J)
(12.6-24)
.. _eq-12.6-24:
.. math:: W_{2}(J) = W_{1}(J) + \Delta W(J)
(12.6-25)
.. _eq-12.6-25:
.. math:: T_{2}(J) = T_{1}(J) + \Delta T(J)
Since
(12.6-26)
.. _eq-12.6-26:
.. math:: \overline{T}_{2}(J) = \frac{1}{2}\left( T_{2}(J) + T_{2}(J + 1) \right)
.. math:: = \frac{1}{2}\left( T_{1}(J) + \Delta T(J) + T_{1}(J + 1) + \Delta T(J + 1) \right)
.. math:: = \overline{T}_{1}(J) + \frac{1}{2}(\Delta T(J) + \Delta T(J + 1))
then
(12.6-27)
.. _eq-12.6-27:
.. math:: \Delta\overline{T}(J) = \frac{1}{2}\left( \Delta T(J) + \Delta T(J + 1) \right)
The heat flux then becomes
(12.6-28)
.. _eq-12.6-28:
.. math:: Q_{2}(J) = \gamma(J)\left[\frac{T_{e_{2}}(J) - {\overline{T}}_{1}(J) - \Delta\overline{T}(J)}{R_{ec_{2}}(J)} + \gamma_{2}(J)\frac{T_{s_{2}}(J) - {\overline{T}}_{1}(J) - \Delta\overline{T}(J)}{R_{{sc}_2}(J)}\right]
.. math:: = \gamma(J)\Bigg[ \frac{T_{e_{2}}(J) - {\overline{T}}_{1}(J)}{R_{ec_{2}}(J)} + \gamma_{2}(J)\frac{T_{s_{2}}(J) - {\overline{T}}_{1}(J)}{R_{{sc}_2}(J)}
.. math:: - \frac{\left( \Delta T(J) + \Delta T(J + 1) \right)}{2}\left\lbrack \frac{1}{R_{ec_{2}}(J)} + \frac{\gamma_{2}(J)}{R_{sc_{2}}(J)} \right\rbrack \Bigg]
The heat flux expression can be written in a more compact form by
defining the variable :math:`Q_{o}(J)`
(12.6-29)
.. _eq-12.6-29:
.. math:: Q_{0}(J) = \gamma(J)\left\lbrack \frac{T_{e_{2}}(J) - \overline{T}(J)}{R_{ec_{2}}(J)} + \gamma_{2}(J)\frac{T_{s_{2}}(J) - {\overline{T}}_{1}(J)}{R_{sc_{2}}(J)} \right\rbrack
The derivative with respect to vapor temperature of :math:`Q_{o}` is then
(12.6-30)
.. _eq-12.6-30:
.. math:: \frac{dQ_{0}}{dT_{J}} = - \gamma(J)\left\lbrack \frac{1}{R_{ec_{2}}(J)} + \frac{\gamma_{2}(J)}{R_{sc_{2}}(J)} \right\rbrack
These can be substituted into the above equation for
:math:`Q_{2}(J)` to give
(12.6-31)
.. _eq-12.6-31:
.. math:: Q_{2}(J) = Q_{0}(J) + \frac{1}{2}\frac{dQ_{0}}{dT_{J}}\left( \Delta T(J) + \Delta T(J + 1) \right)
Substituting Eqs. :ref:`12.6-22 <eq-12.6-22>` through :ref:`12.6-25 <eq-12.6-25>` and
:ref:`12.6-31 <eq-12.6-31>` for the advanced time terms in the differenced continuity
equation, Eq. :ref:`12.6-20 <eq-12.6-20>` gives
(12.6-32)
.. _eq-12.6-32:
.. math:: \frac{W_{1}(J + 1) - W_{1}(J) + \Delta W(J + 1) - \Delta W(J)}{z(J + 1) - z(J)} + A_{c}(J)\frac{\Delta\rho(J) + \Delta\rho(J + 1)}{2\Delta t}
.. math:: = \frac{A_{c}(J)}{\lambda_{1}(J) + \Delta\lambda(J)}\left( Q_{0}(J) + \frac{1}{2}\frac{dQ_{0}}{dT_{J}}\left( \Delta T(J) + \Delta T(J + 1) \right) \right)
If :math:`\Delta \lambda (J) \ll \lambda_{1}(J)`,
(12.6-33)
.. _eq-12.6-33:
.. math:: \frac{1}{\lambda_{1}(J) + \Delta\lambda(J)} \approx \frac{1}{\lambda_{1}(J)}\left[1 - \frac{\Delta\lambda(J)}{\lambda_{1}(J)}\right]
which changes Eq. :ref:`12.6-32 <eq-12.6-32>` to
(12.6-34)
.. _eq-12.6-34:
.. math:: \frac{W_{1}(J + 1) - W_{1}(J) + \Delta W(J + 1) - \Delta W(J)}{z(J + 1) - z(J)} + A_{c}(J)\frac{\Delta\rho(J) + \Delta\rho(J + 1)}{2\Delta t}
.. math:: = A_{c}(J)\frac{Q_{0}(J)}{\lambda_{1}(J)}\left( 1 + \frac{\Delta T(J) + \Delta T(J + 1)}{2Q_{0}(J)}\frac{dQ_{0}}{dT_{J}} \right)\left( 1 - \frac{\Delta\lambda(J)}{\lambda_{1}(J)} \right)
Eliminating second-order terms gives the linearized differenced
continuity equation
(12.6-35)
.. _eq-12.6-35:
.. math:: \frac{W_{1}(J + 1) - W_{1}(J) + \Delta W(J + 1) - \Delta W(J)}{z(J + 1) - z(J)} + A_{c}(J)\frac{\Delta\rho(J) + \Delta\rho(J + 1)}{2\Delta t}
.. math:: = A_{c}(J)\frac{Q_{0}(J)}{\lambda_{1}(J)}\left(1 + \frac{\Delta t(J) + \Delta T(J + 1)}{2Q_{0}(J)}\frac{dQ_{0}}{dT_{J}} - \frac{\Delta\lambda(J)}{\lambda_{1}(J)}\right)
.. _section-12.6.2.1.3:
Utilization of Saturation Conditions to Limit the Independent Variables to :math:`\Delta W` and :math:`\Delta p`
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Since saturation conditions are assumed, :math:`\rho`, :math:`T`, and :math:`\lambda` are all
functions of :math:`p` only. Therefore, the difference terms can be expressed
in terms of :math:`\Delta p` as follows:
(12.6-36)
.. _eq-12.6-36:
.. math:: \Delta T(J) = \Delta p(J)\frac{dT}{dp_{J}}
(12.6-37)
.. _eq-12.6-37:
.. math:: \Delta\rho(J) = \Delta p(J)\frac{d\rho}{dp_{J}}
(12.6-38)
.. _eq-12.6-38:
.. math:: \Delta\lambda(J) = \frac{1}{2}\left( \Delta p(J) + \Delta p(J + 1) \right)\frac{d \lambda}{dp_{J}}
with :math:`\Delta p (J) = p_{2} (J) - p_{1} (J)`.
All derivatives are evaluated at the start of the time step. When Eqs. :ref:`12.6-36 <eq-12.6-36>`
through :ref:`12.6-38 <eq-12.6-38>`. are inserted into Eq. :ref:`12.6-35 <eq-12.6-35>` the result is a differenced
continuity equation which is linear in :math:`\Delta W` and :math:`\Delta p`:
(12.6-39)
.. _eq-12.6-39:
.. math:: \frac{W_{1}(J + 1) - W_{1}(J) + \Delta W(J + 1) - \Delta W(J)}{z(J + 1) - z(J)} + A_{c}(J)\frac{1}{2\Delta t}\left( \Delta p(J)\frac{d\rho}{dp_{J}} + \Delta p(J + 1)\frac{d\rho}{dp}(J + 1) \right)
.. math:: = A_{c}(J)\frac{Q_{0}(J)}{\lambda_{1}(J)}\Bigg[ 1 + \left( \Delta p(J)\frac{dT}{dp_{J}} + \Delta p(J + 1)\frac{dT}{dp_{J + 1}} \right)\frac{1}{2Q_{0}(J)}\frac{dQ_{0}}{dT_{J}}
.. math:: - \left( \Delta p(J) + \Delta p(J + 1) \right)\frac{1}{2\lambda_{1}(J)}\frac{d\lambda}{dp_{J}} \Bigg]
Multiply through by :math:`\left(z(J + 1) - z(J)\right) \Delta t`
and define the coefficients
(12.6-40)
.. _eq-12.6-40:
.. math:: c_{1,j} = \left( z(J + 1) - z(J) \right)A_{c}(J)\left\lbrack\frac{1}{2}\frac{d\rho}{dp_{J}} - \frac{\Delta t}{2\lambda_{1}(J)}\frac{dT}{dp_{J}}\frac{dQ_{0}}{dT_{J}} + \frac{Q_{0}(J)\Delta t}{2\lambda_{1}^{2}(J)}\frac{d\lambda}{dp_{J}}\right\rbrack
(12.6-41)
.. _eq-12.6-41:
.. math:: c_{2,j} = - \Delta t
(12.6-42)
.. _eq-12.6-42:
.. math:: c_{3,j} = \left( z(J + 1) - z(J) \right)A_{c}(J)\left\lbrack \frac{1}{2}\frac{d\rho}{dp_{J + 1}} - \frac{\Delta t}{2\lambda_{1}(J)}\frac{dT}{dp_{J + 1}}\frac{dQ_{0}}{dT_{J}} + \frac{Q_{0}(J)\Delta t}{2\lambda_{1}^{2}(J)}\frac{d\lambda}{dp_{J}} \right\rbrack
(12.6-43)
.. _eq-12.6-43:
.. math:: c_{4,J} = \Delta t
Also, define the source term
(12.6-44)
.. _eq-12.6-44:
.. math:: h_{J} = W_{1}(J) - W_{1}(J + 1) + \left( z(J + 1) - z(J) \right)A_{c}(J)\frac{Q_{0}(J)}{\lambda_{1}(J)}\Delta t
Then the differenced mass equation becomes
(12.6-45)
.. _eq-12.6-45:
.. math:: c_{1,J}\Delta p(J) + c_{2,J}\Delta W(J) + c_{3,J}\Delta p(J + 1) + c_{4,J}\Delta W(J + 1) = h_{j}
.. _section-12.6.2.2:
Finite Differencing of the Continuity Equation for a Mesh Segment Which Contains a Bubble Interface
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. _section-12.6.2.2.1:
Adjustment of the Form of the Continuity Equation to Account for Interface Motion
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
The representation of the partial derivative with respect to time is not
as straightforward at an interface as it is at an interior node.
Referring to :numref:`figure-12.6-3`, in which node boundary :math:`J` is the interface,
the partial derivative of the density :math:`\rho` is still calculated as
(12.6-46)
.. _eq-12.6-46:
.. math:: \left. \frac{\partial\rho}{\partial t} \right|_{J + 1/2} = \frac{1}{2}\left\lbrack \left. \frac{\partial\rho}{\partial t} \right|_{J} + \left. \frac{\partial\rho}{\partial t} \right|_{J + 1} \right\rbrack
The partial derivatives at :math:`J` and :math:`J + 1` are computed just as they are
for an interior node; in particular,
(12.6-47)
.. _eq-12.6-47:
.. math:: \left. \frac{\partial\rho}{\partial t} \right|_{J} = \frac{\rho\left( z\left( J,t \right),t + \Delta t \right) - \rho\left(z\left( J,t \right),t\right)}{\Delta t}
.. _figure-12.6-3:
.. figure:: media/Figure12.6-3.jpeg
:align: center
:figclass: align-center
Evaluation of the total Derivative at a Liquid-Vapor Interface
so the partial derivative at :math:`J` is differenced as though the interface
were stationary at :math:`z\left(J,t\right)`. However, the interface is moving at
a velocity :math:`v^{i}` and moves a distance :math:`v^{i}\Delta t` from
:math:`z\left(J,t\right)` to :math:`z\left(J,t + \Delta t\right)` during the time step, so
the spatial derivative in the continuity equation must be differenced
from :math:`z\left(J + 1,t + \Delta t\right)` to :math:`z\left(J, t + \Delta t\right)`, rather
than to :math:`z\left(J,t\right)`. Thus, the continuity equation for an interface
node will involve advanced time quantities at three spatial points
(:math:`z \left(J,t \right)`, :math:`z \left( J,t + \Delta t \right)`, and
:math:`z\left(J + 1, t + \Delta t\right)`)
rather than at two. This difficulty can be resolved if the partial derivative
with respect to time is expressed in terms of :math:`\rho\left(z\left(J,t + \Delta t\right), t + \Delta t\right)`
rather than :math:`\rho\left(z\left(J,t\right), t + \Delta t\right)`
by introducing the Lagrangian total derivative, :math:`\frac{d\rho}{dt}`.
(12.6-48)
.. _eq-12.6-48:
.. math:: \frac{d\rho}{dt} = \frac{\partial\rho}{\partial t} + v\frac{\partial\rho}{\partial z}
or
(12.6-49)
.. _eq-12.6-49:
.. math:: \frac{\partial\rho}{\partial t} = \frac{d\rho}{dt} - v\frac{\partial\rho}{\partial z}
with
(12.6-50)
.. _eq-12.6-50:
.. math:: \frac{d\rho}{dt} = \frac{\rho\left( z\left( J,t + \Delta t \right),t + \Delta t \right) - \rho\left( z\left( J,t \right),t \right) + \rho\left( z(J + 1),t + \Delta t \right) - \rho(z(J + 1),t)}{2\Delta t}
(12.6-51)
.. _eq-12.6-51:
.. math:: v = \frac{v^{i}}{2}
and
(12.6-52)
.. _eq-12.6-52:
.. math:: \frac{\partial\rho}{\partial z} = \frac{\rho\left( z\left( J,t + \Delta t \right),t + \Delta t \right) - \rho(z(J + 1),t + \Delta t)}{z\left( J,t + \Delta t \right) - z\left(J + 1, t + \Delta t \right)}
The factor of (1/2) appears in the velocity term to represent the average of
the velocity of the interface point :math:`J` (which is :math:`v^{i}`) and that of the
point :math:`J + 1` (which is zero). This formulation is equivalent to
expressing :math:`\rho\left(z\left( J,t\right),t + \Delta t\right)` as an interpolated
value between
:math:`\rho\left(z\left(J,t + \Delta t\right)\right)` and :math:`\rho\left(z(J + 1), t + \Delta t\right)`.
The continuity equation now takes the form
(12.6-53)
.. _eq-12.6-53:
.. math:: \frac{\partial W}{\partial z} + A_{c}\frac{d\rho}{dt} - v\frac{\partial\rho}{\partial z} = A_{c}\frac{Q}{\lambda}
Let a node boundary quantity :math:`x` be represented at the lower interface
by :math:`x_{j}^{i}(J1)` and at the upper interface by :math:`x_{j}^{i}(J2)` where :math:`j = 1` or :math:`2`
to designate the beginning or end of the time step. A node midpoint quantity :math:`y`
is represented at the lower interface by :math:`y_{j}(J1)` (since it is contained
in node :math:`J1`) and at the upper interface by :math:`y_{j}(J2-1)` (since it is
located in node :math:`J2-1`). As shown in Fig. 12.6-4, :math:`J1` is the number of the
fixed mesh node boundary below the lower interface, while :math:`J2` is the fixed
node boundary number above the upper interface. Looking first at the lower
interface, the total time derivative of the density is then
(12.6-54)
.. _eq-12.6-54:
.. math:: \frac{d\rho}{dt} = \frac{\rho_{2}(J1 + 1) - \rho_{1}(J1 + 1) + \rho_{2}^{i}(J1) - \rho_{1}^{i}(J1)}{2\Delta t}
or, using Eqs. :ref:`12.6-22 <eq-12.6-22>` and :ref:`12.6-37 <eq-12.6-37>` to express
:math:`\frac{d\rho}{dt}` in terms of the change in vapor
pressure,
(12.6-55)
.. _eq-12.6-55:
.. math:: \frac{d\rho}{dt} = \frac{1}{2\Delta t}\left(\Delta \rho(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} + \Delta\rho^{i}(J1)\frac{d\rho^{i}}{dp_{J1}} \right)
The spatial derivative term becomes
(12.6-56)
.. _eq-12.6-56:
.. math:: v\frac{\partial\rho}{\partial z} = \frac{v_{2}^{i}(J1)}{2}\frac{(\rho_{2}(J1 + 1) - \rho_{2}^{i}(J1))}{(z(J1 + 1) - z_{2}^{i}(J1))}
The advanced time velocity is expressed in the linearized form of Eq.
:ref:`12.6-21 <eq-12.6-21>` as
(12.6-57)
.. _eq-12.6-57:
.. math:: v_{2}^{i}(J1) = v_{1}^{i}(J1) + \Delta v^{i}(J1)
It is apparent from Eqs. :ref:`12.6-56 <eq-12.6-56>` and :ref:`12.6-57 <eq-12.6-57>`
that the added complication
of treating the bubble interface introduces two extra unknowns beyond
those described in :numref:`section-12.6.2.1`, namely, the change in interface
velocity :math:`\Delta v^{i}` and the advanced time interface position :math:`z_2^i` (the
current time interface velocity, :math:`v_1^i`, is assumed known). As will be shown in :numref:`section-12.6.2.2.2` and :numref:`section-12.6.2.2.3`,
both of these quantities can be expressed as
functions of :math:`\Delta p` and of known terms, so that their effect on the final
differenced equation is to add terms to the coefficients
:math:`c_{1,J}` and :math:`c_{3,J}` and to the source term :math:`h_{J}` in
Eq. :ref:`12.6-45 <eq-12.6-45>`.
.. _figure-12.6-4:
.. figure:: media/Figure12.6-4.jpeg
:align: center
:figclass: align-center
Placement of a Bubble in the SASSYS-1 Axial Coolant Mesh
.. _section-12.6.2.2.2:
Formulation of :math:`\Delta v^i` as a Function of :math:`\Delta p`
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
If the bubble under consideration is labeled bubble :math:`K`, the change in
interface velocity is calculated from the change in mass flow rate
:math:`\Delta W_{l}\left(K\right)` in liquid slug :math:`K`, which is the slug below bubble :math:`K`:
(12.6-58)
.. _eq-12.6-58:
.. math:: \Delta v^{i}(J1) = \frac{\Delta W_{l}\left( K \right)}{\rho_{l}^{i}(J1)A_{c}(J1)}f_{rwt}(K)
where :math:`f_{rwt}\left(K\right)` is the liquid film rewetting factor at the top of
slug :math:`K` and :math:`\rho_{l}^{i}(J1)`
is the liquid sodium density at the interface. The film rewetting factor
accounts for the effect of liquid film on the cladding and structure on the
interface velocity; this has already been discussed in :numref:`section-12.3`.
From Eq. :ref:`12.3-6 <eq-12.3-6>` of :numref:`section-12.3`,
the liquid film rewetting factor is
(12.6-59)
.. _eq-12.6-59:
.. math:: f_{rwt} = \frac{v^{i}}{v_{l}} = \frac{1 - P_{e}\left( w_{fe}v_{fe} + \gamma_{2}w_{fs}v_{fs} \right)/(v_{l}A_{c})}{1 - P_{e}\left( w_{fe} + \gamma_{2}w_{fs} \right)/A_{c}}
The change in mass flow rate is computed from the liquid slug momentum
equation given in Eq. :ref:`12.2-34 <eq-12.2-34>`:
(12.6-60)
.. _eq-12.6-60:
.. math:: \Delta W_{l}\left( K \right) = \frac{\left( AA_{0}\left( K \right) + \theta_{2}\left( \Delta p_{b}\left( K \right) - \Delta p_{t}\left( K \right) \right) \right)\Delta t}{I_{1}\left( K \right) + \theta_{2}\Delta I_{1}\left( K \right) + BB_{0}\left( K \right)\theta_{2}\Delta t}
with :math:`\theta_{1} = \theta_{2} = 1/2` for the semi-implicit method
and :math:`\theta_{1} = 0`, :math:`\theta_{2} = 1` for the fully implicit
formulation. The pressure changes are defined as
:math:`\Delta p_{b}\left(K\right) =` change in pressure at lower end of liquid slug K.
:math:`\Delta p_{t}\left(K\right) =` change in pressure at upper end of liquid slug K.
Since :math:`\Delta p_{t}\left(K\right)` is the pressure change at the bubble interface,
(12.6-61)
.. _eq-12.6-61:
.. math:: \Delta p_{t}\left( K \right) = \Delta p^{i}(J1)
Therefore, :math:`\Delta v^{i}(J1)` can be split into two parts as
(12.6-62)
.. _eq-12.6-62:
.. math:: \Delta v^{i}(J1) = \Delta v_{0}^{i}(J1) + \Delta p^{i}(J1)\frac{dv^{i}}{dp_{J1}}
with
(12.6-63)
.. _eq-12.6-63:
.. math:: \Delta v_{0}^{i}(J1) = \frac{f_{fwt}(K)}{\rho_{l}^{i}(J1)A(J1)}\frac{\left( AA_{0}\left( K \right) + \theta_{2}\Delta p_{b}\left( K \right) \right)\Delta t}{I_{1}\left( K \right) + \theta_{2}\Delta I_{1}\left( K \right) + BB_{0}\left( K \right)\theta_{2}\Delta t}
being the change in interface velocity with no change in interface
pressure and
(12.6-64)
.. _eq-12.6-64:
.. math::
\frac{dv^{i}}{dp_{J1}} = - \frac{f_{rwt}\left( K \right)}{\rho_{l}^{i}(J1)A(J1)}\frac{\theta_{2}\Delta t}{I_{1}\left( K \right) + \theta_{2}\Delta I_{1}\left( K \right) + BB_{0}\left( K \right)\theta_{2}\Delta t}
being the velocity change with respect to the change in interface
pressure. With the change in interface velocity now expressed as a
linear function of the change in interface pressure (Eq. :ref:`12.6-62 <eq-12.6-62>`), the
advanced time interface velocity in the expression for the spatial
derivative term :math:`v\frac{\partial\rho}{\partial z}` (Eq. :ref:`12.6-56 <eq-12.6-56>`)
can be replaced using Eqs. :ref:`12.6-57 <eq-12.6-57>` and :ref:`12.6-62 <eq-12.6-62>`.
In addition, the advanced time densities in Eq. :ref:`12.6-56 <eq-12.6-56>` can be
linearized through Eq. :ref:`12.6-22 <eq-12.6-22>` and written in terms of pressure
changes form Eq. :ref:`12.6-37 <eq-12.6-37>`. This
gives :math:`v\frac{\partial\rho}{\partial z}` as
(12.6-65)
.. _eq-12.6-65:
.. math:: v\frac{\partial\rho}{\partial z} = \frac{1}{2}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) + \Delta p_{i}(J1)\frac{dv^{i}}{dp_{J1}} \right)
.. math:: \times \frac{\rho_{i}(J1 + 1) + \Delta p(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} - \rho_{1}^{i}(J1) - \Delta p^{i}(J1)\frac{d\rho^{i}}{dp_{J1}}}{z(J1+1) - z_2^i(J1)}
Now, all that remains is to express the interface position
:math:`z_{2}^{i}` as a function of the change in pressure; this will be
done in the next subsection.
.. _section-12.6.2.2.3:
Formulation of :math:`z_2^i` as a Function of :math:`\Delta p`
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
The function height can be written as
(12.6-66)
.. _eq-12.6-66:
.. math:: z_{2}^{i}(J1) = z_{1}^{i}(J1) + \frac{\Delta t}{2}\left( v_{1}^{i}(J1) + v_{2}^{i}(J1) \right)
which, substituting Eqs. :ref:`12.6-57 <eq-12.6-57>` and :ref:`12.6-62 <eq-12.6-62>`
for :math:`v_{2}^{i}`,
becomes
(12.6-67)
.. _eq-12.6-67:
.. math:: z_{2}^{i}(J1) = z_{1}^{i}(J1) + \frac{\Delta t}{2}v_{1}^{i}(J1) + \left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J) + \Delta p^{i}(J1)\frac{dv^{i}}{dp_{J1}} \right)
.. math:: = z_{1}^{i}(J1) + v_{1}^{i}(J1)\Delta t + \frac{\Delta t}{2}v_{0}^{i}(J1) + \Delta p^{i}(J1)\frac{\Delta t}{2}\frac{dv^{i}}{dp_{J1}}
Define
(12.6-68)
.. _eq-12.6-68:
.. math:: z_{0}^{i}(J1) = z_{1}^{i}(J1) + v_{1}^{i}(J1)\Delta t + \frac{\Delta t}{2}\Delta v_{0}^{i}(J1)
as the interface height with no pressure change at the interface and
(12.6-69)
.. _eq-12.6-69:
.. math:: \frac{dz^{i}}{dp_{J1}} = \frac{1}{2}\Delta t\frac{dv^{i}}{dp_{J1}}
as the change in interface height with respect to the change in
interface pressure. Then :math:`z_{2}^{i}` is the linear function
(12.6-70)
.. _eq-12.6-70:
.. math:: z_{2}^{i}(J1) = z_{0}^{i}(J1) + \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}
Using Eq. :ref:`12.6-70 <eq-12.6-70>` in Eq. :ref:`12.6-65 <eq-12.6-65>`,
:math:`v\frac{\partial\rho}{\partial z}` can be expressed as a function of
only :math:`\Delta p`,
(12.6-71)
.. _eq-12.6-71:
.. math:: v\frac{\partial \rho}{\partial z} = \frac{1}{2}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) + \Delta p^{i}(J1)\frac{dv^{i}}{dp_{J1}} \right)
.. math:: \times \frac{\rho_{1}(J1 + 1) + \Delta p(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} - \rho_1^i(J1) - \Delta p^i\frac{d\rho^i}{dp_{J1}}}{z(J1 + 1) - z_{0}^{i}(J1) - \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}}
.. _section-12.6.2.2.4:
Final Differenced Form of the Continuity Equation at the Lower Interface as a Function of Only :math:`\Delta p` and :math:`\Delta W`
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
With :math:`\frac{d\rho}{dt}` in the interface continuity
equation, Eq. :ref:`12.6-53 <eq-12.6-53>`, written in the form in Eq. :ref:`12.6-55 <eq-12.6-55>` and
:math:`v\frac{\partial\rho}{\partial z}` given by Eq. :ref:`12.6-71 <eq-12.6-71>` and all
other terms in Eq. :ref:`12.6-53 <eq-12.6-53>` treated as in :numref:`section-12.6.2.1`,
the differenced form of the continuity equation at the lower interface is, after
multiplication by
:math:`\left( z(J1 + 1) - z_{2}^{i}(J1) \right)` and
substitution of Eq. :ref:`12.6-70 <eq-12.6-70>` for :math:`z_{2}^{i}`,
(12.6-72)
.. _eq-12.6-72:
.. math:: \Delta t\left( W_{1}(J1 + 1) - W_{1}^{i}(J1) + \Delta W(J1 + 1) - \Delta W^{i}(J1) \right)
.. math:: + \frac{A_{c}(J1)}{2}\left( \Delta \rho(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} + \Delta \rho^{i}(J1)\frac{d\rho^{i}}{dp_{J1}} \right)\bigg( z(J1 + 1) - z_{0}^{i}(J1)
.. math:: - \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}} \bigg) - \frac{A_{c}(J1)\Delta t}{2}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) + \Delta p^{i}(J1)\frac{dv^{i}}{dp_{J1}} \right)
.. math:: \times \left( \rho_{1}(J1 + 1) + \Delta p(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} - \rho_{1}^{i}(J1) - \Delta p^{i}(J1)\frac{d\rho^{i}}{dp_{J1}} \right)
.. math:: = A_{c}(J1)\Delta t\frac{Q_{0}(J1)}{\lambda_{1}(J1)}\bigg[ 1 + \left( \Delta p^{i}(J1)\frac{dT^{i}}{dp_{J1}} + \Delta p(J1 + 1)\frac{dT}{dp_{J1 + 1}} \right)\frac{1}{2Q_{0}(J1)}\frac{dQ_{0}}{dT_{J1}}
.. math:: - \left( \Delta p^{i}(J1) + \Delta p(J1 + 1)\frac{1}{2\lambda_{1}(J1)}\frac{d\lambda}{dp_{J1}} \right) \bigg]\bigg( z(J1 + 1) - z_{0}^{i}(J1)
.. math:: + \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}} \bigg)
Setting
:math:`\Delta z_{0}(J1) = z(J1 + 1) - z_{0}^{i}(J1)`
and eliminating second-order terms gives a differenced interface
continuity equation which is linear in the pressure and mass flow rate
changes:
(12.6-73)
.. _eq-12.6-73:
.. math:: \Delta t\left( W_{1}(J + 1) - W_{1}^{i}(J1) + \Delta W(J1 + 1) - \Delta W^{i}(J1) \right) \\
.. math:: + \frac{A_{c}(J1)}{2}\Delta z_{0}(J1)\left( \Delta p(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} + \Delta p^{i}(J1)\frac{d\rho^{i}}{dp_{J1}} \right) \\
.. math:: - \frac{A_{c}(J1)\Delta t}{2} \bigg\{ \left\lbrack v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right\rbrack\left\lbrack \rho_{1}(J1 + 1) - \rho_{1}^{i}(J1) \right\rbrack \\
.. math:: + \left[v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right] \left\lbrack \Delta p(J1 + 1)\frac{d\rho}{dp_{J1 + 1}} - \Delta p^{i}(J1)\frac{d\rho^{i}}{dp_{J1}} \right\rbrack \\
.. math:: + \Delta p^{i}(J1)\frac{dv^{i}}{dp_{J1}}\left\lbrack \rho_{1}(J1 + 1) - \rho_{1}^{i}(J1) \right\rbrack \bigg\} - A_{c}(J1)\Delta t\frac{Q_{0}(J1)}{\lambda_{1}(J1)} \bigg\{ \Delta z_{0}(J1) \\
.. math:: - \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}} + \Delta z_{0}(J)\left\lbrack \Delta p^{i}(J1)\frac{dT^{i}}{dp_{J1}} + \Delta p(J1 + 1)\frac{dT}{dp_{J1 + 1}} \right\rbrack\frac{1}{2Q_{0}(J1)}\frac{dQ_{0}}{dT_{J1}} \\
.. math:: - \Delta z_{0}(J1)\frac{1}{2\lambda_{1}(J1)}\left\lbrack \Delta p^{i}(J1) + \Delta p(J1 + 1) \right\rbrack\frac{d\lambda}{dp_{J1}} \bigg\} = 0
Equation :ref:`12.6-73 <eq-12.6-73>` can be simplified by defining the coefficients
(12.6-74)
.. _eq-12.6-74:
.. math:: c_{1,J1} = \frac{A_{c}(J1)}{2}\bigg\{\Delta z_{0}(J1)\frac{d\rho^{i}}{dp_{J1}} + \Delta t\bigg[ \left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right)\frac{d\rho^{i}}{dp_{J1}}
.. math:: - \left( \rho_{1}(J1 + 1) - \rho_{1}^{i}(J1) \right)\frac{dv^{i}}{dp_{J}} \bigg] + \Delta t\bigg[ \frac{2Q_{0}(J1)}{\lambda_{1}(J1)}\frac{dz^{i}}{dp_{J1}}
.. math:: - \frac{\Delta z_{0}(J1)}{\lambda_{1}(J1)}\frac{dQ_{0}}{dT_{J1}}\frac{dT^{i}}{dp_{J1}} + \frac{Q_{0}(J1)\Delta z_{0}(J1)}{\lambda_{1}^{2}(J1)}\frac{d\lambda}{dp_{J1}} \bigg]\bigg\}
(12.6-75)
.. _eq-12.6-75:
.. math:: c_{2,J1} = - \Delta t
(12.6-76)
.. _eq-12.6-76:
.. math:: c_{3,J1} = \frac{A_{c}(J1)}{2}\bigg\{\Delta z_{0}(J1)\frac{d\rho}{dp_{J1 + 1}} - \Delta t \left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right)\frac{d\rho}{dp_{J1 + 1}}
.. math:: - \Delta t\frac{\Delta z_{0}(J1)}{\lambda_{1}(J1)}\left\lbrack \frac{dQ_{0}}{dT_{J1}}\frac{dT}{dp_{J1 + 1}} - \frac{Q_{0}(J1)}{\lambda_{1}(J1)}\frac{d\lambda}{dp_{J1}} \right\rbrack \bigg\}
(12.6-77)
.. _eq-12.6-77:
.. math:: c_{4,J1} = \Delta t
and the source term
(12.6-78)
.. _eq-12.6-78:
.. math:: h_{J1} = \Delta t \bigg\{ W_{1}^{i}(J1) - W_{1}(J1 + 1) + A_{c}(J1) \bigg[\frac{1}{2}\big( v_{1}^{i}(J1)
.. math:: + \Delta v_{0}^{i}(J1) \big)\left( \rho_{1}(J1 + 1) - \rho_{1}^{i}(J1) \right) + \frac{Q_{0}(J1)\Delta z_{0}(J1)}{\lambda_{1}(J1)} \bigg]\bigg\}
Then the differenced continuity equation at the lower interface can be
written as
(12.6-79)
.. _eq-12.6-79:
.. math:: c_{1,J1}\Delta p^{i}(J1) + c_{2,J1}\Delta W^{i}(J1) + c_{3,J1}\Delta p(J1 + 1) + c_{4,J1}\Delta W(J1 + 1) = h_{J1}
Note that with
:math:`v_{1}^{i}`, :math:`\Delta v_{0}^{i}`, :math:`\frac{dv^{i}}{dp}`, and :math:`\frac{dz^{i}}{dp}` set
to zero, the expressions for the coefficients :math:`c_{1,J1}`, and the
source term :math:`h_{J1}` reduce to the forms given at the end of
:numref:`section-12.6.2.1` for the continuity equation for a mesh segment which
lies entirely within the bubble.
.. _section-12.6.2.2.5:
Differenced Form of the Continuity Equation at the Upper Interface
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
The derivation of the differenced continuity equation at the upper
interface is nearly identical to that for the lower interface, and so
just a brief explanation will be given here. The terms in the interface
continuity equation of Eq. :ref:`12.6-53 <eq-12.6-53>` are now expressed as follows:
The time derivative is
(12.6-80)
.. _eq-12.6-80:
.. math:: \frac{d\rho}{dt} = \frac{\rho_{2}^{i}(J2) - \rho_{1}^{i}(J2) + \rho_{2}(J2 - 1) - \rho_{1}(J2 - 1)}{2\Delta t}
.. math:: = \frac{1}{2\Delta t}\left( \Delta p^{i}(J2)\frac{d\rho^{i}}{dp_{J2}} + \Delta p(J2 - 1)\frac{d\rho}{dp_{J2 - 1}} \right)
The :math:`v\frac{\partial\rho}{\partial z}` term is
(12.6-81)
.. _eq-12.6-81:
.. math:: v\frac{\partial\rho}{\partial z} = \frac{v_{1}^{i}(J2) + \Delta v^{i}(J2)}{2}\frac{\rho_{2}^{i}(J2) - \rho_{2}(J2 - 1)}{z_{2}^{i}(J2) - z(J2 - 1)}
where
(12.6-82)
.. _eq-12.6-82:
.. math:: \Delta v^{i}(J2) = \frac{\Delta W_{l}\left( K + 1 \right)}{\rho_{l}^{i}(J2)A_{c}(J2 - 1)}f_{rwb}(K + 1)
and
(12.6-83)
.. _eq-12.6-83:
.. math:: \Delta W_{l}\left( K + 1 \right) = \frac{\left( AA_{0}\left( K + 1 \right) + \theta_{2}\left( \Delta p^{i}(J2) - \Delta p^{i}\left( K + 1 \right) \right) \right)\Delta t}{I_{1}\left( K + 1 \right) + \theta_{2}\Delta I_{1}\left( K + 1 \right) + BB_{0}\left( K + 1 \right)\theta_{2}\Delta t}
The rewetting factor at the bottom of liquid slug :math:`K + 1`, :math:`f_{rwb}\left(K + 1\right)`,
is computed as in Eq. :ref:`12.6-59 <eq-12.6-59>`. The expression for
:math:`\Delta v_{i}(J2)` then reduces to
(12.6-84)
.. _eq-12.6-84:
.. math:: \Delta v^{i}(J2) = \Delta v_{0}^{i}(J2) + \Delta p^{i}(J2)\frac{dv^{i}}{dp_{J2}}
with the definitions of
:math:`\Delta v_{0}^{i}(J2)` and :math:`\frac{dv^{i}}{dp_{J2}}` being
set as were those for
:math:`\Delta v_{0}^{i}(J1)` and :math:`\frac{dv^{i}}{dp_{J1}}`.
Also,
(12.6-85)
.. _eq-12.6-85:
.. math:: z_{2}^{i}(J2) = z_{0}^{i}(J2) + \Delta p^{i}(J2)\frac{dz^{i}}{dp_{J2}}
can be derived in exactly the same fashion as was Eq. :ref:`12.6-70 <eq-12.6-70>`.
Therefore,
(12.6-86)
.. _eq-12.6-86:
.. math:: v\frac{\partial\rho}{\partial z} = \frac{1}{2}\left( v_{1}^{i}(J2) + \Delta v_{0}^{i}(J2) + \Delta p^{i}(J2)\frac{dv^{i}}{dp_{J2}} \right)
.. math:: \times \left( \frac{\rho_{1}^{i}(J2) + \Delta p^{i}(J2)\frac{d\rho^{i}}{dp_{J2}} - \rho_{1}^{i}(J2 - 1) - \Delta p(J2 - 1)\frac{d\rho}{dp_{J2 - 1}}}{z_{0}^{i}(J2) + \Delta p^{i}(J2)\frac{dz^{i}}{dp_{J2}} - z(J2 - 1)} \right)
Substituting Eqs. :ref:`12.6-80 <eq-12.6-80>` and :ref:`12.6-86 <eq-12.6-86>`
into Eq. :ref:`12.6-53 <eq-12.6-53>` and expressing
:math:`\frac{\partial W}{\partial z}` and :math:`A_{c}\frac{Q}{\lambda}` as
in :numref:`section-12.6.2.1` gives the differenced upper interface continuity
equation
(12.6-87)
.. _eq-12.6-87:
.. math:: c_{1,J2 - 1}\Delta p(J2 - 1) + c_{2,J2 - 1}\Delta W(J2 - 1) + c_{3,J2 - 1}\Delta p^{i}(J2)
.. math:: + c_{4,J2 - 1}\Delta W^{i}(J2) = h_{J2 - 1}
where
(12.6-88)
.. _eq-12.6-88:
.. math:: c_{1,J2 - 1} = \frac{A_{c}(J2 - 1)}{2}\bigg\{\Delta z_{0}(J2 - 1)\frac{d\rho}{dp_{J2 - 1}} + \Delta t\bigg[\left( v_{1}^{i}(J2) + \Delta v_{0}^{i}(J2) \right)\frac{d\rho}{dp_{J2 - 1}}
.. math:: + \frac{\Delta z_{0}(J2 - 1)}{\lambda_{1}(J2 - 1)}\left(\frac{dQ_{0}}{dT_{J2 - 1}}\frac{dT}{dp_{J2 - 1}} - \frac{Q_{0}(J2 - 1)}{\lambda(J2 - 1)}\frac{d\lambda}{dp_{J2}} \right) \bigg] \bigg\}
(12.6-89)
.. _eq-12.6-89:
.. math:: c_{2,J2 - 1} = - \Delta t
(12.6-90)
.. _eq-12.6-90:
.. math:: c_{3,J2 - 1} = \frac{A_{c}(J2 - 1)}{2}\Bigg\{\Delta z_{0}(J2 - 1)\frac{d\rho^{i}}{dp_{J2}} - \Delta t\Bigg[\left( v_{1}^{i}(J2) + \Delta v_{0}^{i}(J2) \right)\frac{d\rho^{i}}{dp_{J2}}
.. math:: + \left\lbrack \rho_{1}^{i}(J2) - \rho_{1}(J2 - 1) \right\rbrack\frac{dv^{i}}{dp_{J2}}\Bigg] - \Bigg[\Delta t\frac{2Q_{0}(J2 - 1)}{\lambda_{1}(J2 - 1)}\frac{dz^{i}}{dp_{J2}}
.. math:: + \frac{\Delta z_{0}(J2 - 1)}{\lambda_{1}(J2 - 1)}\left(\frac{dQ_{0}}{dT_{J2 - 1}}\frac{dT^{i}}{dp_{J2}} - \frac{Q_{0}(J2 - 1)}{\lambda_{1}(J2 - 1)}\frac{d\lambda}{dp_{J2}}\right)\Bigg]\Bigg\}
(12.6-91)
.. _eq-12.6-91:
.. math:: c_{4,J2 - 1} = \Delta t
(12.6-92)
.. _eq-12.6-92:
.. math:: h_{J2 - 1} = \Delta t\bigg\{ W_{1}(J2 - 1) - W_{1}^{i}(J2) + A_{c}(J2 - 1)\bigg[\frac{1}{2}( v_{1}^{i}(J2)
.. math:: + \Delta v_{0}^{i}(J2))\left( \rho_{1}^{i}(J2) - \rho_{1}(J2 - 1) \right) + \frac{Q_{0}(J2 - 1)\Delta z_{0}(J2 - 1)}{\lambda_{1}(J2 - 1)}\bigg]\bigg\}
and
(12.6-93)
.. _eq-12.6-93:
.. math:: \Delta z_{0}(J2 - 1) = z_{1}^{i}(J2) - z(J2 - 1)
.. _section-12.6.2.3:
Finite Differencing of the Momentum Equation for a Mesh Segment Contained Entirely in One Bubble
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Now that differenced forms of the continuity equation both within the
bubble and at liquid-vapor interfaces have been derived, a similar
process must be applied in the momentum equation. The differencing of
the momentum equation is accomplished in essentially the same way as
that of the continuity equation. Once again, bubble interface will be
ignored for the sake of simplicity until the equations for an interior
node have been developed.
.. _section-12.6.2.3.1:
Differencing of the Momentum Equation
'''''''''''''''''''''''''''''''''''''
The differenced forms of the terms in the momentum equation, Eq. :ref:`12.6-2 <eq-12.6-2>`,
are a follows. The channel pressure drop is given by
(12.6-94)
.. _eq-12.6-94:
.. math:: \frac{\partial p}{\partial z}\Big|_{J + \frac{1}{2}} = \frac{p_{2}(J + 1) - p_{2}(J)}{z(J + 1) - z(J)}
The momentum convection takes the form
(12.6-95)
.. _eq-12.6-95:
.. math:: \left. \frac{1}{A_{c}}\frac{\partial}{\partial z}\left( \frac{W^{2}}{\rho A_{c}} \right) \right|_{J + \frac{1}{2}} = \frac{1}{A_{c}^{2}(J)\left( z(J + 1) - z(J) \right)}\left( \frac{W_{2}^{2}(J + 1)}{\rho_{2}(J + 1)} - \frac{W_{2}^{2}(J)}{\rho_{2}(J)} \right)
(12.6-96)
.. _eq-12.6-96:
The condensation term is expressed in terms of the contributions from
cladding and structure as
.. math:: F_{c}(J) = F_{ce}(J) + F_{cs}(J)
where the functions for :math:`F_{ce}` and :math:`F_{cs}` are linearly approximated as
follows:
(12.6-97a)
.. _eq-12.6-97a:
.. math:: \text{if }Q_{e_{2}}(J) \geq 0,\ F_{ce}(J) = 0
(12.6-97b)
.. _eq-12.6-97b:
.. math:: \text{if }Q_{e_{2}}(J) < 0,\ F_{ce}(J) = \frac{Q_{e_{2}}(J)}{2\lambda_{2}(J)A_{c}(J)}\left( \frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_{2}(J + 1)} \right)
and
(12.6-98a)
.. _eq-12.6-98a:
.. math:: \text{if }Q_{s_{2}}(J) \geq 0,\ F_{cs}(J) = 0
(12.6-98b)
.. _eq-12.6-98b:
.. math:: \text{if }Q_{s_{2}}(J) \leq 0,\ F_{cs}(J) = \frac{Q_{s_{2}}(J)}{2\lambda_{2}(J)A_{c}(J)}\left( \frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_2(J + 1)} \right)
:math:`Q_{e_2}(J)` and :math:`Q_{s_2}(J)` are defined in Eq.
:ref:`12.6-19 <eq-12.6-19>`. Now, linearize :math:`Q_{e_2}` and :math:`Q_{s_2}` by
defining
(12.6-99)
.. _eq-12.6-99:
.. math:: Q_{oe}(J) = \gamma(J) \frac{T_{e_{2}}(J) - \overline{T}_{1}\left(1\right)}{R_{ec_2}(J)}
(12.6-100)
.. _eq-12.6-100:
.. math:: \frac{dQ_{oe}}{dT_{J}} = -\frac{\gamma(J)}{R_{ec_2}(J)}
(12.6-101)
.. _eq-12.6-101:
.. math:: Q_{os}(J) = \gamma(J)\gamma_{2}(J)\frac{T_{s_{2}}(J) - {\overline{T}}_{1}(J)}{R_{sc_{2}}(J)}
(12.6-102)
.. _eq-12.6-102:
.. math:: \frac{dQ_{os}}{dT_{J}} = - \frac{\gamma(J)\gamma_{2}(J)}{R_{sc_{2}}(J)}
These expressions are related to :math:`Q_{o}(J)`, defined in :numref:`section-12.6.2.1.2`, by
(12.6-103)
.. _eq-12.6-103:
.. math:: Q_{o}(J) = Q_{oe}(J) + Q_{os}(J)
and
(12.6-104)
.. _eq-12.6-104:
.. math:: \frac{dQ_{o}}{dT_{J}} = \frac{dQ_{oe}}{dT_{J}} + \frac{dQ_{os}}{dT_{J}}
Substituting these functions into the above formulations of :math:`F_{ce}` and :math:`F_{cs}` gives
(12.6-105a)
.. _eq-12.6-105a:
.. math:: \text{if }Q_{e_{2}}(J) \geq 0,\ F_{ce}(J) = 0
(12.6-105b)
.. _eq-12.6-105b:
.. math::
\text{if }Q_{e_{2}}(J) < 0,\ F_{ce}(J) = \frac{1}{2\lambda_{2}(J)A_{c}(J)}\left( \frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_{2}(J + 1)} \right) \\
\times \left[ Q_{oe}(J) + \frac{1}{2}\frac{dQ_{oe}}{dT_{J}}\left( \Delta T(J) + \Delta T(J + 1) \right) \right]
(12.6-106a)
.. _eq-12.6-106a:
.. math:: \text{if}\ Q_{s_{2}}(J) \geq 0,\ F_{cs}(J) = 0
(12.6-106b)
.. _eq-12.6-106b:
.. math::
\text{if }Q_{s_{2}}(J) < 0,\ F_{cs}(J) = \frac{1}{2\lambda_{2}(J)A_{c}(J)}\left( \frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_{2}(J + 1)} \right) \\
\times \left[ Q_{os}(J) + \frac{1}{2}\frac{dQ_{os}}{dT_J}\left( \Delta T(J) + \Delta T(J+1) \right) \right]
If the definitions
(12.6-107a-b)
.. _eq-12.6-107a-b:
.. math::
F_{oe}(J) = \Bigg\{ \begin{matrix}
0,\ Q_{e_{2}} \geq 0\ \\
Q_{oe}(J),\ Q_{e_{2}}(J) < 0\ \\
\end{matrix}
(12.6-108a-b)
.. _eq-12.6-108a-b:
.. math::
F_{os}(J) = \Bigg\{ \begin{matrix}
0,\ Q_{s_{2}} \geq 0, \\
Q_{os}(J),\ Q_{s_{2}}(J) < 0 \\
\end{matrix}
(12.6-109)
.. _eq-12.6-109:
.. math:: F_{o}(J) = F_{oe}(J) + F_{os}(J)
are used in the expressions for :math:`F_{ce}` and :math:`F_{cs}`, the condensation
term becomes
(12.6-110)
.. _eq-12.6-110:
.. math:: F_{c}(J) = \frac{F_{o}(J)}{2A_{c}(J)\lambda_{2}(J)}\bigg[ 1 + \left( \frac{F_{oc}(J)}{2Q_{oe}(J)F_{o}(J)}\frac{dQ_{oe}}{dT_{J}} + \frac{F_{os}(J)}{2Q_{os}(J)F_{o}(J)}\frac{dQ_{os}}{dT_{J}} \right)
.. math:: \times(\Delta t(J) + \Delta T(J + 1))\bigg]\left(\frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_{2}(J + 1)} \right)
This expression can be simplified by setting
(12.6-111a-b)
.. _eq-12.6-111a-b:
.. math::
\frac{dF_{oe}}{dT_{J}} = \Bigg\{ \begin{matrix}
0,\ Q_{e_{2}}(J) \geq 0 \\
\frac{dQ_{oe}}{dT_{J}},Q_{e_{2}}(J) < 0 \\
\end{matrix}
which can be stated more compactly as
(12.6-112)
.. _eq-12.6-112:
.. math:: \frac{dF_{oe}}{dT_{J}} = \frac{F_{oe}(J)}{Q_{oe}(J)}\frac{dQ_{oe}}{dT_{J}}
Similarly,
(12.6-113)
.. _eq-12.6-113:
.. math:: \frac{dF_{os}}{dT_{J}} = \frac{F_{os}(J)}{Q_{os}(J)}\frac{dQ_{os}}{dT_{J}}
so that the definition
(12.6-114)
.. _eq-12.6-114:
.. math:: \frac{dF_{o}}{dT_{J}} = \frac{dF_{oe}}{dT_{J}} + \frac{dF_{os}}{dT_{J}}
can be made. This gives
(12.6-115)
.. _eq-12.6-115:
.. math:: F_{c}(J) = \frac{F_{o}(J)}{2A_{c}(J)\lambda_{2}(J)}\left\lbrack 1 + \frac{dF_{o}}{2F_{o}(J)}\frac{dF_{o}}{dT_{J}}\left( \Delta T(J) + \Delta T(J + 1) \right) \right\rbrack
.. math:: \times\left( \frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_{2}(J + 1)} \right)
as the differenced form of the condensation term.
The change in mass flow rate with time is
(12.6-116)
.. _eq-12.6-116:
.. math:: \frac{\partial W}{\partial t} = \frac{W_{2}(J) + W_{2}(J + 1) - W_{1}(J) - W_{1}(J + 1)}{2\Delta t}
The orifice pressure drop becomes
(12.6-117)
.. _eq-12.6-117:
.. math:: \frac{W\left| W \right|}{2\rho_{v}A_{c}^{2}}\frac{\partial k_{or}}{\partial z} = \frac{1}{4A_{c}^{2}(J)}\left( \frac{W_{2}(J)\left| W_{2}(J) \right|}{\rho_{2}(J)} + \frac{W_{2}(J + 1)\left| W_{2}(J + 1) \right|}{\rho_{2}(J + 1)} \right)
.. math:: \times\frac{k_{or}(J)}{z(J + 1) - z(J)}
Finally, if the friction factor
(12.6-118)
.. _eq-12.6-118:
.. math:: f_{v} = A_{fr}\left( \frac{D_{h}W}{\mu A_{c}} \right)^{b_{fr}}
is considered in the differencing, the friction term may be written as
(12.6-119)
.. _eq-12.6-119:
.. math:: \frac{f_{v}W\left| W \right|}{2\rho A_{c}^{2}D_{h}} = \frac{A_{fr}W\left| W \right|^{1 + b_{fr}}}{2\rho\mu^{b_{fr}}D_{h}^{1 - b_{fr}}A_{c}^{2 + b_{fr}}}
.. math:: = \frac{A_{fr}}{4D_h^{1-b_{fr}}(J)A_c^{2+b_{fr}}(J)}\left[ \frac{W_{2}(J+1)\left|W_{2}(J+1)\right|^{1+b_{fr}}}{\rho_{2}(J+1)\mu_2^{b_{fr}}(J+1)} + \frac{W_{2}(J)\left|W_{2}(J)\right|^{1+b_{fr}}}{\rho_{2}(J)\mu_2^{b_{fr}}(J)}\right]
Using Eqs. :ref:`12.6-94 <eq-12.6-94>`, :ref:`12.6-95 <eq-12.6-95>`, :ref:`12.6-115 <eq-12.6-115>`,
:ref:`12.6-116 <eq-12.6-116>`, :ref:`12.6-117 <eq-12.6-117>`, and :ref:`12.6-119 <eq-12.6-119>`
results in a differenced momentum equation as follows:
(12.6-120)
.. _eq-12.6-120:
.. math:: \frac{p_{2}(J + 1) - p_{2}(J)}{z(J + 1) - z(J)} + \frac{1}{A_{c}^{2}(J)\left( z(J + 1) - z(J) \right)}\left( \frac{W_{2}^{2}(J + 1)}{\rho_{2}(J + 1)} - \frac{W_{2}^{2}(J)}{\rho_{2}(J)} \right)
.. math:: + \frac{A_{fr}}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\left\lbrack \frac{W_{2}(J + 1)\left| W_{2}(J + 1) \right|^{1 + b_{fr}}}{\rho_{2}(J + 1)\mu_{2}^{b_{fr}}(J + 1)} + \frac{W_{2}(J)\left| W_{2}(J) \right|^{1 + b_{fr}}}{\rho_{2}(J)\mu_{2}^{b_{fr}}(J)} \right\rbrack
.. math:: - \frac{F_{o}(J)}{2A_{c}(J)\lambda_{2}(J)}\left\lbrack 1 + \frac{1}{2F_{o}(J)}\frac{dF_{o}}{dT_{J}}\left( \Delta T(J) + \Delta T(J + 1) \right) \right\rbrack\left( \frac{W_{2}(J)}{\rho_{2}(J)} + \frac{W_{2}(J + 1)}{\rho_{2}(J + 1)} \right)
.. math:: + \frac{W_{2}(J) + W_{2}(J + 1) - W_{1}(J) - W_{1}(J + 1)}{2A_{c}(J)\Delta t} + \frac{1}{4A_{c}^{2}(J)}\Big( \frac{W_{2}(J)\left| W_{2}(J) \right|}{\rho_{2}(J)}
.. math:: + \frac{W_{2}(J + 1)|W_{2}(J + 1)|}{\rho_{2}(J + 1)} \Big)\frac{k_{or}(J)}{z(J + 1) - z(J)} = 0
The equation must now be linearized and reduced to a linear function of
:math:`\Delta p` and :math:`\Delta W`; this procedure will be discussed in the next subsection.
.. _section-12.6.2.3.2:
Linearization of Advanced Time Quantities and Introduction of Saturation Conditions
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Equation :ref:`12.6-120 <eq-12.6-120>` is made linear by expressing the advanced time
variables in linearized form as in Eq. :ref:`12.6-21 <eq-12.6-21>`. This means that, in
addition to using the expressions in Eqs. :ref:`12.6-22 <eq-12.6-22>` through
:ref:`12.6-25 <eq-12.6-25>` for density, heat of vaporization, mass flow rate, and
temperature, the linearized form of Eq. :ref:`12.6-120 <eq-12.6-120>` must include
the following representations of pressure and viscosity:
(12.6-121)
.. _eq-12.6-121:
.. math:: p_{2}(J) = p_{1}(J) + \Delta p(J)
(12.6-122)
.. _eq-12.6-122:
.. math:: \mu_{2}(J) = \mu_{1}(J) + \Delta\mu(J)
The differenced momentum equation therefore becomes
(12.6-123)
.. _eq-12.6-123:
.. math:: \frac{p_{1}(J + 1) + \Delta p(J + 1) - p_{1}(J) - \Delta p(J)}{z(J + 1) - z(J)} + \frac{1}{A_{c}^{2}(J)\left( z(J + 1) - z(J) \right)}
.. math:: \times\left( \frac{\left( W_{1}(J + 1) + \Delta W(J + 1) \right)^{2}}{\rho_{1}(J + 1) + \Delta\rho(J + 1)} - \frac{\left( W_{1}(J) + \Delta W(J) \right)^{2}}{\rho_{1}(J) + \Delta\rho(J)} \right)
.. math:: + \frac{A_{fr}}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\Bigg( \frac{\left( W_{1}(J + 1) + \Delta W(J + 1) \right)\left| W_{1}(J + 1) + \Delta W(J + 1) \right|^{1 + b_{fr}}}{\left( \rho_{1}(J + 1) + \Delta\rho(J + 1) \right)\left( \mu_{1}(J + 1) + \Delta\mu(J + 1) \right)^{b_{fr}}}
.. math:: + \frac{\left( W_{1}(J) + \Delta W(J) \right) \left| W_{1}(J) + \Delta W(J) \right|^{1 + b_{fr}}}{(\rho_{1}(J) + \Delta\rho(J)\left( \mu_{1}(J) + \Delta\mu(J) \right)^{b_{fr}}} \Bigg) - \frac{F_{o}(J)}{2A_{c}(J){(\lambda}_{1}(J) + \Delta\lambda(J))}
.. math:: \times\left( 1 + \frac{1}{2F_{o}(J)}\frac{dF_{o}}{dT_{J}}\left( \Delta T(J) + \Delta T(J + 1) \right) \right)\left( \frac{W_{1}(J) + \Delta W(J)}{\rho_{1}(J) + \Delta\rho(J)} + \frac{W_{1}(J + 1) + \Delta W(J + 1)}{\rho_{1}(J + 1) + \Delta\rho(J + 1)} \right)
.. math:: + \frac{\Delta W(J) + \Delta W(J + 1)}{2A_{c}(J)\Delta t} + \frac{1}{4A_{c}^{2}(J)}\Bigg( \frac{\left(W_{1} + \Delta W(J)\right)\left| W_{1}(J) + \Delta W(J) \right|}{\rho_{1}(J) + \Delta\rho(J)}
.. math:: + \frac{\left(W_{1}(J + 1) + \Delta W(J + 1)\right)\left|W_{1}(J + 1) + \Delta W(J + 1)\right|}{\rho_{1}(J + 1) + \Delta\rho(J + 1)} \Bigg)\frac{k_{or}(J)}{z(J + 1) - z(J)} = 0
Note that the orifice coefficient :math:`k_{or}` and the hydraulic diameter
:math:`D_{h}` are, like the channel area :math:`A_{c}`, considered to be constant over
the time step. Now, Eq. :ref:`12.6-123 <eq-12.6-123>` can be expressed in terms of only two
independent variables, :math:`\Delta p` and :math:`\Delta W`, by imposing the restriction
that :math:`\rho`, :math:`T`, :math:`\lambda`, and :math:`\mu` be functions of :math:`p` only, i.e., by imposing
saturation conditions. This is accomplished by substituting Eqs. :ref:`12.6-36 <eq-12.6-36>`
through :ref:`12.6-38 <eq-12.6-38>`, as well as the expression
(12.6-124)
.. _eq-12.6-124:
.. math:: \Delta\mu(J) = \Delta p(J)\frac{d\mu}{dp_{J}}
into Eq. :ref:`12.6-123 <eq-12.6-123>`. The resulting equation, after multiplying by
:math:`\Delta t \left( z (J+1) - z (J) \right)` and letting
:math:`\Delta z(J) = z(J+1) - z(J)`, is the very long
expression
(12.6-125)
.. _eq-12.6-125:
.. math:: \left( p_{1}(J + 1) - p_{1}(J) \right)\Delta t + \Delta t\Delta p(J + 1) - \Delta t\Delta p(J)
.. math:: + \frac{\Delta t}{A_{c}^{2}(J)}\Bigg[ \frac{W_{1}^{2}(J + 1) + 2W_{1}(J + 1)\Delta W(J + 1) + \left( \Delta W(J + 1) \right)^{2}}{\rho_{1}(J + 1)}\left(1 - \frac{\Delta p(J + 1)}{\rho_{1}(J + 1)}\frac{d\rho}{dp_{J + 1}} \right)
.. math:: - \frac{W_{1}^{2}(J) + 2W_{1}(J)\Delta W(J) + \left( \Delta W(J) \right)^{2}}{\rho_{1}(J)}\left( 1 - \frac{\Delta p(J)}{\rho_{1}(J)}\frac{d\rho}{dp_{J}} \right) \Bigg]
.. math:: + \frac{A_{fr}\Delta t\Delta z(J)}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\Bigg[ \frac{W_{1}(J + 1)\left| W_{1}(J + 1) \right|^{1 + b_{fr}}\left( \frac{1 + \left( 2 + b_{fr} \right)\Delta W(J + 1)}{W_1(J+1)}\right)}{\rho_{1}(J + 1)\mu_{1}^{b_{fr}}(J + 1)}
.. math:: \times \left( 1 - \frac{\Delta p(J + 1)}{\rho_{1}(J + 1)}\frac{d\rho}{dp_{J + 1}} - b_{fr}\frac{\Delta p(J + 1)}{\mu_{1}(J + 1)}\frac{d\mu}{dp_{J + 1}} \right)
.. math:: + \frac{W_{1}(J)\left| W_{1}(J) \right|^{1 + b_{fr}}\left( 1 + \frac{\left( 2 + b_{fr} \right)\Delta W(J)}{W_{1}(J)} \right)}{\rho_{1}(J)\mu_{1}^{b_{fr}}(J)}\left( 1 - \frac{\Delta p(J)}{\rho_{1}(J)}\frac{d\rho}{dp_{J}} - b_{fr}\frac{\Delta p(J)}{\mu_{1}(J)}\frac{d\mu}{dp_{J}} \right) \Bigg]
.. math:: - \frac{F_{o}(J)\Delta t\Delta z(J)}{2A_{c}(J)\lambda_{1}(J)} \left( 1 - \frac{\Delta p(J) + \Delta p(J + 1)}{2\lambda_{1}(J)}\frac{d\lambda}{dp_{J}} \right)
.. math:: \times\left\lbrack 1 + \frac{1}{2F_{o}(J)}\frac{dF_{o}}{dT_{J}}\left( \Delta p(J)\frac{dT}{dp_{J}} + \Delta p(J + 1)\frac{dT}{dp_{J + 1}} \right) \right\rbrack
.. math:: \times\Bigg[ \frac{W_{1}(J + 1)}{\rho_{1}(J)} \left(1 + \frac{\Delta W (J)}{W_{1}(J)}\right)\left(1 - \frac{\Delta p (J)}{\rho_{1} (J)} \frac{d\rho}{dp_{J}}\right)
.. math:: + \frac{W_{1}(J + 1)}{\rho_{1}(J + 1)}\left( 1 + \frac{\Delta W(J + 1)}{W_{1}(J + 1)} \right)\left(1 - \frac{\Delta p (J+1)}{\rho_{1} (J+1)}\frac{d\rho}{dp_{J+1}} \right) \Bigg]
.. math:: + \frac{\Delta z(J)}{2A_{c}(J)}\Delta W(J) + \frac{\Delta z(J)}{2A_{c}(J)}\Delta W(J + 1)
.. math:: + \frac{\Delta t}{4A_{c}^{2}(J)}\Bigg[ \frac{W_{1}(J)\left| W_{1}(J) \right|\left( 1 + \frac{2\Delta W(J)}{W_{1}(J)} \right)}{\rho_{1}(J)}\left( 1 - \frac{\Delta p(J)}{\rho_{1}(J)}\frac{d\rho}{dp_{J}} \right)
.. math:: + \frac{W_{1}(J + 1)\left| W_{1}(J + 1) \right|\left( 1 + \frac{2\Delta W(J + 1)}{W_{1}(J + 1)} \right)}{\rho_{1}(J + 1)}\left( 1 - \frac{\Delta p(J + 1)}{\rho_{1}(J + 1)}\frac{d\rho}{dp_{J + 1}} \right) \Bigg]
.. math:: \times\,k_{or}(J) = 0
where it is assumed that
:math:`\Delta \mu \ll \mu_{1}` and :math:`\Delta \rho \ll \rho_{1}`, so that
(12.6-126)
.. _eq-12.6-126:
.. math:: \frac{1}{\rho_{1}(J) + \Delta\rho(J)} = \frac{1}{\rho_{1}(J)}\left( 1 - \frac{\Delta\rho(J)}{\rho_{1}(J)} \right)
and
(12.6-127)
.. _eq-12.6-127:
.. math:: \frac{1}{\left( \mu_{1}(J) + \Delta\mu(J) \right)^{b_{fr}}} = \frac{1}{\mu_{1}^{b_{fr}}(J)}\left( 1 - b_{fr}\frac{\Delta\mu(J)}{\mu_{1}(J)} \right)
Eliminating second-order terms produces the final linearized form of the
differenced momentum equation:
(12.6-128)
.. _eq-12.6-128:
.. math:: b_{1,J}\Delta p(J) + b_{2,J}\Delta W(J) + b_{3,J}\Delta p(J + 1) + b_{4,J}\Delta W(J + 1) = g_{J}
with
(12.6-129)
.. _eq-12.6-129:
.. math:: b_{1,J} = - \Delta t + \frac{\Delta t}{A_{c}^{2}(J)}\frac{W_{1}^{2}(J)}{\rho_{1}^{2}(J)}\frac{d\rho}{dp_{J}} - \frac{A_{fr}\Delta t\Delta z(J)}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\frac{W_{1}(J)\left| W_{1}(J) \right|^{1 + b_{fr}}}{\rho_{1}(J)\mu_{1}^{b_{fr}}(J)}
.. math:: \times \left( \frac{1}{\rho_{1}(J)}\frac{d\rho}{dp_{J}} + \frac{b_{fr}}{\mu_{1}(J)}\frac{d\mu}{dp_{J}} \right)
.. math:: + \frac{F_{o}(J)\Delta t\Delta z(J)}{4A_{c}(J)\lambda_{1}(J)}\Bigg[ \frac{1}{\lambda_{1}(J)}\frac{d\lambda}{dp_{J}}\left( \frac{W_{1}(J)}{\rho_{1}(J)} + \frac{W_{1}(J + 1)}{\rho_{1}(J + 1)} \right)
.. math:: - \frac{1}{F_{o}(J)}\frac{dF_{o}}{dT_{J}}\frac{dT}{dp_{J}}\left( \frac{W_{1}(J)}{\rho_{1}(J)} + \frac{W_{1}(J + 1)}{\rho_{1}(J + 1)} \right) + \frac{2W_{1}(J)}{\rho_{1}^{2}(J)}\frac{d\rho}{dp_J} \Bigg]
.. math:: - \frac{W_{1}(J)\left| W_{1}(J) \right|\Delta t}{4A_{c}^{2}(J)\rho_{1}^{2}(J)}\frac{d\rho}{dp_{J}}k_{or}(J)
(12.6-130)
.. _eq-12.6-130:
.. math:: b_{2,J} = - \frac{2W_{1}(J)\Delta t}{\rho_{1}(J)A_{c}^{2}(J)} + \frac{A_{fr}\Delta t\Delta z(J)}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\frac{\left( 2 + b_{fr} \right)\left| W_{1}(J) \right|^{1 + b_{fr}}}{\rho_{1}(J)\mu_{1}^{b_{fr}}(J)}
.. math:: - \frac{F_{o}(J)\Delta t\Delta z(J)}{2A_{c}(J)\lambda_{1}(J)\rho_{1}(J)} + \frac{\Delta z(J)}{2A_{c}(J)} + \frac{\Delta t\left| W_{1}(J) \right|}{2A_{c}^{2}(J)\rho_{1}(J)}k_{or}(J)
(12.6-131)
.. _eq-12.6-131:
.. math:: b_{3,J} = \Delta t - \frac{\Delta t}{A_{c}^{2}(J)}\frac{W_{1}^{2}(J+1)}{\rho_{1}^{2}(J+1)}\frac{d\rho}{dp_{J + 1}} - \frac{A_{fr}\Delta t\Delta z(J)}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\frac{W_{1}(J + 1)\left| W_{1}(J + 1) \right|^{1 + b_{fr}}}{\rho_{1}(J + 1)\mu_{1}^{b_{fr}}(J + 1)}
.. math:: \times \left( \frac{1}{\rho_{1}(J + 1)}\frac{d\rho}{dp_{J + 1}} + \frac{b_{fr}}{\mu_{1}(J + 1)} \right)\frac{d\mu}{dp_{J + 1}} + \frac{F_{o}(J)\Delta t\Delta z(J)}{4A_{c}(J)\lambda_{1}(J)}
.. math:: \times \Bigg[ \frac{1}{\lambda_{1}(J)}\frac{d\lambda}{dp_{J}}\left( \frac{W_{1}(J)}{\rho_{1}(J)} + \frac{W_{1}(J + 1)}{\rho_{1}(J + 1)} \right)
.. math:: - \frac{1}{F_{o}(J)}\frac{dF_{o}}{dT_{J}}\frac{dT}{dp_{J+1}}\left( \frac{W_{1}(J)}{\rho_{1}(J)} + \frac{W_{1}(J + 1)}{\rho_{1}(J + 1)} \right) - \frac{2W_{1}(J + 1)}{\rho_{1}^{2}(J + 1)}\frac{d\rho}{dp_{J+1}} \Bigg]
.. math:: - \frac{W_{1}(J + 1)\left| W_{1}(J + 1) \right|\Delta t}{4A_{c}^{2}(J)\rho_{1}^{2}(J + 1)}\frac{d\rho}{dp_{J + 1}}k_{or}(J)
(12.6-132)
.. _eq-12.6-132:
.. math:: b_{4,J} = - \frac{2W_{1}(J + 1)\Delta t}{\rho_{1}(J + 1)A_{c}^{2}(J)} + \frac{A_{fr}\Delta t\Delta z(J)}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\frac{\left( 2 + b_{fr} \right)\left| W_{1}(J + 1) \right|^{1 + b_{fr}}}{\rho_{1}(J + 1)\mu_{1}^{b_{fr}}(J + 1)}
.. math:: - \frac{F_{o}(J)\Delta t\Delta z(J)}{2A_{c}(J)\lambda_{1}(J)\rho_{1}(J + 1)} + \frac{\Delta z(J)}{2A_{c}(J)} + \frac{\Delta t\left| W_{1}(J + 1) \right|}{2A_{c}^{2}(J)\rho_{1}(J + 1)}k_{or}(J)
(12.6-133)
.. _eq-12.6-133:
.. math:: g_{J} = \left( p_{1}(J) - p_{1}(J + 1) \right)\Delta t - \frac{\Delta t}{A_{c}^{2}(J)}\left( \frac{W_{1}^{2}(J + 1)}{\rho_{1}(J + 1)} - \frac{W_{1}^{2}(J)}{\rho_{1}(J)} \right)
.. math:: - \frac{A_{fr}\Delta t\Delta z(J)}{4D_{h}^{1 - b_{fr}}(J)A_{c}^{2 + b_{fr}}(J)}\left( \frac{W_{1}(J + 1)\left| W_{1}(J + 1) \right|^{1 + b_{fr}}}{\rho_{1}(J + 1)\mu_{1}^{b_{fr}}(J + 1)} + \frac{W_{1}(J)\left| W_{1}(J) \right|^{1 + b_{fr}}}{\rho_{1}(J1)\mu_{1}^{b_{fr}}(J1)} \right)
.. math:: + \frac{F_{o}(J)\Delta t\Delta z(J)}{2A_{c}(J)\lambda_{1}(J)}\left( \frac{W_{1}(J)}{\rho_{1}(J)} + \frac{W_{1}(J + 1)}{\rho_{1}(J + 1)} \right)
.. math:: - \frac{\Delta tk_{or}(J)}{4A_{c}^{2}(J)}\left( \frac{W_{1}(J)\left| W_{1}(J) \right|}{\rho_{1}(J)} + \frac{W_{1}(J + 1)\left| W_{1}(J + 1) \right|}{\rho_{1}(J + 1)} \right)
.. _section-12.6.2.4:
Finite Differencing of the Momentum Equation for a Mesh Segment Which Contains a Bubble Interface
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The inclusion of a bubble interface in node :math:`J` introduces the Lagrangian
total derivative, as described in :numref:`section-12.6.2.2.1`. The momentum
equation then becomes
(12.6-134)
.. _eq-12.6-134:
.. math:: \frac{\partial p}{\partial z} + \frac{1}{A_{c}^{2}}\frac{\partial}{\partial z}\frac{W^{2}}{\rho_{v}} + \frac{f_{v}W|W|}{2\rho_{v}D_{h}A_{c}^{2}} - F_{c} + \frac{1}{A_{c}}\frac{dW}{dt} - v\frac{\partial W}{\partial z} + \frac{W\left| W \right|}{2\rho_{v}A_{c}^{2}}\frac{\partial k_{or}}{\partial z} = 0
As in the case of the continuity equation, the impact of the bubble
interface on the differenced equation is transmitted via the interface
velocity :math:`v^{i}` and the interface position :math:`z^{i}`. Therefore, only those
terms in Eq. :ref:`12.6-134 <eq-12.6-134>` which include either :math:`v^{i}` or :math:`z^{i}` will alter
the coefficient :math:`b_{i}` or the source term :math:`g` in Eq. :ref:`12.6-128 <eq-12.6-128>`.
.. _section-12.6.2.4.1:
Differenced Form of the Momentum equation at the Lower Interface
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
The terms of Eq. :ref:`12.6-134 <eq-12.6-134>` can be differenced as in :numref:`section-12.6.2.3`; in addition, the term
:math:`v\frac{\partial W}{\partial z}` can be expressed as
(12.6-135)
.. _eq-12.6-135:
.. math:: v\frac{\partial W}{\partial z} = \frac{v_{2}^{i}(J1)}{2}\frac{W_{2}(J1 + 1) - W_{2}^{i}(J1)}{\Delta z^{i}(J1)}
where :math:`\Delta z^{i}(J1) = z (J1 + 1) - z_{2}^{i}(J)`. Therefore, the differenced form of Eq. :ref:`12.6-134 <eq-12.6-134>` is
(12.6-136)
.. _eq-12.6-136:
.. math:: \frac{p_{2}(J1 + 1) - p_{2}^{i}(J1)}{\Delta z^{i}(J1)} + \frac{1}{A_{c}^{2}(J1)\Delta z^i(J1)}\left( \frac{W_{2}^{2}(J1 + 1)}{\rho_{2}(J1 + 1)} - \frac{W_{2}^{2}(J1)}{\rho_{2}(J1)} \right)
.. math:: + \frac{A_{fr}}{4D_{h}^{1 - b_{fr}}(J1)A_{c}^{2 + b_{fr}}(J1)}\left( \frac{W_{2}(J1 + 1)\left| W_{2}(J1 + 1) \right|^{1 + b_{fr}}}{\rho_{2}(J1 + 1)\mu_{2}^{b_{fr}}(J1 + 1)} + \frac{W_{2}^{i}(J1)\left| W_{2}^{i}(J1)^{1 + b_{fr}} \right|}{\rho_{2}(J1)\mu_{2}^{b_{fr}}(J1)} \right) - F_{e}(J1)
.. math:: + \frac{W_{2}(J1 + 1) + W_{2}^{i}(J1) - W_{1}(J1 + 1) - W_{1}^{i}(J1)}{2A_{c}(J1)\Delta t} - \frac{v_{2}^{i}(J1)}{2A_{c}(J1)}\frac{W_{2}^{i}(J1 + 1) - W_{2}^{i}(J1)}{\Delta z^{i}(J1)}
.. math:: + \frac{1}{4A_{c}^{2}(J1)}\left( \frac{W_{2}^{i}(J1)\left| W_{2}^{i}(J1) \right|}{\rho_{2}(J1)} + \frac{W_{2}(J1 + 1)\left| W_{2}(J1 + 1) \right|}{\rho_{2}(J1 + 1)} \right)\frac{k_{or}^{i}(J1)}{\Delta z(J1)} = 0
Multiplying Eq. :ref:`12.6-136 <eq-12.6-136>` by :math:`\Delta t\Delta z^{i}(J1)` eliminates
:math:`z_{2}^{i}(J1)` from all but the third, fourth, and fifth terms. The
sixth term is the only one which includes :math:`v_{2}^{i}(J1)`. Thus,
only these four terms need to be examined to determine the contribution
of the bubble interface to the final differenced momentum equation.
.. _section-12.6.2.4.2:
Specification of Bubble Interface Terms
'''''''''''''''''''''''''''''''''''''''
Of the four terms, all but the one involving :math:`v_{2}^{i} (J1)` are present in the momentum
equation as derived in :numref:`section-12.6.2.3`. They affect the interface
contribution only because they involve the factor :math:`z\Delta^{i}(J1)`. For example,
consider the condensation term:
(12.6-137)
.. _eq-12.6-137:
.. math:: \Delta t\Delta z^{i}(J1)F_{c}(J1) = \Delta tF_{c}(J1)\Delta z_{0}^{i}(J1) - \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\Delta tF_{c}(J1)
The first term on the right-hand side is identical to the condensation
term as expressed in the differenced momentum equation, Eq. :ref:`12.6-128 <eq-12.6-128>`.
Therefore, the condensation term contribution to the momentum equation
is modified to account for the bubble interface simply by adding the
second term on the right-hand side of Eq. :ref:`12.6-136 <eq-12.6-136>`
to Eq. :ref:`12.6-128 <eq-12.6-128>`. Written out in full, this
term is, after linearization of advanced time quantities and
substitution of the saturation condition relations,
(12.6-138)
.. _eq-12.6-138:
.. math:: \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\Delta tF_{c}(J1) = \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\frac{F_{o}(J1)\Delta t}{2A_{c}(J1)\lambda_{1}(J1)}
.. math:: \times\left( 1 - \frac{\Delta p^{i}(J1) + \Delta p(J1 + 1)}{2\lambda_{i}(J1)}\frac{d\lambda}{dp_{J1}} \right)
.. math:: \times\left( 1 + \frac{2}{F_{o}(J1)}\frac{dF_{o}}{dT_{J1}}\left( \Delta p(J1)\frac{dT^{i}}{dp_{J1}} + \Delta p(J1 + 1)\frac{dT}{dp_{J1 + 1}} \right) \right)
.. math:: \times \Bigg( \frac{W_{1}(J1)}{\rho_{1}(J1)} \left( 1 + \frac{\Delta W^{i}(J1)}{W_{1}^{i}(J1)} \right)\left( 1 - \frac{\Delta p^{i}(J1)}{\rho_{1}(J1 + 1)}\frac{d\rho}{dp_{J1 + 1}} \right)
.. math:: + \frac{W_{1}(J1 + 1)}{\rho_{1}(J1 + 1)}\left( 1 + \frac{\Delta W(J1 + 1)}{W_{1}(J1 + 1)} \right)\left( 1 - \frac{\Delta p(J1 + 1)}{\rho_{1}(J1 + 1)}\frac{d\rho}{dp_{J1 + 1}} \right) \Bigg)
Neglecting second-order terms, this reduces to
(12.6-139)
.. _eq-12.6-139:
.. math:: \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\Delta tF_{c}(J1) = \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\frac{F_{o}(J1)\Delta t}{2A_{c}(J1)\lambda_{1}(J1)}\left( \frac{W_{1}^{i}(J1)}{\rho_{1}(J1)} + \frac{W_{1}(J1 + 1)}{\rho_{1}(J1 + 1)} \right)
Similarly, the friction term contribution comes from
(12.6-140)
.. _eq-12.6-140:
.. math:: \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\Delta t\frac{A_{fr}}{4D_{h}^{1 - b_{fr}}(J1)A_{c}^{2 + b_{fr}}(J1)}\Bigg( \frac{W_{2}(J1 + 1)\left| W_{2}(J1 + 1) \right|^{1 + b_{fr}}}{\rho_{2}(J1 + 1)\mu_{2}^{b_{fr}}(J1 + 1)}
.. math:: + \frac{W_{2}^{i}(J1)\left| W_{2}^{i}(J1) \right|^{1 + b_{fr}}}{\rho_{1}(J1)\mu_{2}^{b_{fr}(J1)}} \Bigg) = \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}}\frac{A_{fr}\Delta t}{4D_{h}^{1 - b_{fr}}(J1)A_{c}^{2 + b_{fr}}(J1)}
.. math:: \times \left( \frac{W_{1}(J1 + 1)\left| W_{1}(J1 + 1) \right|^{1 + b_{fr}}}{\rho_{1}(J1 + 1)\mu_{2}^{b_{fr}}(J1 + 1)} + \frac{W_{1}^{i}(J1)\left| W_{1}^{i}(J1) \right|^{1 + b_{fr}}}{\rho_{1}(J1)\mu_{1}^{b_{fr}}(J1)} \right)
if second-order terms are dropped. The time derivative of the mass flow
rate is
(12.6-141)
.. _eq-12.6-141:
.. math:: \Delta t\Delta z^{i}(J1)\frac{\Delta W(J1 + 1) - \Delta W^{i}(J1)}{2A_{c}(J1)\Delta t}
.. math:: = \left( \Delta z_{o}^{i}(J1) - \Delta p^{i}(J1)\frac{dz^{i}}{dp_{J1}} \right)\frac{\Delta W(J1 + 1) + \Delta W^{i}(J1)}{2A_{c}(J1)}
.. math:: = \frac{\Delta z_{o}^{i}(J1)}{2A_{c}(J1)}\left( \Delta W(J1 + 1) + \Delta W^{i}(J1) \right)
if second-order terms are dropped. This is identical to the differenced
form of :math:`\frac{1}{A_{c}}\frac{\partial W}{\partial t}` for mesh
segments which are contained within the bubble and so this term actually
makes no contribution with respect to the bubble interface.
The last term to analyze is :math:`v\frac{\partial W}{\partial z}`.
This reduces to
(12.6-142)
.. _eq-12.6-142:
.. math:: \frac{v_{2}^{i}(J1)\Delta t}{2A_{c}(J1)}\left( W_{2}^{i}(J1 + 1) - W_{2}^{i}(J1) \right)
.. math:: = \frac{\Delta t}{2A_{c}(J1)}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) + \Delta p^{i}(J1)\frac{dv_{i}}{dp_{J1}} \right)
.. math:: \times \left( W_{1}(J1 + 1) + \Delta W(J1 + 1) - W_{1}^{i}(J1) - \Delta W^{i}(J1) \right)
.. math:: = \frac{\Delta t}{2A_{c}(J1)}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right)\left( W_{1}(J1 + 1) - W_{1}^{i}(J1) \right)
.. math:: + {(v}_{1}^{i}(J1) + \Delta v_{0}^{i}(J1))\Delta W(J1 + 1) - (v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1))\Delta W^{i}(J1)
.. math:: + \Delta p^{i}(J1)\frac{dv^{i}}{dp_{J1}}{(W}_{1}(J1 + 1) - W_{1}^{i}(J1))
Adding Eqs. :ref:`12.6-139 <eq-12.6-139>`, :ref:`12.6-140 <eq-12.6-140>` and :ref:`12.6-142 <eq-12.6-142>`
to Eqs. :ref:`12.6-128 <eq-12.6-128>` produces the final differenced momentum equation at the lower interface:
(12.6-143)
.. _eq-12.6-143:
.. math:: {\widetilde{b}}_{1,J1}\Delta p^{i}(J1) + {\widetilde{b}}_{2,J1}\Delta W^{i}(J1) + {\widetilde{b}}_{3,J1}\Delta p(J1 + 1) + {\widetilde{b}}_{4,J1}\Delta W(J1 + 1) = {\widetilde{g}}_{J1}
where
(12.6-144)
.. _eq-12.6-144:
.. math:: {\widetilde{b}}_{1,J1} = b_{1,J1} + \frac{dz^{i}}{dp_{J1}}\Bigg[ \frac{F_{o}(J1)\Delta t}{2A_{c}(J1)\lambda_{1}(J1)}\left( \frac{W_{1}^{i}(J1)}{\rho_{1}(J1)} + \frac{W_{1}(J1 + 1)}{\rho_{1}(J1 + 1)} \right)
.. math:: - \frac{A_{fr}\Delta t}{4D_{h}^{1 - b_{fr}}(J1)A_{c}^{2 + b_{fr}}(J1)}\Bigg(\frac{W_{1}(J1 + 1)\left| W_{1}(J1 + 1) \right|^{1 + b_{fr}}}{\rho_{1}(J1 + 1)\mu_{1}^{b_{fr}}(J1 + 1)}
.. math:: + \frac{W_{1}^{i}(J1)\left| W_{1}^{i}(J1) \right|^{1 + b_{fr}}}{\rho_{1}(J1)\mu_{1}^{b_{fr}}(J1)} \Bigg)\Bigg]
.. math:: - \frac{dv^{i}}{dp_{J1}}\frac{\Delta t}{2A_{c}(J1)}\left(W_{1}(J1 + 1) - W_{1}^{i}(J1) \right)
(12.6-145)
.. _eq-12.6-145:
.. math:: {\widetilde{b}}_{2,J1} = b_{2,J1} + \frac{\Delta t}{2A_{c}(J1)}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right)
(12.6-146)
.. _eq-12.6-146:
.. math:: {\widetilde{b}}_{3,J1} = b_{3,J1}
(12.6-147)
.. _eq-12.6-147:
.. math:: {\widetilde{b}}_{4,J1} = b_{4,J1} - \frac{\Delta t}{2A_{c}(J1)}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right)
(12.6-148)
.. _eq-12.6-148:
.. math:: {\widetilde{g}}_{J1} = g_{J1} + \frac{\Delta t}{2A_{c}(J1)}\left( v_{1}^{i}(J1) + \Delta v_{0}^{i}(J1) \right)\left( W_{1}(J1 + 1) - W_{1}^{i}(J1) \right)
with :math:`b_{1,J1}`, :math:`b_{2,J1}`, :math:`b_{3,J1}`, :math:`b_{4,J1}`, and :math:`g_{J1}`
being defined as in Eq. :ref:`12.6-128 <eq-12.6-128>`.
.. _section-12.6.2.4.3:
Differenced Form of the Momentum Equation at the Upper Interface
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
The development of the momentum equation at the upper interface is
identical to that explained in :numref:`section-12.6.2.4.1` and
:numref:`section-12.6.2.4.2` for the lower interface, and so only the final equation will be given here:
(12.6-149)
.. _eq-12.6-149:
.. math:: {\widetilde{b}}_{1,J2 - 1}\Delta p(J2 - 1) + {\widetilde{b}}_{2,J2 - 1}\Delta W(J2 - 1) + {\widetilde{b}}_{3,J2 - 1}\Delta p^{i}(J2)
.. math:: {} + {\widetilde{b}}_{4,J2 - 1}\Delta W^{i}(J2) = {\widetilde{g}}_{J2 - 1}
with
(12.6-150)
.. _eq-12.6-150:
.. math:: {\widetilde{b}}_{1,J2 - 1} = b_{1,J2 - 1}
(12.6-151)
.. _eq-12.6-151:
.. math:: {\widetilde{b}}_{2,J2 - 1} = b_{2,J2 - 1} + \frac{\Delta t}{2A_{c}(J2 - 1)}\left( v_{1}^{i}(J2) + \Delta v_{0}^{i}(J2) \right)
(12.6-152)
.. _eq-12.6-152:
.. math:: {\widetilde{b}}_{3,J2 - 1} = b_{3,J2 - 1} - \frac{\Delta t}{2A_{c}(J2 - 1)}\frac{dv^{i}}{dp_{J2}}\left( W_{1}^{i}(J2) - W_{1}(J2 - 1) \right)
.. math:: + \frac{A_{fr}\Delta t}{4D_{h}^{1 - b_{fr}}(J2 - 1)A_{c}^{2 + b_{fr}}(J2 - 1)}\frac{dz^{i}}{dp_{J2}}\Bigg( \frac{W_{1}^{i}(J2)\left| W_{1}^{i}(J2) \right|^{1 + b_{fr}}}{\rho_{1}(J2)\mu_{1}^{b_{fr}}(J2)}
.. math:: + \frac{W_{1}(J2 - 1)\left| W_{1}(J2 - 1) \right|^{1 + b_{fr}}}{\rho_{1}(J2)\mu_{1}^{b_{fr}}(J2)} \Bigg)
.. math:: - \frac{\Delta t}{2A_{c}(J2 - 1)}\frac{dz^{i}}{dp_{J2}}\frac{F_{o}(J2 - 1)}{\lambda_{1}(J2 - 1)}\left( \frac{W_{1}^{i}(J2)}{\rho_{1}(J2)} + \frac{W_{1}(J2 - 1)}{\rho_{1}(J2 - 1)} \right)
(12.6-153)
.. _eq-12.6-153:
.. math:: {\widetilde{b}}_{4,J2 - 1} = b_{4,J2 - 1} - \frac{\Delta t}{2A_{c}(J2 - 1)}v_{1}^{i}(J2) + \Delta v_{0}^{i}(J2)
(12.6-154)
.. _eq-12.6-154:
.. math:: {\widetilde{g}}_{J2 - 1} = g_{J2 - 1} + \frac{\Delta t}{2A_{c}(J2 - 1)}\left( v_{1}^{i}(J2) + \Delta v_{0}^{i}(J2) \right)\left( W_{1}^{i}(J2) - W_{1}^{i}\left( J - 2 \right) \right)
.. _section-12.6.2.5:
Additional Details Concerning Interface Nodes
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Because at least one of the boundaries of an interface mesh segment is
not usually aligned with the fixed segment boundaries of the axial mesh,
variables defined at the interface segment midpoint must be handled
somewhat differently from what they are for a segment located in the
interior of a bubble. Midpoint quantities such as :math:`A_c` and :math:`D_{h}`,
which are considered to be piecewise constant over the fixed mesh
segments, are assigned the values of the fixes mesh segment in which the
bubble interface lies. Thus, for example, the flow area at the upper
interface is taken as :math:`A_{c}(J2 - 1)` because the upper interface
lies in axial mesh segment :math:`J2 - 1`. This convention has been used
throughout the equations in :numref:`section-12.6.2.2` and
:numref:`section-12.6.2.4`. However, three midpoint variables are exceptions to this rule: the cladding
temperature :math:`T_{e}`, the structure temperature :math:`T_{s}`, and the average
coolant temperature :math:`\overline{T}`. :math:`T_{e}` and :math:`T_{s}` are calculated outside the
two-phase flow model once every heat-transfer time step. As used in the
boiling calculation, both temperatures are extrapolated values
calculated using the temperatures at the beginning and the end of the
last heat-transfer time step and extrapolated to the end of the current
coolant time step. All three temperatures are taken to be located at the
midpoints of the Eulerian grid intervals, whereas the boiling
calculation needs them at the midpoints of the interface vapor mesh
segments (for example, :math:`T_{e}`, :math:`T_{s}`, and :math:`\overline{T}` are needed at a point halfway
between :math:`z (J2 - 1)` and the upper interface). Therefore, :math:`T_{e}` and :math:`T_{s}`
are interpolated in space to the center of the interface segments
using the values in the interface fixed mesh segments and the adjacent
vapor segments, while :math:`\overline{T}` is interpolated to the
points using the temperatures at the interfaces and at the nearest fixed
mesh boundaries contained in the bubble.
Special considerations must be given to any mesh boundary that is
crossed by a liquid-vapor interface during a time step. When setting up
the interior segments for the pressure-drop calculation, any segment
which is outside the bubble at either the beginning or the end of the
time step is bypassed and treated as part of the interface segment.
Also, if at the end of the step the interface is within :math:`\varepsilon` of a segment
boundary (:math:`\varepsilon = .002`) then the segment is also bypassed. These
tests are all made on the basis of :math:`z_{1}^{i}(J)` and
:math:`z_{0}^{i}(J)` as given by Eq. :ref:`12.6-68 <eq-12.6-68>`. Therefore, a
segment :math:`JC` is included as an interior segment only if
.. math:: z_{1}^{i}(J1) < z(JC)
.. math:: z_{0}^{i}(J1) < z\left( JC \right) - \varepsilon
.. math:: z_{1}^{i}(J2) > z(JC + 1)
and
.. math:: z_{0}^{i}(J2) > z\left( JC + 1 \right) + \varepsilon
After the calculations of :math:`\Delta p` and :math:`\Delta W` for segments inside the
bubble are completed, :math:`p_{2}` and :math:`W_{2}` for bypassed
segments are obtained by linear interpolation between the interface and
the first segment inside the interface.
If an interface segment, :math:`JT = J1` or :math:`JT = J2`, represents parts of
more than one fixed Eulerian mesh segment, then :math:`A_{c}(JT)` and :math:`D_{h}(JT)`
for the interface segment are obtained as weighted averages over the Eulerian
segment values. The weighting factors are proportional to the :math:`\Delta z` within the
bubble. Also, :math:`1/D_{h}` is averaged:
(12.6-155)
.. _eq-12.6-155:
.. math:: \frac{1}{D_{h}(JT)} = \frac{\sum_{\text{JC}}^{}\frac{\Delta z(JC)}{D_{h}(JC)}}{\sum_{\text{JC}}^{}{\Delta z(JC)}}
The :math:`\Delta z(JC)`\ s are based on :math:`z_0^i(JT)`.
.. _section-12.6.3:
Simultaneous Solution of the Differenced, Linearized Mass and Momentum Equations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
If the discretized mass and momentum equations derived in :numref:`section-12.6.2`
(Eqs. :ref:`12.6-45 <eq-12.6-45>` and :ref:`12.6-128 <eq-12.6-128>`) are applied to
each axial segment in the bubble, the result is a set of simultaneous linear equations which have as independent
variables only the changes in vapor pressure and in mass flow rate. To this set
must be added the equations at the upper and lower bubble interfaces (Eqs.
:ref:`12.6-79 <eq-12.6-79>`, :ref:`12.6-87 <eq-12.6-87>`, :ref:`12.6-143 <eq-12.6-43>`, and :ref:`12.6-149 <eq-12.6-149>`),
giving a set of :math:`2 \left( J2 - J1 \right)` equations in :math:`2 \left( J2 - J1 + 1 \right)`
unknowns. To complete the set, a boundary condition must be imposed at both
interfaces. This requirement is satisfied by expressing the interface mass flow
rate as
(12.6-156)
.. _eq-12.6-156:
.. math:: W_{1}^{i}\left( JT \right) = v_{1}^{i}\left( JT \right)\rho_{1}^{i}\left( JT \right)A_{c}(JX)
which can be differenced to give
(12.6-157)
.. _eq-12.6-157:
.. math:: \Delta W^{i}\left( JT \right) = \left( v_{1}^{i}\left( JT \right)\Delta\rho^{i}\left( JT \right) + \rho_1^{i}\left( JT \right)\Delta v^{i}\left( JT \right) \right)A_{c}(JX)
The index :math:`JT` is :math:`J1` for the lower interface and :math:`J2` for the upper
interface, while :math:`JX` is :math:`J1` for the lower interface and :math:`J2-1` at the upper
one. Using the expressions derived earlier for
:math:`\Delta \rho^{i}\left( JT \right)` and :math:`\Delta v^{i}\left( JT \right)`
(Eqs. :ref:`12.6-37 <eq-12.6-37>` and :ref:`12.6-62 <eq-12.6-62>`) produces
(12.6-158)
.. _eq-12.6-158:
.. math:: \Delta W^{i}\left( JT \right) = \left( v_{1}^{i}\left( JT \right)\frac{d\rho^{i}}{dp_{JT}} + \rho_{1}^{i}\left( JT \right)\frac{dv^{i}}{dp_{JT}} \right)A_{c}\left( JX \right)\Delta p^{i}\left( JT \right)
.. math:: + \Delta v_{0}^{i}\left( JT \right)A_{c}\left( JX \right)\rho_{1}^{i}(JT)
which can be rearranged to give
(12.6-159)
.. _eq-12.6-159:
.. math:: - \left( v_{1}^{i}\left( JT \right)\frac{d\rho^{i}}{dp_{JT}} + \rho_{1}^{i}\left( JT \right)\frac{dv^{i}}{dp_{JT}} \right)A_{c}\left( JX \right)\Delta p^{i}\left( JT \right) + \Delta W^{i}\left( JT \right)
.. math:: = \rho_{1}^{i}\left( JT \right)\Delta v_{0}^{i}\left( JT \right)A_{c}(JX)
Define
(12.6-160)
.. _eq-12.6-160:
.. math:: a_{1,JX} = - \left( v_{1}^{i}\left( JT \right)\frac{d\rho^{i}}{dp_{JT}} + \rho_{1}^{i}\left( JT \right)\frac{dv^{i}}{dp_{JT}} \right)A_{c}(JX)
(12.6-161)
.. _eq-12.6-161:
.. math:: a_{2,JX} = 1.0
and
(12.6-162)
.. _eq-12.6-162:
.. math:: d_{JX} = \rho_{1}^{i}\left( JT \right)\Delta v_{0}^{i}\left( JT \right)A_{c}(JX)
Then Eq. 12.6-159 becomes
(12.6-163)
.. _eq-12.6-163:
.. math:: a_{1,JX}\Delta p^{i}\left( JT \right) + a_{2,JX}\Delta W^{i}\left( JT \right) = d_{JX}
If Eq. :ref:`12.6-163 <eq-12.6-163>` is applied to both upper and lower interfaces and the resulting
two equations are added to the set of mass and momentum equations, a set of
:math:`2 \left( J2 - J1 + 1 \right)` equations in :math:`2 \left( J2 - J1 + 1 \right)`
unknowns will result. This set can be written in matrix form as
(12.6-164a)
.. _eq-12.6-164a:
.. math:: \mathbf{AX} = \mathbf{B}
where the form of matrix :math:`\mathbf{A}` is given by Eq. :ref:`12.6-164b<eq-12.6-164b>`
and the :math:`\mathbf{X}` and :math:`\mathbf{B}` vectors are given by Eqs. :ref:`12.6-164c<eq-12.6-164c>` and
:ref:`12.6-164d<eq-12.6-164d>`.
(12.6-164b)
.. _eq-12.6-164b:
.. math::
\mathbf{A} = \begin{bmatrix}
a_{1,J1} & a_{2,J1} & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\
\widetilde{b}_{1,J1} & \widetilde{b}_{2,J1} & \widetilde{b}_{3,J1} & \widetilde{b}_{4,J1} & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\
c_{1,J1} & c_{2,J1} & c_{3,J1} & c_{4,J1} & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\
0 & 0 & b_{1,J1+1} & b_{2,J1+1} & b_{3,J1+1} & b_{4,J1+1} & \cdots & 0 & 0 & 0 & 0 \\
0 & 0 & c_{1,J1+1} & c_{2,J1+1} & c_{3,J1+1} & c_{4,J1+1} & \cdots & 0 & 0 & 0 & 0 \\
& & & & \ddots & \ddots & \ddots & \ddots & \ddots & & \\
0 & 0 & 0 & 0 & 0 & 0 & \cdots & \widetilde{b}_{1,J2-1} & \widetilde{b}_{2,J2-1} & \widetilde{b}_{3,J2-1} & \widetilde{b}_{4,J2-1} \\
0 & 0 & 0 & 0 & 0 & 0 & \cdots & c_{1,J2-1} & c_{2,J2-1} & c_{3,J2-1} & c_{4,J2-1} \\
0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & a_{1,J2-1} & a_{2,J2-1} \\
\end{bmatrix}
(12.6-164c)
.. _eq-12.6-164c:
.. math::
\mathbf{X} = \begin{bmatrix}
\begin{matrix}
\Delta p^{i}(J1) \\
\Delta W^{i}(J1) \\
\Delta p(J1 + 1) \\
\end{matrix} \\
\Delta W(J1 + 1) \\
\vdots \\
\Delta p(J2 - 1) \\
\Delta W(J2 - 1) \\
\Delta p^{i}(J2) \\
\Delta W^{i}(J2) \\
\end{bmatrix}
and
(12.6-164d)
.. _eq-12.6-164d:
.. math::
B = \begin{bmatrix}
\begin{matrix}
d_{J1} \\
{\widetilde{g}}_{J1} \\
h_{J1} \\
\end{matrix} \\
g_{J1 + 1} \\
h_{J1 + 1} \\
\vdots \\
{\widetilde{g}}_{J2 - 1} \\
h_{J2 - 1} \\
d_{J2 - 1} \\
\end{bmatrix}
This matrix equation can be solved using Gaussian elimination to give
the changes in vapor pressure and mass flow rate at all nodes in the
bubble. :numref:`section-12.appendices` presents a step-by-step description of how Eq.
:ref:`12.6-164 <eq-12.6-164a>` is solved using Gaussian elimination. Once the new vapor
pressures are known, the vapor temperatures and all thermodynamic sodium
properties can be calculated.