Importance of mass transfer fluctuations on reaction outcome in multiphase reactors.

This notebook has been written in Mathematica by 

Mark J. McCready
Professor and Chair of Chemical 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 whatever laws pertain to the World Wide Web and the Internet.

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

Version:  5/30/00
A more recent version of this notebook may be available at
http: // www . nd . 
                            edu/~mjm/
                        Reaction_Selectivity  _Mult . nb

Summary

This notebook is intended to show how periodic fluctuations in the mass transfer coefficient can cause not only differences in the total conversion but differences in the selectivity of the products if the fluctuation frequency is of the correct order of magnitude.  The intended application is gas-liquid reactions involving a solid phase catalyst such as trickle bed reactors.   

A CSTR (Continuous flow Stirred Tank Reactor) model is used.  First the model is demonstrated for a single phase reaction.  Next a simple liquid-solid model is used.  When fluctuations in the liquid-solid mass transfer coefficient are included, there is a change in selectivity with frequency.  Finally the complete gas-liquid-solid model of Wu, McCready and Varma (1995) is done.  This shows a maximum in selectivity at an intermediate frequency.

Mathematica aside

In Mathematica, it is convenient to give all expressions a "name".  I try to pick ones that are consistent with what is being done (but sometimes "temp" is used).  This assignment is done with an "=" sign.  To make an equation, a "==" is used.  This distinction is very useful in computer algebra and is employed in all of the packages with which I am familiar .

Input notation

I would enter the dynamic boundary condition from a key pad as

\(dc = \(-\ γ\)\ D[η[x, t], {x, 
              2}]\  - 
		\(ρ\_1\) \((\(-D[ϕ\_1[x, y, t], t]\) - 
              U1\ D[ϕ\_1[x, y, t], x] - 
              g\ η[x, 
                  t])\) + 
			\(ρ\_2\) \((\(-D[ϕ\_2[x, y, t], 
                  t]\) - 
					U2\ D[ϕ\_2[x, y, t], x] - 
              g\ η[x, t])\)\)

In Standard form, which is a little shorter, you would need the typeset window to make this practical.

\(dc = \(-\((γ\ ∂\_{x, 2}η[x, 
                  t])\)\) - ρ\_1\ \((\(-∂\_t ϕ\_1[x, y, 
                    t]\) - U1\ ∂\_x ϕ\_1[x, y, t] - 
              g\ η[x, 
                  t])\) + ρ\_2\ \((\(-∂\_t ϕ\_2[x, y, 
                    t]\) - U2\ ∂\_x ϕ\_2[x, y, t] - 
              g\ η[x, t])\)\)

Finally it may be more easily read in Traditional form as shown here.  Note that this particular cell is inactive and will not be evaluated.  You can change this, if you really want to, by selecting the cell, going into the "preferences" under the Edit menu and under Cell Options, Evaluation Options, set Evaluatable to True.

dc=-(γ (((∂^2) η(x, 
                      t))/(∂(x^2))))-(Subscript[ρ, 1]) (-((∂(Subscript[ϕ, 1])[x, y, 
                        t])/(∂t))-U1 ((∂(Subscript[ϕ, 1])[x, y, 
                        t])/(∂x))-g η(x, t))+(Subscript[ρ, 2]) (-((∂(Subscript[ϕ, 2])[x, y, 
                        t])/(∂t))-U2 ((∂(Subscript[ϕ, 2])[x, y, 
                        t])/(∂x))-g η(x, t))

Back to start of calculations.

Note that you can convert any "cryptic" input expression by selecting the cell, going up to the cell menu and selecting "convert to traditional form".  This is "shift curly t" on the keypad.

Keywords

Chemical reaction engineering, reaction selectivity, gas-liquid reaction, multiphase reactor, trickle-bed reactor, reactor dynamics

References

R. Wu, M. J. McCready and A. Varma (1995), Chemical Engineering Science, 50, pp3333-3334.   T. W. F. Russell and M. M. Denn (1972), Introduction to Chemical Engineering Analysis, Wiley.  R. W. Felder and R. W. Rousseau (1986) Elementary Principles of Chemical Processes, 2nd Edition. Wiley.

Development of the model equations and solution of the model

Basic problem set-up

Here is convenient CSTR model that is solved by Mathematica.

We have a CSTR (Continuous flow Stirred Tank Reactor) in which the reaction, A + B --> 2 D is carried out.  Second order reaction kinetics are assumed given by
(rate =  k ca cb).  The procedure for solving these kinds of problems is given in the references Russell and Denn or Felder and Rouseau.  It involves writing mass balances around the reactor for the chemical species and the total mass.  It is critical to realize that it is mass that is conserved.  However, for the examples shown here, there are no volume changes associated with reaction.  Thus  we will write equations that appear to conserve volume.  Note that this is not a general way of solving CSTR problems and you will have errors if you use it for gas phase reactions or liquid reactions with significant volume changes.

