Guidance on Derivatives for Concentration-Based Density Equation

[burritotukituki]

I'm currently in the process of developing my own PorousFlow material and aiming to incorporate a concentration-based density equation for water. Initially, I'm considering a straightforward relation:

ρ=1000+0.001×C

Where 1000 represents the density of pure water, and C denotes concentration. However, I'm encountering difficulties in accurately determining the derivatives, I am using as a base PorousFlowSingleComponentFluid.C, and looks like there the derivatives are computed by an external function... Could someone provide guidance on how the derivatives should be computed?

to the code

Thanks,

[1runer]

Hi @burritotukituki
PorousFlowSingleComponentFluid doesn't propagate the derivatives towards the kernel.

You should have a look into PorousFlowMulticomponentFluid
and into PorousFlowBrine

Those objects call FluidProperties including a concentration variable. You can follow their structure and the structure in BrineFluidProperties to create your own FluidProperties object.

[burritotukituki]

wops... looks more complicated than what I had in mind. Okay, thanks. I will take a look and come back to you.

[burritotukituki]

Hi @1runer,

I've been delving into how these fluid properties function. From what I've gathered, it seems I need to create a C++ object similar to BrineFluidProperties.C. Within this object, there's a function called rho_from_p_T_X, which appears to be what I need...

Now, I can take this function as a starting point to develop my own. Water density is determined by rho_from_p_T, while the brine density is calculated by scaling the water density with the ratio of brine molar mass to pure water molar mass.

So, by modifying the brine density function and tailoring it to my requirements, it should work, right?

[1runer]

Yes. This should work.

But take care to also calculate the derivatives w.r.t. to the mass fraction.
You can do this either by your self or follow the approach in BrineFluidProperties using the AD formulation.

[burritotukituki]

Okay great, now can you please guide me and tell me where the derivatives are computed in BrineFluidProperties? Is it here?

  FPDualReal water_density;
  if (_water_fp_derivs)
  {
    Real rho, drho_dp, drho_dT;
    _water_fp->rho_from_p_T(pressure.value(), Tv.value() + _T_c2k, rho, drho_dp, drho_dT);
    water_density = rho;

    water_density.derivatives() = pressure.derivatives() * drho_dp + Tv.derivatives() * drho_dT;
  }
  else
    water_density = _water_fp->rho_from_p_T(pressure.value(), Tv.value() + _T_c2k);

  // The brine density is given by the water density scaled by the ratio of
  // brine molar mass to pure water molar mass
  return water_density * molarMass(xnacl) / _Mh2o;

I am sorry I am bit lost on this one...
Thanks again

[1runer]

Please see this function in BrineFluidProperties

void
BrineFluidProperties::rho_from_p_T_X(Real pressure,
                                     Real temperature,
                                     Real xnacl,
                                     Real & rho,
                                     Real & drho_dp,
                                     Real & drho_dT,
                                     Real & drho_dx) const
{
  // Initialise the AD value and derivatives
  FPDualReal p = pressure;
  Moose::derivInsert(p.derivatives(), 0, 1.0);
  FPDualReal T = temperature;
  Moose::derivInsert(T.derivatives(), 1, 1.0);
  FPDualReal x = xnacl;
  Moose::derivInsert(x.derivatives(), 2, 1.0);

  _water_fp_derivs = true;
  FPDualReal ad_rho = this->rho_from_p_T_X(p, T, x);

  rho = ad_rho.value();
  drho_dp = ad_rho.derivatives()[0];
  drho_dT = ad_rho.derivatives()[1];
  drho_dx = ad_rho.derivatives()[2];
}

The derivatives are calculated in the AD FluidProperties function:
FPDualReal
BrineFluidProperties::rho_from_p_T_X

and in the "void function" extracted w.r.t. pressure, temperature and mass_fraction to fill the material properties the PorousFlow module needs.