function THsolver_simple()
global Af Ac;
global Tin Pressure w;
global nr nz dz length;
global Rf0 Rf1 Rc0 Rc1;
close all;
% --- System parameters
Tin = 538; % Inlet temperature (K)
w = 1.0*7100*6*1000/60/60/313/126; % Flow rate (kg/s) based on total flow
% rate
%w = 1500e6/(35*5300)/313/126; % Flow rate (kg/s) based on heat
% transfer
w = 0.5*w;
Pressure = 12.3e6; % System pressure (Pa)
% --- Rod geometry
Rf0 = 0.06e-2; % Fuel inner radius (m)
Rf1 = 0.38e-2; % Fuel outer radius (m)
Rc0 = 0.3865e-2; % Cladding inner radius (m)
Rc1 = 0.4535e-2; % Cladding outer radius (m)
length = 248e-2; % Total length (m)
pitch = 1.23e-2; % Pin pitch (m)
% --- Fuel and coolant area
%Ac = pitch^2-Rc1^2*pi;
Af = Rf1^2*pi;
% --- Number of axial and radial zones in fuel rod
nz = 100;
nr = 50;
dz = length/nz;
Z = linspace(-length/2, length/2, nz+1)';
R = [0 linspace(Rf0, Rf1, nr-3+1) Rc0 Rc1]';
% --- Calculate the area of each annulus
A = pi*(R(2:end).*R(2:end) - R(1:end-1).*R(1:end-1));
% --- Get power density
Pden = GetPowerDensity(Z, R);
% --- Solve the coolant temperature distribution
Tc = SolveCoolantTemperature(Z, Pden, A);
% --- Solve the fuel temperature distribution
Tf = SolveFuelTemperature(R, Pden, A, Tc);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display(GetFuelK(1000))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display(GetGasK(700))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display(GetCladK(560))
figure
%plot(0.5*(Z(1:end-1)+Z(2:end)), Pden,'rx')
hold on
plot(Z, Tc,'b-')
plot(0.5*(Z(1:end-1)+Z(2:end)), Tf(:,end),'k-')
hold off
grid on
xlabel('Axial coordinate')
ylabel('Temperature (K)')
legend('Coolant', 'Cladding','location','northwest')
figure
plot(0.5*(Z(1:end-1)+Z(2:end)), sum(Pden.*repmat(A',nz,1),2)*dz,'k.')
hold on
plot(0.5*(Z(1:end-1)+Z(2:end)), Tf(:,1),'r-')
hold off
grid on
xlabel('Axial coordinate')
ylabel('Fuel centerline temperature / power level')
figure
%plot(0.5*(R(1:end-1)+R(2:end)),Pden(ceil(nz/2),:),'rx')
hold on
plot(R*100,Tf(ceil(nz/2),:),'bx')
hold off
grid on
xlabel('Radius (cm)')
ylabel('Fuel temperature (K)')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Tf = SolveFuelTemperature(R, P, A, Tc)
global nz nr Tin w dz Rf0 Rf1 Rc0 Rc1;
% --- R is a vertical vector containing the r-limits of the annular
% zones
% --- P is a matrix containing the power density in each of the
% fuel rod annular zones
% --- Tc is a vector containing the coolant temperatures at each of
% the axial elevations (nz + 1)
% --- We'll start with an initial guess for the fuel temperature
% distribtion
Tf = zeros(nz,nr+1)+100;
% --- Loop over axial zones
for i = 1:1:nz
% --- Get bulk coolant temperature for this axial zone
Tb = 0.5*(Tc(i) + Tc(i+1));
% --- We will start at the cladding outer surface
% --- Calculate the heat flux through the cladding surface
Hflux = sum(P(i,1:end).*A'*dz)/(dz*2*pi*R(end));
% --- Get heat transfer coefficient between cladding and coolant
% (now constant)
hclad = GetCladH(Tb);
% --- Calculate and store cladding outer temperature
Tf(i, nr + 1) = Tb + Hflux/hclad;
Tout = Tf(i, nr + 1);
for j = nr:-1:1
% --- Store zone inner and outer radius
Ri0 = R(j);
Ri1 = R(j+1);
% --- Calculate the heat flux through the zone inner surface
if (j > 2)
Hflux = sum(P(i,1:j-1).*A(1:j-1)'*dz)/(dz*2*pi*Ri0);
else
Hflux = 0;
end
% --- Calculate the power produced in this zone
Pdeni = P(i,j);
% --- Calculate increase in temperature over this segment
epsilon = 1.0;
Told = Tf(i,j);
while (epsilon > 0.01)
% --- Get the heat conductivity of this zone
Teff = 0.5*(Tf(i,j+1) + Tf(i,j));
if (Ri0 >= Rc0)
k = GetCladK(Teff);
elseif ((Ri0 >= Rf1) || (Ri1 <= Rf0))
k = GetGasK(Teff);
else
k = GetFuelK(Teff);
end
% --- Calculate the constants C0, C1 and C2
C0 = -Pdeni/(4*k);
if (j > 1)
C1 = Pdeni*Ri0^2/(2*k) - Hflux*Ri0/k;
else
C1 = 0;
end
C2 = Tout + Pdeni/(4*k)*Ri1^2 - C1*log(Ri1);
% --- Calculate and store zone inner temperature
if (j > 1)
Tf(i, j) = C0*Ri0^2 + C1*log(Ri0) + C2;
else
Tf(i, j) = C0*Ri0^2 + C2;
end
epsilon = (Tf(i,j) - Told)/Told;
Told = Tf(i,j);
end
% --- This inner temperature will be the outer temperature for
% the next zone
Tout = Tf(i, j);
end
% --- Calculate heat transfer in each of the fuel rings
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Tc = SolveCoolantTemperature(Z, Pden, A)
global nz Tin w dz Af;
% --- Z is a vertical vector containing the z-limits of the axial
% nodes
% --- Pden is a matrix containing the power density in each of the
% fuel rod annular zones
% --- A is a matrix containing the cross sectional area for each of the
% fuel rod annular zones
% --- We'll start with an initial guess for the coolant temperature
% distribtion
Tc = ones(nz + 1,1)*Tin;
% --- Loop over the axial length of the rod and solve the temperature
% of each axial zone
for i = 2:1:nz+1
% --- Calculate amount of power flowing into this segment
Pin = sum(Pden(i-1,:).*A'*dz);
% --- Calculate increase in temperature over this segment
epsilon = 1.0;
Told = Tc(i);
while (epsilon > 0.01)
% --- Get specific heat capacity based on current temperature
cp = GetSpecificHeat(0.5*(Tc(i)+Tc(i-1)));
% --- Calculate new upper temperature
Tc(i) = Tc(i-1) + Pin/(w*cp);
epsilon = (Tc(i) - Told)/Told;
Told = Tc(i);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cp = GetSpecificHeat(T)
cptbl = [4926 5061 5229 5443 5729 6136 6781 6976 6976]; % J/(kg K)
Ttbl = [538 548 558 568 578 588 598 608 1000]; % K
% --- Interpolate data to current temperature
cp = interp1(Ttbl, cptbl, T);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function k = GetFuelK(T)
if (1 == 2)
% --- Use constant heat conductivity
k = 3.5; % W / (m K)
else
% --- Use modified NFI model (FRAPCON/FRAPTRAN)
% For fresh fuel, 0% gadolinium, UO2
A = 0.0452;
B = 2.46e-4;
E = 3.5e9;
F = 16361;
k = 1/(A + B*T) + E/T^2*exp(-F/T);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function k = GetGasK(T)
if (1 == 2)
% --- Use constant heat conductivity
k = 0.3; % W / (m K)
else
% --- Use power law (FRAPCON/FRAPTRAN)
k = 2.531e-3*T^0.7146;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function k = GetCladK(T)
if (1 == 2)
% --- Use constant heat conductivity
k = 17.0; % W / (m K)
else
% --- Use polynomial relation (FRAPCON/FRAPTRAN)
k = 7.51 + 2.09e-2*T - 1.45e-5*T^2 + 7.67e-9*T^3;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = GetCladH(Tb)
% --- Use constant heat transfer coefficient
h = 2e6/35.0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Pden = GetPowerDensity(Z,R)
global nz nr Rf0 Rf1 dz length Af;
PRod = 1500e6/313/126;
C = PRod/(2*Af);
Pden = zeros(nz, nr);
for i = 1:nz
for j = 1:nr
if (R(j) >= Rf0) && (R(j+1) <= Rf1)
Pden(i,j) = C*cos(0.5*(Z(i) + Z(i+1))/(length/2)*(pi/2));
else
Pden(i,j) = 0.0;
end
end
end
end