Diffusion and reaction in a simulated biofilm

This notebook has been written in  Mathematica by 

Mark J. McCready
Professor of Chemical and Biomolecular Engineering
University of Notre Dame
Notre Dame, IN 46556
USA

mjm@nd.edu
http://www.nd.edu/~mjm/


It is copyrighted to the extent allowed by which ever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version: 12/27/19
More recent versions of this notebook may be available at:
http://www.nd.edu/~mjm/tumor_calcs_1.nb

Use of the diffusion-reaction equation for transport through the shell and within the IL core

Diffusion through the film. of thickness b

b=1;

c0=1;

Pure diffusion

k=0;
α=1;

eqs={ D[c[x,t],t]  ==   α D[ D[c[x,t],x],x] - k   c[x,t], 
(D[c[x,t],x]/.x->b) ==0  ,
c[0,t] == c0 ,c[x,0] ==c0  UnitStep[-x]}

sol = NDSolve[{eqs},{c},{x,0,b},{t,0,20}]

Plot[Evaluate[((c[x,t]/.sol[[1]])/.t->{.01,.1,1})],{x,0,b},PlotRange->All,AxesLabel->{"x","c"},PlotLabel->"pure diffusion, α=1",PlotLegends->{"t=.01","t=.1","t=1"}]

c is getting stuck along the way by first order binding.

k=1;
α=1;

eqs={ D[c  [x,t],t]  ==   α D[ D[c[x,t],x],x] - k   c[x,t],
          D[cb[x,t],t]  ==                                                   + k   c[x,t], 
(D[c[x,t],x]/.x->b) ==0  ,
c[0,t] == c0 ,c[x,0] ==c0  UnitStep[-x],
                             cb[x,0] ==0}

sol = NDSolve[{eqs},{c,cb},{x,0,b},{t,0,20}]

Plot[Evaluate[((c[x,t]/.sol[[1]])/.t->{.01,.1,1})],{x,0,b},PlotRange->All,AxesLabel->{"x","c"},PlotLabel->"diffusion+ binding, α=1, k=1",PlotLegends->{"t=.01","t=.1","t=1"}]

Plot[Evaluate[((cb[x,t]/.sol[[1]])/.t->{.01,.1,1})],{x,0,b},PlotRange->All,AxesLabel->{"x","cb"},PlotLabel->"bound c, α=1 ,k=1",PlotLegends->{"t=.01","t=.1","t=1"}]

Plot[Evaluate[(((cb[x,t]+c[x,t])/.sol[[1]])/.t->{.01,.1,1})],{x,0,b},PlotRange->All,AxesLabel->{"x","cb"},PlotLabel->"total c, α=1 ,k=1",PlotLegends->{"t=.01","t=.1","t=1"}]

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.  Say 800 microns in diameter

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.  This gives a thickness of 20 microns

r1 =380 10^-4;

Length of time of calc.

t1 =50

The interface concentration of  CO2 in moles/cm^3    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Here are the boundary conditions for CO2.  No flow/symmetry at r=0 and a constant value on the outside.

bc1 = {Ab[r1,t] ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.01;
 (Subscript[α, 1]) = (0.000021752053691275165`  + 0.00001835 ( UnitStep[ r1-r]))*10;

Diffusion is faster in the IL

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r]       }

inits = {Abinit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab},{r,r1,R},{t,0,t1},MaxSteps->50000]

Plot[((Ab[r,t]/.sol[[1]])/.t->{.01,.02,.1,.5}),{r,r1,R},PlotRange->All]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,100,500}]

So it works!

Let’s assume that the process will be “shell-controlled”

((5.321696428571429`*^-7 moleCO2)/((cm^2) s))

In 1 cm^3 of processing equipment (for an idealized packed bed) we would have

1cm^3 (1-.37)/cm3bed /(4/3 π (200 10^-4 cm)^3/particle)

% 4 π (200 10^-4 cm)^2/particle

% ((5.321696428571429`*^-7 molCO2)/((cm^2) s))

Compare to a packed bed absorber with MEA

CO2 into MEA from Imperial

The mass transfer coefficient is about

.0013 kmole/hr/m^2/kPa

% 101 kPa

%/.{kmole->1000 mole, hr-> 3600 s,m->100 cm}

The prescribed contact area350 m^2/m^3

350 m^2/(cubicmeter)

%/.{m-> 100 cm}

%/.cubicmeter-> (100  cubiccm)^3

((3.6472222222222224`*^-6 mole)/((cm^2) s)) %

So for an idealized packed bed of 400 micron particles, the particles give about a factor of 5 on a volumetric basis increase over a gas-liquid MEA absorber.

