updateTotalFields

void pRveUnsIncrTotalLagrangianStress::updateTotalFields()
{
    unsIncrTotalLagrangianStress::updateTotalFields();

    avgEpsilon_ += avgDEpsilon_;

    totPointD_.internalField() +=
        pointDD().internalField()
      + (avgDEpsilon_ & mesh().points());

    totD_ += DD() + (avgDEpsilon_ & mesh().C());
  volSymmTensorField avgDSigma = 
        2*mu()*avgDEpsilon_ + lambda()*tr(avgDEpsilon_)*I;

    totSigma_ += DSigma() + avgDSigma;

    surfaceSymmTensorField avgDSigmaf = 
        2*muf()*avgDEpsilon_ + lambdaf()*tr(avgDEpsilon_)*I;

    totSigmaf_ += DSigmaf() + avgDSigmaf;


这个是更新固体所有受到的力的函数,这里的具体语句代表是么呢?

主程序在子迭代结束以后,出现了固体模型下的这个函数,fsi.stress().updateTotalFields();我在固体模型文件夹 pRveUnsIncrTotalLagrangianStress 中将 updateTotalField() 找到,为了画出此语句的流程图必须看懂这个程序,但是里面有些变量我不知道代表什么

哪些变量?都做了什么操作?

void unsIncrTotalLagrangianStress::updateTotalFields()
{
    Info << "Update total fields" << endl;

    D_ += DD_;
    pointD_ += pointDD_;
    gradD_ += gradDD_;
    gradDf_ += gradDDf_;
    sigma_ += DSigma_;
    sigmaf_ += DSigmaf_;
    epsilon_ += DEpsilon_;
    epsilonf_ += DEpsilonf_;

    U_ = fvc::ddt(D_);

    if (rheology_.plasticityActive())
    {
        epsilonP_ += rheology_.DEpsilonP();
//         epsilonPf_ += rheology_.DEpsilonPf();
        rheology_.updateYieldStress();
    }
}

这个文件夹下面也有这个updatetoalfields?

bool unsIncrTotalLagrangianStress::evolve()
{
    Info << "Evolving stress model unsIncrTotalLagrangianStress" << 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"));

    componentReferenceList cr
    (
        stressProperties().lookup("componentReference"),
        componentReference::iNew(mesh())
    );

    dimensionedScalar K("K", dimless/dimTime, 0);
    if (stressProperties().found("K"))
    {
        K = dimensionedScalar(stressProperties().lookup("K"));
    }

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

    lduMatrix::debug = debug;

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

    // Activate old time fields
    fvc::ddt(D_);

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

        DD_.storePrevIter();

        fvVectorMatrix DDEqn
        (
            rho_*fvm::d2dt2(DD_)
          - fvm::laplacian(2*muf_ + lambdaf_, DD_, "laplacian(DDD,DD)")
         == fvc::div
            (
                mesh().Sf()
              & (
                  - (muf_ + lambdaf_)*gradDDf_
                  + muf_*gradDDf_.T() + lambdaf_*(I*tr(gradDDf_))
                )
            )
        );

        // Add damping if not zero
        if (K.value() > SMALL)
        {
            DDEqn += K*rho_*fvm::ddt(DD_);
        }

        // Update strain increment
        DEpsilonf_ = symm(gradDDf_);
        if (nonLinear && !enforceLinear)
        {
            DEpsilonf_ += 0.5*symm(gradDDf_ & gradDDf_.T());
            DEpsilonf_ += 0.5*symm(gradDDf_ & gradDf_.T());
            DEpsilonf_ += 0.5*symm(gradDf_ & gradDDf_.T());
        }

        DSigmaf_ = 2*muf_*DEpsilonf_ + I*(lambdaf_*tr(DEpsilonf_));
        
        if (rheology_.plasticityActive())
        {
            DSigmaf_ -= 2*muf_*fvc::interpolate(rheology_.DEpsilonP());
//             DSigmaf_ -= 2*muf_*rheology_.DEpsilonPf();
        }

        if (nonLinear && !enforceLinear)
        {

            DDEqn -=
                fvc::div
                (
                    muf_*(mesh().Sf() & (gradDDf_ & gradDf_.T()))
                  + muf_*(mesh().Sf() & (gradDf_ & gradDDf_.T()))
                  + muf_*(mesh().Sf() & (gradDDf_ & gradDDf_.T()))
                  + 0.5*lambdaf_*tr(gradDDf_ & gradDf_.T())*mesh().Sf()
                  + 0.5*lambdaf_*tr(gradDf_ & gradDDf_.T())*mesh().Sf()
                  + 0.5*lambdaf_*tr(gradDDf_ & gradDDf_.T())*mesh().Sf()
                )
              + fvc::div(mesh().Sf() & (DSigmaf_ & gradDf_))
              + fvc::div(mesh().Sf() & ((sigmaf_ + DSigmaf_) & gradDDf_));
        }

        if (rheology_.plasticityActive())
        {
            DDEqn += 
                fvc::div
                (
                    2*muf_
                   *(
                        mesh().Sf() 
                      & fvc::interpolate(rheology_.DEpsilonP())
                    )
                );

//             DDEqn += 
//                 fvc::div(2*muf_*(mesh().Sf() & rheology_.DEpsilonPf()));
        }

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

        forAll (cr, crI)
        {
            DDEqn.setComponentReference
            (
                cr[crI].patchIndex(),
                cr[crI].faceIndex(),
                cr[crI].dir(),
                cr[crI].value()
            );
        }

        solverPerf = DDEqn.solve();

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

        DD_.relax();

        if (interface().valid())
        {
            interface()->updateDisplacementIncrement(pointDD_);
            interface()->updateDisplacementIncrementGradient
            (
                gradDD_, 
                gradDDf_
            );
        }
        else
        {
            volToPoint_.interpolate(DD_, pointDD_);
            gradDD_ = fvc::grad(DD_, pointDD_);
            gradDDf_ = fvc::fGrad(DD_, pointDD_);
        }

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

            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 momentu residual
        res = residual();

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

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

        if (lduMatrix::debug)
        {
            Info << "Relative residual = " << res << endl;
        }

        // Calculate strain increment
        {
            DEpsilon_ = symm(gradDD_);

            if(nonLinear && !enforceLinear)
            {
                DEpsilon_ += 0.5*symm(gradDD_ & gradDD_.T());
                DEpsilon_ += 0.5*symm(gradDD_ & gradD_.T());
                DEpsilon_ += 0.5*symm(gradD_ & gradDD_.T());
            }
        }

        // Correct plasticity term
        rheology_.correct();
    }
    while
    (
        (res > convergenceTolerance) 
     && (++iCorr < nCorr)
    );

    fvc::ddt(D_);

    // Calculate second Piola-Kirchhoff stress increment
    {
        DSigma_ = 2*mu_*DEpsilon_ + I*(lambda_*tr(DEpsilon_));

        if (rheology_.plasticityActive())
        {
            DSigma_ -= 2*mu_*rheology_.DEpsilonP();
        }
    }

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

    lduMatrix::debug = 1;

    if (nonLinear && enforceLinear)
    {
        return false;
    }

    return true;
}

