13.2. Mathematical Model

13.2.1. Sodium Vapor Flow Model

The sodium vapor flow generates the axial pressure gradients and interfacial shear forces, which influence the motion of molten cladding. The cladding motion in turn changes flow areas, hydraulic diameters, and surface roughness which all affect the sodium vapor flow. The governing equations for the sodium vapor, while presented elsewhere (Section 12.6.1) in more detail, are reproduced here to show the intimate coupling between the sodium vapor dynamics and cladding motion. The conservation equations for the sodium vapor are

(13.2-1)t(Avρv)+wz=˙mv
(13.2-2)wt+z(w2Avρv)+Av[pz+Fv] ={˙mvwAvρv,if condensing(13.2‑2a)0,if vaporizing(13.2‑2b)

where

Av = vapor flow area

ρv = vapor density

w = vapor mass flowrate

˙mv = the rate of vapor generation per unit length of channel

pz = channel axial pressure gradient

Fv = cladding/vapor interfacial friction force per unit volume of vapor.

The friction force term is evaluated from the following equation:

(13.2-3)Fv=fsfMw|w|2ρvA2vDv

where

fsf = single-phase friction factor

M = flooding multiplier

Dv = hydraulic diameter for the vapor.

The hydraulic diameter is estimated by

(13.2-4)Dv=Dhα

where

Dh = hydraulic diameter for bare fuel or fuel pin (depending upon the zone)

α = vapor fraction based on area available for molten steel plus vapor.

This equation assumes that the hydraulic diameter Dv varies as the square root of the area v (which would be exact for a circular vapor flow passage) and hence scales as the square root of vapor fraction (proportional to the flow area).

The smooth wall (single-phase) friction factor is given by

(13.2-5)fsf=AFRV(Re)BFRVv
(13.2-6)(Re)v=wDvAvμv

where

AFRV, BFRV = input constants

μv = vapor viscosity.

For consistency, the inputed constants AFRV and BFRV used in CLAP are identical to those used in the sodium voiding model.

The flooding multiplier is flow-regime dependent. The transition from smooth wall friction (M=1) to a chaotic two-phase flow structure (M=M2Φ, where M2Φ is defined later) occurs nominally when the vapor velocity exceeds the flooding velocity. The current version of CLAP uses the correlation of Fauske, et al. [13-4]:

(13.2-7)vflood=[323ρcgΔrcρvfsfM2Φ]1/2

where

Δrc = one-half the molten cladding thickness

ρc = molten cladding density.

The CLAP model has a hysteresis effect in the treatment of the flooding multiplier whereby flooding and the inverse, deflooding, occur at slightly different velocities as follows:

(13.2-8)set Mn=M2Φ,if Mn1=1and |wAvρv|>1.1vflood
(13.2-9)set Mn=1,if Mn1=M2Φand |wAvρv|<0.9vflood

where n is the current time-step number. If the conditions in either Eq. (13.2-8) or Eq. (13.2-9) are not met, then the previous value of M is retained.

The two-phase friction multiplier is given by

(13.2-10)M2Φ=I[1+ε(1α)]for α>αcrit
(13.2-11)M2Φ=I[1+ε(1αcrit)]for ααcrit

where ε and αcrit are input and I is an incoherence multiplier. For values of ε=75 and αcrit=0, Eq. (13.2-10) is equivalent to the Wallis correlation [13-5].

Radial incoherency in cladding melting, which occurs during the early portion of cladding motion, results in a lower bundle pressure drop and lower cladding velocity than predicted by one-dimensional models that ignore this effect [13-6]. The incoherence multiplier allows for this incoherency effect by temporarily reducing the two-phase friction multiplier during early cladding motion in the following manner:

(13.2-12)I=[fps(fps)0]xfor fps<(fps)0
(13.2-13)I=1for fps(fps)0

where fps are full-power seconds since cladding first moved in a channel and (fps)0 and x are input constants.

The code sets M2ϕ equal to unity, equivalent to a smooth wall friction factor, if the value of M2ϕ computed by Eq. (13.2-10) is less than unity. Consequently, the effect of the multiplier “I” is to delay flooding by (fps)0 full-power seconds. For cases where the vapor velocity is below that necessary for flooding, neither (fps)0 nor ε will have any effect on the cladding motion. Values for (fps)0 and x of 0.3 and 3 are provisionally recommended; better values will be determined later by calibration of the CLAP model using experimental data.

13.2.2. Moving Cladding Basic Equations

The mass, momentum, and energy conservation equations for moving (molten) cladding are

(13.2-14)t(ρcAc)+z(ρcAcvc)=˙mc
(13.2-15)t(ρcAccc)+z(ρcAcv2c)+Acpz +AcFpAvFv+Acρcg={˙mcvcif freezing(13.2‑13a)0,if melting(13.2‑13b)
(13.2-16)t(ρcAcec)+z(ρcAcvcec)=ϕcPr+˙mcec

Ac = moving cladding cross-sectional area

vc = moving cladding velocity

˙mc = mass rate of melting of substrate per unit length of channel (or rate of freezing if negative)

Fp = pin friction-force per unit volume of molten cladding

ec = moving cladding internal energy

ϕc = flux of sensible heat from the solid/liquid interface into the molten cladding

Pr = perimeter of solid/liquid interface.

The pin or bare-fuel friction force is evaluated by the following:

(13.2-17)Fp=cfρcvc|vc|2Dc

where

cf = pin friction factor

Dc = molten cladding hydraulic diameter.

The CLAZAS [13-3] formulation for the friction factor has been incorporated into CLAP:

(13.2-18)cf=64Refor Re<(Re)break
(13.2-19)cf=bffor Re(Re)break

where Re is the molten-cladding Reynolds number defined by:

(13.2-20)Re=Dcρc|vc|μc

(Re)break is the critical Reynolds number for the transition from laminar to turbulent friction and bf is the turbulent Moody friction factor. The CLAZAS derived values for (Re)break and bf of 2.1×103 and 2.0×102 are currently recommended as input values.

The cladding viscosity depends upon the cladding temperature Tc and melt-fraction f as follows

(13.2-21)μc=μmexp(aTcaTm)for Tc>Tm
(13.2-22)μc=(μtμm)(1f)q+μmfor Tc=Tm
(13.2-23)μc=μsfor Tc<Tm

where Tm is the cladding melting temperature, μs, μt, μm are the user-supplied viscosities for solid cladding, at the solidus temperature, and cladding at the liquidus temperature, and a and q are input constants. Values for μm and a of 6.42×103 Pa-s and 5492 K are recommended [13-7]. The input parameters μt,μs, and q are not properties per se, but are constants in a phenomenological model and hence need to be determined by calibration with cladding motion test data. In general, the values selected for μt and μs should be sufficiently large to effectively arrest cladding motion for solid cladding. Provisional values for both μt and μs of 104 Pa-s (105 poise), from the CLAZAS sample problem [13-3], and for q of 0.5 are recommended.

The momentum and energy Eq. (13.2-15) and Eq. (13.2-16) are converted to non-conservative form using Eq. (13.2-14), giving the following:

(13.2-24)ρc[vct+vcvct]+pz+FpAvAcFv+gρc={0if freezing(13.2‑19a)mc˙vcAcif melting(13.2‑19b)
(13.2-25)ect+vcecz=ϕcPrAcρc

13.2.3. Refrozen Cladding Basic Equations

The conservation equations for the refrozen steel layer are

(13.2-26)t(ρsAs)=˙mc
(13.2-27)t(ρsAses)=ϕrPeϕPr˙mees

where

ρs = refrozen steel density

As = refrozen steel cross-sectional area

es = refrozen steel internal energy

ϕr = heat flux at the interface between intact cladding and refrozen cladding

Pe = outer perimeter of intact cladding

ϕ = flux of sensible heat from the refrozen cladding to the melt (solid/liquid) interface.

Recall that the quantity ˙mc was the rate of melting (negative for freezing) and hence appears as a sink (negative sign) in the refrozen cladding Eq. (13.2-26) and Eq. (13.2-27).

By manipulating the preceding equations, we obtain the nonconservative form of the energy equation:

(13.2-28)est=(ϕrPeϕPr)ρsAs

Note that the sensible heat fluxes on either side of the melt layer are not the same (see Figure 13.2.1). The difference ϕhf is the heat flux for fusion, defined by

(13.2-29)ϕhf=ϕϕc

The rate of melting (or freezing, if negative) is directly related to the fusion heat flux:

(13.2-30)˙mc=Prϕhfλ

where λ is an effective heat of fusion of the cladding defined as

(13.2-31)λ=eces

13.2.4. Heat-transfer Relationships

13.2.4.1. Moving Cladding on Bare Fuel

The configuration in the molten zone is shown in Figure 13.2.1. There is no solid cladding in this region and, therefore, the heat flux for melting, ϕhf, is zero. The sensible heat flux is given by

(13.2-32)ϕc=h(TfTc)

where

Tf = fuel pellet surface temperature

h = coefficient of heat transfer from the fuel to molten cladding.

The heat-transfer coefficient for the moving cladding is defined from a correlation for liquid-metal heat transport as

(13.2-33)hDckc=C1(Dcρccpc|vc|kc)C2+C3

where

kc = molten cladding thermal conductivity

cpc = molten cladding specific heat capacity

C1,C2,C3 = input constants.

../../_images/image416.png

Figure 13.2.1 Schematic Showing the Different Heat-Transfer Configuration in CLAP

The same correlation and constants are used for computing heat-transfer coefficients for single-phase sodium flow in the channel, (Eq. (3.3-9) of Chapter 3). The correlation of Maresca and Dwyer [13-8] for in-line flow in rod bundles is recommended, for which the constants are:

(13.2-34)C1=0.0155(¯ψ)0.86
(13.2-35)C1=0.86
(13.2-36)C3=6.66+3.126(PD)+1.184(PD)2

where ¯ψ = mean ratio of the thermal and momentum diffusivities, and PD = rod pitch-to-diameter ratio. Since the code does not allow for flow-rate-dependent values for C1, it is suggested that ¯ψ be set equal to unity in Eq. (13.2-34) with the understanding that the resulting heat-transfer coefficients will be slightly overestimated.

13.2.4.2. Moving Cladding on Intact and Refrozen Cladding

Various configurations can exist in the region outside the molten zone, depending upon whether or not moving cladding and/or refrozen cladding are present in addition to the original intact solid cladding. In this subsection, we consider the case where all three components are present (see Figure 13.2.1). The various heat fluxes shown in the figure are given by

(13.2-37)ϕr=kc(TiTs)(Δri+Δrs)
(13.2-38)ϕr=kc(TsTm)Δrs
(13.2-39)ϕc=h(TmTc)

where

Ti = temperature of the intact solid cladding

Ts = temperature of the refrozen cladding

Tm = cladding melting temperature

Δri = half-thickness of the intact solid cladding

Δrs = half-thickness of the refrozen cladding.

The film coefficient formula is identical to Eq. (13.2-33) and the fusion heat flux ϕhf is given by Eq. (13.2-29).

13.2.4.3. Intact and Refrozen Cladding

In the absence of moving cladding (see Figure 13.2.1), we set ϕ=0 and ϕhf=0. The heat transfer between the initial and refrozen cladding is computed from

(13.2-40)ϕr=kc(TiTs)(Δri+Δrs)

which is identical to Eq. (13.2-37).

13.2.4.4. Intact and Moving Cladding

The moving cladding may be sufficiently hot to melt through the layer of refrozen cladding; this results in the configuration (see Figure 13.2.1) in which just the intact cladding and moving cladding are present. In the numerical model this situation also arises when moving cladding first contacts the intact cladding in an axial segment, in which case refrozen cladding is formed and will be present in the next time step.

The heat fluxes ϕ1 and ϕ2 and trial heat flux are calculated as follows

(13.2-41)ϕ1=kc(TiTm)Δri
(13.2-42)ϕ2=h(TmTc)
(13.2-43)ϕtrial=ϕ1ϕ2

Recalling that ϕr is the outer surface heat flux for the initial cladding, ϕc the sensible heat flux for moving cladding, and ϕhf the latent heat flux associated with refrozen cladding, we have

(13.2-44)ϕhf=0,ϕr=ϕ2, and ϕc=ϕ2for ϕtrial0
(13.2-45)ϕhf=ϕtrial,ϕr=ϕ1, and ϕc=ϕ2for ϕtrial<0

13.2.5. Cladding Physical Properties Relationships

The density of the solid and liquid cladding are assumed to be linear functions of temperature:

(13.2-46)ρs=ρ0s[13β(TTref)]
(13.2-47)ρc=ρ0c[1C(TTm)]

where

ρ0s = density of solid cladding at reference temperature Tref

ρ0c = density of cladding at the liquidus temperature

β = linear coefficient of thermal expansion for solid

C = volumetric coefficient of thermal expansion for liquid.

For partly-molten cladding, the density is evaluated using

(13.2-48)ρ=[(1f)ρs(Tm)+fρc]1

where f is the melt fraction as defined by Eq. (13.2-51).

The internal energies of the solid and molten cladding are also assumed to be linear functions of temperature

(13.2-49)es=cpsT
(13.2-50)ec=eoc+cpcT

where eoc=λo+(CpsCpc)Tm, and math:lambda^{o} = heat of fusion. For mixtures, the mass melt fraction is given by

(13.2-51)f=(ecpsTm)λo