exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Solid models update: clean-up

Browse source
This commit is contained in:
Hrvoje Jasak 2014-02-23 10:54:20 +00:00
parent 3810b29f37
commit bb03a935fa
11 changed files with 224 additions and 211 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateRelativeResidual.H
View file

@ -16,7 +16,6 @@
(
DU.internalField()
- DU.prevIter().internalField()
)
/magDU
)/magDU
);
}

2
applications/solvers/solidMechanics/elasticIncrSolidFoam/createFields.H
View file

@ -24,7 +24,7 @@
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector("zero", dimLength, vector::zero)
dimensionedVector("zero", dimLength, vector::zero)
);
volSymmTensorField DEpsilon

5
applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C
View file

@ -79,7 +79,7 @@ int main(int argc, char *argv[])
# include "calculateDivDSigmaExp.H"
//- linear momentum equation
// Linear momentum equation
fvVectorMatrix DUEqn
(
fvm::d2dt2(rho, DU)
@ -127,7 +127,8 @@ int main(int argc, char *argv[])
&& ++iCorr < nCorr
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
Info<< nl << "Time " << runTime.value()
<< ", Solving for " << DU.name()
<< ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual()
<< ", Relative residual = " << relativeResidual

31
applications/solvers/solidMechanics/elasticIncrSolidFoam/readDivDSigmaExpMethod.H
View file

@ -1,9 +1,22 @@
//- how explicit component of sigma is to be calculated
word divDSigmaExpMethod(mesh.solutionDict().subDict("solidMechanics").lookup("divSigmaExp"));
Info << "Selecting divSigmaExp calculation method " << divDSigmaExpMethod << endl;
if(divDSigmaExpMethod != "standard" && divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose" && divDSigmaExpMethod != "laplacian")
{
FatalError << "divSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
<< "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
<< exit(FatalError);
}
const dictionary& stressControl =
mesh.solutionDict().subDict("solidMechanics");
// Read how explicit component of sigma is to be calculated
word divDSigmaExpMethod(stressControl.lookup("divSigmaExp"));
Info<< "Selecting divSigmaExp calculation method: "
<< divDSigmaExpMethod << endl;
if
(
divDSigmaExpMethod != "standard"
&& divDSigmaExpMethod != "surface"
&& divDSigmaExpMethod != "decompose"
&& divDSigmaExpMethod != "laplacian"
)
{
FatalErrorIn(args.executable())
<< "divSigmaExp method " << divDSigmaExpMethod << " not found! "
<< "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
<< exit(FatalError);
}

58
applications/solvers/solidMechanics/elasticIncrSolidFoam/writeFields.H
View file

@ -1,36 +1,34 @@
if (runTime.outputTime())
{
{
volScalarField epsilonEq
(
IOobject
(
"epsilonEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
(
IOobject
(
"epsilonEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField sigmaEq
(
IOobject
(
"sigmaEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
(
IOobject
(
"sigmaEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
runTime.write();
}
}

6
applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/createFields.H
View file

@ -25,11 +25,11 @@
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedSymmTensor("zero", dimless, symmTensor::zero)
mesh,
dimensionedSymmTensor("zero", dimless, symmTensor::zero)
);
//- second Piola-Kirchhoff stress tensor
//- Second Piola-Kirchhoff stress tensor
volSymmTensorField sigma
(
IOobject

80
applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/elasticNonLinTLSolidFoam.C
View file

@ -44,21 +44,21 @@ Author
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "createSolidInterfaceNonLin.H"
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "createSolidInterfaceNonLin.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
Info<< "\nStarting time loop\n" << endl;
while(runTime.loop())
while(runTime.loop())
{
Info<< "Time: " << runTime.timeName() << nl << endl;
Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
# include "readSolidMechanicsControls.H"
int iCorr = 0;
scalar initialResidual = 0;
@ -70,36 +70,40 @@ int main(int argc, char *argv[])
{
U.storePrevIter();
surfaceTensorField shearGradU =
((I - n*n)&fvc::interpolate(gradU));
surfaceTensorField shearGradU
(
"shearGradU",
(I - sqr(n)) & fvc::interpolate(gradU)
);
fvVectorMatrix UEqn
(
fvm::d2dt2(rho, U)
==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
// + fvc::div
// (
// -(mu + lambda)*gradU
// + mu*gradU.T()
// + mu*(gradU & gradU.T())
// + lambda*tr(gradU)*I
// + 0.5*lambda*tr(gradU & gradU.T())*I
// + (sigma & gradU),
// "div(sigma)"
// )
+ fvc::div
(
fvm::d2dt2(rho, U)
==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
// + fvc::div(
// -( (mu + lambda) * gradU )
// + ( mu * gradU.T() )
// + ( mu * (gradU & gradU.T()) )
// + ( lambda * tr(gradU) * I )
// + ( 0.5 * lambda * tr(gradU & gradU.T()) * I )
// + ( sigma & gradU ),
// "div(sigma)"
// )
+ fvc::div(
mesh.magSf()
*(
- (muf + lambdaf)*(fvc::snGrad(U)&(I - n*n))
+ lambdaf*tr(shearGradU&(I - n*n))*n
+ muf*(shearGradU&n)
+ muf * (n & fvc::interpolate(gradU & gradU.T()))
+ 0.5*lambdaf
*(n * tr(fvc::interpolate(gradU & gradU.T())))
+ (n & fvc::interpolate( sigma & gradU ))
)
)
);
mesh.magSf()*
(
- (muf + lambdaf)*(fvc::snGrad(U) & (I - n*n))
+ lambdaf*tr(shearGradU & (I - n*n))*n
+ muf*(shearGradU & n)
+ muf*(n & fvc::interpolate(gradU & gradU.T()))
+ 0.5*lambdaf*(n*tr(fvc::interpolate(gradU & gradU.T())))
+ (n & fvc::interpolate(sigma & gradU))
)
)
);
if (solidInterfaceCorr)
{
@ -153,9 +157,9 @@ int main(int argc, char *argv[])
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
Info<< "End\n" << endl;
return(0);
return(0);
}

19
applications/solvers/solidMechanics/elasticNonLinULSolidFoam/moveMesh.H
View file

@ -1,15 +1,16 @@
if(moveMeshMethod == "inverseDistance")
{
{
# include "moveMeshInverseDistance.H"
}
else if(moveMeshMethod == "leastSquares")
{
}
else if(moveMeshMethod == "leastSquares")
{
# include "moveMeshLeastSquares.H"
}
else
{
FatalError << "move mesh method " << moveMeshMethod << " not recognised" << nl
}
else
{
FatalErrorIn(args.executable())
<< "move mesh method " << moveMeshMethod << " not recognised" << nl
<< "available methods are:" << nl
<< "inverseDistance" << nl
<< "leastSquares" << exit(FatalError);
}
}

40
applications/solvers/solidMechanics/elasticSolidFoam/aitkenRelaxation.H
View file

@ -1,4 +1,4 @@
// aitken acceleration
// aitken acceleration
aitkenDelta.storePrevIter();
@ -7,26 +7,26 @@ aitkenDelta = (U - U.prevIter()) / aitkenInitialRes;
// update relaxation factor
if(iCorr == 0)
{
{
aitkenTheta = 0.1;
}
else
{
vectorField b = aitkenDelta.internalField() - aitkenDelta.prevIter().internalField();
// scalar sumMagB = gSum(mag(b));
scalar sumMagB = gSum(magSqr(b));
if(sumMagB < SMALL)
{
// Warning << "Aitken under-relaxation: denominator less then SMALL"
// << endl;
sumMagB += SMALL;
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)
/
sumMagB;
}
}
else
{
vectorField b = aitkenDelta.internalField()
- aitkenDelta.prevIter().internalField();
scalar sumMagB = gSum(magSqr(b));
if(sumMagB < SMALL)
{
// Warning << "Aitken under-relaxation: denominator less then SMALL"
// << endl;
sumMagB += SMALL;
}
aitkenTheta = -aitkenTheta*
gSum(aitkenDelta.prevIter().internalField() & b)/
sumMagB;
}
// correction to the latest U
U += aitkenTheta*aitkenDelta*aitkenInitialRes;

63
applications/solvers/solidMechanics/elasticSolidFoam/createFields.H
View file

@ -12,21 +12,20 @@
mesh
);
// volTensorField gradU = fvc::grad(U);
volTensorField gradU
(
IOobject
(
"grad(U)",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
"grad(U)",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedTensor("zero", dimless, tensor::zero)
);
//surfaceVectorField snGradU = fvc::snGrad(U);
);
surfaceVectorField snGradU
(
IOobject
@ -38,7 +37,7 @@
IOobject::NO_WRITE
),
mesh,
dimensionedVector("zero", dimless, vector::zero)
dimensionedVector("zero", dimless, vector::zero)
);
volVectorField V
@ -85,19 +84,19 @@
dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
);
volVectorField divSigmaExp
(
volVectorField divSigmaExp
(
IOobject
(
"divSigmaExp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
"divSigmaExp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedVector("zero", dimForce/dimVolume, vector::zero)
);
);
constitutiveModel rheology(sigma, U);
@ -127,18 +126,20 @@
// for aitken relaxation
volVectorField aitkenDelta
(
IOobject
(
"aitkenDelta",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
(
IOobject
(
"aitkenDelta",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
// aitken relaxation factor
scalar aitkenInitialRes = 1.0;
scalar aitkenTheta = 0.1;

128
applications/solvers/solidMechanics/elasticSolidFoam/elasticSolidFoam.C
View file

@ -62,7 +62,7 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
while(runTime.loop())
{
{
Info<< "Time: " << runTime.timeName() << nl << endl;
# include "readSolidMechanicsControls.H"
@ -74,84 +74,80 @@ int main(int argc, char *argv[])
lduMatrix::debug = 0;
if (predictor)
{
Info << "\nPredicting U, gradU and snGradU based on V,"
<< "gradV and snGradV\n" << endl;
U += V*runTime.deltaT();
gradU += gradV*runTime.deltaT();
snGradU += snGradV*runTime.deltaT();
}
{
Info<< "\nPredicting U, gradU and snGradU based on V,"
<< "gradV and snGradV\n" << endl;
U += V*runTime.deltaT();
gradU += gradV*runTime.deltaT();
snGradU += snGradV*runTime.deltaT();
}
do
{
U.storePrevIter();
{
U.storePrevIter();
# include "calculateDivSigmaExp.H"
// linear momentum equation
fvVectorMatrix UEqn
(
rho*fvm::d2dt2(U)
==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+ divSigmaExp
);
// linear momentum equation
fvVectorMatrix UEqn
(
rho*fvm::d2dt2(U)
==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+ divSigmaExp
);
if (solidInterfaceCorr)
{
solidInterfacePtr->correct(UEqn);
}
if (solidInterfaceCorr)
{
solidInterfacePtr->correct(UEqn);
}
// if (relaxEqn)
// {
// UEqn.relax();
// }
solverPerf = UEqn.solve();
solverPerf = UEqn.solve();
if (iCorr == 0)
{
initialResidual = solverPerf.initialResidual();
aitkenInitialRes = gMax(mag(U.internalField()));
}
if (iCorr == 0)
{
initialResidual = solverPerf.initialResidual();
aitkenInitialRes = gMax(mag(U.internalField()));
}
if (aitkenRelax)
{
if (aitkenRelax)
{
# include "aitkenRelaxation.H"
}
else
{
U.relax();
}
}
else
{
U.relax();
}
gradU = fvc::grad(U);
gradU = fvc::grad(U);
# include "calculateRelativeResidual.H"
if (iCorr % infoFrequency == 0)
{
Info<< "\tTime " << runTime.value()
<< ", Corrector " << iCorr
<< ", Solving for " << U.name()
<< " using " << solverPerf.solverName()
<< ", res = " << solverPerf.initialResidual()
<< ", rel res = " << relativeResidual;
if (aitkenRelax)
{
Info<< ", aitken = " << aitkenTheta;
}
Info<< ", inner iters = " << solverPerf.nIterations() << endl;
}
}
if (iCorr % infoFrequency == 0)
{
Info<< "\tTime " << runTime.value()
<< ", Corrector " << iCorr
<< ", Solving for " << U.name()
<< " using " << solverPerf.solverName()
<< ", res = " << solverPerf.initialResidual()
<< ", rel res = " << relativeResidual;
if (aitkenRelax)
{
Info<< ", aitken = " << aitkenTheta;
}
Info<< ", inner iters = " << solverPerf.nIterations() << endl;
}
}
while
(
iCorr++ == 0
||
(solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance
&&
iCorr < nCorr)
);
(
iCorr++ == 0
|| (
solverPerf.initialResidual() > convergenceTolerance
//relativeResidual > convergenceTolerance
&& iCorr < nCorr
)
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << U.name()
<< ", Initial residual = " << initialResidual
@ -178,9 +174,9 @@ int main(int argc, char *argv[])
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
Info<< "End\n" << endl;
return(0);
return(0);
}