Prev: replace object properties
Next: approximating the derivative of an array
From: Dave on 28 Jun 2010 15:17
I'm trying to get a numerical solution for the mass diffusion from a sphere into a fluid of finite volume, analagous to thermal diffusion in the cooling of a sphere in a well-stirred solution. The problem I am having is that I cannot apply the right boundary condition in pdepe to represent flux from the sphere into the surrounding fluid, which is initially free of diffusing material. I have a dC/dt term which would need updating at each time step. The equations I am using are:
dc/dt = (2/x)* dc/dx + d2c/dx2 - fick's 2nd law in spherical, dimensionless form
i.c. c = 0.5 @t=0 - concentration of sphere at time = 0
b.c. dc/dy = 0 @y=0 - b.c. at the centre of the sphere
dc/dy = alpha*dc/dt @y=1 - b.c. at the sphere surface,
where alpha=(D,surf area,Vol of solution,dt,dr - obtained by mass balance of
flux through the sphere into the surrounding fluid at the sphere surface)
As dc/dt would have to be obtained semi-explicitly it can be approximated by:
c(n+1)-c(n)/dt - where c(n/n+1) is the concentration of the surrounding fluid, c, at time n/n+1 respectively.
Also, dc/dy @y=1 can be approx by (cr - c(n))/dr - where cr = c@ y=1, ie the concentration of the sphere at the surface.

The terms I have used in pdepe are as follows:
function [c,f,s] = pdex2pde(x,t,u,DuDx)
c = 1;
f = DuDx;
s = 0;
____________________________________
function u0 = pdex2ic(x)
u0 = 0.5;
____________________________________
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = 0;
ql = 1;
for h = 1:length(n)
pr = (V/(D*S*dt))*(csol(h+1) -(csol(h)));
qr = 1;
end;

I dont know how to represent to pr/qr terms, please help.
From: Torsten Hennig on 28 Jun 2010 23:00
> I'm trying to get a numerical solution for the mass
> diffusion from a sphere into a fluid of finite
> volume, analagous to thermal diffusion in the cooling
> of a sphere in a well-stirred solution. The problem I
> am having is that I cannot apply the right boundary
> condition in pdepe to represent flux from the sphere
> into the surrounding fluid, which is initially free
> of diffusing material.

This flux is usually represented by a boundary condition
of the kind
-D*dc/dx = beta*(c-c_inf)
(D: diffusion coefficient,
beta: mass transfer coefficient,
c_inf : concentration in the surrounding fluid)

Why can't you apply this condition ?

>I have a dC/dt term which
> would need updating at each time step. The equations
> I am using are:
> dc/dt = (2/x)* dc/dx + d2c/dx2 - fick's 2nd
> law in spherical, dimensionless form
> i.c. c = 0.5 @t=0 - concentration of sphere
> at time = 0
> b.c. dc/dy = 0 @y=0 - b.c. at the centre of
> the sphere
> dc/dy = alpha*dc/dt @y=1 - b.c. at the
> . at the sphere surface,
> where alpha=(D,surf area,Vol of
> a,Vol of solution,dt,dr - obtained by mass balance
> of
> flux through the sphere into the surrounding
> rounding fluid at the sphere surface)
> As dc/dt would have to be obtained semi-explicitly it
> can be approximated by:
> c(n+1)-c(n)/dt - where c(n/n+1) is the concentration
> of the surrounding fluid, c, at time n/n+1
> respectively.
> Also, dc/dy @y=1 can be approx by (cr - c(n))/dr -
> where cr = c@ y=1, ie the concentration of the
> e sphere at the surface.
>
> The terms I have used in pdepe are as follows:
> function [c,f,s] = pdex2pde(x,t,u,DuDx)
> c = 1;
> f = DuDx;
> s = 0;
> ____________________________________
> function u0 = pdex2ic(x)
> u0 = 0.5;
> ____________________________________
> function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
> pl = 0;
> ql = 1;
> for h = 1:length(n)
> pr = (V/(D*S*dt))*(csol(h+1) -(csol(h)));
> qr = 1;
> end;
>
> I dont know how to represent to pr/qr terms, please
> help.

A boundary condition of the kind
dc/dt = alpha*dc/dx at x=1 can not be set in pdepe.

Alternatively you could discretize the diffusion
equation in space and use ODE15s to solve the
resulting system of ordinary differential equations
(method of lines - approach).

Best wishes
Torsten.
From: Dave Esc on 30 Jun 2010 14:09
Thanks Torsten

I dont know how to apply the c_inf term as this concentration outside of the sphere is not constant and changes with time.

