[tomdylan]

Hi foamers,

I'd like to expand the capabilities of a solver for coupled flow-deformation in porous media in order to consider unsaturated flow conditions. The problem consist of a diffusion-type equation for pore pressure p and an momentum equation for the displacement U that are handled in a segregated manner (see the code sniplet).

Under unsaturated conditions the coefficients Dp and Dp2 depend on the pressure, so the whole system gets a bit more non-linear. The "clean" way would be to evaluate p in a inner iteration loop (Picard iteration?) and then after convergence to pass over grad(p) to the momentum equation.

Since pEqn receives a kind of "lagged" contribution fvc::div(fvc::ddt(Dp2,U)) from U.Eqn and exibits therefore some linearity which is handled in the iteration loop I wonder if just updating the Dp and Dp2 coefficients within the existing iteration loop would be a good idea. For this I need the values of the underrelaxed p-field. Is there a way to access the under-relaxed p-field values after p.relax() in order to evaluate the p-dependent coefficients according to the corresponding relaxed p-values?

do
        {
            p.storePrevIter();

        fvScalarMatrix pEqn
        (
            fvm::ddt(p) == fvm::laplacian(Dp, p) - fvc::div(fvc::ddt(Dp2,U))
        );

        pResidual = pEqn.solve().initialResidual();
        p.relax(); 
    
            U.storePrevIter();
        
            fvVectorMatrix UEqn
            (
                fvm::d2dt2(rho, U)
             ==
                fvm::laplacian(K, U, "laplacian(DU,U)")

              + fvc::div(sigma)

              - fvc::div(K & gradU)

              - fvc::grad(p)
            );

            UResidual = UEqn.solve().initialResidual();
            U.relax();
            gradU = fvc::grad(U);
 
    #       include "calculateSigmaEpsilonOrtho.H"  
         
            initialResidual = max(pResidual,UResidual);
            if (iCorr % 1 == 0){ Info << "\tTime " << runTime.value()
             << ", Corrector " << iCorr
             << ", Solving for " << U.name()
            // << " using " << solverPerf.solverName()
             << ", residual = " << initialResidual << endl;   
              }   
        } while (initialResidual > convergenceTolerance && ++iCorr < nCorr);