WilkAndy
commented
10 years ago I have decided against building this at present. However, i will push up changes to RichardsMaterial and associated kernels that: (1) Make RichardsMaterial much nicer to read (2) Prepare RichardsMaterial and the kernels, etc, for easy inclusion of a PSS material, since all the derivatives in RichardsMaterial are now wrt general variables rather than just porepressures, and these are used as such in the kernels (exception noted below).
In the future, if i want to create a PSS material it should be fairly straightforward but lengthy: (1) The exception mentioned above involves RichardsFlux: the SUPG kernel must be generalised in the way that i've noted in comments. This is not a big deal, but not appropriate to do it now. (2) The Seff and Sat UserObjects must be generalised to provide, for instance, pressures as functions of saturations. Again this is not a big deal, but it won't be a small task either partly because various tests and doco will be needed. (3) RichardsMaterial will have to be enhanced to generate the porepressures and saturations, their gradients and derivatives in the same way as the PPPP case does presently. (4) LOTS of tests!
For posterity, I want to describe why i wanted to create a PSS material, and why now I don't. There were two reasons for wanting to create a PSS material: (1) Because it might be advantageous to use different formulations in different problems, to minimise runtimes (2) To get around the problem of singular Jacobian for saturation=0=residual_saturation in >=2 phase problems.
The reason number (2) is why i don't want to create PSS material now. In a PPP formulation, for a 2-phase problem with variables P0 and P1, the mass term has derivative
( dR0/dP0 dR0/dP1 ) = (r0'*S0 + r0*S0' -r0*S0' )
( dR1/dP0 dR1/dP1 ) (r1*S1' r1'*S1 - r1*S1' )Here r0 is the density of "0", r1 the density of "1", S0 the saturation of "0", S1 the saturation of "1", r0' = dr0/dP0, S0' = dS0/dP0, etc. Now consider the case when S0 is at its residual. This means, when using the standard vanGenuchten curve, that S0'=0. Then the top row becomes
(r0'*S0 0)If the residual saturation is zero then this is
(0 0)and the Jacobian is singular. Therefore, we must only use residual>0 in problems where one phase might disappear. However, this no good when porosity is time-varying, since the saturation can physically be pushed below the residual by large positive increases in porosity. Then Darcy-Richards flow is no longer physically valid. This manifests itself numerically by a singular Jacobian for S<S_res, once again. Therefore, for time-varying porosity we either: (1) use residual=0 and get a singular jacobian when one phase disappear; (2) use residual>0 and get physically meaningless situation with singular jacobian when saturation goes less than residual.
I initially thought that this problem went away for a PS formulation. The Jacobian is
( dR0/dP0 dR0/dS1 ) = (r0'*S0 -r0 )
( dR1/dP0 dR1/dS1 ) (-r1'*S1 r1'*(dp/dS)*S1 + r1 )Notice that the determinant is not obviously zero even for S0=0 and S1=1. However, another problem arises: dp/dS=infinity for S0=0 when using the standard van-Genuchten expression. So we are stuck with a crappy Jacobian even in this formulation when S0=0.
So, I've decided not to pursue this PSSS formulation for now. What i've decided to pursue is modifying the van Genuchten expression near S0=0. See issue #3350 .
Build it. This should allow users to use a (pressure, saturation, saturation, saturation, ...) formulation for multiphase fluid flow instead of a (pressure, pressure, pressure, ...) formulation.