Stuck on a hideous differential equation

I'm stuck on this PDE: $$\rho \cdot C \cdot \dfrac{\partial T}{\partial \tau} + u \cdot \rho \cdot C \cdot \dfrac{\partial T}{\partial r} = \dfrac{\lambda}{r^{2}} \cdot \dfrac{\partial}{\partial r}\left(r^{2} \cdot \dfrac{\partial T}{\partial r}\right) + K - \dfrac{\rho \cdot C \cdot T}{V} \cdot \dfrac{\partial V}{\partial \tau}$$ This should describe the variation of enthalpy during a flash pyrolysis. I know all the parameters ($\rho, C, \lambda, K$) and I can obtain $u$ and $V$ from analytic expression.

Assuming I know the initial condition and two boundary conditions, I'd like to solve it numerically using matlab.

Can anybody help me out? Thanks in advance.

IC and BC: $$T_{\tau = 0} = T_{0} \hspace{2cm} 0 \leq r \leq R$$ $$\dfrac{\partial T}{\partial r}\bigg|_{r=0} = 0$$ $$\dfrac{\partial T}{\partial r}\bigg|_{r=R} = C1 \cdot (T_{m} - T_{R})$$

Yes, $\tau$ is a time variable and $r$ is the space coordinate. I'm expecting a profile of $T$ with respect to both variables.

I'm sorry if I can't provide so much useful information, but I was given a scientific paper to study where these equations (there is a lot of equations, but this one is the trickiest) are displayed but there is no mention of how they solved and got to the results.

• Can you provide a few more details, i.e. what variables you are solving for? Are you solving for $T$ over some spatial coordinate $r$? If so, what is $\tau$? And what do your boundary conditions look like? – Pedro Oct 25 '13 at 11:34
• This Question was originally posted at Math.SE, where because of its computational nature this forum was recommended. I surmise that $\tau$ is a time-like variable and $r$ is a 1D space variable. Note however the singularity of coefficients at $r=0$. – hardmath Oct 25 '13 at 12:26
• Hi @user103012, and welcome to scicomp! Are you familiar with finite difference methods for PDE's? – Paul Oct 25 '13 at 13:03
• You can separate this even further with $g(\tau)h(r)=T(\tau,r)$, since $\tau$ and $r$ are not coupled. Likely, the decoupled equation may have an analytical solution. – Bort Oct 28 '13 at 0:31

It's just an advection diffusion equation in one space dimension. I don't know how to do it in matlab, but there are tutorial program in most of the open source finite element libraries (deal.II, fenics, libmesh, ...) that demonstrate how to solve such problems. Solving this equation with a finite element code isn't going to be particularly difficult as long as either (i) the diffusion dominates the advection, or (ii) you properly stabilize the equation, or (iii) you can afford a sufficiently fine mesh, which you likely can because your model is only 1d.

(Disclaimer: I'm one of the authors of deal.II.)

• I'm pretty sure that your library would solve this and I'd gladly give it a try but I don't think I have the competence to write down this equation in a C++ code. – user103012 Oct 25 '13 at 12:51
• There is a singularity in $r=0$ in the coefficients of the equations, so it's not exactly the classical AD equation. – Dr_Sam Oct 25 '13 at 13:17
• @Dr_Sam: The BC $\partial_rT=0|_{r=0}$ should prevent that singularity, no? – Kyle Kanos Oct 26 '13 at 3:15
• The singularity disappears if you take the appropriate inner product for spherically symmetric problems, which is $(u,v)=\int_0^R u(r) v(r) r^2 \; dr$. – Wolfgang Bangerth Oct 26 '13 at 10:40

Given the way they the boundary conditions are given I would suggest a finite difference code as hinted by Paul up above. It seems that you can set up a grid for $\tau$ and $r$ and pick a forward centered scheme for $\tau$ and possibly a centered difference scheme for $r$. I can go into more detail if you know finite differencing, of if you eventually want to use this method. Just as a first step though you may want to get all of the derivative with respect to to $r$ on one side and those w.r.t to $\tau$ on the other.

I know that Matlab has built in ODE solvers but I'm not sure about PDE solvers. I think you will learn the most by making a first attempt by yourself, by writing some code in a .m file.

• Thanks @Paul. You see, I'm not familiar with PDE or FEA. I searched for some tutorial of how to solve PDEs in matlab but I've never been able to find an example that suits my equation. – user103012 Oct 26 '13 at 8:51

I'd also suggest a finite difference code (as suggested by Paul and AshikIdrisy), but rather than discretizing in both $r$ and $\tau$, I would only discretize in $r$. Semi-discretizing the differential equations in $r$ will yield a system of ODEs in $\tau$, which can then be solved in MATLAB using one of MATLAB's built-in ODE solvers.

As far as PDE solvers go, you could try using MATLAB's pdepe solver, since your equation looks parabolic.