Carbon distribution in an open system iterative solve block.xmcd

[Joeboy]

Here I have used Mathcad to solve a set of equiibrium chemistry equations which compute the various carbonate species as function of the partial pressure of carbon dioxide.  Note the vast range of the solution values, approximately 10 orders of magnitude.  I find this unique to chemistry problems.

[ValeryOchkov]

OK!

But constants is not constants - functions of temperature:

[ValeryOchkov]

We can think about

Ca2+ + CO32- = CaCO3

in real water systems.

[ValeryOchkov]

I use for the solution of this task not the Given-Find but the programming:

[ValeryOchkov]

And the last point!

It will be good to see UNITS in this calculation: Pco2 := 10^-2 atm etc.

I have seen a lot of same Mathcad task but all was without units!

[Joeboy]

I originally wrote this file  some time in the late 1990s as a teaching tool for Civil/Environmental engineers.  I suspect I have 20 versions of it.

The dissociation constants are functions of temperature but for a given problem or design a specific temperature is normally assumed

Look at your programming solution ....all equations (jumbled together) very little text....your solution essentially assumes the reader already knows how to use Mathcad and Mathacd programming and how to solve the problem.  I am generally dealing with engineers not water chemists or programmers, this solution would be essentially useless to them.  No client is going to look at this.  My students expected and my clients pay for solutions with a degree of understandability (is that a word ?)  I started using Mathcad 1..5 around 1990 and one of the reasons I did (and still do) is because of the free form format which gives me the ability  to place text when and where its needed. I find the solve block algorithm far easier to document and explain than code, keep in mind very few engineering students and none of my clients actually write code.  A useful solution always contains far more text than calculations and MUST be understandable to the USER, not simply  a demonstration of how many equations can be crammmed into how little space................

I have several versions of this file with units....somewhere.  I did not use units because equilibrium constants as well as solubility product values are rarely if ever given with units in any text or reference book, simply not done in chemistry.  While I use units the vast majority of the time I feel that, in this case it simply confuses the issue.  Although I did not do it here I generally include text explaining why units were not used.

P.S.  If one of my students had turned in your solution I would have returned it with a note saying I would not grade it until I had complete documentation .  The same is true in industry, if an enginer submiited this to his/her supervisor for approval as part of a project I can tell you that approval will not be given until both the supervisor and the client can understand what is on the page  (I have any number of former students that will attest to that). If this is the best one can do with Mathcad they may as well be using Matlab or Mathematica.

[AlanStevens]

I found this to be an interesting problem Dixie, and fully agree with your comments on documentation (I come from an industrial engineering background).

One potential problem with your formulation of the solution method though, is that the results are potentially very sensitive to exactly how the equations are written.  If, for example, your students decided to write the ion product of water as OH = Kw/H instead of H = Kw/OH (as they well might, given that OH occurs on its own in the charge balance equation) then the results they got would be nonsense as a result of numerical errors.  (Hence the importance of checking the results against the original equations of course). 

In fact by using a few simple substitutions your equations can be reduced to a cubic in H (as a function of PCO2) which is easily solved for (positive) H by a Given ... Find solve block.  The other ion concentrations can then be derived in a straightforward manner.  The results are also slightly more accurate (you don't get the carbonic acid inaccuracies).

However, I guess many of your engineers would be uncomfortable with this approach too!

Alan

[Joeboy]

You hit the nail on the head.  I find that when solving chemistry problems involving more than 2 or 3 equations, writing them in "ratio" form often improves the odds of convergence.  I also may increase the number of decimal places before Mathcad sees a number as "zero".  Even so the governing equations are not always satisfied under all conditions.  I find chemistry problems somewhat unique in this respect in that it is not uncommon to have a solution vector spanning 10-15 orders of magnitude.  If you examine chemistry texts you will see that such problem solutions are (were) approximated using graphical approaches because chemists often lacked the math tools to solve them but theyhad enough chemistry knowledge to make intelligent approximations.  My students often lacked sufficient chemistry knowledge but they did understand that a system of say 5 independent equations will result in a 5 element solution vector, if they can be solved.

I (and my students) are certainly aware that when solving small numbers of equations it may just as easy to combine them into a single equation.  I sometimes have them do this; however more realistic problems involves larger numbers of equations and such problems are often part of a numerical methods course where I am also trying to demonstrate to students actual applications of various methods of solving linear and nonlinear equations.  Keep in mind that a single applied problem can be used to illustrate a variety of issues, a problem involving chemistry is not JUST a chemistry problem

Intersting bit of knowledge:  One of the first applications of numerical techniques was solving the equations governing the reactions leading to the explosion of the first atomic bomb. 

[ValeryOchkov]

AlanStevens wrote:

...
  If, for example, your students decided to write the ion product of water as OH = Kw/H instead of H = Kw/OH (as they well might, given that OH occurs on its own in the charge balance equation) then the results they got would be nonsense as a result of numerical errors.

Alan

Sorry but real equation is a H+ * a OH- = Kw where a is an activity of ions - a function of temperature, ionic force etc

This system is not for the Given-Find Block but for programming by the bisection method for example...

[ValeryOchkov]

AlanStevens wrote:

...

In fact by using a few simple substitutions your equations can be reduced to a cubic in H (as a function of PCO2) which is easily solved for (positive) H by a Given ... Find solve block.  The other ion concentrations can then be derived in a straightforward manner.  The results are also slightly more accurate (you don't get the carbonic acid inaccuracies).
...

Alan

[Joeboy]

Nice use of Mathcad's symbolic solving capabilities.

[Joeboy]

Come on !  Any chemistry book points out that equilibruim coefficients, ion products, and solubility products are functions of temperature and ionic strength but  then hastens to point out that for dilute solutions the activity is equal to the molar concentration, that is, the activity coefficient is essentially 1.0.  The only time I have found it necessary to incorporate activity is with brine solutions from oil wells

[ValeryOchkov]

But we must think not about water - about sea water with ionic activity <<1

[Joeboy]

I like the use of the symbolic solve for the carbonate ion, however that solution does not account for the charge balance equation nor the ion product of water which, I believe, would have to be solved simultaneously with the other equations,  thus you would need to reduce 5 equations to 1.  The carbonate inaccuracies are minimal I believe, about a 10% variation in concentrations that run around 10-7 moles per liter.  You can't measure carbonate ion concentrations that accrately using standard techniques.