kill(all)$

heat : equation 1D with volumetric heat generation R^n:(R-direction,n=0,1,2),version 1.0
peter.vlasschaert@gmail.com,10/07/2019

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

load(lrats);

see : use p10b and below , lratsubst([],expr)

a)  taylor :expansion
----------------------------

p1:taylor(f(x),x,x0,2);

p1a: x-x0=Δx;

p1b:subst([p1a],p1);

b)  f(x+Δx)
---------------

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

zero order approximation
-----------------------------------

p2a:m(f,x,Δx,0);

first order :approximation.  
-----------------------------------

p2b:m(f,x,Δx,1);

rem :  Qt = C*m*dT/dt  (Watt = J/s)  = >  Qv =ρ*dT/dt*Δx*Δy*Δz ==ρ*dT/dt*Δv (= Watt/m^3)
=====================================================================
T : temperature    (K)
t  : time                  (s)
C : heat capacity  (J/kg*K), specific heat of material
m : mass               (Kg)          
ρ  : density            (Kg/m^3)
Δx = Δy = Δz         (m)
Δv                           (m^3)
=> m = ρ*A*dx    ' x-direction'
====================================================================

Qnet : Net rate of heat gain by conduction
================================

p3:Qnet = m(Q(x),x,Δx,0)-m(Q(x),x,Δx,1);

Qgen:  heat generation
=======================

p4:Qgen = g(x)*A*Δx;

rem : g(x) : volumetric heat generation (Watt/m^3)

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

p5:Qv = ρ*C*A*Δx*'diff(T,t);

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

p6a:Qnet+Qgen = Qv;

p6b:rhs(p3+p4) = rhs(p5);

definition : fourier law : Q(x)
=====================

p7:Q(x) = -k*A*'diff(T,x);

k : function of direction (W/m*K)
k:thermal conductivity

use : fourier law  into energy balance
--------------------------------------------------

declare(k,constant);

p6c:subst([p7],p6b);

p6d:expand(p6c/(A*Δx*k));

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

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

use :α into 'p6d'

p6e:lhs(p6d) =lratsubst([p8],rhs(p6d));

because : complex conversion no dependencies define.
=> use trick

a)  A=A(x)
========

p6f:part(lhs(p6e),1) +subst([k=1],part(lhs(p6e),2) )= rhs(p6e);

rem :p6f : A=A(x)

b) A = constant
============

p6g:part(lhs(p6f),1) +subst([A=1],part(lhs(p6f),2) )= rhs(p6f);

p9:-'diff((-'diff(T,x,1)),x,1) = 'diff(T,x,2);

p6h:lratsubst([p9],p6g);

p6i:expand(p6h);

cylindrical coordinate system
======================
cylinder (A) = 2*π*r*L

use : trick from above

p10: A= 2*π*r*L;

use : p6f

p10a:lratsubst([p10,x=r],p6f);

p10a1:expand(p10a);

p10b::part(lhs(p10a1),1) +lratsubst([2*π*L=1],part(lhs(p10a1),2) )= rhs(p10a1);

spherical coordinate system
======================
cylinder (A) =4*π*r^2

p11:volume= 4/3*π*r^3 ; /*volume sphere*/

p11a: 'diff(volume,r) = diff(rhs(p11),r);  /*area full sphere*/

p11b:A= rhs(p11a);

use : p6f

p11c:subst([p11b,x=r],p6f);

p11d:part(lhs(p11c),1) +lratsubst([4*π=1],part(lhs(p11c),2) )= rhs(p11c);

general equation
=============

trick: r^2 = R^n

p11e:subst([r^2=R^n],p11d);

cartesian coordinates.
=================
R=x , n=0

p11e1:subst([R=x,n=0,r=x],p11e);

g(r) = volumetric generation.

cylindrical coordinates..
==================
R=r, n=1

p11e2:subst([R=r,n=1],p11e);

spherical coordinates.
=================
R=r, n=2

p11e3:subst([R=r,n=2],p11e);

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,x,1))*k,   A =A(x,z) = Δx*Δz
Q(z)=-A*('diff(T,x,1))*k,   A =A(x,y) = Δx*Δy
m = ρ*Δx*Δy*Δz
g(x) = g  (= W/m^3)
=================================
rem : A perpendicular to the direction.
=================================