Let’s look at equilibrium

Langmuir

y3eq = p3 Exp[-(dHrxn-T3 dSrxn)/R/T3] /(1+ p3 Exp[-(dHrxn-T3 dSrxn)/R/T3])

Table[Log[y3eq/.{dHrxn->-41100, T3->298, dSrxn->-125, R->8.314}],{p3, .01,.3,.01}]

y3eq/((1-y3eq))

FullSimplify[%]

%/.{dHrxn->-41100, T3->298, dSrxn->-125, R->8.314}

So the conclusion is that the rate on is very much bigger than the rate off!

What is a reasonable value for kon?

Burcu measured ~ 50000 l/mol/s

50000 l/mol/s

%/.l->1000 cm^3

% (1/400 moleIL/cm^3)  Ab0 moles CO2/cm^3

%/400.

2816./(1/400) /(mole/cm^3)

6 molesCO2/cm^3/s

% 4/3 π (200 10^-4 cm)^3/( 4 π  (200 10^-4 cm)^2)

Just Diffusion and reaction in the core.

If the viscosity is 50 cP

100 .001

1.14 10^-10 ((.05)^-.47) m^2/s

%/.m->100 cm

IL diffusivity == 3.4 10^-6, kon = 10000

So the time scale for the particle if just diffusion is

(180 10^-4 cm)^2 /(4.7 10^-6 cm^2/s)

Time scale of reaction  1/( k CIL)  If k=

50000. (l/mole/s)

%/.l->1000 cm^3

1/(%  (1/400. mole/cm^3))

%/.l->1000 cm^3

1000 g/l/(400 g/mole)

50000. (l/mole/s) %

1/%

1/250.

So we are not limited by reaction.

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag 

 

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 300

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-r1]

Ag is the concentration of ionic liquid.  If the density is about 1 g/cm^3 and the molecular weight is ~400 g/mol.

Aginit=Ag[r,0]==(1/400)

The total moles of IL and hence the total moles of CO2 that could be absorbed

Integrate[ 4 π r^2 1/400,{r,0, 180. 10^-4}]

ILcapacity = 1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

So in 300 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[r1,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=1000;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = 4.7 10^-6  (* the IL diffusivity *)

Diffusion is faster in the shell than the IL!... which is distressing.

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,r1},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->t1/100,(Ab[r,t]/.sol[[1]])/.t->t1/50,(Ab[r,t]/.sol[[1]])/.t->t1/10,(Ab[r,t]/.sol[[1]])/.t->t1/5,(Ab[r,t]/.sol[[1]])/.t->t1/2,(Ab[r,t]/.sol[[1]])/.t->t1},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

NIntegrate[4π  r^2 (Ab[r,t]/.sol[[1]])/.t->t1,{r,0,r1}]

totalmolesCO2in=NIntegrate[4.7 10^-6  D[(Ab[r,t]/.sol[[1]]),r]/.r->r1,{t,10,t1}] 4 π r1^2/(t1-10) t1

Plot[{(Ab[r,t]/.sol[[1]])/.t->t1/100,(Ab[r,t]/.sol[[1]])/.t->t1/50,(Ab[r,t]/.sol[[1]])/.t->t1/10,(Ab[r,t]/.sol[[1]])/.t->t1/5,(Ab[r,t]/.sol[[1]])/.t->t1/2,(Ab[r,t]/.sol[[1]])/.t->t1},{r,.8 r1,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->t1/10,(Ag[r,t]/.sol[[1]])/.t->t1/5,(Ag[r,t]/.sol[[1]])/.t->t1/2,(Ag[r,t]/.sol[[1]])/.t->t1},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->t1/2,(B[r,t]/.sol[[1]])/.t->t1/10,(B[r,t]/.sol[[1]])/.t->t1/5,(B[r,t]/.sol[[1]])/.t->t1/2,(B[r,t]/.sol[[1]])/.t->t1},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

FractionalILreacted = NIntegrate[4 π  r^2 (B[r,t]/.sol[[1]])/.t->t1,{r,0,r1}] moleIL/ILcapacity

Plot[ 4 10^-6 (D[Ab[r,t]/.sol[[1]],r]/.r->r1),{t,0,t1}]

The average flux is:

NIntegrate[3.4 10^-6 D[(Ab[r,t]/.sol[[1]]),r]/.r->r1,{t,200,t1}]/(t1-200)

1.44 10^-8 mole/cm^2/s

.37 m^3/(4 π /3 (.01 cm)^3) sphere

%/.m->100 cm

%/. sphere -> 4 π (.01 cm)^2

%4 %

%/.cm->.01 m

Just Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 100

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-r1]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[r1,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = 1 10^-5  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,r1},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%892

