Solving bioheat PDE numerically

[eprparadox]

Hello!

My goal here is to plot the solution to the bioheat equation for a tumor as a function of time. I'm plotting this for a fixed radius at r = 0 (the very center of the tumor).

The equation to solve is this:

$$\rho_1c_1 \frac{\partial T}{\partial t} = 3(\frac{\partial }{\partial r}(\frac{\partial T}{\partial r}) + \frac{w_{b1}c_b}{\lambda_1}((T_b - T) + \frac{P\lambda_1}{w_{b1}c_b})$$

The numerical solution to this equation is given in a paper I'm reading. Using the central finite difference method, they get:

$$T_0^{n+1} = (1 - 6 \frac{\lambda_1}{\rho_1c_1}(\frac{\Delta t}{\Delta r^2}))T_0^n + 6(\frac{\lambda_1}{\rho_1c_1}\frac{\Delta t}{\Delta r^2}) + (\frac{\Delta t}{\rho_1c_1})w_{b1}c_b((T_b - T_0^n) + \frac{P\lambda_1}{w_{b1}c_b})$$Note: the lambda, rho, c constants are thermal conductivity, density and specific heat of the tumor and anything with a subscript b refers to blood.

The graph presented in the paper for the above discrete solution (for 0 < t < 600 seconds) has the Temperature linear for about the first 100 seconds and then it begins to taper off after that (see attached photo).

But when I plot this solution in mathematica, it's just linearly grows without bound. And I don't know why.

I was thinking that maybe I needed to solve the PDE for some sort of steady state solution and then plot the sum of this steady state + the transient solution above. But the paper doesn't include a steady state solution for the PDE so that got me confused.

Any insight into this would be great.

Thanks so much.

 

[Simon Bridge]

$$T_0^{n+1} = (1 - 6 \frac{\lambda_1}{\rho_1c_1}(\frac{\Delta t}{\Delta r^2}))T_0^n + 6(\frac{\lambda_1}{\rho_1c_1}\frac{\Delta t}{\Delta r^2}) + (\frac{\Delta t}{\rho_1c_1})w_{b1}c_b((T_b - T_0^n) + \frac{P\lambda_1}{w_{b1}c_b})$$

You have an equation of form:
$T_{n+1}=AT_n + B + C - DT_n + E$

$$A=1 - 6 \frac{\lambda_1}{\rho_1c_1}(\frac{\Delta t}{\Delta r^2})\\
B=6(\frac{\lambda_1}{\rho_1c_1}\frac{\Delta t}{\Delta r^2})\\
C=\frac{\Delta t}{\rho_1c_1})w_{b1}c_b T_b\\
D=\frac{\Delta t}{\rho_1c_1})w_{b1}c_b\\
E=\frac{\Delta t}{\rho_1c_1})w_{b1}c_b\frac{P\lambda_1}{w_{b1}c_b}$$

rearranging: $T_{n+1}=(A-D)T_n+(B+C+E) = \lambda T_n+\nu$

If all A through E are constants, then the curve can be characterized by the value of A-D.
i.e. if A-D=0, the curve is a horizontal line.

To get the type of curve in the thumbnail, you need $0<A-D<1$

 

[Chestermiller]

The finite difference equation you presented does not look like the correct finite difference approximation to the partial differential equation.

Chet

 

[eprparadox]

Hello,

First, thanks Simon for your input. I looked at my data and it seems like Chestermiller, you're right. For the constant data provided, I just can't seem to get the output the author gets.

But I guess my next question is how do you know it doesn't look correct? And in order to take the PDE and find the finite difference equation, do I simply substitute in the approximations for the first and second derivatives with the central finite difference equation and then simplify?

Thanks for the help.

 

[eprparadox]

Hello!

Alright so as you can see in the attached image, I'm working solving the provided PDE using the central finite difference method.

I"ve substituted in the approximations for derivatives.