I have tried your suggested way of discretising the spherically-symmetric diffusion equation in space:
From:
dc/dt = (d^2c/dx^2)+((2/y)*(dc/dx))
to
dc(j)/dt = ((c(j+1)-2*c(j)+c(j-1))/dx^2)+((2/x(j)*(c(j+1)-c(j))/dx))
with the same Neumann b.c.s:
dc/dx=0 at x=0, discretised: c(2)=c(1), and dc/dx=alpha*dc/dt for sphere sphere surface. I have tried to use the ODE solver to solve the discretised equation but am again having trouble with the boundary conditions and having to update c_inf term (given as csol).
The error I get with ODE45 is:
??? Attempted to access c(1,2); index out of bounds because size(c)=[0,0].

Error in ==> fd_rhs at 6
c(k,1) = c(k,2);

and the script is:

global nt dt c delx delx2 csol V dr D S
% set up the parameters
time = 10000;
D = 3.95E-10;
R0 = 1.50E-03;
tau = D*time/(R0^2);
S = 1.26E+01;
V = 296.2962963;
dr = 0.01010101;
D = 3.95E-10;
N = 100;
dt = tau/(N-1);
nt = 50;
F = D*S*dt/(V*dr);
delx = 1/(nt-1);
delx2 = delx*delx;
x(1) = 0.;
for i=2:nt-1;
c0(i-1) = 0.;
x(i) = x(i-1)+delx;
end;
x(nt) = x(nt-1)+delx;

tbreak = linspace(0,tau,N);
for k=1:100;
tspan = [tbreak(k) tbreak(k+1)];
[t cc] = ode45('fd_rhs',tspan,c0);
nn=size(t);
c(k,:) = cc(nn(1),:);
c0(:) = c(k,:);
end;
csol(k,1)=zeros;
% add in the first and last point and i.c.
for k=1:100;
d(k,1) = d(k,2);
csol(k) = (csol(k-1))*(1+F*((c(k-1,nt))-1));
for i=2:nt-1;
d(k,i)=c(k,i-1);
d(1,:)=0.0005;
end;
end;
===============================================
function ydot=fd_rhs(time,y)
global nt dt c delx delx2 csol V dr D S
f = (V*dr)/(D*S*dt);
% add in the boundary conditions
for k=1:100;
c(k,1) = c(k,2);
c(nt) = (c(nt-1))+(f*(csol(k)-csol(k-1)));
end;
% transfer the y to c
for i=1:nt-2
c(i+1) = y(i);
end;

% calculate the right-hand side
for i=1:nt-2;
dffd(i) = c(i) - 2*c(i+1) + c(i+2);
dfd(i) = 2*(c(i+1)-c(i))/i;
end;
dffd = dfd/delx2;
dfd = dfd/delx;
sol = dfd + dffd;
% put result in the ydot
ydot = sol';

Could you tell me where I am going wrong and how you suggest I resolve the problem?
Many thanks
From: Dave Esc on 12 Jul 2010 13:38
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1443157399.21287.1277794843549.JavaMail.root(a)gallium.mathforum.org>...
> > I'm trying to get a numerical solution for the mass
> > diffusion from a sphere into a fluid of finite
> > volume, analagous to thermal diffusion in the cooling
> > of a sphere in a well-stirred solution. The problem I
> > am having is that I cannot apply the right boundary
> > condition in pdepe to represent flux from the sphere
> > into the surrounding fluid, which is initially free
> > of diffusing material.
>
> This flux is usually represented by a boundary condition
> of the kind
> -D*dc/dx = beta*(c-c_inf)
> (D: diffusion coefficient,
> beta: mass transfer coefficient,
> c_inf : concentration in the surrounding fluid)
>
> Why can't you apply this condition ?
>

Hi Torsten

I'm trying to apply the above condition using pdepe solver but c_inf has to be updated at each time step. I have tried to form a loop where at first c_inf = 0 and then the pde solver is used over a time step dt subject to the boundary condition you prescribe, then c_inf is updated subject to dc_inf/dt = J*S/V where J is flux across sphere surface area S into surrounding fluid volume V. This value of c_inf is then to be used in the next step of the pdepe and continues for the required time span. The trouble I am having is using the latest solution values of the pde (ie c distribution within the sphere). Using this method I guess I would effectively have to update the inital conditions each time step?
Here is the code that I have attempted to use but will only give solutions for the first time step:

