I am trying to solve the 1-D heat equation numerically with a variable source term. The system is basically a tank containing styrene in which it polymerizes to liberate heat. I have assumed that the cross-section of the tank is much smaller compared to the length so any variations in the y and z directions can be neglected. I am starting with $T_0 = 313.15 \text{K}$. The domain space I am looking at is $0 \leq x\leq 5\text{m}$ and $t \geq 0$. The equations are as follows. I write the heat equation as- $$\frac{\partial T}{\partial t} = \frac{k}{\rho C_p}\frac{\partial^2 T}{\partial x^2} + \frac{1}{\rho C_p}\left(\frac{\rho\Delta H}{MW}\frac{\partial m}{\partial t}\right)$$ The $\frac{\partial m}{\partial t}$ is basically the rate of reaction expressed in terms of the mole fraction $m$. The rate of reaction is- $$\frac{\partial m}{\partial t} = -A\left(\frac{\rho}{MW}\right)^{\frac{3}{2}}m^{\frac{5}{2}}$$ Here: $$k = 0.03 \hspace{8pt} \text{[cal/(s-m-K)]}$$ $$A = A_0\exp(A_1m_p + A_2m_p^2 + A_3m_p^3)$$ where: $$m_p = 1-m$$ $$A_0 = 1.964\times10^5\exp\left(\frac{-10,040}{T}\right)$$ $$A_1 = 2.57-5.05\times10^{-3}T$$ $$A_2 = 9.56-1.76\times10^{-2}T$$ $$A_3 = -3.03+7.85\times10^{-3}T$$ The rest of the values are constants and they are- $$\rho = 906 \hspace{6pt} \text {[g/L]}$$ $$MW = 104.15 \hspace{6pt} \text{[g/mol]}$$ $$C_p = 0.4365 \hspace{6pt} \text{[cal/g-K]}$$ $$\Delta H = -17,300 \hspace{6pt}\text{[cal/mole]}$$ The initial conditions are- $$T = T_0$$ $$m = 1$$ The boundary conditions are- $@ x= 0$ $$ \frac{\partial T}{\partial x} = 0, \hspace{5pt} \frac{\partial m}{\partial x} = 0$$ $@ x= 5m$ $$ \frac{\partial T}{\partial x} = 0, \hspace{5pt} \frac{\partial m}{\partial x} = 0$$ I am assuming that there is no flux at either boundary. I used the 'pdepe' command in MATLAB to solve this system of coupled PDEs. Here's the code I wrote-

clear all
global rho Cp MW dH k
tspan = 60*(1:500:30000);
x = 0:0.25:5;
m = 0;
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, tspan);

u1 = sol(:,:,1);
u2 = sol(:,:,2);
xlabel('Distance x')
ylabel('Time t')

function [c,f,s] = pdefun(x,t,u,dudx)
    rho = 906;
    Cp = 0.4365;
    MW = 104.15;
    dH = -17800;
    k = 0.03;
    y2 = 1-u(2);
    A0 = 1.964*(10^5)*exp(-10040/u(1));
    A1 = 2.57-5.05*u(1)*(10^(-3));
    A2 = 9.56-1.76*u(1)*(10^(-2));
    A3 = -3.03+7.85*u(1)*(10^(-3));
    A = A0*exp(A1*(y2) + A2*(y2^2) + A3*(y2^3));
    F = -A*((rho/MW)^(3/2))*((u(2))^(5/2));
    c = [1; 1];
    f = [k*rho/Cp; 0].*dudx;
    s = [(dH/(MW*Cp))*F; F];

function u0 = pdeic(x)
    u0 = [313.15; 1];

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
    pl = [0; 0];
    ql = [1; 1];
    pr = [0; 0];
    qr = [1; 1];

The problem is that on drawing the surface plot of the temperature [u(1)], I don't see any variation with the length of the tank. This is what the surface plot looks like- enter image description here

Another funny thing is that even if I change the function 'f' in 'pdefun' to garbage values like f = [65675; 767].*dudx (these numbers are meaningless, I just typed out something randomly), I get the same plot regardless of the function f. I even cleared all variables before running the code but the same thing happens. I don't understand what's going on here. Any help is appreciated.


You are starting from a uniform temperature and you have insulated boundary conditions; so there is no heat conduction occurring. Likewise, your initial mole fraction is also constant in $x$ so that the heat input is uniform along the length. So the fact that the temperature and mole fraction don't change as a function of $x$ is exactly what you should expect.

Another way of looking at this is: As you have defined it, this problem could be modeled as two coupled ODE with time as the only independent variable.

| cite | improve this answer | |
  • $\begingroup$ That makes sense. I guess this is the reason why changing 'f' doesn't affect the solution? $\endgroup$ – Nishant Kumar Gupta Aug 9 at 18:23
  • $\begingroup$ Correct. $dT/dx=0$ so it doesn't matter what the conductivity is, the heat flux is always zero. $\endgroup$ – Bill Greene Aug 9 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.