Issue with the computation of solution volume.

[SaiP]

Hello David,

I am using PhreeqcRM for coupling PHREEQC (Version 3.6) with a CFD solver for 2 phase. After testing for the flow, I can confirm my 2Phase flow solver works well. So, I am using GAS_PHASE_MODIFY and considering a non-isothermal system.

I have an issue with the mass of the solution.

From my understanding PHREEQC always considers 1 kilogram of solution. So, if my solution density is 0.95947 g/cm3 (Sol. at 100C and pressure 24.375 atm) then my solution volume is 1.042 lt to have 1 kg of water. To maintain 1 kg of solution, if solution density falls to 0.9 g/cm3, then solution volume has to proportionally increase to 1/0.9 = 1.11 lt.

Now coming to my case considering coupling of PHREEQC with CFD code. I tried 2 cases. Note that for the pressure and temperature I use, I always have only liquid phase (Pressure is inbetween 25 atm - 15 atm and Temperature is inbetween 150C - 25C).

Case 1:
Lets say, my porosity, representative volume and saturation are all initially equal to 1. And I start the run.
In one of my cell, pressure is 24.75 atm and temperature is 25C. Density of solution is 0.998112 g/cm3 and the solution volume is 1lt. So, this means at this moment we have 0.998112 kg of water.
 
After a few time steps, the temperature in this cell increases to 93.01C and pressure is 24.778 atm. Here, I have solution density of 0.964382 g/cm3 and solution volume of 1.04464 lt. That equals 1.00743 kg of water. So, you can already notice that I gained several grams of water into the system. 

Further, my temperature increases to 139.84C and pressure 24.8106 atm and here the solution density is 0.92746 g/cm3 and solution volume is 1.15968 lt. I have a mass of 1.0755 kg. So, here its clear that I am gaining substantial quantities of solution.

Case 2: In the above case, I restrict the representative volume to be 1 at all times. So, as explained above, if my solution volume is 1.15 lt, my saturation is 1.15 (note my porosity is 1 at all times). Saturation over 1 is unphysical. So, I update the representative volume at all time steps to be equal to the solution volume. So, my saturation is at most 1.
For such a scenario, I start with temperature 25C, pressure 24.75 atm and have solution density 0.998112 kg ans solution volume of 1. So, I have 0.998112 kg of water. Over time, the cell has temperature of 141.71C and pressure 24.75 atm. The solution density is 0.925778 g/cm3 and solution volume is 1.14746 lt. The mass being 1.0623 kg. And at the end, my temperature and pressure doesnt change anymore but the solution volume goes on increasing.
---
From my observation, the solution volume goes on increasing and is not really proportional to the change in the fluid density. This is eventually resulting in gaining additional mass of the solution which is unphysical. Is there something obviously wrong/ missing in my calculation? I tried switching UsingSolutionDensityVolume but I didnt see any impact (the mass keeps increasing).

Thanks!!

[dlparkhurst]

Haven't we been through this before?

I didn't notice what unit you were transporting, but let's say it was mol/L. Assuming representative volume, porosity, and saturation are all 1.0, PhreeqcRM calculates moles in the cell for each element to be mol/L * 1.0 * 1.0 * 1.0. Assuming no other reactions, PhreeqcRM will calculate the new distribution of species and new density and solution volume. If concentrations have changed during the transport step, then the density and the solution volume will change, in particular, as you demonstrate, the volume will probably no longer be 1.0 L.

What needs to happen is that the change in density needs to be incorporated in the transport equations. In my simple minded approach, you would take the PhreeqcRM calculated densities and iterate with the transport equation to arrive at concentrations, that when submitted to PhreeqcRM produce produce a solution volume of 1.0. Multiphase flow is way too complicated and it is the reason I stick to the chemistry.

Although PhreeqcRM concentrations are often reported as mol/kgw, cells do not usually contain exactly 1 kg of water. The mass of water will vary as transport proceeds and molality concentrations are calculated by dividing moles of elements by water mass. The molarity concentration values returned by GetConcentrations are moles divided by solution volume. (There is another option in converting units for GetConcentrations by using volume defined by rv, porosity, saturation, and density defined by the user with SetDensity.)