Consider a tank reactor that has inflows, qa and qb and outflow, q.  The overall mass (actually volume) balance is the differential equation:

veq = v'[t]  -  qa  - qb +  q  == 0

The change in volume per time is v'[t] and is nonzero if q does not equal qa + qb.

We can do a component balance for "ca".  The inflow is: (qa caf), the outflow is: q ca[t] and second order reaction that removes "a" is: k ca[t] cb[t] v[t].

caeq = 
      v[t] ca'[t]  - qa caf +  q ca[t] +  
          k ca[t] cb[t] v[t]  == 0

We also need a similar component balance for "cb" with inflow, (qb cbf), outflow (q cb[t]) and second order reaction k ca[t] cb[t] v[t], terms.

cbeq = 
      v[t] cb'[t]  - qb cbf +  q cb[t] +  
          k ca[t] cb[t] v[t]  == 0

We can do a component balance for "cd".  There is no inflow,  but it does have outflow (q cd[t]) and second order reaction 2 k ca[t] cb[t] v[t], terms.

cdeq = 
      v[t] cd'[t]   +  q cd[t] - 
          2  k ca[t] cb[t] v[t]  == 0

To recap and make a Mathematica list the equations are:

eqs = {veq, caeq, cbeq, cdeq, v[0] == 100, ca[0] == 0, 
        cb[0] == 0, cd[0] == 0}

Note that we have include some initial conditions in the reactor, v=100 and the concentrations are initially 0.

Transient reactor example

Here we make some substitutions for a transient reaction example.  At time =0, the feed compositions of qa and qb are set to caf= .1, and cbf = .2.  The flow rates are qa=10 and qb=10 and the exit is set to q = 20 so no volume change occurs in the reactor.  We expect that this example will have a transient period leading to a steady state.

eqsolv1 = 
      eqs /. { caf ->  .1, cbf ->  .2, qa -> 10, 
          qb -> 10, 
          q -> 20, 
k ->  .1}

Here the numerical ODE package solves them.

ans1 = 
      NDSolve[eqsolv1, {ca[t], cb[t], cd[t], v[t]}, {t, 0, 100}]

The command to plot an "interpolating function" requires the "Evaluate" command.  In this figure we are plotting the concentrations  a,  b, and  d as a function of time from the solution "ans1".

We see that the concentrations reach a steady state.  The highest one is "B", "A" is next and the product is lowest.

Plot[Evaluate[{ca[t], cb[t], cd[t]} /. ans1], {t, 0, 100}, 
        PlotStyle -> {Dashing[{ .02,  .02}], Dashing[{ .05,  .02}], 
            Dashing[{1, 0}]},  
        AxesLabel -> {"time", "concentration"}];

return to conclusions

Additional calculations for this part

As was mentioned above the behavior is a simple transient leading to a steady state.  For your interest you could:

1.  Vary some of the numbers to see how the transient time varies with the flow rates, reactor volume and reaction rate.

2. Change the form of the reaction expression to see its effect.

3.  Nondimensionalize the expression to find the dimensionless groups and use the values of the groups to note the relative importance of the reaction on the time scale of the transient.

Problem set up.

In this section we consider that there is a liquid and a solid phase.  The solid phase is a catalyst and it is mass transfer to this catalyst influences the reaction.  First this mass transfer will be steady.  Later it will be oscillatory. 

The reaction is now a + b --> d+ e.  The reaction to make e occurs only on the surface of the catalyst.  The reaction to make d occurs in the bulk.  We can now consider the selectivity of the two products by looking at the homogeneous versus heterogeneous reaction rates.  Note that mass transfer will play an important role as well in the synthesis of e.    

First we have an equation for a.  The new term is the transfer of a from the liquid phase to the solid phase as governed by concentration and a mass transfer coefficient ka.

aeq = 
      v ca'[t]  - qa caf +  q ca[t] +  k ca[t] cb[t] v  +  
          ka (ca[t] - cas[t])  == 0

Since there are two phases, we need to keep track of the "a" on the solid phase.  "cas" is the concentration of a on the catalyst.  The reaction term depletes a from the solid and the mass transfer term enhances it.

aseq =  
      vs cas'[t]  + 
 k2 cas[t] cbs[t] vs  -  
          ka (ca[t] - cas[t])  == 0

Here are the analogous equations for b.

beq = 
      v cb'[t]  - qb cbf +  q cb[t] +  
          k ca[t] cb[t] v  +  
 
          ka (cb[t] - cbs[t]) == 0

bseq = 
      vs cbs'[t]   +  
          k2 cas[t] cbs[t] v  -  
 
          ka (cb[t] - cbs[t]) == 0

Here are the two equations for the product, e.

eeq = 
      v ce'[t]   +  q ce[t]  + ka (ce[t] - ces[t]) == 
        0

eseq = 
      vs ces'[t]   -     k2 cas[t] cbs[t] vs  - 
          ka (ce[t] - ces[t]) == 0

