Solving diffusion- reversible chemical reaction equation: How to solve both terms numerically?

If you know the form of the kinetic equation for this chemical reaction, then you put it on the right side of the partial differential equation (accumulation = mass transport - kinetics). In the simplest situation: kinetics = r = k1cAcB-k2cC. Obviously, cB and cC are derived from stoichiomertic relations. Mass transport is not separated from the chemical reaction (the chemical reaction runs parallel to the internal mass transport).

Mirosław Grzesik Well, Lets say that I am solving for A, then the chemical kinetic term will be in term of K1CACB-K2CC which is nonlinear. I tried to use Euler explicit method to solve for accumulation and mass transport and Newtons method for reaction term (for simultaneous nonlinear equation), but I am not sure since it became extremely complicated with many undefined variables.

Matt Mj, Let us specify whether you solve the parabolic equation once for the given parameter values, or you solve this equation multiple times, trying to fit this equation to the experimental data set and thus determining the value of its parameters (inverse problem). The nonlinearity of the kinetic term, especially when taking the temperature into account, can cause difficulties (stiff equation). The diffusion model other than Fick's model can also pose problems. The kinetic term, however, does not require the Newton's method. What "simultaneous nonlinear equation" you are writing about. Also, specify what the equation is about: process in catalyst grain or process in a reactor.

Mirosław Grzesik I am trying to simulate the spatiotemporal evolution of species' concentrations which occur within a polymer assuming A is a gas around the polymer and equilibrate on the surface via Henry's law and B is a reactive site within the polymer. The temperature is constant and yes I am considering Fick's second law for mass transport.

@Matt you can simplify the kinetic reaction equation by considering solid reactant B as the excess reactant and express the kinetic rate law as dependent on reactant A (Pseudo reaction rate), then Langmiur-Hinselwood kinetic model expression . To add the mass transfer resistance (external film, molecular diffusion, surface diffusion) determines the concentration of reactant A available to take part in the reaction. The significance of internal mass transfer (pore) resistance can be determine by calculating the Thiele modulus for the system in question. Also the Fick's second law considered makes the diffusion-reaction model non-linear because the differential of the concentration variable is second order multiplied by the diffusion coefficient

Matt Mj, And such a complicated, probably multi-parameter problem can be described by a single parabolic equation? I wonder under what assumptions? After all, there are two phases in adsorption equilibrium (Henry's law). Well, so you want to treat adsorption as a chemical reaction with the scheme A + S AS. It can be both linear and non-linear kinetics. However, I still do not know where this system of nonlinear equations came from, that you want to solve with Newton's method.

Mirosław Grzesik I meant the reaction term as a nonlinear since it is second order for forward reaction. I will look into multi-species diffusion reaction. Thank you for your generous explanation.

Karinate Okiy Do you mean considering k'=k1CB? but still we have Cc as the kinetic reaction equation= -k'CA+k2Cc.

Matt Mj, The reaction/pseudo-reaction heading towards AS may or may not be second order. It all depends on the assumptions (e.g. compliance with mass action law) and subsequent model verification. I would like to add that you will not get the Henry isotherm based on the adsorption scheme: A + S AS, only for the scheme: A AS. So the second order goes against Henry's law. Well, What I wrote above was just a guess. You have not provided information on gas flow geometry, it is not known whether the polymer is porous, and there is no information on hydrodynamics. Finally, perhaps your process can be described by the Bohart-Adams equation, but without knowing the details, it's still a guess.

@Matt there is need for determining the reaction mechanism or is it just an equilibrium reaction. I think you should consider adsorption and desorption schemes separately. Surface coverage of the polymer substrate by the adsorbed A (Oa), so in the langmuir hinshelwood kinetic model R=E*k'Oa where initial concentration of adsorbed A (Oa) is zero for adsorption and 1 for the desorption schemes respectively in the diffusion-reaction model.

Oa = Ka*Ca/1+Ka*Ca

Where E (epsilon) is the effectiveness factor (inverse of the Thiele modulus)

Karinate Okiy I am trying to simulate temporal evolution of D-R equation of the reactive species in a single bead. I am trying to understand the behavior of the species within the polymer So I am not looking at Langumuir model or lets say adsorption/desorption. But thanks for the suggestions.

Matt Mj your model will consist of diffusion terms (Fick's second law) and chemical reaction terms-parabolic partial differential equations. This will require numerical solution.

An overall mass conservation equation, individual species mass conservation equations and the initial conditions and boundary conditions will complete the system of PDE's.

Finite volume method or Matlab PDE solver are possible approaches to finding the concentration of the species as a function of time and position. Just my suggestions, hope this will help

Yes, I have already identified all IC's, BC's and assumptions for each specie I have and used finite difference method to discretize the domain for mass conservation equation. The only cliff which I see prevents me from completing the equation to be ready to simulate is kinetic reversible reaction equation which is coupled with multiple species concentrations and as we know all of them are changing. For example, if I want to evaluate A, using FTCS method: Ca(Xi,Tj+1)=Ca(Xi,Tj)+D*dt/(dx)^2(Ca(Xi+1,Tj)-2Ca(Xi,Tj)+Ca(Xi-1,Tj)) - dt(-k1Ca(Xi,Tj)*Cb(Xi,Tj)+k2*Cc(Xi,Tj)), i and j are intervals. As you see, we have B and C and they are not constants.

Thanks for your help.