10.3 Heat Conduction in a Sphere with a Nuclear Heat Source > restart; > A:=r-> 4*Pi*r^2; area of conduction in a spherical shell > Sn := r->Sn0*(1+b*(r/RF)^2); define the volumetric nuclear energy source: Sn as in eq. 10.3-1 > Q:= r->qr(r)*A(r); eq. 10.3-2 giving the thermal energy flowing across the surface of a sphere of rarius r > LHS:=limit((Q(r+dr)-Q(r))/dr,dr=0); Left hand side of eq. 10.3-5 > de1:=simplify((LHS=Sn(r)*A(r))/(4*Pi)); Equivalent to eq. 10.3-6 > de2:=simplify((LHS=0)/(4*Pi)); The DE for the cladding > s:=dsolve({de1,qr(0)=finite},qr(r)); solving the ODE and using the BC 10.3-10 > assign(s); We must assign the result given by dsolve and > qF:=unapply(qr(r),r); use unapply to make the flux in the fissionable region: qF a function. This agrees with 10.3-12 > qr:='qr'; We must get rid of the value assigned to qr to use it in the next step. > s:=dsolve({de2,qr(RF)=qF(RF)},qr(r)); Now we can find qr in the cladding using 10.3-7 and the BC 10.3-11 > assign(s);qC:=unapply(qr(r),r); eq. 10.3-13 > de3:= -kF*D(TF)(r)=qF(r); de4:=-kC*D(TC)(r)=qC(r); equations 10.3-14 and 10.3-15: applications of Fourier's Law > s:=dsolve({de4,TC(RC)=T0},TC(r)); assign(s); solve the ODE for the cladding with BC 10.3-19 > TC:=unapply(TC(r),r); > dTC:=r->simplify(TC(r)-T0); This agrees with eq. 10.3-21 as seen by: > dTCbook:=Sn0*RF^2/(3*kC)*(1+3*b/5)*((RF/r)-(RF/RC)); > simplify(dTC(r)-dTCbook); > s:=dsolve({de3,TF(RF)=TC(RF)},TF(r)); assign(s); solve 10.3-14 with BC 10.3-18 > TF := unapply(TF(r),r); > dTF:=r->simplify(TF(r)-T0);dTFbook:=Sn0*RF^2/(6*kF)*((1-(r/RF)^2)+3*b/10*(1-(r/RF)^4))+Sn0*RF^2/(3*kC)*(1+3*b/5)*(1-RF/RC);differ:=simplify(dTF(r)-dTFbook); The difference TF(r) - T0 we found compared to the result in eq. 10.3-20. > with(plots): > tauC:=rho->simplify(dTC(RC*rho)/(Sn0*RF^2/3/kC)); rho is the dimesionless variable: r/RC, tauC is (TC-T0)*3*kc/(Sn0*RF^2) > tauC(rho); > RF:=alpha*RC;tauC:=unapply(tauC(rho),rho,alpha,b); alpha=RF/RC > tauF:=rho->simplify(dTF(RC*rho)/(Sn0*RF^2/3/kC)); tauF is (TF-T0)*3*kc/(Sn0*RF^2) > tauF(rho); > kF:=beta*kC; beta=kF/kC > tauF:=unapply(tauF(rho),rho,alpha,beta,b); > p1:=plot(tauF(rho,.75,1.5,.2),rho=0...0.75,title=`(T-T0)*3*kc/(Sn0*RF^2) vs r/RC`): alpha=0.75, beta=1.5, b=0.2, rho goes from 0 to 0.75 > p2:=plot(tauC(rho,.75,.2),rho=0.75....1.0): alpha=0.75, b=0.2, rho goes from .75 to 1 > display({p1,p2}); Note the change in slope at rho=0.75
> ?display >