MATLAB BioHeat Equation solution in MATLAB using pdepe

  • Thread starter Thread starter m.r.fouladi
  • Start date Start date
  • Tags Tags
    Matlab Pde
AI Thread Summary
The discussion focuses on solving the bioheat equation using MATLAB's pdepe function, specifically for a tumor located at the center of healthy tissue. The equation incorporates various constants related to thermal properties and metabolic energy release rates. The user seeks assistance in correcting their MATLAB code, which they believe does not align with the physical principles of the problem. Key assumptions include uniform initial temperatures and zero flux at the tumor's center and the outer surface of the healthy tissue. The conversation emphasizes the need for accurate coding to reflect the underlying physics of the bioheat transfer scenario.
m.r.fouladi
Messages
2
Reaction score
0
We have this Equation as bioheat equation:
∂T/∂t = α ∇2T + 1/ρc[S+Sp+Sm]
and also this:
Sp=mbcb(Tab-T)
that all α,ρ,c,S,Sm,mb,cb,Tab are constants, now I want to solve this equation in conditions below with pdepe in MATLAB:
There is a Tumor as a sphere with radius 1 cm exactly in center of a Normal Tissue with radius of 5 cm, an electrode at t=0 gives an Energy to the center of Tumor for 400 seconds that T will vary and we want to have this variation by solving this equation, also below assumptions is considered:
We assume that Tumor has the radius of 1 cm and that is perfused at the same level as healthy tissue. The tumor is located at the center of a sphere of healthy tissue that has a radius of 5 cm and all of the tissue has a metabolic energy release rate of 145 W/m3.
The power is delivered locally to the tumor by positioning the electrode in the center of the tumor. the initial temperature of tumor and surrounding tissue is uniform at 37 °C. in addition, It is assumed that because of symmetry there is zero flux at the center of the tumor and that at the outer surface of healthy tissue the flux is also zero. the following numerical values for the parameters are used:
mb = 0.18
cb = 3300
Sm = 145
Tab = 37 °C
α = 10∧(-7)
c=3800
S= 4 * 10∧(5)
ρ= 850.I coded it in MATLAB such below, but my code isn't right as my knowledge, please correct it my masters:
Matlab:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pdex1

m = 2;
x = 0.001:0.001:05;
t = 0:1:400;

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);

% A surface plot is often a good way to study a solution.
surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,400)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1e7;
f = DuDx;
if (x<0.01)
s = (18.089*10^(-2)-(18*10^(-5)*u))*(1e7);
else
s=594*(310-u)/(850*3800)
end
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 310;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul-exp((-18)*10^(-5)*t)-1004.94;
ql = 0;
pr = ul-exp((-18)*10^(-5)*t)-1004.94;
qr = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
why It's failed?
 
Last edited by a moderator:
Physics news on Phys.org
m.r.fouladi said:
why It's failed?
Why do you say it failed?

Also, I have added
Code:
tags around your code to make it more readable. Please do that every time you post code.
 
DrClaude said:
Why do you say it failed?

Also, I have added
Code:
tags around your code to make it more readable. Please do that every time you post code.
thanks for your attention.
because It doesn't agree with my physics knowledge of this problem
 

Similar threads

Replies
8
Views
3K
Replies
6
Views
4K
Replies
1
Views
4K
Replies
3
Views
4K
Replies
7
Views
4K
Replies
6
Views
2K
Back
Top