Can Matlab be used to solve a moving boundary problem

I'm trying to determine the neutralization time (tn) of strong acid and strong base. Strong acid, A, is initially in a slab of thickness -2δ and strong base, B, is initially in a slab of thicknss 2δ. Slabs alternate with a slab of A next to one of B next to one of A next to one of B and so on. Because of symmetry, only the half slab thicknesses of each slab of reactant needs to be considered. The slabs only have diffusion of the reactants and when the reactants meet they react instantaneously in a stoichiometry of 1 to 1. So at the reaction plane, the concentration of each reactant is zero. The problem is a moving boundary problem and is illustrated in the attachment.

The boundary conditions at -δ and +δ are no flux conditions. The other boundary is the reaction plane which moves with time. At this plane the concentration of either A or B is zero.

Fick's 2nd law applies. Because the reaction is instantaneous and irreversible, the effect of the reaction is to keep the concentration of each reactant at zero at the reaction plane. Therefore, no reaction term is need nor a convection term; it is just Fick's 2nd law involving only diffusion. I have not been able to come up with an analytical solution so I want to try a numerical method. Could I use MatLab to solve this problem or is the moving boundary aspect of the problem beyond its capabilities?

[Torsten]

Matlab's "pdepe" can solve partial differential equations in one space dimension.

I think the main problem in your case is that you need to solve additionally an ordinary differential equation for the position of the reaction plane. This ordinary differential equation can't be solved simultaneously in "pdepe".

So my guess is that you will need to discretize your PDE in space and solve the resulting system of ordinary differential equations together with the equation for the position of the reaction plane using ODE15S.

But maybe you find something in the internet. Google "Stefan problem & Matlab".

[Robert Demyanovich]

Typically moving boundary problems with constant concentrations at the boundaries have a typical relationship between the reaction plan and time, as follows:

where is the reaction plane which changes with time and β is a constant. So this might simplify things, but on the other hand, β is unknown. I don't know enough about numerical simulations to know if I can run a simulation where β is a guess that is updated with each iteration.

[Torsten]

You can update beta with time, but not with iteration in a PDE solver.

Further, thinking about your problem, it seems that at the position of the reaction plane, you have an inner boundary condition. PDEPE won't be able to handle this.

[Bill Greene]

I'm still confused about how the position of the front is determined. You show an equation in your comment above but say that β is unknown. So how do you determine β? Can you please provide a complete mathematical description of your problem?

In Stefan-type problems that Torsten mentioned, there is an ODE that determines the front position. But, for this class of problem, solution approaches have been developed that do not require explicit tracking of the front.

[Robert Demyanovich]

I've attached an updated (although imperfect) diagram. Because of symmetry we only need to consider diffusion within the slab bounded by . As you can see in the diagram, because there is initially more A (acid) than B (base) on a stoichiometric basis, the reaction plane moves towards . When it reaches , neutralization of B has occurred and the neutralization time on the drawing is indicated by . After infinite time, A equilibrizes over the entire region (). At infinite time, the concentration of A is simply (.

The mass transfer is governed by Fick's law of diffusion:

with the following initial conditions:

at t = 0 and

and at t = 0 and

Now at , A diffuses in the region bounded by and the reaction plane [let's call it ].

At , B diffuses in the opposite direction of A within the region bounded by and .

For diffusion of A one boundary condition at x = is:

For diffusion of B one boundary condition at x = is:

The other boundary condition for both A and B is at the reaction plane.

For A: :

and for B: :

Now in similar diffusion problems (but with different boundary conditions), the following is often used to relate the location of the reaction plane to time:

Although I don't know much about the Stefan problem, I believe that a similar equation results from the melting of ice.

At this point in time, I'm not really sure about the value of the constant, β. I need to do some more research, and it's likely it can be calculated.