But the key thing here is that the PDE in question is for r = 0 (at the very center of the tumor) and we have a condition that states:

$$T_{i-1}^n = T_1^n$$

But since we're dealing with r = 0, that essentially means i = 0.

This leaves me stuck on how to deal with the second derivative w.r.t r on the right side of the equation. If I plug in i = 0, then the last term becomes

$$T_{-1}^n$$

And I don't understand what I'm suppose to do with that.
Any ideas?

 

[the_wolfman]

You're finite difference equation is wrong. For starters your equation is missing terms that account for spatial variation in temperature.

My guess is that is should be

$T_0^{n+1} = (1 - 6 \frac{\lambda_1}{\rho_1c_1}(\frac{\Delta t}{\Delta r^2}))T_0^n + 6(\frac{\lambda_1}{\rho_1c_1}\frac{\Delta t}{\Delta r^2})T_1^n + (\frac{\Delta t}{\rho_1c_1})w_{b1}c_b((T_b - T_0^n) + \frac{P\lambda_1}{w_{b1}c_b})$

note that I've inserted a term $T_1^n$. This is just an educated guess and you should check it.


First start by constructing a general difference equation that applies to points in the interior.
You can then use this equation to construct a new equation for T_0 by assuming that a ghost point exists at T_{-1}. You are assuming that T'=0 at r=0. Using 2nd order central difference this gives you the equation for T_{-1}...T_{-1}=T_{1}. Use this in you original the equation for T_0 and see what you get.

If you want to find the steady state solution let $\partial T /\partial t=0$ in the initial PDE. This will give you a solvable second order ode for T that depends only on r.

 

[Chestermiller]

Not only does the finite difference equation look wrong, but so also does the PDE. Where is the thermal conductivity on the first term on the right hand side of the equation? Also, where does the factor of 3 come from?

You alluded to the parameter r being the "radial" location. If this is in cylindrical coordinates or in spherical coordinates, then the first term on the right hand side of the equation is incorrect (geometrically). Before we start getting into how to handle the boundary condition at r = 0, we need to ascertain which of these (if either) is what you intend. And, before we do that, we need to give the PDE a clean bill of health. Please show us your derivation of the PDE. Also please tell us about the boundary condition at the outer radius.

Getting the PDE straightened out shouldn't be too much trouble. But it might also be worthwhile for you to take a crash course in numerical analysis. Even when the spatial finite difference approximations are substituted into the equation, there will still be numerical solution issues to be considered with respect to the use of implicit or explicit solutions in time, and the whole issue of the numerical stiffness of the associated set of coupled ODEs.

Chet

 

[eprparadox]

Hey Chet,

Thanks so much for your response. So I thought I'd back up. All of this data is taken from a thesis paper. I figured it would be correct, but let me start from the beginning. We have the spherical pennes bioheat equation for inside of a tumor of radius R. So looking at 0 < r < R:

$$\rho_1 c_1 \frac{\partial T_1}{\partial t} = \frac{\lambda_1}{r^2}\frac{\partial}{\partial r}(r^2\frac{\partial T_1}{\partial r}) + w_{b1}\rho_b c_b(T_b - T_1) + P$$

They go on to state that: a central differencing scheme is used to discretize this equation at every node. The boundary condition of zero temperature gradient at r = 0 is imposed by using a fictitious node and prescribing

$$T_{i-1}^n = T_1^n$$

in the discretized governing equation at r = 0. In order to avoid a singularity at r = 0, this equation is replaced at r = 0 with