function [c,f,s] = pde1(x,t,u,DuDx)
c = 1;
f = DuDx;
s = 0;
%---------------------------------------------------
function u0 = pde1ic(x)
u0 = 0.0005;
%---------------------------------------------------
function [pl,ql,pr,qr] = pde1bc(xl,ul,xr,ur,t)
global c_inf
betaD = 0.02;
pl = 0;
ql = 1;
pr = betaD*(ur - c_inf);
qr = 1;
%---------------------------------------------------
global c_inf
R0 = 0.0015;
S = 2.83E-05;
V = 0.000001;
A = R0*S/V;
dx = 0.02;
tau = 1.7556;
dt = tau/99;
tstep = (tau/dt)+1;
c_inf = zeros(1,tstep);
for k = 1:tstep;
m = 2;
x = 0:.02:1;
t = [0 dt/2 dt];
sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
u(:,:,k) = sol(:,:,end);
c_inf(1,:) = c_inf(1,k)-A*dt*((u(3,51,k)-u(3,50,k))/dx);
end;
%---------------------------------------------------

The problem seems to be getting the pde to use updated u, for each step of the loop the pde reverts back to the ic. u0=0.0005 when I want it to use the latest values of u obtained.
Any idea how to do this or how to resolve the problem?
Regards
From: Torsten Hennig on 12 Jul 2010 22:53
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <1443157399.21287.1277794843549.JavaMail.root(a)gallium.
> mathforum.org>...
> > > I'm trying to get a numerical solution for the
> mass
> > > diffusion from a sphere into a fluid of finite
> > > volume, analagous to thermal diffusion in the
> cooling
> > > of a sphere in a well-stirred solution. The
> problem I
> > > am having is that I cannot apply the right
> boundary
> > > condition in pdepe to represent flux from the
> sphere
> > > into the surrounding fluid, which is initially
> free
> > > of diffusing material.
> >
> > This flux is usually represented by a boundary
> condition
> > of the kind
> > -D*dc/dx = beta*(c-c_inf)
> > (D: diffusion coefficient,
> > beta: mass transfer coefficient,
> > c_inf : concentration in the surrounding fluid)
> >
> > Why can't you apply this condition ?
> >
>
> Hi Torsten
>
> I'm trying to apply the above condition using pdepe
> solver but c_inf has to be updated at each time step.
> I have tried to form a loop where at first c_inf = 0
> and then the pde solver is used over a time step dt
> subject to the boundary condition you prescribe, then
> c_inf is updated subject to dc_inf/dt = J*S/V where J
> is flux across sphere surface area S into surrounding
> fluid volume V. This value of c_inf is then to be
> used in the next step of the pdepe and continues for
> the required time span. The trouble I am having is
> using the latest solution values of the pde (ie c
> distribution within the sphere). Using this method I
> guess I would effectively have to update the inital
> conditions each time step?
> Here is the code that I have attempted to use but
> will only give solutions for the first time step:
>
> function [c,f,s] = pde1(x,t,u,DuDx)
> c = 1;
> f = DuDx;
> s = 0;
> %---------------------------------------------------
> function u0 = pde1ic(x)
> u0 = 0.0005;
> %---------------------------------------------------
> function [pl,ql,pr,qr] = pde1bc(xl,ul,xr,ur,t)
> global c_inf
> betaD = 0.02;
> pl = 0;
> ql = 1;
> pr = betaD*(ur - c_inf);
> qr = 1;
> %---------------------------------------------------
> global c_inf
> R0 = 0.0015;
> S = 2.83E-05;
> V = 0.000001;
> A = R0*S/V;
> dx = 0.02;
> tau = 1.7556;
> dt = tau/99;
> tstep = (tau/dt)+1;
> c_inf = zeros(1,tstep);
> for k = 1:tstep;
> m = 2;
> x = 0:.02:1;
> t = [0 dt/2 dt];
> sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);

Use here
u(:,:,k) = sol(:,:,1);
c_inf = c_inf - A*dt*(u(3,51,k)-u(3,50,k))/dx;
end;

c_inf should be a scalar variable for use in pde1bc.

> u(:,:,k) = sol(:,:,end);
> c_inf(1,:) =
> :) = c_inf(1,k)-A*dt*((u(3,51,k)-u(3,50,k))/dx);
> end;
> %---------------------------------------------------
>
> The problem seems to be getting the pde to use
> updated u, for each step of the loop the pde reverts
> back to the ic. u0=0.0005 when I want it to use the
> latest values of u obtained.
> Any idea how to do this or how to resolve the
> problem?
> Regards

Make _u_ and _k_ global variables, initialize u = 0.0005 before the first call to pdepe and use u(3,:,k) in pde1ic at time step k to set up u0.

Best wishes
Torsten.
 |  Next  |  Last
Pages: 1 2 3
Prev: replace object properties
Next: approximating the derivative of an array