kill(all)$

heat : initial and boundary conditions  (ic,bc),version 1a, final
peter.vlasschaert@gmail.com,08/07/2019

ic : initial condition        ( t=0,begin heat conduction: mathematical form)
bc: boundary condition ( t≠ 0, physical effects:mathematical form)

depends(T,r);

declare(n,integer);

info : T : temperature
          k:thermal conductivity (also function of r)
          r: radius
          n: integer  (0≤n≤2)
          g(r): function of heat generation.

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

r: x  one dimensional,n=0 

steady : independent of time
======================

p2:ev(p1,r=x,n=0);

kill(T);

dependencies;

depends(T,[x,t]);

unsteady : dependent of time
======================

p2a:g(x)/k+'diff(T,x,2)= z*diff(T,t);

z= ? need find from unit check on both sides.

p2b:expand(k*p2a);

dependencies;

rem : we need use always 'p2a,p2b',because
of the dependencies.

1e) steady : only lhs (p2a)
-----------------------------------

p3:lhs(p2a) = 0;

a)  prescribed temperature bc
=======================

example wall : thickness = L
GIVEN :  x=0 , left      =>  T  =TL
                x=L , right    =>  T = TR
g(x) = a*x+b (heat generation inside)

p4:ode2(p3,T,x);

1e)  x=0 , left      =>  T  =TL
-------------------------------------

p4a1:subst([g(x)=a*x+b,T=TL],p4);

p4a2:ev(p4a1,nouns);

p4a3:subst([x=0],p4a2);

2e)  x=L , left      =>  T  =TR
-------------------------------------

p4b1:subst([g(x)=a*x+b,T=TR],p4);

p4b2:ev(p4b1,nouns);

p4b3:subst([x=L],p4b2);

3e)  find %k1,%k2
--------------------------

p5a:solve([p4a3,p4b3],[%k1,%k2]);

p5a1:flatten(p5a);

4e)  find  T (x) ?
---------------------

p6a1:subst([g(x)=a*x+b],p4);

p6a2:ev(p6a1,nouns);

p6a3:subst(p5a1,p6a2);

b) prescribed heat flux boundary conditions
=================================
gradient  @ boundary.

definition qx  (heat flux)
==========

p7:q[x] = - k *'diff(T,x);

left face of the wall   :q[0],inside the wall direction. 'use p7'
right face of the wall :q[L],outside wall direction.    'use p7'

q(x) = a*x+b
use : from above. 'p6a2'

1e) x=0 , left      =>  q[0]
---------------------------------

p7:p6a2;

  => (dT/dx)  

p7a1:diff(p7,x);

  => -k*(dT/dx)

p7a2:-k*p7a1;

p7a3:q[0] =expand(subst([x=0],rhs(p7a2)));

p7a4:solve(p7a3,%k2);

p7a5:subst(p7a4,p7);

p7a6:ratsimp(p7a5);

2e) x=L , left      =>  q[L]
---------------------------------

=> (dT/dx)  

p7b1:diff(p7a5,x);

  => -k*(dT/dx)

p7b2:-k*p7b1;

p7b3:expand(p7b2);

p7b4:q[L] =subst([x=L],rhs(p7b3));

how to use bc2 for 'p3'
use : g(x) =a*x+b
=================

p7c:subst([g(x) = a*x+b],p3);

p7c1:ode2(p7c,T,x);

find : %k1,%k2
============

p7c2a:bc2( p7c1,x=0,-k*(diff(T,x)) =q0,x=L,-k*(diff(T,x)) =qL);

p7c2b:bc2( p7c1,x=0,-k*('diff(T,x)) =q0,x=L,-k*('diff(T,x)) =qL);

rem : q[0]  => q0,q[L] => qL

test : p7c2a ? p7c2b
================

p7c3:is(rhs(p7c2a)=rhs(p7c2b));

find : T(x)
=======

p7c4:solve(p7c2,T);

=>   no solution found for  T=T(x)

dT/dx = 0 ,insulated boundary

Thermal symmetry
==============
centre plane      : dT/dx |  x= L/2 ,'plane symmetry'  => insulated boundary
centre cylinder   :dT/dr  |  r=0

c) convection bc : surface energy balance
================================
condition on both surfaces
hL = heat transfer coef. on the left
hR = heat transfer coef. on the right
--------------------------------------------------------=
convection flux : fluid to surface

@ x=0  surface (left)

p8a1:hL*(TinfL-Tx0) = -k*'diff(T,x);

@ x=L surface  (right)

p8a2:hR*(TxL-TinfL) = -k*'diff(T,x);

c) radiation bc 
================================
condition on both surfaces
εL = emissivities coef. on the left
εR =  emissivities coef. on the right
σ   = stefan boltzmann constant = 5.67 W/(m^2*K^4)
inf = ∞
L = left
R = right
--------------------------------------------------------=
radiation boundary : both sides

@ x=0  surface (left)
left : radiation => conduction
'conduction : insde'

 -k*'diff(T,x) =σ*εL*(TinfL^4-Tx0^4);

@ x=L  surface  (right)
left : conduction =>radiation
'conduction : insde'

 -k*'diff(T,x) =σ*εR*(TxL^4-TinfR^4);