Finally we need only 1 equation for d as it is formed in the liquid phase.  (We could consider d adsorbing onto the catalyst and doing nothing -- which it might.   But then again we have not considered an adsorption onto the walls of the reactor for a, b, e and this occurs also.

deq = 
      v cd'[t]   +  q cd[t] - 2  k ca[t] cb[t] v  == 
        0

Here is the complete two phase equation set.

eqstwophase = {aeq, aseq, beq, bseq, eeq, eseq, deq, ca[0] == 0, 
        cas[0] == 0,  cbs[0] == 0, cb[0] == 0, 
        ces[0] == 0, ce[0]  == 0, cd[0] == 0}

Solution to the transient problem with a steady mass transfer coefficient.

We set the initial volume of the reactor and the feed rates and concentrations.

eqssolv2 = 
      eqstwophase /. {v -> 100,  caf ->  .1, cbf ->  .2, 
          qa -> 10, qb -> 10, 
          q -> 20, 
k ->  .1, k2 -> 100, 
          ka -> 1000, vs -> 300}

Now solve these numerically.

ansTP = 
      NDSolve[eqssolv2, {ca[t], cas[t], cb[t], cbs[t], cd[t], ce[t], 
          ces[t]}, {t, 0, 100}]

Here is a plot of the answers for the liquid phase concentrations.  a is the solid line, b  is the short dash, d cannot be seen and e  is the long-short dash. s

Plot[Evaluate[{ca[t], cb[t], cd[t], ce[t]} /. ansTP], {t, 0, 100}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .05,  .05}], Dashing[{ .06,  .02,  .02,  .02}]}, 
        AxesLabel -> {"time",  "concentration"}];

return to conclusions

For these conditions we do not make much d.

Plot[Evaluate[{ cd[t]} /. ansTP], {t, 0, 100} ];

Note that the a concentration initially increases until the reaction rate is fast enough to reduce it.

Plot[Evaluate[{ ca[t]} /. ansTP],{t, 0, 100}, PlotLabel->{"\!\(\*
StyleBox[\"a\",\nFontSlant->\"Italic\"]\) transient concentration"} ];

Here are the concentrations on the catalyst.

Plot[Evaluate[{cas[t], cbs[t], ces[t]} /. ansTP], {t, 0, 100}, 
        PlotRange -> All, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .06,  .02,  .02,  .02}]}];

Additional calculations for this part

1.  The relative amounts of d and e can be varied by changing the mass transfer coefficient and the catalyst loading.  Do these work exactly together.    

2. Change the relative values of the reaction rate constants.

3.  Calculate the average concentrations across the transient period.

4.  Look at the effect of different initial concentrations in the reactor.  Consider a reactor startup in which you wish to minimize the formation of d, how do you set the initial concentrations in the reactor?

5.  Make the equations dimensionless and see the size of the relative effects by calculating the values of the dimensionless groups.  How many groups are there now?

Problem set up

We now wish to let the mass transfer rate vary in time.  Choosing a single frequency may not be completely realistic, but it is simple and will show the effect.  

Note that a mass transfer coefficient never is negative so we need something a little more complex than a sine wave.  Here we have a background minimum value of ka, a frequency ω and a maximum value, kb.

masst[ω_, t_, ka_,  kb_] :=  
      If[Sin[ω t] ≥ 0, ka  +  kb Sin[ω t], 
        ka]

Plot[
        masst[ω, t, ka,  kb] /. {ka -> 3, 
            kb -> 2, ω -> 5}, {t, 0, 10}, 
        AxesLabel -> {"time", "mass transfer \
coefficient"}];

Here we change the values.

Plot[
        masst[ω, t, ka,  kb] /. {ka -> 0, 
            kb -> 9, ω -> 2}, {t, 0, 10}];

This seems to work well for the mass transfer coefficient.  We can now just incorporate it into the model we have been using.

The two-phase fluctuating equations are:

aeqfluc = 
      v ca'[t]  - qa caf +  q ca[t] +  k ca[t] cb[t] v  +  
          masst[ω, t, ka,  kb] (ca[t] - cas[t])  == 
        0

beqfluc = 
      v cb'[t]  - qb cbf +  q cb[t] +  
          k ca[t] cb[t] v  +  
 
          masst[ω, t, ka,  kb] (cb[t] - cbs[t]) == 
        0

aseqfluc =  
      vs cas'[t]  + 
 k2 cas[t] cbs[t] vs  -  
          masst[ω, t, ka,  kb] (ca[t] - cas[t])  == 
        0

bseqfluc = 
      vs cbs'[t]   +  
          k2 cas[t] cbs[t] vs  -  
 
          masst[ω, t, ka,  kb] (cb[t] - cbs[t]) == 
        0

eeqfluc = 
      v ce'[t]   +  q ce[t]  + 
          masst[ω, t, ka,  kb] (ce[t] - ces[t]) == 
        0

