fsi.stress().evolve()

[xye111]

bool unsTotalLagrangianStress::evolve()
{
    Info << "Evolving stress model" << endl;

    int nCorr
    (
        readInt(stressProperties().lookup("nCorrectors"))
    );

    scalar convergenceTolerance
    (
        readScalar(stressProperties().lookup("convergenceTolerance"))
    );

    scalar relConvergenceTolerance = 0;
    if (stressProperties().found("relConvergenceTolerance"))
    {
        relConvergenceTolerance =
            readScalar(stressProperties().lookup("relConvergenceTolerance"));
    }

    Switch nonLinear(stressProperties().lookup("nonLinear"));

    Switch debug(stressProperties().lookup("debug"));

    dimensionedScalar K("K", dimless/dimTime, 0);
    if (stressProperties().found("K"))
    {
        K = dimensionedScalar(stressProperties().lookup("K"));
    }
        
    dimensionedVector g("g", dimVelocity/dimTime, vector::zero);
    if (stressProperties().found("g"))
    {
        g = dimensionedVector(stressProperties().lookup("g"));
    }

    bool enforceLinear = false;
    stressProperties().set("enforceLinear", enforceLinear);

    if (thermalStress())
    {
        threeK();
        alpha();
        threeKf();
        alphaf();
    }

    int iCorr = 0;
    scalar initialResidual = 0;
    lduMatrix::solverPerformance solverPerf;
    scalar res = 1;
    scalar maxRes = 0;
    scalar curConvergenceTolerance = convergenceTolerance;

    lduMatrix::debug = debug;


    do
    {
        if (lduMatrix::debug)
        {
            Info<< "Time: " << runTime().timeName() 
                << ", outer iteration: " << iCorr << endl;
        }

        D_.storePrevIter();

        fvVectorMatrix DEqn
        (
            rho_*fvm::d2dt2(D_)
         == fvm::laplacian(2*muf_ + lambdaf_, D_, "laplacian(DD,D)")
          + fvc::div
            (
                mesh().Sf()
              & (
                  - (muf_ + lambdaf_)*gradDf_
                  + muf_*gradDf_.T() + lambdaf_*(I*tr(gradDf_))
                )
            )
          + rho_*g
        );

        // Add damping 
        if (K.value() > SMALL)
        {
            DEqn += K*rho_*fvm::ddt(D_);
        }

        if (nonLinear && !enforceLinear)
        {
            surfaceSymmTensorField Ef = 
                symm(gradDf_) + 0.5*symm(gradDf_ & gradDf_.T());

            surfaceSymmTensorField sigmaf = 2*muf_*Ef  + I*(lambdaf_*tr(Ef));

            DEqn -= 
                fvc::div
                (
                    muf_*(mesh().Sf() & (gradDf_ & gradDf_.T()))
                  + 0.5*lambdaf_*tr(gradDf_ & gradDf_.T())*mesh().Sf()
                )
              + fvc::div(mesh().Sf() & (sigmaf & gradDf_));
        }

        if (thermalStress())
        {
            DEqn += 
                fvc::div
                (
                    mesh().Sf()*threeKf()*alphaf()*DTf()
                );
        }

        if (interface().valid())
        {
            interface()->correct(DEqn);
        }

        solverPerf = DEqn.solve();

        D_.relax();

        if(iCorr == 0)
        {
            initialResidual = solverPerf.initialResidual();
        }

        if (interface().valid())
        {
            interface()->updateDisplacement(pointD_);            
            interface()->updateDisplacementGradient(gradD_, gradDf_);
        }
        else
        {
            volToPoint_.interpolate(D_, pointD_);
            gradD_ = fvc::grad(D_, pointD_);
            gradDf_ = fvc::fGrad(D_, pointD_);
        }

        if (nonLinear && !enforceLinear)
        {
            surfaceScalarField Det = det(I+gradDf_);

            scalar minDetFf = min(Det).value();
            reduce(minDetFf, minOp<scalar>());

            scalar maxDetFf = max(Det).value();
            reduce(maxDetFf, maxOp<scalar>());

            if ( (minDetFf<0.01) || (maxDetFf>100) )
            {
                Info << minDetFf << ", " << maxDetFf << endl;

                enforceLinear = true;
                stressProperties().set("enforceLinear", enforceLinear);
            }
        }

        // Calculate relative momentum residual
        res = residual();

        if (res > maxRes)
        {
            maxRes = res;
        }

        curConvergenceTolerance = maxRes*relConvergenceTolerance;
        if (curConvergenceTolerance < convergenceTolerance)
        {
            curConvergenceTolerance = convergenceTolerance;
        }

        if (lduMatrix::debug)
        {
            Info << "Relative residual = " << res << endl;
        }
    }
    while
    (
        res > curConvergenceTolerance 
     && ++iCorr < nCorr
    );

    U_ = fvc::ddt(D_);

    // Calculate second Piola-Kirchhoff stress
    {
        epsilon_ = symm(gradD_);

        if(nonLinear)
        {
            epsilon_ += 0.5*symm(gradD_ & gradD_.T());
        }

        sigma_ = 2*mu_*epsilon_ + I*(lambda_*tr(epsilon_));

        if (thermalStress())
        {
            sigma_ -= symmTensor(1,0,0,1,0,1)*threeK()*alpha()*DT();
        }
    }

    Info << solverPerf.solverName() << ": Solving for " << D_.name()
        << ", Initial residula = " << initialResidual
        << ", Final residual = " << solverPerf.initialResidual()
        << ", No outer iterations = " << iCorr
        << "\nMax relative residual = " << maxRes 
        << ", Relative residual = " << res 
        << ", enforceLinear = " << enforceLinear << endl;

    lduMatrix::debug = 1;

    if (nonLinear && enforceLinear)
    {
        return false;

    }

    return true;
}


void unsTotalLagrangianStress::predict()
{
    Info << "Predicting stress model" << endl;

    D_ = D_ + U_*runTime().deltaT();

    if (interface().valid())
    {
        interface()->updateDisplacement(pointD_);
        interface()->updateDisplacementGradient(gradD_, gradDf_);
    }
    else
    {
        volToPoint_.interpolate(D_, pointD_);
        gradD_ = fvc::grad(D_, pointD_);
        gradDf_ = fvc::fGrad(D_, pointD_);
    }

    D_.correctBoundaryConditions();
}

这是fsi对象下全拉格朗日固体模型下的evolve函数，对应主程序里的fsi.stress().evolve()，这个的流程图需要画出来吗