sandbox/easystab/graetz.m

    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);

    \displaystyle 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 \displaystyle 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 \displaystyle 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 \displaystyle 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 \displaystyle \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.