How to implement temperature-dependent pressure calcs in Thermal and Flow Engines

Hey Zoheir,

I am assuming the relationship between your gas temperature and the pore pressure is non-linear? If so, and you want a strong coupling, you'd need to formulate the fully coupled discretization and solve the two unknowns (pore temperature and pore pressure) simultaneously using an iterative approach. We recently did something analogous using a newton-raphson optimization for volume and pressure in partially saturated media.

If the fluid flux is linearly related to the pressure field (and still stoke's flow...big assumption), you could use a weak coupling and it would be relatively trivial to implement. In this case, you would solve pressure and temperature as is already implemented in ThermalEngine. However, you'd then need to use the new temperature to compute mid-step pore pressures with your non-linear relationship between temperature and pressure <- this is the piece that would be a trivial implementation in python or C++. All of this would be followed at the next time-step by using the existing implementation of solving new pressure field to solve the fluid fluxes and pressures (and temperatures) at the next timestep. This weak coupling means you are never truly at equilibrium and you run the risk of losing stability/accuracy on your solution.

If you wanted a strong coupling, I could advise you or implement it for you - but it would require significant development hours, and therefore, I would have little choice but to charge for that kind of software development.

I hope it helps,

Robert

Hey Robert,

Sorry for the late response. I was trying to make a weak coupling between pressure and temperature. I ended up slightly modifying [1] in thermal.cpp to allow changes in dV rates due to temperature changes. I know this function is for considering fluid expansion in the pressure calculations, but I set beta=1 and then divide [2] by the current temperature of the cell (cell->info().temp()). This way I assumed I have an ideal gas with linear relationship between temperature and volume changes. I recompiled and it seams to be working. I need to find a way to validate the results though. What do you think?

Here is the modified function:
void ThermalEngine::computeCellVolumeChangeFromDeltaTemp(CellHandle& cell, Real /*cavDens*/)
{
 Real beta;
 if (tempDependentFluidBeta) beta = 7.5e-6 * cell->info().temp() + 5.7e-5; // linear model for thermal expansion
 else
  beta = fluidBeta;
 Real poreVolume;
 if (porosityFactor > 0) poreVolume = cell->info().volume() * porosityFactor; // allows us to simulate low porosity matrices
 else
  poreVolume = (1. / cell->info().invVoidVolume());
 // else K = flow->fluidBulkModulus;
 if (!cell->info().isCavity) {
  cell->info().dv() += -poreVolume * cell->info().dtemp()/ cell->info().temp() / thermalDT;
 } else if (cell->info().isCavity) {
  cell->info().dv() += -(cell->info().volume()) * beta * cell->info().dtemp()
          / thermalDT; // ignore the particles used for fluid discretization in the pore (i.e. use volume())
 }
}

>>If you wanted a strong coupling, I could advise you or implement it for you - but it would require significant development hours, and therefore, I would have little choice but to charge for that kind of software development.
Thanks for the offer. I shared this with Ryan. We might end up going that route at the end but for preliminary results a weak coupling should do.

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L739
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L751

Hey Zoheir,

Yes, this is why I started by asking if your pressure and temperature are non-linearly related (are they?). If they are linearly related, then the existing weak coupling in ThermalEngine is perfectly adequate.

>> What do you think?

In the case that they are linearly related, I do not understand why you would modify anything. ThermalEngine is already doing what you need, as you found. Beta is the coefficient of thermal expansion which has units of (1/v) * (dv/dT). Dividing the entire expression by temperature as you are doing makes little sense to me. Perhaps you can explain your rationale?

-rc

Hi Robert,

>>Yes, this is why I started by asking if your pressure and temperature are non-linearly related (are they?).
Yes, I am assuming ideal gas so pressure can be linearly related to temperature under constant volume assumption. I should look back at this assumption later though!

>>Perhaps you can explain your rationale?

I make beta=1. Assuming ideal gas law, we have dv=v dT/T. So, by dividing [1] (which is v dT) by T, I try to get dv due to change in T. Am I missing anything?

Thanks
Zoheir