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-
clc
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);
surf(x,tspan,u1)
title('u_1(x,t)')
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];
end
function u0 = pdeic(x)
u0 = [313.15; 1];
end
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
pl = [0; 0];
ql = [1; 1];
pr = [0; 0];
qr = [1; 1];
end
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-
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.