eseqfluc = 
      vs ces'[t]   -     k2 cas[t] cbs[t] vs  - 
          masst[ω, t, ka,  kb] (ce[t] - ces[t]) == 
        0

deqfluc = 
      v cd'[t]   +  q cd[t] - 2  k ca[t] cb[t] v  == 
        0

eqstwophasefluc = {aeqfluc, aseqfluc, beqfluc, bseqfluc, eeqfluc, 
        eseqfluc, deqfluc, ca[0] == 0, cas[0] == 0,  
        cbs[0] == 0, cb[0] == 0, ces[0] == 0, 
        ce[0]  == 0, cd[0] == 0}

This example shows the selectivity enhancement with frequency.

Here is the calculation of the solution for ω =1.

ans1TPfluc = 
      NDSolve[eqstwophasefluc /. {v -> 100,  caf ->  .2, 
            cbf ->  .22, qa -> 1, qb -> 1, 
            q -> 2, 
k -> 3, k2 -> 10, 
            ka -> 10, kb -> 200, 
            vs -> 100, ω -> 1}, {ca[t], cas[t], cb[t], 
          cbs[t], cd[t], ce[t], ces[t]}, {t, 0, 1000}, 
        MaxSteps -> 10000]

Here is the plot of the answers.  Note that d and e are about the same for this frequency, ω=1.

Plot[
        Evaluate[{ca[t], cb[t], cd[t], ce[t]} /. ans1TPfluc], {t, 0, 600}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .05,  .05}], Dashing[{ .06,  .02,  .02,  .02}]}, 
        AxesLabel -> {"time",  "concentration"}];

Here is just d and e.

Plot[Evaluate[{cd[t], ce[t]} /. ans1TPfluc], {t, 0, 600}, 
        PlotRange -> All, 
        PlotStyle -> { Dashing[{ .05,  .05}], 
            Dashing[{ .06,  .02,  .02,  .02}]}, 
        AxesLabel -> {"time",  "concentration"}];

\(\(Plot[Evaluate[{cas[t], cbs[t], ces[t]} /. ans1TPfluc], {t, 0, 600}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .06,  .02,  .02,  .02}]}, 
        AxesLabel -> {"time", \ "concentration on catalyst"}];\)\)

We see that e is highest on the catalyst also.

Let's see what the steady state concentrations of d and e are.

NIntegrate[
      Evaluate[(cd[t] /. ans1TPfluc[[1]])]/200, {t, 800, 
        1000}]

NIntegrate[
      Evaluate[(ce[t] /. ans1TPfluc[[1]])]/200, {t, 800, 
        1000}]

We see that d is slightly favored.

Now change the frequency to ω=0.5,

ans2TPfluc = 
      NDSolve[eqstwophasefluc /. {v -> 100,  caf ->  .2, 
            cbf ->  .22, qa -> 1, qb -> 1, 
            q -> 2, 
k -> 3, k2 -> 10, 
            ka -> 10, kb -> 200, 
            vs -> 100, ω ->  .5}, {ca[t], cas[t], cb[t], 
          cbs[t], cd[t], ce[t], ces[t]}, {t, 0, 1000}, 
        MaxSteps -> 10000]

Here is the plot of the answers.  Note that d and e are about the same for this frequency, ω=1.

Plot[
        Evaluate[{ca[t], cb[t], cd[t], ce[t]} /. ans2TPfluc], {t, 0, 600}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .05,  .05}], Dashing[{ .06,  .02,  .02,  .02}]}, 
        AxesLabel -> {"time",  "concentration"}];

Here is just d and e.

Plot[Evaluate[{cd[t], ce[t]} /. ans2TPfluc], {t, 0, 600}, 
        PlotRange -> All, 
        PlotStyle -> { Dashing[{ .05,  .05}], 
            Dashing[{ .06,  .02,  .02,  .02}]}, 
        AxesLabel -> {"time",  "concentration"}];

Let's see what the steady state concentrations of d and e are.

NIntegrate[
      Evaluate[(cd[t] /. ans2TPfluc[[1]])]/200, {t, 800, 
        1000}]

NIntegrate[
      Evaluate[(ce[t] /. ans2TPfluc[[1]])]/200, {t, 800, 
        1000}]

We see that d has become more favored.  Now let's increase the frequency.

Now change the frequency to ω=5, this takes a while to calculate and the maximum steps needs to be increased because of the higher fluctuation frequency.

ans3TPfluc = NDSolve[eqstwophasefluc /. {v -> 100, caf -> 0.2,cbf -> 0.22, qa -> 1, qb -> 1, q -> 2, k -> 3, k2 -> 10, ka -> 10, kb -> 200, vs -> 100, ω -> 5}, 
   {ca[t], cas[t], cb[t], cbs[t], cd[t], ce[t], ces[t]}, 
   {t, 0, 1000}, MaxSteps -> 30000]

