0

I'm looking to model the precipitation of a mineral through time at rate, r. As this mineral precipitates, it draws down the concentration of Fe(II) in solution, and this reduces the saturation state of the solution with respect to the mineral, slowing the rate of mineral precipitation.

cla;
clf;

%%
%variables and rate calculation
Fe(1)=0.01*0.26/1000; %mM, * activity coefficient, convert to M
Si=1.8/1000; %mM, convert to M
pH=8;
H=10^(-pH);
k=1.11*10^(-9);
b=1.81;
Ksp=4.68*10^(22);
IAP=(Fe^(3)*Si^(2))/H^(6);
logIAP=log10(IAP);
logKsp=log10(Ksp); 

rate=k*exp(b*(logIAP-logKsp)); %mm Fe Kg-1 hour-1

%% 
%Fe vs time 
dt=0.01;
t=0:dt:10;

for ind=2:length(t)
Fe=Fe(1)-(k*exp(b*(logIAP-logKsp))/1000*t);
end

I eventually want a plot of Fe vs time. The above code doesn't work because Fe depends on logIAP, which depends on Fe from the previous timestep... I need Matlab to iteratively calculate the two for each time step.

Any help appreciated, thanks!

Geochemist1
  • 3
  • 1
  • can you please write down the dependencies of logLAP on FE – Irreducible Oct 12 '17 at 14:03
  • logIAP is just the log of IAP. IAP is the ion activity product, which is a function of pH and Si (both fixed), and Fe, which varies through time. Therefore, the rate of change of Fe depends on Fe, and so IAP, rate and Fe need to recalculated after each time step. – Geochemist1 Oct 12 '17 at 14:26

2 Answers2

0

As the others pointed out you need to make Fe an array. Also you need to index t correctly in your equation as this is a vector.

Fe=zeros(size(t)); % now Fe is an array
Fe(1)=0.01*0.26/1000; % first value assigned
...
Ksp=4.68*10^(22);
IAP=(Fe(1)^(3)*Si^(2))/H^(6);
logIAP=log10(IAP);
logKsp=log10(Ksp); 
....
for ind=2:length(t)

    logIAP=log10((Fe(ind-1)^(3)*Si^(2))/H^(6));
    Fe(ind)=Fe(ind-1)-(k*exp(b*(logIAP-logKsp))/1000*t(ind));
end
Irreducible
  • 864
  • 11
  • 24
-1

Use Fe as an array :

Fe=[0.01*0.26/1000]; % Fe is now an array

....

for ind=2:length(t)
    Fe(ind)=Fe(ind-1)-(k*exp(b*(logIAP-logKsp))/1000*t); % We fill Fe
end

Then plot Fe in function of t

justanidea
  • 1
  • Thanks! I'm getting this error: "In an assignment A(:) = B, the number of elements in A and B must be the same.". Should I be replacing Fe in the equation for IAP With Fe(ind)? And should I be plotting Fe(ind) or Fe? – Geochemist1 Oct 12 '17 at 14:28