Show[%892,%647,%912]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.015,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%937

Show[%937,%967,%998]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[ 4 10^-6 (D[Ab[r,t]/.sol[[1]],r]/.r->r1),{t,0,t1}]

The average flux is:

NIntegrate[ 1 10^-5 D[(Ab[r,t]/.sol[[1]]),r]/.r->r1,{t,100,500}]/400

Just Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 100

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 300

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-r1]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[r1,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=10;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = 1 10^-5  (*The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,r1},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->200},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%892

Show[%892,%647,%912]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->300},{r,.015,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%937

Show[%937,%967,%998]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->200,(Ag[r,t]/.sol[[1]])/.t->300},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->200,(B[r,t]/.sol[[1]])/.t->300},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[ 4 10^-6 (D[Ab[r,t]/.sol[[1]],r]/.r->r1),{t,0,t1}]

The average flux is:

NIntegrate[ 1 10^-5 D[(Ab[r,t]/.sol[[1]]),r]/.r->r1,{t,50,300}]/250

Just Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 100

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 300

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-r1]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[r1,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = 1 10^-5  (*The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,r1},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%892

Show[%892,%647,%912]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.015,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%937

Show[%937,%967,%998]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->200,(Ag[r,t]/.sol[[1]])/.t->300},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->200,(B[r,t]/.sol[[1]])/.t->300},{r,r0,r1},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[ 4 10^-6 (D[Ab[r,t]/.sol[[1]],r]/.r->r1),{t,0,t1}]

The average flux is:

NIntegrate[ 1 10^-5 D[(Ab[r,t]/.sol[[1]]),r]/.r->r1,{t,50,300}]/250

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 100

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

6.10725611857856`*^-8 moleCO2 /(500 s) /( 4 π 200. 10^-4 cm)^2

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = (0.000021752053691275165`*1  - (0.000021752053691275165 *1-4.1612050095996035`*^-6) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%892

Show[%892,%647,%912]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.019,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

%937

Show[%937,%967,%998]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .95,{t,100,500}]/400

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 100

Change shell diff

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

6.10725611857856`*^-8 moleCO2 /(500 s) /( 4 π 200. 10^-4 cm)^2

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = (0.000021752053691275165` *10  - (0.000021752053691275165 *10-4.1612050095996035`*^-6) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.019,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165 10` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165 10`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .95,{t,100,500}]/400

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 1000

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

6.10725611857856`*^-8 moleCO2 /(500 s) /( 4 π 200. 10^-4 cm)^2

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=1000;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = (0.000021752053691275165`  - (0.000021752053691275165-4.1612050095996035`*^-6) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.019,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .95,{t,10,400}]/390

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 10000

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

6.10725611857856`*^-8 moleCO2 /(500 s) /( 4 π 200. 10^-4 cm)^2

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=10000;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = (0.000021752053691275165`  - (0.000021752053691275165-4.1612050095996035`*^-6) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.019,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .95,{t,10,400}]/390

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-6, kon = 100000

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

1/400 moleIL/cm^3 4/3 π (180. 10^-4 cm)^3

6.10725611857856`*^-8 moleCO2 /(500 s) /( 4 π 200. 10^-4 cm)^2

So in 500 sec if we react all the IL

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100000;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = (0.000021752053691275165`  - (0.000021752053691275165-4.1612050095996035`*^-6) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->3,(Ab[r,t]/.sol[[1]])/.t->20,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->200,(Ab[r,t]/.sol[[1]])/.t->500},{r,.019,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .95,{t,10,400}]/390

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-5, kon = 100

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

diffIL = 4.1612050095996035`*^-5;

(Subscript[α, shell]) = 0.000021752053691275165;

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = ((Subscript[α, shell])  - ((Subscript[α, shell])-diffIL) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500}]/500

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-5, kon = 100

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

diffIL = 4.1612050095996035`*^-5;

(Subscript[α, shell]) = 0.000021752053691275165;

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=100;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = ((Subscript[α, shell])  - ((Subscript[α, shell])-diffIL) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500}]/500

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-5, kon = 1000

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

diffIL = 4.1612050095996035`*^-5;

(Subscript[α, shell]) = 0.000021752053691275165;

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=1000;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = ((Subscript[α, shell])  - ((Subscript[α, shell])-diffIL) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

Plot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,100}]/100

Diffusion through the shell.  Diffusion and reaction in the core.

IL diffusivity == 4.2 10^-5, kon = 10000

((∂Ab)/(∂t)) = α(1/(r^2))((((∂^( )) )/(∂r))(r^2)((∂ Ab)/(∂r))) - (Subscript[k, on]) Ab Ag 