$$\rho_1c_1 \frac{\partial T}{\partial t} = 3(\frac{\partial }{\partial r}(\frac{\partial T}{\partial r}) + \frac{w_{b1}c_b}{\lambda_1}((T_b - T) + \frac{P\lambda_1}{w_{b1}c_b})$$

And I apologize for all of these questions (I do need to ultimately go back and take a crash course on PDE's and numerical analysis), but how do I go about dealing with this singularity and simplifying this PDE at r = 0?

Do I need to study more about PDE's before attempting this? I think just starting from the top is probably the best plan of attack here.

Thanks so much.

 

[bigfooted]

For simplicity, let's assume that all timesteps are the same and that the mesh is equidistant. The simplest discretization is like this:

$\frac{\partial T}{\partial t} \approx \frac{\Delta T}{\Delta t} = \frac{T_i^{n+1} -T_i^n}{\Delta t}$
$\frac{\partial T}{\partial r} \approx \frac{\Delta T}{\Delta x} = \frac{T_{i+1}^n -T_{i}^n}{\Delta r}$
and for second derivatives:
$\frac{\partial^2 T}{\partial r^2} \approx \frac{T_{i+1}^n -T_{i}^n}{\Delta r} - \frac{T_{i}^n -T_{i-1}^n}{\Delta r} = \frac{T_{i+1}^n -2T_i^n + T_{i-1}^n}{(\Delta r)^2}$

Let's say your equation is $\frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial r^2}$
I will treat the singularity later.
The discretized version is
$\frac{T_i^{n+1} -T_i^n}{\Delta t} = \frac{T_{i+1}^n -2T_i^n + T_{i-1}^n}{(\Delta r)^2}$
rearranged we have:

$T_i^{n+1} = T_i^n + \left(T_{i+1}^n -2T_i^n + T_{i-1}^n \right)\frac{\Delta t}{(\Delta r)^2}$

The domain is from r=0..R, divided into N equidistant cells, and we number the nodes as i=0...N+1
For simplicity I will remove the superscript indicating the current timestep 'n'. The new timestep we will simply call timestep 'new':

$T_i^{new} = T_i + \left(T_{i+1} -2T_i + T_{i-1} \right)\frac{\Delta t}{(\Delta r)^2}$

At r=0, we have:
$T_0^{new} = T_0 + \left(T_{1} -2T_0 + T_{-1} \right)\frac{\Delta t}{(\Delta r)^2}$

The symmetry condition says that $\frac{\partial T}{\partial r}|_{r=0}=0$
A second order accurate discretization would be:
$\frac{T_1-T_{-1}}{2\Delta r}=0$
rearrange:
$T_{-1}=T_1$
So: conveniently, the value at the 'extra' node is equal to the value at node i=1
So at r=0 we can now write:
$T_0^{new} = T_0 + \left(T_{1} -2T_0 + T_{1} \right)\frac{\Delta t}{(\Delta r)^2}$
This problem can now be solved.
Now for the singularity.
We now consider the actual equation in spherical coordinates $\frac{\partial T}{\partial t} = \frac{1}{r^2}\frac{\partial}{\partial r} (r^2\frac{ \partial T}{\partial r})$
The usual (?) approach is to write down the discretized equation for the center node of the pde in Cartesian coordinates : $\frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}$
$\frac{T_{0,0}^{new} -T_{0.0}}{\Delta t} =\frac{T_{1,0} -2T_{0,0} + T_{-1,0}}{(\Delta x)^2} +\frac{T_{0,1} -2T_{0,0} + T_{0,-1}}{(\Delta y)^2}$
Because of symmetry, the nodes around (0,0) are all the same:
$T_{1,0}= T_{-1,0} = T_{0,-1} = T_{0,1}$, so let's replace them all with $T_{1,0}$:
$\frac{T_{0,0}^{new} -T_{0.0}}{\Delta t} =\frac{T_{1,0} -2T_{0,0} + T_{1,0}}{(\Delta x)^2} +\frac{T_{1,0} -2T_{0,0} + T_{1,0}}{(\Delta y)^2}$
$\frac{T_{0,0}^{new} -T_{0.0}}{\Delta t} =\frac{2T_{1,0} -2T_{0,0}}{(\Delta r)^2} +\frac{2T_{1,0} -2T_{0,0}}{(\Delta r)^2}$
$\frac{T_{0,0}^{new} -T_{0.0}}{\Delta t} =\frac{4T_{1,0} -4T_{0,0}}{(\Delta r)^2}$
and go back to radial coordinates to get the boundary condition at r=0:
$\frac{T_{0}^{new} -T_{0}}{\Delta t} =\frac{4T_{1} -4T_{0}}{(\Delta r)^2}$

So when you have the initial solution at time t=0, the problem can be solved.

Do I need to study more about PDE's before attempting this? I think just starting from the top is probably the best plan of attack here.

You should definitely read something about PDEs, heat transfer (heat or conduction equation), and numerical methods (finite difference methods).

 

[Chestermiller]

Hey Chet,

Thanks so much for your response. So I thought I'd back up. All of this data is taken from a thesis paper. I figured it would be correct, but let me start from the beginning. We have the spherical pennes bioheat equation for inside of a tumor of radius R. So looking at 0 < r < R:

$$\rho_1 c_1 \frac{\partial T_1}{\partial t} = \frac{\lambda_1}{r^2}\frac{\partial}{\partial r}(r^2\frac{\partial T_1}{\partial r}) + w_{b1}\rho_b c_b(T_b - T_1) + P$$

They go on to state that: a central differencing scheme is used to discretize this equation at every node. The boundary condition of zero temperature gradient at r = 0 is imposed by using a fictitious node and prescribing

$$T_{i-1}^n = T_1^n$$

in the discretized governing equation at r = 0. In order to avoid a singularity at r = 0, this equation is replaced at r = 0 with

$$\rho_1c_1 \frac{\partial T}{\partial t} = 3(\frac{\partial }{\partial r}(\frac{\partial T}{\partial r}) + \frac{w_{b1}c_b}{\lambda_1}((T_b - T) + \frac{P\lambda_1}{w_{b1}c_b})$$

And I apologize for all of these questions (I do need to ultimately go back and take a crash course on PDE's and numerical analysis), but how do I go about dealing with this singularity and simplifying this PDE at r = 0?

Do I need to study more about PDE's before attempting this? I think just starting from the top is probably the best plan of attack here.

Thanks so much.

OK. So you're working in spherical coordinates.

There are two ways of dealing with the boundary condition at r=0 to get the proper second order spatial finite difference approximation for the first term on the right hand side of the equation. Both methods give exactly the same answer. The following is the simpler method.

Since T is maximum or minimum at r = 0, we can expand T in a power series about r = 0 (at any time t):

$$T=T(0) + Br^2 +...$$
So:
$$\frac{∂T}{∂r}=2Br$$
$$r^2\frac{∂T}{∂r}=2Br^3$$
$$\frac{∂}{∂r}\left(r^2\frac{∂T}{∂r}\right)=6Br^2$$
$$\frac{1}{r^2}\frac{∂}{∂r}\left(r^2\frac{∂T}{∂r}\right)=6B$$
But, from the first equation,
$$B≈\frac{T_1-T_0}{(Δr)^2}$$
Therefore,
$$\frac{1}{r^2}\frac{∂}{∂r}\left(r^2\frac{∂T}{∂r}\right)≈6\left(\frac{T_1-T_0}{(Δr)^2}\right)$$
This is second order accurate, the same degree of accuracy as the central difference approximations at the other grid points. This is where the 6 in the_wolfman's response came from. That modified PDE with the factor of 3 is incorrect, so please get rid of it.

For the finite difference approximation to the spatial derivative term at the other grid points, I strongly recommend the following:

$$\frac{1}{r^2}\frac{∂}{∂r}\left(r^2\frac{∂T}{∂r}\right)=\frac{1}{r_i^2(Δr)^2} \left((r_i+\frac{Δr}{2})^2(T_{i+1}-T_i) -(r_i-\frac{Δr}{2})^2 (T_{i}-T_{i-1})\right)$$

Chet