evolve函数也在这个文件里面出现了,这些应该就是固体方程求解的过程

        fvVectorMatrix DDEqn
        (
            rho_*fvm::d2dt2(DD_)
          - fvm::laplacian(2*muf_ + lambdaf_, DD_, "laplacian(DDD,DD)")
         == fvc::div
            (
                mesh().Sf()
              & (
                  - (muf_ + lambdaf_)*gradDDf_
                  + muf_*gradDDf_.T() + lambdaf_*(I*tr(gradDDf_))
                )
            )
        );

这个就是固体的方程了

固体模型跟updatetotalfields函数大概会涉及到什么样的改动呢

这个函数在每个模型函数下面都必须有。根据模型的不同,函数也有所区别

没有对模型进行修改,我想这里应该不会有修改吧。

这里有两个 updateTotal 函数,我到底是用哪一个呢?这个 unsIncrTotalLagrangianStress 和pRveUnsIncrTotalLagrangianStress 有什么区别呢

    D_ += DD_;
    pointD_ += pointDD_;
    gradD_ += gradDD_;
    gradDf_ += gradDDf_;
    sigma_ += DSigma_;
    sigmaf_ += DSigmaf_;
    epsilon_ += DEpsilon_;
    epsilonf_ += DEpsilonf_;

这些参数都分别代表固体的哪些属性呢

  • unsTotalLagrangianStress: Large strain elastic stress analysis solver based on total Lagrangian displacement formulation 全拉格朗日法
  • unsIncrTotalLagrangianStress: Large strain elastic stress analysis solver based on total Lagrangian displacement increment formulation 全拉格朗日增量法
  • pRveUnsTotalLagrangianStress: 修正拉格朗日法
  • pRveUnsIncrTotalLagrangianStress: 修正拉格朗日增量法

在拉格朗日有限变形的框架下,若以初始构形为参考构形,则称全拉格朗日法 (total Lagrangian, T.L.),若以前一个前一计算结果的构形为参考构形,则称为修正的拉格朗日法 (updated Lagrangian, U.L.)。

1 Like

在流固耦合中,具体是应该使用的哪一种呢

小变形的话,用 unsTotalLagrangianStress 就可以了。

可是在全拉格朗日程序下,没有 updateTotalField() 函数