PhreeqcRM gas phase pressure

[carlbirger]

Hello!

I am trying to use PhreqcRM in a two phase (liquid/gas) CFD code. My problem is that the gas phase pressure stays fixed to that specified in the input file. It is not affected by RM_SetPressure(id, pressure).

How can I update the gas phase pressure with the pressure from the CFD code?

In a desperate attempt, I tried to delete the gas phase and create a new one. The gas phase is deleted, but newly created gas phase is not taken into account. What am I doing wrong here?

Please have a look at the input file and fortran lines below. This is for a water / CO2 system

Any help or advice will be greatly appreciated.

Carl Birger

input file:

TITLE test input.
SOLUTION 1  Inlet 1
        units               mg/kgw
        temp               20                              # C
        pressure         1.                               # atm
        pH                  7.00
        water              1.0                              # kg
END
GAS_PHASE 1 Fixed pressure gas phase
        -fixed_pressure
        -pressure       1.              # atm
        -temperature    25.         # C
        -volume         1.              # liter
        CO2(g)          1.               # partial pressure (atm)
END
SELECTED_OUTPUT
        -file                       test.txt
        -selected_out       true
        -reset                   false
        -ph                        true
        -saturation_indices CO2(g)
        -gases                  CO2(g)         # amount in the gas phase (moles)
        -totals                  C                  # mol/kgw
USER_PUNCH
  -heading Viscosity Compressibility
  10 PUNCH VISCOS, KAPPA
END

fortran lines:

    input = ""
    input = trim(input) // "DELETE; -gas_phase" // new_line(ch)
    input = trim(input) // "GAS_PHASE 1" // new_line(ch)
    input = trim(input) // "  -fixed_pressure "// new_line(ch)
    input = trim(input) // "  -pressure       1. "// new_line(ch)
    input = trim(input) // "  -temperature    25. "// new_line(ch)
    input = trim(input) // "  -volume         1.e0  "// new_line(ch)
    input = trim(input) // "  CO2(g)          1. "// new_line(ch)
    input = trim(input) // "END"// new_line(ch)
    status = RM_RunString(id, 1, 1, 1, input)

   ! Transfer data to PhreeqcRM for reactions
   status = RM_UseSolutionDensityVolume(id, 0)     ! 1 = false
   status = RM_SetTemperature(id, temperature)     ! If temperature changes
   status = RM_SetPressure(id, pressure)           ! If pressure changes
   status = RM_SetConcentrations(id, c)            ! Transported concentrations
   status = RM_SetTimeStep(id, time_step)          ! Time step for kinetic reactions
   time = time + time_step
   status = RM_SetTime(id, time)                   ! Current time

   ! Run cells with transported conditions
   status = RM_RunCells(id)

[dlparkhurst]

If you assume a fixed_volume gas phase, the pressure changes to achieve equilibrium between the solution and the gas phase. This would apply to a gas pocket in a cell.

I can see the problem with your case. When a fixed_pressure gas phase is defined, that is the pressure remains constant for the whole simulation. This makes sense for a system with a fixed pressure distribution.  For PHREEQC, it is possible to modify the pressure with a REACTION_PRESSURE n data block, or a GAS_PHASE_MODIFY n data block.

You are transporting the gas phase, so, I think you would need to change the GAS_PHASE composition as well as the solution composition. Unfortunately, PhreeqcRM does not have an easy way to do this. In the simplest case of no parallelization, there are three IPhreeqc instances, a worker, an initialization, and a utility instance. To change only the pressure, you could define GAS_PHASE_MODIFY blocks for each cell and run them through the worker instance with RM_RunString. I think you could alternatively use REACTION_PRESSURE data blocks. I also think you may need to use GAS_PHASE instead to define the new gas composition from your transport calculations. Note that the cell number is for chemistry cells, which may be different from the numbering system in your model. The mapping can be found with RM_GetBackwardMapping.

For OpenMP parallelization, there will be multiple workers with groups of chemistry cells in sequence. The methods RM_GetStartCell and GetEndCell provide the beginning and ending cells in each worker (if load balancing is used, the distribution of cells may change), and use RM_GetIPhreeqcid to apply the definitions to the appropriate worker. For simplicity, I think you could define all GAS_MODIFY blocks for all workers and use RM_RunString for the workeers without causing a problem.

For MPI, you may need to run all of the GAS_PHASE_MODIFY (or REACTION_PRESSURE or GAS_PHASE) for all of the processes with RM_RunString for workers. The definitions for cells not defined in a process should be ignored.

[dlparkhurst]

I'm going to go so far as to say this may be, could be, sort of a bug. I guess I can't think of a case where one would want the solution to be at one pressure and a GAS_PHASE at a different fixed pressure in a transport environment. To be honest, I wasn't thinking too much about gas transport in the development of PhreeqcRM.

I am still not sure how you will deal with the gas phase composition that you transport. Seems like you need to redefine GAS_PHASE of each cell with the transported gases, or if you put all of the mass in solution (so that PHREEQC creates the equilibrium gas phase at the specified pressure), you would still need to redefine the gas phases to have no moles at the start of the reaction calculation. It probably means that there should be a SetGasConcentrations, but that is more work.

If you do want only to set the GAS_PHASE pressures at the same time that the solution pressures are set with RM_SetPressure, you can use the attached code. Replace the method SetPressure in PhreeqcRM.cpp at about line 9586 with the attached code; then recompile the library. The pressure of the GAS_PHASEs will have the pressures specified in the array of pressures used as the argument to RM_SetPressure.