kill(all)$

maxima: heat:heat conduction (1D),heat generation 1.0 (g(T)),part2a(solid cylinder)
peter.vlasschaert@gmail.com,01/08/2021

properties:

k: thermalconductivity ('fixed')
g : heat generation  (=g(T))
H: fixed temperature
T :temperature
r : location in the solid cylinder(with  length   =  L)

see :appendix 'special functions'

model :solid cylinder (r), 0≤ r ≤R
=====

assume(k>0);

→ ode :

pa:1/r*'diff(k*r*'diff(T,r),r)+g=0;

pb:g=b*T;

depends(T,r);

pc:subst([pb],pa/k);

≥ final ode 

pc1:ev(pc,diff);

pc1:ratexpand(pc1);

→ bc : (b)oundary (c)onditions (L =length of the  solid cylinder)

≥ r=0 (begin : center of cylinder)

bc1:diff(T(0),r)=0;

≥ r=R (end : to surface of the solid cylinder)

assume(H≥0);

bc2:T(R)=H;

→ heat flow 'through the surface solid cylinder'

pd:q=-k*A*diff(T,r);

pd1:subst([pb],pd);

solve : model ode solid cylinder .
=========================

load :packages: bessel functions

load(contrib_ode)$

m:odelin(pc1,T,r);

m1:setp(m);/*check type of m of set*/ 

ma:contrib_ode(pc1,T,r);

m1a:listp(ma);/*check type of m of list*/ 

m2:ma[1];

→ differentation T with respect r

md:diff(T,r)=subst([m2],diff(T,r));

md1:ev(md,nouns);

md2:expand(md1);

listofvars(md2);

use :{bc1,bc2}

→ bc1 (r=0,dT/dr=0), 'problem with infinity'

bessel_y: bessel_y(1,0) is undefined.
=============================
" an error. To debug this try: debugmode(true)"  
what to do : %k2 = 0

bc1:subst([%k2=0],md2);

bc1a:lhs(bc1)=subst([r=0],rhs(bc1));

bc1b:lhs(bc1)='limit(rhs(bc1),r,0);

bc1b:ev(bc1b,nouns);

m2a:subst([%k2=0],m2);

→ bc2(r=R, T(R)=H), H=fixed

bc2:subst([r=R],rhs(m2a))=H;;

→ find

p1:solve([bc2],%k1);

fnal solution : ode+bc
================
use :m2a

p2:subst(p1,m2a);

T0:temperature along axis (r=0)

p2a: limit(T0,r,0)= subst([r=0],rhs(p2));

TR:temperature on the surface solid cylinder (r=R)

p2b: TR= subst([r=R],rhs(p2));

q: heat transfer loss through the surface(r=R)
=========================================

first step :differentation p2

p3:'diff(T,r)=diff(rhs(p2),r);

second step: heat flow 'through the wall'

p4:ev(pd1);

p4a:subst([p3],p4);

 :A=2*π*R*L, 'surface area: cylinder length L'

A:2*%pi*R*L;

p4b:fullratsimp(ev(p4a));

p4c:subst([r=R],p4b);

what about :check condition bc1:
---------------------------------------------
use :p3

p5:'limit(lhs(p3),r,0)=limit(rhs(p3),r,0);

what about :check condition bc2:
---------------------------------------------
use :p2

p5a:'limit(T(r),r,R)=limit(rhs(p2),r,R);

rem : T → T(r) ,otherwise symbol limit doesn't exit on the lhs

-----------------------------------------------------------------------------------------------------------------------------------------------

? What about H = 0 when  T(r=H)=0 