kill(all)$

heat equation : cylindrical coordinate system , version 1.0
peter.vlasschaert@gmail.com,10/07/2019

load packages
============

load(lrats);

elementary cylindrical volume element
=============================
Δv = r*Δθ*Δr*Δz
=============================

rem :3d, fourier law , cartesian coordinates
=================================

Q(x)=-A*('diff(T,x,1))*k,   A= A(z,y) = Δz*Δy
Q(y)=-A*('diff(T,y,1))*k,   A =A(x,z) = Δx*Δz
Q(z)=-A*('diff(T,z,1))*k,   A =A(x,y) = Δx*Δy
m = ρ*Δx*Δy*Δz
g(x) = g  (= W/m^3)
=================================
rem : A perpendicular to the direction.
=================================
error : corrected below

see: website
===================================================================
heat : equation 1D with volumetric heat generation R^n:(R-direction,n=0,1,2),version 1.0
n=0,cartesian
n=1,cylindrical
n=2,spherical
===================================================================

definition : taylor series
==================
'.'  => non commutative multiplication for taylor expansion.

p1:m(f,x,Δx,n) :=sum('diff(f,x,i).Δx^i/i !,i,0,n);

first order: approximation.
-----------------------------------
x=r 

p2:m(f,r,Δr,1);

fourier law : cylindrical coordinates
===========================
(x,y,z) => (r,θ,z)
rem : θ     :no dimension.
          r.θ  : dimension.

p3a:Q(r)= -k*r*Δθ*Δz*'diff(T,r);          /*θ-z plane*/
p3b:Q(θ) = -k*Δr*Δz*1/r*'diff(T,θ);    /*r-z plane*/
p3c:Q(z) = -k*r*Δθ*Δr*'diff(T,z);   /*r-θ plane*/

p4a:Q(r+Δr) = m(Q(r),r,Δr,1);
p4b:Q(θ+Δθ) = m(Q(θ),θ*r,Δθ*r,1);
p4c:Q(z+Δz) = m(Q(z),z,Δz,1);

rem :Q(θ+Δθ)=(1/r*'diff(Q(θ),θ))*r*Δθ+Q(θ)

we need modified 'p4b',see : rem 'above'

p4b1:lhs(p4b) = lratsubst([r*θ =v2],part(rhs(p4b),1))+part(rhs(p4b),2);

p4b2:['diff(Q(θ),v2)=1/r*diff(Q(θ),θ)];

p4b3:lratsubst(p4b2,p4b1);

net rate of heat conduction.= Qnet
=====================

p5a:Qnet = lhs(p3a+p3b+p3c)-lhs(p4a+p4b+p4c);

p5b:Qnet=-(part(rhs(p4a),1)+part(rhs(p4b3),1)+part(rhs(p4c),1));

p5c:factor(p5b);

use :fourier law ,{p3a,p3b,p3c}

p5d:subst([p3a,p3b,p3c],p5c);

trick:;
Δv = r*Δθ*Δr*Δz

p5d1:[Δr=1,Δθ=1,Δz=1];

p5d2:lratsubst(p5d1,p5d);

p5d3:lhs(p5d)=subst([r=1],part(rhs(p5d2),1))*Δv+1/r*part(rhs(p5d2),2)*Δv+part(rhs(p5d2),3,1,1)*Δv;

volumetric heat generation = Qgen
====================

p6a: r*Δθ*Δr*Δz =Δv;

p6b:Qgen=g*rhs(p6a);

net rate change internal energy = Qv
=======================

p7a:Qv = ρ*Δv*C*'diff(T,t);

Energy balance on element
=====================

p8a:lhs(p7a) = lhs(p5d3+p6b);

p8b:rhs(p7a) =rhs(p5d3+p6b);

p8c:expand(p8b/Δv);

by definition : α 'thermal diffusivity'
==========================

p9: k/(ρ*C)=α;

use :p9, k = constant

p8d:lhs(p8c)=k*subst([k=1],part(rhs(p8c),1)+part(rhs(p8c),2) +part(rhs(p8c),3))+part(rhs(p8c),4);

p8e:expand(p8d/k);

heat conduction equation in cylindrical coordinates
=======================================

p8f:lratsubst([p9],lhs(p8e)) = rhs(p8e);