MATLAB: Pdepe the problem of negative concentrations

MATLABnegative concentrationspdepe

pdepe the problem of negative concentrations.

May be a bit stupid question,

Assume that we have the simplest form of the diffusion equation

dC(x,t)/dt = D * d2C(x,t)/dx2

with initial conditions C(x,0) = C0

and Neumann boundary conditions (representing for example, consumption of the ions from the electrode): – left side (x=0), no-flux: dC(0,t)/dx = 0 ==> pl = 0; ql = 1 – right side (x=L), constant current flux: dC(L,t)/dx = I/(F*A) ==> pr = I/(F*A); qr = 1 where I is applied constant current; F is Faraday number; A is the electrode surface. So, I/(F*A) = const.

When I run the simulation at some moment I've got the negative values of the concentrations C(x,t) at the whole distance of the electrode. This solution is definitely non-physical.

How can I solve the problem of negative concentrations within pdepe? Is it possible? (may be there are some methods like changing of the boundary conditions when C(x,t) reaches 0 at the right boundary, assign the "negative concentrations" 0 (like C(x,t)=0) after the simulation, etc.)

Thank you very much for your answer and advice beforehand,

Sincerely,

Alexander Ledovskikh

Best Answer

In this case, pr and qr must be defined as

qr = 1 
pr = -D_Li*I/(F*A)

Best wishes

Torsten.

Related Solutions

MATLAB: Error using pdepe (line 293) Spatial discretization has failed.

Try

    function Kronstadt_HW4_problem3_pdepe
     %This is a dimensional formulation 
      close all
    Nx=201;    %number of points in x direction
    Nt=101;    %number of points in time
    m=0;       %cartesian
    D=2.4E-9;  %m^2/s, diffusivity
    Na=0.025;  %mol/m^3-s, sink
    Ci=0.04;   %mol/m^3, initial condition
    Cw = 0.04;  %mol/m^3, this is the set concentration at the right boundary
    Xmax=.1;   %m
    tmax=1000000;  %sec
    dx=Xmax/(Nx-1);
    dt=tmax/(Nt-1);
    t=linspace(0,tmax,Nt);
    x=linspace(0,Xmax,Nx);
    options=odeset; %('MaxStep',100);
     %anonymous functions to pass parameters
      pbBC=@(yl,Cl,yr,Cr,t)pbBC1(yl,Cl,yr,Cr,t,Cw);
      pbIC = @(x)pbIC1(x,Xmax,Cw);  
      pdepb=@(x,t,C,DCDx)pdepb1(x,t,C,DCDx,Na,D);
      sol=pdepe(m,pdepb,pbIC,pbBC,x,t,options);  %solution      
      C=sol(:,:,1);                               %extract solution 
    x_value1=6;
    x_value2=16;   
    x_value3=36; 
    x_value4=51;
    plot(t,C(:,x_value1),t,C(:,x_value2),t,C(:,x_value3),t,C(:,x_value4))   %Plot vs t at specified x locations
    xlabel('Time (s)')
    ylabel('Concentration (mol/m^3)')
    legend([num2str(x_value1*dx-dx) ' m'],[num2str(x_value2*dx-dx) ' m'],[num2str(x_value3*dx-dx) ' m'],[num2str(x_value4*dx-dx) ' m'])
    figure(2)   %Plot vs x at specified times
    t1=2;
    t2=11;
    t3=31;
    t4=101;
    plot(x,C(t1,:),x,C(t2,:),x,C(t3,:),x,C(t4,:))
    xlabel('x (m)')
    ylabel('Concentration (mol/m^3)')
    ylim([0 Cw])
    xlim([0 Xmax])
    legend([num2str(t1*dt-dt) ' sec'],[num2str(t2*dt-dt) ' sec'],[num2str(t3*dt-dt) ' sec'], [num2str(t4*dt-dt) ' sec'])
    figure(3)   %Flux at both wall surfaces
    flux_l=-D*(C(:,2)-C(:,1))/dx;    %mol/s-m

    flux_r=-D*(C(:,Nx)-C(:,Nx-1))/dx;    %mol/s-m
    plot(t,flux_l,t,flux_r)
    legend('Flux (left wall)','Flux (right wall)')
    xlabel('Time (sec)')
    ylabel('Mass Flux (mol/s-m^2)')
    ylim([0 .5E-8])
    figure(4)   %Flux at both left surface
    semilogy(t,flux_l)
    legend('Flux (left wall)')
    xlabel('Time (sec)')
    ylabel('Mass Flux (mol/s-m^2)')
    grid on
    figure(5)
    Ns=5;
    ta=linspace(0,tmax,(Nt-1)/Ns);
    xa=linspace(0,Xmax,(Nx-1)/Ns);
    [X,T]=meshgrid(xa,ta);
    for i=1:length(ta)
        for j=1:length(xa)
            C1(i,j)=C(1+i*Ns,1+j*Ns);
        end
    end
    surface(X,T,C1)
    ylabel('Time (s)')
    xlabel('Position (m)')
    zlabel('Concentration (mol/m^3)')
    function [c,f,s]=pdepb1(x,t,C,DCDx,Na,D)  %define pde system – see Matlab help for pdepe
    c=1;           %coefficient on time derivative
    f=D*DCDx;      %general flux   mol/s-m^3 (negative of diffusive flux)
    s=Na;         %sink   mol/s-m^3
    function C0=pbIC1(x,Xmax,Cw)   %initial condition
    if x==Xmax
      C0=Cw;
    else  
      C0=0.04;   %mol/m^3
    end  
    function [pl,ql,pr,qr]=pbBC1(yl,Cl,yr,Cr,t,Cw)  %boundary conditions
    pl=0;   %left side bc. This one sets flux to zero since the left side is impermeable
    ql=1;       %Form is pl + ql*f = 0 where f is flux from definition of pde system in function pdepb
    pr=Cr-Cw;       %right side bc. This one sets concentration to Cw
    qr=0;

MATLAB: Unsteady difussion in spherical coordinates boundary conditions pdepe

To define dC/dt = 0 at r=0 does not make sense. Maybe you mean dC/dr = 0 ?

If this is the case, try

pl=0;
ql=1;
pr=u;
qr=0;

Best wishes

Torsten.

Related Question