Here is the plot of the answers.  Note that d and e are about the same for this frequency, ω=1.

Plot[Evaluate[{ca[t], cb[t], cd[t], ce[t]} /. ans3TPfluc], 
   {t, 0, 600}, PlotStyle -> {Dashing[{0., 0.}], 
     Dashing[{0.01, 0.01}], Dashing[{0.05, 0.05}], 
     Dashing[{0.06, 0.02, 0.02, 0.02}]}, 
   AxesLabel -> {"time", "concentration"}];

Here is just d and e.  Now we see that e is favored.

Plot[Evaluate[{cd[t], ce[t]} /. ans3TPfluc], {t, 0, 600}, 
   PlotRange -> All, PlotStyle -> {Dashing[{0.05, 0.05}], 
     Dashing[{0.06, 0.02, 0.02, 0.02}]}, 
   AxesLabel -> {"time", "concentration"}];

return to conclusions

Let's see what the steady state concentrations of d and e are.

NIntegrate[Evaluate[cd[t] /. ans3TPfluc[[1]]]/200, {t, 800, 1000}]

NIntegrate[Evaluate[ce[t] /. ans3TPfluc[[1]]]/200, {t, 800, 1000}]

We see that e now has become favored.Changing the frequency of fluctuation of the mass transfer coefficient can affect the reaction selectivity!

We should be sure that the frequency change does not affect the average mass transfer coefficient.

NIntegrate[
        masst[ω, t, ka, kb] /. {ka -> 10, 
            kb -> 200, ω -> 1}, {t,  0, 100}]/
      100

NIntegrate[
        masst[ω, t, ka, kb] /. {ka -> 10, 
            kb -> 200, ω ->  .5}, {t,  0, 100}]/
      100

NIntegrate[
        masst[ω, t, ka, kb] /. {ka -> 10, 
            kb -> 200, ω ->  .5}, {t,  0, 100}, 
        PrecisionGoal -> 20]/100

The numbers are not exactly the same, but they should be.  The differences are in the accuracy of the numerical integration.  We can fix this by doing the integral analytically.

Integrate[
        masst[ω, t, ka, kb] /. {ka -> 10, 
            kb -> 200, ω -> 1}, {t,  0, 100} ]/
      100

Now get the number.

N[%]

Now do another

Integrate[
        masst[ω, t, ka, kb] /. {ka -> 10, 
            kb -> 200, ω -> 1/2}, {t,  0, 100} ]/
      100

N[%]

Integrate[
        masst[ω, t, ka, kb] /. {ka -> 10, 
            kb -> 200, ω -> 5}, {t,  0, 100}]/
      100

N[%]

We see that the frequency does not affect the average value of the mass transfer coefficient.

Additional calculations for this part

We expect that the fluctuation frequency will not always be having an influence on the reaction outcome.  How can we explore this.

1.  At very low frequency we might find that the reaction is always at a pseudo state at the current value of the mass transfer coefficient.  Check to see if this is true.

2.  At a very high frequency, we might find that the reaction system is not capable of responding to any fluctuations and thus behaves as if the value of the transfer coefficient is just its average.  Check this limit.

3.  Make the equations dimensionless, or just figure it out, and tell how the "high" or "low" frequency is determined.

Here is the first example (Case I) from Wu, McCready and Varma, (1995).  In this example, A goes to either B or C depending upon which catalyst site, S or S* is involved.

The application of their research is to cocurrent gas-liquid, packed bed ("trickle bed") reactors.  These reactors exhibit a range of different flow regimes.  A flow regime is the term for the static location and dynamic behavior of the individual phases.  An important regime for these reactors is the "pulse" flow regime where strong disturbances travel through the bed.  These have a well-defined frequency and they cause strong increases in the mass transfer as they flow past a certain point.  Thus incorporation of a fluctuating mass transfer coefficient in the reaction model should help to describe the behavior seen in these reactors.

Reaction Mechanism

Reaction Mechanism involving a gaseous reactant A, two catalyst sites, S and S*, and two products B and C.  

A(g)  +  S  (Subscript[(Subscript[(Underscript[⇌ , k]), -1]^k), 1]) AS
AS + A(g) (Overscript[⟶, (Subscript[k, 2])])2B + S

