diffusion from point source boundary conditions

Question

Juliana 2022-7-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1764705-diffusion-from-point-source-boundary-conditions

编辑： Juliana 2022-7-25

I am using MATLAB's pdepe partial diffeq solver to solve diffusion of a solute. I have a point source at r=0 where there should be intake of solute (so a sink, or negative source). if i use the "s" variable for the solver, does this input it at r=0, or for every value of r? if the latter, then should i incorporate the source term in my BC somehow?

any help is appreciated. thanks.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-7-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1764705-diffusion-from-point-source-boundary-conditions#answer_1012075

编辑：Torsten 2022-7-20

在 MATLAB Online 中打开

if i use the "s" variable for the solver, does this input it at r=0, or for every value of r?

For every value of r (at least if you don't restrict s to special r-values by an if-statement).

if the latter, then should i incorporate the source term in my BC somehow?

The boundary condition at r=0 must be dc/dr = 0 - all other boundary conditions will be ignored.

In my opinion, you have two options:

Define a sphere with a hole in the middle of a certain radius and set the sink term in the boundary conditions part of the code as a flux over the inner surface (mol/(m^2*s)) or

Define a volume source in the s-term (mol/m^3*s) and restrict it to a specified inner radius:

rmin = ...;
if r <= rmin
  s = ...
end 

In either case: These are difficult error-prone and numerically challenging settings and you should verify your results by making a global molar or mass balance.

12 个评论
显示 10更早的评论隐藏 10更早的评论

Torsten 2022-7-21

编辑：Torsten 2022-7-21

在 MATLAB Online 中打开

If you change rmin to 1/37, nothing happens in the code below.

The function "heatsphere" is not called with the boundary values for r (thus r=0 and r=1).

The minimum r value with which it is called is 1/36. I don't know how MATLAB comes up with this value. It's neither part of the grid vector r nor a cell center. Maybe because m=2 instead of m=0, MATLAB includes a certain propagation factor in the r-values.

Since you also don't know the volume of the cells the source term is connected with, it is nearly impossible to dose the source according to what you want to set.

I suggest to include the source as a boundary source at a value r=rmin > 0.

Then you have the area A = 4*pi*rmin^2 over which the substance diffuses and you can set the flux accordingly as D*dn/dr = ndot/A.

r = linspace(0,1,25);

t = linspace(0,1000,25);

m = 2;

sol = pdepe(m,@heatsphere,@heatic,@heatbc,r,t);

u = sol(:,:,1);

surf(r,t,u)

xlabel('r')

ylabel('t')

zlabel('u(r,t)')

view([150 25])

plot(t,sol(:,1))

xlabel('Time')

ylabel('Temperature u(0,t)')

title('Temperature change at center of sphere')

function [c,f,s] = heatsphere(r,t,u,dudr)

D = 1e-4;

c = 1;

f = D*dudr;

s = 0;

rmin = 1/36;

if r <= rmin

s = -12;

end

%----------------------------------------------

function u0 = heatic(r)

u0 = 300;

end

%----------------------------------------------

function [pl,ql,pr,qr] = heatbc(rl,ul,rr,ur,t)

pl = 0.0;

ql = 1.0;

pr = ur - 300;

qr = 0;

end

%----------------------------------------------

Juliana 2022-7-22

编辑：Juliana 2022-7-22

在 MATLAB Online 中打开

Thanks again, the value of 1/36 worked well.

i am attempting to update the source value for every few time values, since the pdepe needs at least 3 values of time to eval for

i am having trouble with the syntax of it. not sure where the loop should go, because right now i get an error that the variable s1 is not recognized, but if i were to put the loop within the function, then i cant do multiple time step vectors. my code is pasted below.

%prompt = "What is concentration of nutrient in medium? ";    %user inputs
medium = 100; %prompt;
%prompt = "What is diffusivity of nutrient in medium? ";
D = .0001; %prompt;
%prompt = "What is the maximum radial distance for the system? ";
rinput = .5; %prompt;
%prompt = "What is maximum time span for the system? ";
tinput = 100; %prompt;
r = linspace(0,rinput,1000);        %space vector
u=linspace(100,0,100);
tn=1;
tm=4;
     
 for tm=4:tinput
    t0=linspace(tn,tm,4);   %make vector with small t values to do steps of diffusion
    rmin=1/36;              %smallest value that matlab computes
    if r<=rmin              %specifies source at approxr=0
    s1=-2*u(tm,rmin);     %rbitraty function for source that depends on u
    else                    %no sink anywhere else
    s1=0;
    end
Unew=diffusio(medium,D,r,t0);      %new diffusion calculation each time
% u= find some way to store new values of matrix with previous iterations
     tn=tn+4;               %add 4 to time vector so we get next terms
     tm=tm+4;
    end     
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (1)                          %surfplot
surf(r,t,u,'EdgeColor', 'none')     %so grid doesnt make image all black
colorbar
clim ([0 100])
xlabel('r')
ylabel('t')
zlabel('u(r,t)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (2)
colormap cool           %aesthetic choice
imagesc(r,t,u)
colorbar 
clim ([0 100])           %fixed colorbar range
xlabel('Radial Distance r','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Diffusion','interpreter','latex')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (3)
plot(t,u(:,1))      %not sure if this is showing at surface of cell?
xlabel('Time')
ylabel('Concentration')
title('Concentration change at center of sphere')
%%%%%%%%%%%%%%%%%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%
function [u]=diffusio(medium,D,r,t)
%medium=100;               %for if we dont wanna input each time
%D=.0001;
%rinput=.5;
%tinput=500;
m = 2;                          %indicates spherical sym
sol = pdepe(m,@diffpde,@icfun,@bcfun,r,t);
u = sol(:,:,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c,f,s] = diffpde(r,t,u,dudr)     %pde definition
c = 1;
f = D*dudr;                 %diffusivity times gradient
s=s1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u0 = icfun(r)     %initial condition
  u0 = medium;              %at t=0, all r, conc=medium
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pl,ql,pr,qr] = bcfun(rl,ul,rr,ur,t)   %boundary conditions
pl = 0;                 %ignored bc spherical   
ql = 0;                 %ignored bc spherical   
pr = 0;                 %indicates no flux on boundary
qr = 1;                 %indicates no flux on boundary
end
end

Torsten 2022-7-22

编辑：Torsten 2022-7-22

在 MATLAB Online 中打开

I don't understand why you make things that complicated.

You can directly set the sink term in the function itself depending on u.

You also have access to the time t in the function if you want to set something time-dependent.

medium = 100; %for if we dont wanna input each time

D = .0001;

rend = .5;

nr = 1000;

r = linspace(0,rend,nr);

rmin = 1/36;

tend = 500;

nt = 100;

t = linspace(0,tend,100);

u = diffusio(medium,D,rmin,r,t);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure (1) %surfplot

surf(r,t,u,'EdgeColor', 'none') %so grid doesnt make image all black

colorbar

clim ([0 100])

xlabel('r')

ylabel('t')

zlabel('u(r,t)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure (2)

colormap cool %aesthetic choice

imagesc(r,t,u)

colorbar

clim ([0 100]) %fixed colorbar range

xlabel('Radial Distance r','interpreter','latex')

ylabel('Time t','interpreter','latex')

title('Diffusion','interpreter','latex')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure (3)

plot(t,u(:,1)) %not sure if this is showing at surface of cell?

xlabel('Time')

ylabel('Concentration')

title('Concentration change at center of sphere')

function u = diffusio(medium,D,rmin,r,t);

m = 2; %indicates spherical sym

sol = pdepe(m,@(r,t,u,dudr)diffpde(r,t,u,dudr,D,rmin),@(r)icfun(r,medium),@bcfun,r,t);

u = sol(:,:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

function [c,f,s] = diffpde(r,t,u,dudr,D,rmin) %pde definition

c = 1;

f = D*dudr; %diffusivity times gradient

s = 0.0;

if r <=rmin

s = -2e-1*u(1);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u0 = icfun(r,medium) %initial condition

u0 = medium; %at t=0, all r, conc=medium

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pl,ql,pr,qr] = bcfun(rl,ul,rr,ur,t) %boundary conditions

pl = 0; %ignored bc spherical

ql = 1; %ignored bc spherical

pr = 0; %indicates no flux on boundary

qr = 1; %indicates no flux on boundary

end

Juliana 2022-7-25

let me explain a bit better:

i am monitoring diffusion of, say, glucose into a cell.

the cell's rate of intake (source term of pde) depends on the surrounding concentration (so, concentration at r=rmin)

therefore, for each new value of concentration at each time t, the source term will change dependent on the new concentration.

i am trying to find a way to solve the pde multiple times over and over, updating the source term for each new concentration value.

so say i want to solve the pde from t=linspace(0,10,1000) so delta t = 0.01 s

so for t=0, i will get concentration=c1 at rmin

so for t=.01, i want to take c1 and input it into my source term function, lets say s=-20*(concentration at r=rmin). so now s=-20*c1. i will get concentration=c2

then for t=0.02, i want the source to be s=-20*c2

and so on.

but i cannot do this in a number of ways for a number of reasons. I am first trying to solve the pde at each time t, but the pde needs a time span of at least three values.

so now i will try and solve the pde in the script with a for loop, creating a time span of 3 values per my delta t of 0.01. and run the pde solver over and over.

but this does not work because the syntax is wrong. i cannot update the source term outside the function that defines it. and i cannot update the source term within the dunction because it does not recognize the output matrix, obviously.

how can I solve this?

Torsten 2022-7-25

编辑：Torsten 2022-7-25

Setting s = -20*u(1) means that a 1st order reaction takes places that consumes u at a rate -20*u.

I don't understand what you mean by "u(last value)". The "last value" is always the "actual value", and this actual value is taken when you set s = -20*u(1). If the source term depends on a concentration in the past (u(t-tau) for some tau > 0), you have a delay PDE. Such a construct cannot be solved with pdepe.

What doesn't look right to you ? The concentration profile ?

Juliana 2022-7-25

编辑：Juliana 2022-7-25

i see, youre right its the same as u(1), but yeah i would expect the concentration profile to look much different, but i guess im wrong idk

请先登录，再进行评论。

diffusion from point source boundary conditions

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

diffusion from point source boundary conditions

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论