Last updated on March 1, 2013
We have discussed the time behavior of the neutron flux and reactivity feedbacks. Now it is time for the the thermal side of things. The point kinetics model describes how much energy is produced in the fuel, but we also need a model for how the energy is transported within the fuel, through the fuel cladding and into the coolant. To figure it out we need models for heat conduction through the fuel and cladding, heat transfer to the fluid through convection, properties of the fluid at different temperatures and pressures and so on. This blog post will deal with the heat conduction in the pellet and through the cladding.
Heat conduction
Heat conduction is the transport of heat in a solid medium, like a fuel rod and the cladding around it. Lets start by making some simplifications, firstly a fuel rod is rather long compared to its diameter. If we split the rod into a bunch of axial zones then it is safe to say that the heat conduction axially in one zone is pretty constant and we can ignore the axial dimension altogether. Secondly power production in a pellet and the composition of the pellet is pretty constant wherever in the pellet you are so you can ignore any angular dependence on the heat conduction. We can assume power production is pretty flat all over the pellet as well. All that is left then is the radial dependence of the temperature.
The temperature distribution is given by the heat equation with internal heat generation and in cylindrical geometry:
$$\large C_p\rho\frac{dT}{dt} = k\nabla^2T+Q=\frac{k}{r}\frac{d}{dr}(r\frac{dT}{dr})+Q$$
$$\large \frac{dU}{dt}$$ describes how temperature changes with time, $$\large Q$$ is power density(watts per cubic meter), $$T$$ is temperature,$$\rho$$ is material density, $$C_p$$ is specific heat capacity and $$\large k$$ is heat conductivity.
This equation is of exactly the same form as the neutron diffusion equation (heat diffuses after all). But for the heat we cant neglect the $$\large \nabla^2$$ term because we want both the time and space dependencies, not only the time dependence like with the point kinetics model.
Analytic treatment
This equation has a nice analytical solution if steady state is assumed, i.e $$\large \frac{dT}{dt}=0$$.
Then we have:
$$\large\frac{d}{dr}(r\frac{dT}{dr})=-\frac{Qr}{k}$$
Lets say we have a fuel rod of radius a within a cladding tube of outer radius b. We can call the pellet for zone 1 and the cladding for zone 2 (lets ignore the gas gap that usually exists between the fuel and cladding because the pellet swell anyway during operation, closing the gap). We then have to solve the above equation for each zone. Zone 2 is a bit simpler because there is no power produced in the cladding which means $$Q=0$$. We integrate the equations twice and end up with these two equations:
for $$r<a$$:
$$\large T_1(r)=\int\int\frac{d}{dr}(r\frac{dT_{1}}{dr})drdr = -\frac{Qr^2}{2k_{f}}+C_{1,1}ln(r)+C_{1,2}$$
for $$a<r<b$$:
$$\large T_2(r)=\int\int\frac{d}{dr}(r\frac{dT_{2}}{dr})drdr =C_{2,1}ln(r)+C_{2,2}$$
To solve it we need some boundary conditions. We can figure out these three boundary conditions.
B.C 1 $$\large -k_f\frac{dT_1}{dr}|_{r=0}=0$$ i.e heat flux at the center is zero.
B.C 2 $$\large -k_f\frac{dT_1}{dr}|_{r=a} = -k_c\frac{dT_2}{dr}$$ i.e heat flux in the pellet-cladding boundary is continous
B.C 3 $$\large -k_c\frac{dT_2}{dr}|_{r=b}=q”(b)=h(T_b-T_m)$$, this relates the heat flux at the cladding surface $$t_c$$, with the temperature of the coolant, $$T_m$$. $$h$$ is the heat transfer coefficient (we will return to it when we discuss the coolant model).
The equations of the form $$\large -k\frac{dT}{dr}$$ above is Fouriers law that relates heat flux to temperature gradient.
If we do some gymnastics with these boundary conditions we end up with these two solutions:
for $$r<a$$:
$$\large T_1(r) =\frac{Q}{4k_f}(a^2-r^2)+\frac{Qa^2}{2}[\frac{1}{k_c}ln(\frac{b}{a})+\frac{1}{hb}]+T_m$$
for $$a<r<b$$:
$$\large T_2(r)=\frac{Qa^2}{2}[\frac{1}{k_c}ln(\frac{b}{r})+\frac{1}{hb}]+T_m$$
If we put in some sane numbers for the heat conductivity and heat transfer coefficients we get a result that looks like this:
About 400 degrees drop in temperature between the center of the fuel pellet and the surface of the cladding. A very pretty plot indeed, but this blog series is about transient behavior, not steady state! Why did I just devote a page to an analytic steady state solution? Well it turns out that the transient numeric solution requires an initial condition and the analytic solution is perfect for that. We always start form a steady state when we do transient calculations.
Numeric treatment
There might exist some analytic solution to the time dependent heat equation in this kind of geometry and with the simplifications I made, but I am not aware of it (and I can’t even remember the more advanced solution methods). I want to make the heat conductivity and specific heat capacity coefficients temperature dependent anyway and that forces us to look at numerical methods instead. Imagine that we divide our fuel pin into radial zones and as before there is no angular dependence. Inner zone is number 1 and outer zone is number N. The zones are sized such that a zone border coincides with the border between pellet and cladding. The zones within the pellet all have the same radial size and so does the zones in the cladding (but cladding zone and pellet zone can be different from eachother).
We want to formulate an equation to calculate what is going on in zone i (yes it is n, not i, in the drawing above but I am to lazy to change it so live with it!), and three equations that contains the boundary conditions for the center , the pellet-cladding interface and the cladding-coolant interface. Lets start with a generic zone that is not at one of the boundaries, we will use the method described in Boundary value problems of heat conduction.
First lets look at the heat equation again and write out the laplacian.
$$\large \frac{C_p\rho}{k} \frac{dT}{dt} =\nabla^2T+\frac{Q}{k}=\frac{1}{r}\frac{d}{dr}(r\frac{dT}{dr})+\frac{Q}{r} = \frac{1}{r}\frac{dT}{dr}+\frac{d^2T}{dr^2}+\frac{Q}{r}$$
If some teacher ever told you that the Taylor series is the solution to every problem he/she was quite correct (if you did not hear that the teacher clearly was guilty of negligence) and we will start with Taylor.
Lets say we have a function that describes temperature in a point in time and space, $$\large T_{i,j}$$, where i denotes space and j time. We can do a Taylor expansion around that point in positive and negative direction:
$$\large T_{i+1,j} \approx T_i + \Delta r\frac{dT_{i,j}}{dr}+\frac{(\Delta r)^2}{2}\frac{d^2T_{i,j}}{dr} …$$
$$\large T_{i-1,j} \approx T_i – \Delta r\frac{dT_{i,j}}{dr}+\frac{(\Delta r)^2}{2}\frac{d^2T_{i,j}}{dr} …$$
We neglect terms in the Taylor expansion of order higher than 2. Lets multiply the first expression by A on both sides and the second expression by B. Lets then add the expression together and bunch the derivatives on one side and the non derivatives on the other. Then we end up with:
$$\large (A-B)\Delta r\frac{dT_{i,j}}{dr} + (A+B)\frac{(\Delta r)^2}{2}\frac{d^2T_{i,j}}{dr^2} = A(T_{i+1,j}-T_i)+B(T_{i-1,j}+T_n) $$
The left side of the expression is an expression describing the change of temperature over space, but that is what the laplacian does in the heat equation! We do a trick and assume that the left hand side of the above expression is approximately identical to the laplacian in the heat equation.
$$\large (A-B)\Delta r\frac{dT_{i,j}}{dr} + (A+B)\frac{(\Delta r)^2}{2}\frac{d^2T_{i,j}}{dr^2} \approx\frac{1}{r}\frac{dT}{dr}+\frac{d^2T}{dr^2}$$
By identification we then see that:
$$\large (A-B)\Delta r = \frac{1}{r}$$
and
$$\large(A+B)\frac{(\Delta r)^2}{2} = 1$$
From those two relations we can easily find A and B. We now replace the laplacian in the heat equation with the right hand side of the Taylor expansion expression(with A and B replaced by their values).
$$\large \frac{C_p\rho}{k} \frac{dT}{dt} =\frac{T_{i-1,j}+T_{i+1,j}+2T_n}{(\Delta r^2)}+ \frac{1}{r_n}\frac{T_{i+1,j}-T_{i,j}}{2\Delta r}+\frac{Q}{k}$$
If we also simply Taylor expand the time derivative:
$$\large\frac{dT}{dt} \approx \frac{T_{i,j+1}-T_{i,j}}{dt} $$
Plug that into the heat equation and then put $$\large T_{i,j+1}$$ on the left hand side and everything else on the right hand side we get this equation.
$$\large T_{i,j+1} = T_{i,j} + \frac{k\Delta T}{\rho C_p}[\frac{T_{i-1,j}+T_{i+1,j}+2T_n}{(\Delta r^2)}+ \frac{1}{r_n}\frac{T_{i+1,j}-T_{i,j}}{2\Delta r}] + Q\frac{\Delta t}{\rho C_p}$$
There we have it, a way to find the temperature distribution at the next time step based on the distribution at the current time step. Now lets look at the boundary conditions, starting at the cladding-coolant boundary. Lets say the boundary is on the surface node b.
B.C 1: $$\large -k_c\frac{dT_2}{dr}|_{r=b}=h(T_b-T_m)$$
We want to use the Taylor expansion again, but this time we cant expand to $$T_{N+1,j}$$ because the node N+1 is outside of the cladding(N is the surface of the cladding). Instead we only have the expansion of $$T_{N-1,j}$$. If we do the same procedure as above but this time multiply the boundary condition by A and the expansion by B and sort all derivatives on one side we get:
$$\frac{B(\Delta r)^2}{2}\frac{d^2T}{dr^2}-(A+B(\Delta r)\frac{dT}{dr})=\large A(\frac{h}{k}(T_{i,j}-T_m)+B(T_{N-1,j}-T_{N,j}) $$
Then identify the left hand side with the laplacian, solve A and B and Taylor expand the time derivative. We then end up with this:
$$\large T_{N,j+1} = T_{N,j} + \frac{k\Delta t}{\rho C_p}[\frac{h}{k}(T_{N,j}-T_m)+\frac{2}{(\Delta r)^2}(T_{N-1,j}-T_{N,j})] $$
There we have a way to calculate the temperature at the cladding-coolant boundary. Now we have the pellet-cladding boundary to sort out. Here we can straight away fiddle together a energy balance. The energy balance is:
the energy stored in the outermost radial zone of the pellet = the energy flow into the outermost zone of the pellet – the energy flow out of the pellet into the cladding + the energy produced in the outermost radial zone of the pellet .
Expressed mathematically we have:
$$\large \rho_p C_{p}\frac{T_{a,j+1}-T_{a,j}}{\Delta T}\pi(r^2_a-r^2_{a-1}) = -2k_p\pi r_{a-1}l\frac{T_{a,j}-T_{a-1,j}}{\Delta r} + 2k_c\pi r_al\frac{T_{a+1,j}-T{a,j}}{\Delta r} \\+ Q\pi(r^2_a-r^2_{a-1}) $$
From this we can right away find what we want:
$$\large T_{a,j+1} = T_{a,j} +\frac{2\Delta t}{\rho_p C_p\Delta r (r^2_a-r^2_{a-1})}[k_cr_a(T_{a+1,j}-T_{a,j})-k_pr_{a-1}(T_{a,j}-T_{a-1,j})]\\ +Q\frac{\Delta t}{\rho_p C_p}$$
Now we are almost finished, we only need to figure out the boundary condition for the center of the pellet. Here we have the boundary condition:
$$\large -k_p\frac{dT_{o,j}}{dr} = 0$$
and the taylor expansion
$$\large T_{1,j} = T_{0,j} +\Delta r\frac{dT_{0,j}}{dr}+\frac{(\Delta r)^2}{2}\frac{d^2T_{0,j}}{dr^2}$$
Do the previous trick of multiplying the boundary condition by A and the expansion by B and sort the derivatives on the left hand side, then identify the left hand side with the laplacian.
$$\large (\Delta r B-Ak_p)\frac{dT_{0,j}}{dr}+B\frac{(\Delta r)^2}{2} = B(T_{1,j}-T_{0,j})$$
We are talking about the center, r=0. If anyone remember there is a rule, L’Hospitals rule, with which one can show that:
$$\large \lim_{r \to 0}\frac{1}{r}\frac{dT}{dr} = \frac{d^2T}{dr^2} $$
and our laplacian simply turns into $$\large 2\frac{d^2T}{dr^2}$$ which means $$\large B=\frac{4}{(\Delta r)^2}$$.
Fiddle around with it and expand the time derivative of the heat equation and we get.
$$\large T_{0,j+1} = T_{0,j}+4\frac{k_p}{\rho_p C_p}\frac{T_{1,j}-T_{0,j}}{(\Delta r)^2}+\frac{Q}{\rho_p C_p}$$
Well that was a proper handful! Now its time to connect the above with the point kinetics equation from the last blog post. I have done all that in a small program that can be downloaded here, the executable is called transa.exe. The rest of the files are the license terms and a readme file that describes the input options and how to run the file on linux. In the folder source one can find the source code for the executable. If we plug in 1400 pcm of reactivity insertion and the rest as “typical pwr” parameters we get this nice animation. The top window shows the fuel temperature in the pellet and cladding for each time step and the bottom window the power level. We see that even though the power jumps dramatically the fuel temperature just goes up by about 400 degrees.
Now we have the neutron interactions and the heat transfer in the fuel fixed. Next blog post will be about the water and how it effects the $$\large h$$ parameter and the moderator feedbacks.
/Johan