kill(all)$

maxima: heat:heat conduction (1D),heat generation 1.0 (g(T)),part3a(solid sphere)
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 sphere

see :appendix 'special functions'

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

assume(k>0,b>0);

→ ode :

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

pb:g=T*b;

depends(T,r);

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

≥ final ode 

pc1:ev(pc,diff);

pc1:expandwrt(lhs(pc1),T);

→ bc : (b)oundary (c)onditions ( solid sphere)

≥ r=0 (center: solid sphere)

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

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

assume(H≥0);

bc2:T(R)=H;

→ heat flow 'through the surface solid sphere'

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

pd1:subst([pb],pd);

solve : model ode solid sphere .
=========================

load(contrib_ode)$

when needed : ?

besselexpand:true;

ma:contrib_ode(pc1,T,r);

mb:rhs(ma[1]);

mb1:diff(T,r)=realpart(diff(mb,r));

Maxima :  %k1,%k2
===============
→ (%k1+%k2) ?   we can overcome the problem : use assume
---------------------------------------------------------------------------------------
result : (%k1+%k2)  = pos → limit (mb1,r→0) = - ∞
            (%k1+%k2)  = neg → limit (mb1,r→0) =  ∞
---------------------------------------------------------------------------------------
result :%k1+%k2 = 0   →  limit (mb1,r→0) =  0
---------------------------------------------------------------------------------------
rem : assume can be used =  ' assume(%k1+%k2 =0), is not allowed

→ solution: extra equation

mb2:%k1+%k2=0;

rem :mb1a:limit(rhs(mb1),r,0);  ' to overcome assume' = mb2

solve : remove : %k2 or %k1 
======================
to find : T(r)  use 'ma'

→ use mb2

pa1:solve(mb2,%k1);

→ remove %k1 ;ma.

pa2:subst(pa1,ma)[1];

→ convert : complex exponential => trig functions  , 'use euler equation=exp(i*θ)'
---------------------------------------------------------------------------------------------------------------
euler equations :            exp(i*θ) = cos(θ)+i*sin(θ)
---------------------------------------------------------------------------------------------------------------

pa3:lhs(pa2)=demoivre(rhs(pa2));

→ general simplify

pa3:ratsimp(%);

→ to find %k2 we need bc2, ' see above' : T(R)=H

pa4:subst([r=R],rhs(pa3)) =H;

pa5:solve(pa4,%k2);

→ subst pa5 into pa3

pa6:subst(pa5,pa3);

rem : H go to zero : bc2

pa7:'limit(lhs(pa6),r,R) = limit(rhs(pa6),H,0);

check ? (b)oundary (c)onditions bci ,i=1,2
================================

→  bc1

pa9:'diff(rhs(pa6),r)=diff(rhs(pa6),r);;/*bc1*/

pa91:'limit(lhs(pa9),r,0)=limit(rhs(pa9),r,0);

→  bc2

pa8:limit(T(r),r,R) = limit(rhs(pa6),r,R);/*bc2*/

heat flow 'through the surface solid sphere' @ r =R
======================================
use :pd1

pa10:pd1;

first step :use  pa9, dT/dr. @ r=R
-----------------------------------

pa101:diff(T,r)=ratsimp(ev(rhs(pa9),diff));

pa102:lhs(pa101) =combine(ev(rhs(pa101),r=R));

second step :   A (solid sphere)
-------------------------------------------

pa103:A=diff(4/3*%pi*r^3,r); /* differentation volume of sphere */

pa1031:subst([r=R],pa103);

final step : 
--------------

pa104:subst([pa102,pa103],pa10);

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