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 个赞
小变形的话，用 unsTotalLagrangianStress 就可以了。
可是在全拉格朗日程序下，没有 updateTotalField() 函数