A(g) + S* \(TraditionalForm\`\( ⇌ \+\((k')\)\_\(-1\)\%\((k')\)\_2\)\) AS*

AS* (Overscript[⟶, (Subscript[k, 2])])C + S*

Problem set up

The eight dimensionless equations are given below along with the definitions of the variables.  We see that these comprise a set that is much more complex than the simple example given above, but we will be able to solve it just as quickly as the previous examples.  

The first three equations are standard CSRT equations for A, B and C.  The next three equations are the A, B and C equations on the solid catalyst.  The final two equations are for the number of empty sites of S and S*.

Here is the page from Wu et al. (1995) with the reaction equations.

Here are the equations in a better form for printing

Define the fluctuating mass transfer coefficient

masst[ω_, t_, ka_,  kb_] :=  
      If[Sin[ω t] ≥ 0, ka  +  kb Sin[ω t], 
        ka]

Plot[
        masst[ω, t, ka,  kb] /. {ka ->  .01, 
            kb ->  .06, ω ->  .01}, {t, 0, 
          4000}];

masstave = 
      Integrate[
          masst[ω, t, ka,  kb] /. {ka ->  .01, 
              kb ->  .01*6, ω ->  .1}, {t, 0, 3000}]/
        3000

The component equations translated from the picture above are

eq1 = -x1'[t] + κ (1 - x1[t])  - 
            masst[ω, t, ka,  kb] (x1[t] - x4[t])  == 
          0;

eq2 = -x2'[t]  - κ x2[t]  + 
            masst[ω, t, ka,  kb] (x5[t] - x2[t]) == 
          0;

eq3 = -x3'[t]  -  κ x3[t] - 
            masst[ω, t, ka,  kb] (x3[t] - x6[t]) == 
          0;

eq4 = -x4'[
                t] + ϵ/(1 - ϵ)/ϵp masst[ω, t, ka,  kb] (x1[t] - x4[t])  + 
            m1 (-k1 x4[t] x7[
                      t]  + 
(km1 - k2 x4[t]) (1 - 
                        x7[t])) + 
         m2 (-kp1 x4[t] x8[t]  + kpm1 (1 - x8[t]))  == 
          0;

eq5 = -x5'[
                t]  - ϵ/(1 - ϵ)/ϵp masst[ω, t, ka,  kb]  (x5[t] - x2[t])  + 
            2 m1 k2 x4[t] (1 - x7[t])  == 0;

eq6 = -x6'[
                t]  + ϵ/(1 - ϵ)/ϵp masst[ω, t, ka,  kb]  (x3[t] - x6[t])  + 
            m2 kp2  (1 - x8[t])  == 0;

eq7 = -x7'[t] - 
            k1 x4[t] x7[
                t]  + (km1 + k2 x4[t]) (1 - x7[t]) == 
          0;

eq8 = -x8'[t]  - 
            kp1 x4[t] x8[t]  + (kpm1 + kp2 ) (1 - x8[t]) == 
          0;

Here is one to check

eq7

Here are all of the equations and initial conditions.

eqsx={eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, 
        x1[0] ==  .00, 
x2[0] == 0.00, 
        x3[0] == 0.00, x4[0] ==  .00, x5[0] == 0.00, 
        x6[0] == 0.00, 
x7[0] == 0.00,  
        x8[0] == 0.00}

Here is the base case from Wu et al.

We save all of the intermediate answers for ω = 0,0.01, etc. and make a final selectivity plot.

As we vary the frequency we need to adjust the time range and maximum steps to get converged answers for each of the frequencies.  If you don't like my answers, make some adjustments to make them better.

First solve for no fluctuations, ω =0.

eqsz1 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05,  
          k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0}

ansz1 = 
      NDSolve[eqsz1, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 10000]

Plot[Evaluate[{x1[t]} /. ansz1], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

Plot[Evaluate[{x2[t]} /. ansz1], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans1 = 
      NIntegrate[(x2[t] /. ansz1)[[1]] , {t, 1200, 1500}]/
        300

Plot[Evaluate[{x3[t]} /. ansz1], {t, 0, 1000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x3ans1 = 
      NIntegrate[(x3[t] /. ansz1)[[1]] , {t, 900, 1000}]/
        100

select[1] =  x2ans1/x3ans1

Now choose ω = 0.01

eqsz2 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.01}

ansz2 = 
      NDSolve[eqsz2, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 10000}, 
        MaxSteps -> 50000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz2], {t, 0, 10000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans2 = 
      NIntegrate[(x2[t] /. ansz2)[[1]], {t, 9000, 10000}]/
        1000

x3ans2 = 
      NIntegrate[(x3[t] /. ansz2)[[1]], {t, 9000, 10000}]/
        1000

select[2] =  x2ans2/x3ans2

Now choose ω = 0.02

eqsz3 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.02}

ansz3 = 
      NDSolve[eqsz3, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 4000}, 
        MaxSteps -> 15000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz3], {t, 0, 4000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans3 = 
      NIntegrate[(x2[t] /. ansz3)[[1]], {t, 2000, 3000}]/
        1000

x3ans3 = 
      NIntegrate[(x3[t] /. ansz3)[[1]], {t, 2000, 3000}]/
        1000

select[3] =  x2ans3/x3ans3

Now choose ω = 0.05

eqsz4 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.05}

ansz4 = 
      NDSolve[eqsz4, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 4000}, 
        MaxSteps -> 15000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz4], {t, 0, 4000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans4 = 
      NIntegrate[(x2[t] /. ansz4)[[1]], {t, 2000, 3000}]/
        1000

x3ans4 = 
      NIntegrate[(x3[t] /. ansz4)[[1]], {t, 2000, 3000}]/
        1000

select[4] =  x2ans4/x3ans4

Now choose ω = 0.1

eqsz5 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.1}

ansz5 = 
      NDSolve[eqsz5, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1400}, 
        MaxSteps -> 20000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz5], {t, 0, 1400}, 
        PlotStyle -> {GrayLevel[0], Dashing[{ .02}]}];

x2ans5 = 
      NIntegrate[(x2[t] /. ansz5)[[1]], {t, 1200, 1400}]/
        200

x3ans5 = 
      NIntegrate[(x3[t] /. ansz5)[[1]], {t, 1200, 1400}]/
        200

select[5] = x2ans5/x3ans5

Now choose ω = 0.2

eqsz6 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω ->  .2}

ansz6 = 
      NDSolve[eqsz6, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 26000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz6], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans6 = 
      NIntegrate[(x2[t] /. ansz6)[[1]], {t, 1300, 1500}]/
        200

x3ans6 = 
      NIntegrate[(x3[t] /. ansz6)[[1]], {t, 1300, 1500}]/
        200

select[6] =  x2ans6/x3ans6

Now choose ω = 0.5

eqsz7 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω ->  .5}

ansz7 = 
      NDSolve[eqsz7, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 26000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz7], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans7 = 
      NIntegrate[(x2[t] /. ansz7)[[1]], {t, 1300, 1500}]/
        200

x3ans7 = 
      NIntegrate[(x3[t] /. ansz7)[[1]], {t, 1300, 1500}]/
        200

select[7] =  x2ans7/x3ans7

Now choose ω = 1

eqsz8 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 10, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 1}

ansz8 = 
      NDSolve[eqsz8, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 26000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz8], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans8 = 
      NIntegrate[(x2[t] /. ansz8)[[1]], {t, 1300, 1500}]/
        200

x3ans8 = 
      NIntegrate[(x3[t] /. ansz8)[[1]], {t, 1300, 1500}]/
        200

select[8] =  x2ans8/x3ans8

Now compare these as a function of frequency

Table[select[i], {i, 2, 8}]

selectivityenhance = 
      Table[100*(select[i]/select[1]  - 1), {i, 2, 8}]

omegas = {  .01,  .02,  .05,  .1,  .2,  .5, 1}

Needs["Graphics`Graphics`"]

\(\(LogLinearListPlot[
        Transpose[{omegas, 
            selectivityenhance}], 
AxesLabel -> {"Frequency", \ " % Selectivity Enhancement"}, 
        PlotRange -> {0, 50}, 
        PlotStyle -> {PointSize[ .03]}];\)\)

return to conclusions

We see that there is an intermediate frequency where the enhancement is a maximum.

Here is a second case from Wu et al.

This is case (2) from their paper.  The parameter (Subscript[m, 1]) is reduced to 1 from 10.  This corresponds to a reduction in the number of catalytic cites necessary for the desired reaction.  This reduction in mass transfer area for the desired reaction makes the effect of fluctuation more effective.

We save all of the intermediate answers for ω = 0,0.01, etc. and make a final selectivity plot.

First solve for no fluctuations, ω =0.

eqsz1 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05,  
          k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0}

ansz1 = 
      NDSolve[eqsz1, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 10000]

Plot[Evaluate[{x1[t]} /. ansz1], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

Plot[Evaluate[{x2[t]} /. ansz1], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans1 = 
      NIntegrate[(x2[t] /. ansz1)[[1]] , {t, 1200, 1500}]/
        300

Plot[Evaluate[{x3[t]} /. ansz1], {t, 0, 1000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x3ans1 = 
      NIntegrate[(x3[t] /. ansz1)[[1]] , {t, 900, 1000}]/
        100

select[1] =  x2ans1/x3ans1

Now choose ω = 0.01

eqsz2 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.01}

ansz2 = 
      NDSolve[eqsz2, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 10000}, 
        MaxSteps -> 50000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz2], {t, 0, 10000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans2 = 
      NIntegrate[(x2[t] /. ansz2)[[1]], {t, 9000, 10000}]/
        1000

x3ans2 = 
      NIntegrate[(x3[t] /. ansz2)[[1]], {t, 9000, 10000}]/
        1000

select[2] =  x2ans2/x3ans2

Now choose ω = 0.02

eqsz3 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.02}

ansz3 = 
      NDSolve[eqsz3, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 4000}, 
        MaxSteps -> 15000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz3], {t, 0, 4000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans3 = 
      NIntegrate[(x2[t] /. ansz3)[[1]], {t, 2000, 3000}]/
        1000

x3ans3 = 
      NIntegrate[(x3[t] /. ansz3)[[1]], {t, 2000, 3000}]/
        1000

select[3] =  x2ans3/x3ans3

Now choose ω = 0.05

eqsz4 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.05}

ansz4 = 
      NDSolve[eqsz4, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 4000}, 
        MaxSteps -> 15000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz4], {t, 0, 4000}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans4 = 
      NIntegrate[(x2[t] /. ansz4)[[1]], {t, 2000, 3000}]/
        1000

x3ans4 = 
      NIntegrate[(x3[t] /. ansz4)[[1]], {t, 2000, 3000}]/
        1000

select[4] =  x2ans4/x3ans4

Now choose ω = 0.1

eqsz5 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 0.1}

ansz5 = 
      NDSolve[eqsz5, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1400}, 
        MaxSteps -> 20000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz5], {t, 0, 1400}, 
        PlotStyle -> {GrayLevel[0], Dashing[{ .02}]}];

x2ans5 = 
      NIntegrate[(x2[t] /. ansz5)[[1]], {t, 1200, 1400}]/
        200

x3ans5 = 
      NIntegrate[(x3[t] /. ansz5)[[1]], {t, 1200, 1400}]/
        200

select[5] = x2ans5/x3ans5

Now choose ω = 0.2

eqsz6 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω ->  .2}

ansz6 = 
      NDSolve[eqsz6, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 26000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz6], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans6 = 
      NIntegrate[(x2[t] /. ansz6)[[1]], {t, 1300, 1500}]/
        200

x3ans6 = 
      NIntegrate[(x3[t] /. ansz6)[[1]], {t, 1300, 1500}]/
        200

select[6] =  x2ans6/x3ans6

Now choose ω = 0.5

eqsz7 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω ->  .5}

ansz7 = 
      NDSolve[eqsz7, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 26000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz7], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans7 = 
      NIntegrate[(x2[t] /. ansz7)[[1]], {t, 1300, 1500}]/
        200

x3ans7 = 
      NIntegrate[(x3[t] /. ansz7)[[1]], {t, 1300, 1500}]/
        200

select[7] =  x2ans7/x3ans7

Now choose ω = 1

eqsz8 = 
      eqsx /. {κ ->  .01, ka ->  .01, kb ->  .05, 
          kc ->  .01,  k1 -> 1, km1 -> 1, 
          k2 -> 1, 
kp1 -> 1, k1p -> 1, 
          kpm1 -> 1, kp2 -> 1, m1 -> 1, 
          m2 -> 
            1, ϵ ->  .5,  ϵp ->  .3, ω -> 1}

ansz8 = 
      NDSolve[eqsz8, {x1[t], x2[t], x3[t], x4[t], 
          x5[t], 
x6[t], x7[t], x8[t]}, {t, 0, 1500}, 
        MaxSteps -> 26000]

Plot[Evaluate[{x2[t], x3[t]} /. ansz8], {t, 0, 1500}, 
        PlotStyle -> {Dashing[{ .00,  .00}], Dashing[{ .01,  .01}], 
            Dashing[{ .03,  .03}], Dashing[{ .04,  .04}]}];

x2ans8 = 
      NIntegrate[(x2[t] /. ansz8)[[1]], {t, 1300, 1500}]/
        200

x3ans8 = 
      NIntegrate[(x3[t] /. ansz8)[[1]], {t, 1300, 1500}]/
        200

select[8] =  x2ans8/x3ans8

Now compare these as a function of frequency

Table[select[i], {i, 2, 8}]

selectivityenhance = 
      Table[100*(select[i]/select[1]  - 1), {i, 2, 8}]

omegas = {  .01,  .02,  .05,  .1,  .2,  .5, 1}

Needs["Graphics`Graphics`"]

\(\(LogLinearListPlot[
        Transpose[{omegas, 
            selectivityenhance}], 
AxesLabel -> {"Frequency", \ " % Selectivity Enhancement"}, 
        PlotRange -> {0, 80}, 
        PlotStyle -> {PointSize[ .03]}];\)\)

For this case, where the available cites for the desired reaction is reduced, the intermediate maximum is even larger.

Now what else can we do?

The Wu et al. paper did about 10 cases for this example so most of the obvious effects are show.  Two big questions are not explicitly answered in their paper however.

1.  What is responsible for the occurrence of a maximum in selectivity at an intermediate frequency?

2.  How come the maximum occurs at ω = 0.1 rather than ω=1?

Conclusions

1.  If you are just interested in transient and steady state CSTR calculations that are just a little too complex to solve analytically, we show how to easily do this.  The transient CSTR calculations for a single phase reactor show a simple transient leading to a steady state.  

2.  If the model includes a second phase, there can be a competition between products that are produced in different phases.   

3.  For a two phase system, there can be a change in selectivity with frequency.    

4.  For a gas-liquid reaction with a solid catalyst.  The reaction selectivity can exhibit a maximum at an intermediate frequency.  This effect can have important consequences for reactor scale-up as the dominant frequency for pulses decreases significantly as the column diameter is increased.