%{
# Graetz solution
This is the thermal entrance problem. A classical problem in thermal flow: internal foreced convection. We have a Poiseuille flow of a given temperature in a pipe which wall is at that temperature, and at $x=0$ there is a sudden change of the pipe wall to another temperature. We are interested in the description of the way the temperature of the fluid will adjust
to that of the pipe. We assume the solution in the form of a series of functions of separated variables, and applying the equation to each term of this series gives a set of eigenvalue/eigenvector solutions. This is what we solve here.
%}
format compact
clear all; clf
% if octave
%setenv("GNUTERM","X11")
clf
%%%% parameters
Ny=100; %number of grid nodes in y
Ly=1; % domain size in y
%%%% chebychev in y
scale=-Ly/2;
[y,DM] = chebdif(Ny,4);
y=(y-1)*scale;
dy=DM(:,:,1)/scale;
dyy=DM(:,:,2)/scale^2;
I=eye(Ny);
%{
$$ u_{Pois}\frac{\partial T}{\partial x} = (
\frac{ 1}{ r}\frac{\partial }{\partial r}r\frac{\partial T}{\partial r})
$$
with $T(x,1)=0$ and $\partial_r T(x,0)=0$.
in $x=0$, $T(0,r)=1$
We look at the solution in term of a series of a function of $x$ times a function of $r$; it is straightforward to see that
by the $x$ function is an exponential, which must decrease
$$T =\Sigma e^{-\alpha_i^2 x } c_i\theta_i(r)$$
and $u_{Pois}=(1-r^2)$,
so coefficients $c_i$ will be obtained latter by $1= \Sigma_i c_i e^{-\alpha_i^2 x } \theta_i(r)$
The problem to solve is
$$
s (1-r^2)\theta_i(r)=(\frac{\partial^2 }{\partial r^2} + \frac{ 1}{ r}\frac{\partial }{\partial r})\theta_i(r)
$$
were $s=-\alpha_i^2$.
%}
% System matrices
A=dyy+diag(1./y)*dy;
E=-diag(1-y.^2);
%{
We can rewrite this system of equations in a matrix form
$$
s_i E\theta_i =A \theta_i
$$
$(\frac{\partial^2 }{\partial r^2} + \frac{ 1}{ r}\frac{\partial }{\partial r})$
is A, and $-(1-r^2) I $ is E
attention
$y(0)=0$ and $y(Ny)=1$
$\theta_i(r=1)=0$ and $\theta_i'( r=0)=0$
r=0 is y=0 so position 1 and r=1 is y=1 position Ny
and with the boundary conditions $\theta_i(1)=0$ and $\partial_r \theta_i(0)=0$.
%}
% boundary conditions
loc=[1,Ny];
C=[dy(1,:); I(Ny,:)];
E(loc,:)=0;
A(loc,:)=C;
%{
# compute eigenvalues
%}
[U,S]=eig(A,E);
s=diag(S);
rem=abs(s)>1e5; s(rem)=[];U(:,rem)=[];
[t,o]=sort(real(s)); s=s(o);U=U(:,o);
%{
Note that the eigen value have been sorted.
# show Eigen functions
%}
plot(y,U(:,1:5),'-');legend('1','2','3','4','5');
%figure
%plot(y,dy*U(:,1:5),'-');legend('1','2','3','4','5');
% figure
% plot(y,dyy*U(:,1:5),'-');legend('1','2','3','4','5');
%{
#Results
Here the five first values of $\alpha_i$
2.7044
6.6790
10.6734
14.6711
18.6699
%}
disp('the first eigen values ');
sqrt(s(1:5))
%{
# Links
# Exercises
* Do the same in plane 2D
* Do the same for the unsteady heat equation
$$ \frac{\partial T}{\partial t} = (
\frac{ 1}{ r}\frac{\partial }{\partial r}r\frac{\partial T}{\partial r})
$$
with $T(t,1)=0$ and $\partial_r T(t,0)=0$.
in $t=0$, $T(0,r)=1$
#Bibliography
* A. Leontiev (1985) ”Théorie des échanges de chaleur et de masse”, ed. MIR.
%}