is anyone familiar with implementing a sorption isotherm under unsaturated conditions? I got some weird results here (amount of sorbed species increaseing more than amount of solved species decreases, causing increase in total amount) and looking into the code in rfmat_cp.cpp it seems like saturation is not considered in computation of the retardation factor.
I’m not familiar with C++ but how I understood the code, the retardation factor seems to be calculated as
(rfmat_cp.cpp, line 1564 ff) which assumes full saturation by using porosity instead of volumetric water content or did I miss something there? Porosity is said to be a constant value ($POROSITY 1 0.392) and I use a Freundlich isotherm ($ISOTHERM 4 KF nF).
How can I correctly account for sorption under unsaturated conditions?
I think there are no unconfined benchmark examples included but to my understanding you could consider to use a different porosity model.
Now you use $POROSITY 1 0.392, i.e. porosity model = 1. In other words the porosity is constant. Instead you could use e.g. a curve input with case 0. Case 13 look also interesting and it my takes saturation already into account.
Here a list of the possible methods for porosity models:
switch (porosity_model)
{
case 0: // n = f(x)
porosity = GetCurveValue(fct_number, 0, primary_variable[0], &gueltig);
break;
case 1: // n = const
porosity = porosity_model_values[0];
break;
case 2: // n = f(sigma_eff), Stress dependance
porosity = PorosityEffectiveStress(number, primary_variable[0]);
break;
case 3: // n = f(S), Free chemical swelling
porosity
= PorosityVolumetricFreeSwellingConstantIonicstrength(number, primary_variable[0], primary_variable[1]);
break;
case 4: // n = f(S), Constrained chemical swelling
porosity = PorosityEffectiveConstrainedSwellingConstantIonicStrength(number, primary_variable[0],
primary_variable[1], &porosity_sw);
break;
case 5: // n = f(S), Free chemical swelling, I const
porosity = PorosityVolumetricFreeSwelling(number, primary_variable[0], primary_variable[1]);
break;
case 6: // n = f(S), Constrained chemical swelling, I const
porosity
= PorosityEffectiveConstrainedSwelling(number, primary_variable[0], primary_variable[1], &porosity_sw);
break;
case 7: // n = f(mean stress) WW
gval = ele_value_dm[number];
primary_variable[0] = -gval->MeanStress(assem->gp) / 3.0;
porosity = GetCurveValue(porosity_curve, 0, primary_variable[0], &gueltig);
break;
case 10:
/* porosity change through dissolution/precipitation */
porosity = PorosityVolumetricChemicalReaction(number);
break;
case 11: // n = temp const, but spatially distributed CB
// porosity = porosity_model_values[0];
porosity = _mesh->ele_vector[number]->mat_vector(por_index);
break;
case 12: // n = n0 + vol_strain, WX: 03.2011
porosity = PorosityVolStrain(number, porosity_model_values[0], assem);
break;
case 13:
{
CRFProcess* m_pcs_flow = PCSGetFlow();
const int idx_n = m_pcs_flow->GetElementValueIndex("POROSITY");
// porosity change through dissolution/precipitation
// Here, you should access porosity from the element value vector of the flow process
// so you have to get the index of porosity above, if porosity model = 13
porosity = m_pcs_flow->GetElementValue(number, idx_n + 1);
break;
}
I’m not sure if I’m on the right track, but with assignig a curve in case 0 I would give the water content values in place of the porosity to correctly account for the undersaturation?
Assigning it that way would make porosity (resp. water cont.) time dependent, right? But for the transient flow simulation I want to implement I wouldn’t know the temporal course, plus I would need the assigned values to vary in space…
Case 13 looks interesting, how is it assigned? Simply
$POROSITY
13
(it doesn’t seem to call porosity_model_values[0])
or rather
$POROSITY
13 0.392
?
Do you know how saturation would be considered in this case?