((∂Ag)/(∂t))  =  - (Subscript[k, on]) Ab Ag  

We will leave this dimensional for the moment

The particle radius in cm.

R =200 10^-4;

We need to go to r=0, but this is a numerical singularity, so we pick a small number, r0

r0=.00001

This is the inner core radius.. So we pick ~20 microns for the thickness

r1 =180 10^-4;

Length of time of calc.

t1 = 500

The interface concentration of  CO2.    This is 1 ATM of pure CO2

Ab0 = 1/(22400*273./298);

Abnot[t] =Ab0  ;

We need to make the initial profile of CO2.  This 0 except for the surface concentration.

Use the unitstep function to get all 0 except for x=1.

Abinit = Ab[r,0] == Ab0  UnitStep[r-R]

Ag is the concentration of ionic liquid.

Aginit=Ag[r,0]==(1/400) UnitStep[r1-r]

The total capacity of the core is hence

(1/400 moles/cm^3) 4/3 π (180. 10^-4 cm)^3

So an average flux could be

6.10725611857856`*^-8 moles/(4 π (200 10^-4 cm)^2)/(100 s)

This is complex IL and CO2.  (Doesn’t really need to be tracked)

Binit = B[r,0] == 0

Here are the boundary conditions for CO2.  No flux/symmetry at r=0 and a constant value on the outside.

bc1 = {( D[Ab[r,t],r]/.r->r0) ==0,Ab[R,t] == Ab0}

The values of the diffusivity and rate constants are

diffIL = 4.1612050095996035`*^-5;

(Subscript[α, shell]) = 0.000021752053691275165;

(*(Subscript[k, on])=    10000 .1(1-1 UnitStep[r-r1]);*)
(Subscript[k, on])=10000;
(Subscript[k, off])=.0001; (* allow for a little reversibility, we could adjust using real delta H*)
 (Subscript[α, 1]) = ((Subscript[α, shell])  - ((Subscript[α, shell])-diffIL) ( UnitStep[ r1-r]));  (*The first number is the shell.  The second number + first is the IL*)

Diffusion is faster in the shell than the IL!... which is distressing.

Plot[(Subscript[α, 1]),{r,0,R}]

The PDE's are

eqs={ D[Ab[r,t],t]  ==  ((Subscript[α, 1])/(r^2)) D[r^2  D[Ab[r,t],r],r] - (Subscript[k, on])   Ab[r,t]  Ag[r,t]+ (Subscript[k, off]) B[r,t]   ,
           D[Ag[r,t],t]  ==   -  (Subscript[k, on])   Ab[r,t]  Ag[r,t]  + (Subscript[k, off]) B[r,t]    ,
              D[B[r,t],t]  ==      (Subscript[k, on])   Ab[r,t]  Ag[r,t]  - (Subscript[k, off]) B[r,t]       }

inits = {Abinit,Aginit,Binit}

bcs = {bc1}

{eqs,inits,bcs}

sol = NDSolve[{eqs,inits,bcs},{Ab, Ag,B},{r,r0,R},{t,0,t1},MaxSteps->50000]

Plot[{(Ab[r,t]/.sol[[1]])/.t->1,(Ab[r,t]/.sol[[1]])/.t->10,(Ab[r,t]/.sol[[1]])/.t->50,(Ab[r,t]/.sol[[1]])/.t->100,(Ab[r,t]/.sol[[1]])/.t->300,(Ab[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "CO2"]

Plot[{(Ag[r,t]/.sol[[1]])/.t->1,(Ag[r,t]/.sol[[1]])/.t->10,(Ag[r,t]/.sol[[1]])/.t->50,(Ag[r,t]/.sol[[1]])/.t->100,(Ag[r,t]/.sol[[1]])/.t->300,(Ag[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Remaining IL"]

Plot[{(B[r,t]/.sol[[1]])/.t->1,(B[r,t]/.sol[[1]])/.t->10,(B[r,t]/.sol[[1]])/.t->50,(B[r,t]/.sol[[1]])/.t->100,(B[r,t]/.sol[[1]])/.t->300,(B[r,t]/.sol[[1]])/.t->500},{r,r0,R},PlotRange->All, PlotLegends->"Expressions",PlotLabel-> "Complexed IL"]

LogPlot[0.000021752053691275165` (D[Ab[r,t]/.sol[[1]],r]/.r->R),{t,0,t1}, PlotRange->All]

LogPlot[0.000021752053691275165` D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98,{t,0,500},PlotRange->All]

The average flux is:

NIntegrate[Abs[.000021752053691275165`D[(Ab[r,t]/.sol[[1]]),r]/.r->1/25 .98],{t,10,90}]/90