3 views (last 30 days)

Show older comments

I have the following code for the reaction interface between acid and base (A-acid, B-base,S-salt):

dS=k1*cA(i+1,:).*cB(i+1,:)*dt;

IDX = (cA(i+1,:)<=NC | cB(i+1,:)<=NC); % Logical indexing

%if in a column for this time iteration either cA or cB is less than NC

%then set the value of dS equal to zero (i.e. no reaction occurred at

%that location

dS(IDX)=0;

where cA and cB are arrays, dS is a vector with the same number of columns as cA and cB and dt is the time increment.

Now I want to analyze each element of the current row (i+1, the time iteration) in cA and cB based on what the corresponding value is in the column (space increment) of the vector dS. If dS equals zero, I want to set cA and cB to one set of values (no reaction has occurred). If dS is greater than zero, then reaction occurred and a different line of code is used to set the value of the concentration in that cell in either cA or cB.

If dS > 0

cA(i+1,:)=cA(i+1,:)-dS;

cB(i+1,:)=cB(i+1,:)-dS;

if dS = 0, seems to be more tricky. In this case the acid and base concentrations are set by the water dissociation constant, Kw. The code needs to determine which array (cA or cB) has the higher value for the concentration. Then the concentration of the one with the lower value is calculated based on Kw and the higher concentration. So if cA = .02 and cB = .0000001 in equivalent cells in cA and cB (equivalent based on row/column), then cA = 0.02 and cB = kW/cA.

Coming from a visual basic (VB) background, I have a hard time understanding how to handle many cells at once as MatLab can do. Clearly, it's more efficient and I don't want to use VB logic to code this issue.

Joel Lynch
on 20 Jun 2021

Edited: Joel Lynch
on 20 Jun 2021

The trick is to recognize that you can use & / | (and/or) to chain multiple logical operations together. In this case, you can get logical inidices that match both dS=0 and cA>cB.

% indicies of columns where dS==0 AND cA>cB

idx = dS==0 & cA(i+1,:)>cB(i+1,:);

cB(i+1,idx) = kW/cA(i+1,idx);

% inidices of colums where dS==0 AND cB>cA

idx = dS==0 & cA(i+1,:)<cB(i+1,:);

cA(i+1,idx) = kW/cB(i+1,idx);

The various operators used in matlab are summarized here, and can be useful when learning the language.

Joel Lynch
on 20 Jun 2021

ODE solvers should be included in the base software. MATLAB has some good documentation on how to choose and implement an ode scheme. I would be vary careful using euler integration, unless you know the system is non-stiff.

Questions on chemical kinetics are probably outside the scope of this forum, but a standard way to test stiff ODE's is to continually reduce your time-step size. If your value converges to a constant result, stiffness probably isn't the most immediate problem.

Good luck!

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!