How to code an instantaneous reaction in Matlab

I have two slabs adjacent to one another. One contains acid and one contains base. The reaction at the interface in instantaneous and irreversible, which means that the concentration of both acid and base at the interface is zero. Acid tries to diffuse into the base and vice versa. How would I code the "reaction" into Matlab. So at the interface node for either acid or base --> in - out = consumed. I don't know how to code the amount counsumed at any time.

Answer 1

[William Rose]

@Robert Demyanovich,

Here is a simulation. Figure 1 shows the concentration profiles at specific time points. The simulation should conserve mass. Figure 2 shows that it does. Figure 2 also shows that the base is fully consumed by about t=0.55.

The time step is dt=1e-5. This value was chosen because negative concentrations occurred when dt=1e-4, which indicates dt=1e-4 is not small enough.

For each time step

• estimate A & B diffusion at all bins with a discretized version of the diffusion equation
• react any A + B => S that are in the same bin, at all bins
• update A, B, S concentrations to reflect reactions at all bins

Repeat

Salt does not diffuse from its site of formation.

[Robert Demyanovich]

Hello William: Thanks for taking the time to consider my problem and setup simulations. Could you provide me with the Matlab code so that I can run the simulations? The basic idea looks good, but I'm not sure about the actual profiles.

[William Rose]

@Robert Demyanovich

Code attached. A zip file with an animation (.mp4) is also attached. The movie is fun to watch. It was created in the script, as an AVI file. I converted it to mp4 to make it smaller. Then I zipped it because you cannot post .mp4 files on this site, but you can post zip files. If you publish this work, I would appreciate acknowlegement (or co-authorship - I see your recent https://doi.org/10.1063/5.0040336), because the algorithm I created to solve this problem numerically is non-trivial and non-obvious. Thank you. William Rose, [e-mail]

Answer 2

[William Rose]

@Robert Demyanovich,

Do you expect to end up with a partial differential equation, as expected for diffusion? Of course the 1D equation is

where u=concentration and α=diffusion constant.

Or do you expect to get an ordinary differential equation, which can result if you can make enough simplifying assumptions about the geometry?

Are the acid slab and the base slab porous solids in which acid and base diffuse, or are the slabs well-stirred liquids, or something else? As the product forms at the interface, does it accumulate to make a growing third slab in between the other two? In that case, the acid and base on each side must diffuse through the growing product slab in order to reach one another. And then we must ask whether reactions can occur within the product slab, or only on its outer edges.

A simple version of this problem is to assume u=0 at x=0, at all times, and that the acid and base sides are independent and have the same diffusion constant and initial value, and the concentrations are initally uniform through each slab. These are boundary conditions. Then it is like a 1-D heat problem where one end of the bar is at T=0.

Create a 2D array u(i,j) whose elements are the acid (or base; they're the same) concentrations at . In other words, row 1 elements are the concentrations at t=0, for various positions. Column 1 elements are the concentrations at x=0, for various times. Column 1 will be zero for the whole simulation, because the reactant reaching x=0 is immediately consumed. Row 1 is a row of constants (except for a 0 in column 1), beacuse we assume a uniform concentration in the slab at t=0. The rest of the rows are unknown. Compute row 2 by applying the discretized diffusion equation to row 1, and think about how you want to treat the edges. Repeat for successive rows, until done, i.e. until completing the final time step. Will it reach a steady state? Do the simulation to find out.

[Star Strider]

Exactly my concerns, however I didn’t have the patience to enunciate all of them.

+1

[Robert Demyanovich]

Attached is an image of the diffusion. These are dilute solutions so no need to worry about the salt product impacting the diffusion.

Acid starts out in the left hand slab and base in the right hand slab. On the diagram is the neutralization time. Acid is always in excess (i.e. ). The acid and base diffuse towards each other, but since there is more acid than base, the acid winds up penetrating into the base slab. At infinite time the excess acid has equilibrated across both halves of the slab.

I have an analytical solution, but I want to check it using Matlab. Since acid and base cannot coexist, we setup the analytical solution by defining a concentration where:

where A is for acid and B for base. This allowed for use of separation of variables in a single region as opposed to two regions and the moving boundary within the total slab did not even factor into the derivation of the final solution. The analytical solution looks pretty good, but I wanted to check some results that are not obvious; hence the desire to use Matlab. I should say I have never used Matlab before and am only learning about it.

The diffusion equation is:

with the following boundary conditions:

at

or in terms of the defined variable c,

at

or in terms of the defined variable c,

The initial conditions are:

for

for

Diffusion constants for hydrogen and hydroxyl ion are assumed to be equal.

It is a moving boundary problem. Initially, the reaction plane starts at x = 0, but then moves towards . Unlike the analytical derivation, I'm not sure that I can setup Matlab code to ignore the moving boundary - yet have the Matlab output actually show the boundary moving (the analytical expression does show the boundary moving). What I have right now, doesn't really show any change within the half-slabs, probably because there is no term that models the consumption of acid and base. In other words, the initial quantities are not getting any lower as a result of the reaction.