Fixes: derivative of krb and bounded saturation equation

[eugenioidaea]

• derivative of krb in krBrooksAndCorey.C was the same as the one defined for kra;
• saturation Sb drops below Sbmin in heterogeneous K fields and code crashes. Testing alternatives for saturation equation in SEqn.H: fvc::div(phib) -> fvm::div(phib, Sb)

[phorgue]

Thank for the fix in krBrooksAndCorey.C

Concerning the modifications of SEqn.H :

• fvm::div(phib, Sb): the phase flux b is directly included in phib, the equation is incorrect if Sb is added.
• I'm not sure if the call to correctboundaryconditions() is necessary, but I'll look into it.

[eugenioidaea]

Thanks for the reply!

You're right about SEqn.H, if Sb is added it becomes redundant.

About saturation dropping below Sbmin & code crashing after hours of simulations, do you have any recommendations to prevent it? I tried lowering Courant to 0.2 but it does not help :(

[phorgue]

This is not exactly the good place to discuss about this.

Usually, this kind of problem occurs more with the van Genuchten model, which exhibits sharp fronts, than with the Brooks and Corey model. I would say that it is rather a problem of configuration (BC or flow properties).