Update and clean-up of solid mechanics

This commit is contained in:
Hrvoje Jasak 2014-04-10 17:45:26 +01:00
parent 158e6ccebf
commit a3b9ee8fac
32 changed files with 483 additions and 359 deletions

View file

@ -20,5 +20,7 @@ wmake elasticThermalSolidFoam
wmake icoFsiElasticNonLinULSolidFoam
wmake viscoElasticSolidFoam
wmake stressFemFoam
(cd utilities ; wmake all)
(cd deprecatedSolvers ; ./Allwmake)

View file

@ -11,8 +11,6 @@ wmake contactStressFoam
wmake newStressedFoam
wmake newContactStressFoam
wmake stressFemFoam
wmake icoFsiFoam
wmake solidDisplacementFoam

View file

@ -1,53 +0,0 @@
if (runTime.outputTime())
{
// Displacement gradient
tetPointTensorField gradU = tetFec::grad(U);
// Stress tensor
tetPointSymmTensorField sigma =
rho*(2.0*mu*symm(gradU) + lambda*I*tr(gradU));
// sigmaXX
tetPointScalarField sigmaXX
(
IOobject
(
"sigmaXX",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sigma.component(symmTensor::XX)
);
// sigmaYY
tetPointScalarField sigmaYY
(
IOobject
(
"sigmaYY",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sigma.component(symmTensor::YY)
);
// sigmaXY
tetPointScalarField sigmaXY
(
IOobject
(
"sigmaXY",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sigma.component(symmTensor::XY)
);
runTime.write();
}

View file

@ -1,5 +1,5 @@
//- Increment of Green finite strain tensor
// Increment of Green finite strain tensor
DEpsilon = symm(gradDU) + 0.5*symm((gradDU & gradU.T()) + (gradU & gradDU.T()) + (gradDU & gradDU.T()));
//- Increment of second Piola-Kirchhoff stress tensor
// Increment of second Piola-Kirchhoff stress tensor
DSigma = 2*mu*(DEpsilon - DEpsilonP) + lambda*(I*tr(DEpsilon));

View file

@ -1,5 +1,5 @@
if(divDSigmaExpMethod == "standard")
{
{
divDSigmaExp = fvc::div
(
(mu*gradDU.T())
@ -7,39 +7,40 @@ if(divDSigmaExpMethod == "standard")
- ((mu + lambda)*gradDU),
"div(sigma)"
);
}
else if(divDSigmaExpMethod == "surface")
{
divDSigmaExp =
fvc::div(
mesh.magSf()
*(
muf*(n&fvc::interpolate(gradDU.T()))
}
else if(divDSigmaExpMethod == "surface")
{
divDSigmaExp = fvc::div
(
mesh.magSf()*
(
muf*(n & fvc::interpolate(gradDU.T()))
+ lambdaf*tr(fvc::interpolate(gradDU))*n
- (muf + lambdaf)*(n&fvc::interpolate(gradDU))
- (muf + lambdaf)*(n & fvc::interpolate(gradDU))
)
);
// divDSigmaExp = fvc::div
// (
// muf*(mesh.Sf() & fvc::interpolate(gradDU.T()))
// + lambdaf*(mesh.Sf() & I*fvc::interpolate(tr(gradDU)))
// - (muf + lambdaf)*(mesh.Sf() & fvc::interpolate(gradDU))
// );
}
else if(divDSigmaExpMethod == "decompose")
{
surfaceTensorField shearGradDU =
((I - n*n)&fvc::interpolate(gradDU));
}
else if(divDSigmaExpMethod == "decompose")
{
surfaceTensorField shearGradDU((I - n*n) & fvc::interpolate(gradDU));
divDSigmaExp =
fvc::div(
mesh.magSf()
*(
- (muf + lambdaf)*(fvc::snGrad(DU)&(I - n*n))
+ lambdaf*tr(shearGradDU&(I - n*n))*n
+ muf*(shearGradDU&n)
divDSigmaExp = fvc::div
(
mesh.magSf()*
(
- (muf + lambdaf)*(fvc::snGrad(DU) & (I - n*n))
+ lambdaf*tr(shearGradDU & (I - n*n))*n
+ muf*(shearGradDU & n)
)
);
// divDSigmaExp = fvc::div
// (
// mesh.magSf()
@ -49,9 +50,9 @@ if(divDSigmaExpMethod == "standard")
// + muf*(shearGradDU&n)
// )
// );
}
else if(divDSigmaExpMethod == "laplacian")
{
}
else if(divDSigmaExpMethod == "laplacian")
{
divDSigmaExp =
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
@ -60,8 +61,10 @@ if(divDSigmaExpMethod == "standard")
+ lambda*(I*tr(gradDU)),
"div(sigma)"
);
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
}
}
else
{
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -1,30 +1,67 @@
if(divDSigmaNonLinExpMethod == "standard")
{
{
divDSigmaNonLinExp = fvc::div
(
( mu * (
(
mu*
(
(gradDU & gradU.T())
+ (gradU & gradDU.T())
+ (gradDU & gradDU.T())
) )
+ ( lambda * 0.5 * tr( (gradDU & gradU.T())
)
)
+ (
0.5*lambda*
tr
(
(gradDU & gradU.T())
+ (gradU & gradDU.T())
+ (gradDU & gradDU.T())
) * I )
+ ( DSigma & gradU )
+ ( (sigma + DSigma) & gradDU ),
)*I
)
+ (DSigma & gradU)
+ ((sigma + DSigma) & gradDU),
"div(sigma)"
);
}
else if(divDSigmaNonLinExpMethod == "surface")
{
divDSigmaNonLinExp =
fvc::div(
mesh.magSf()
*(
( muf * (n & fvc::interpolate( (gradU & gradDU.T()) + (gradDU & gradU.T()) + (gradDU & gradDU.T()) )) )
+ ( 0.5*lambdaf * (n * tr(fvc::interpolate( (gradU & gradDU.T()) + (gradDU & gradU.T()) + (gradDU & gradDU.T()) ))) )
+ (n & fvc::interpolate( (DSigma & gradU) + ((sigma + DSigma) & gradDU) ))
}
else if(divDSigmaNonLinExpMethod == "surface")
{
divDSigmaNonLinExp = fvc::div
(
mesh.magSf()*
(
(
muf*
(
n & fvc::interpolate
(
(gradU & gradDU.T())
+ (gradDU & gradU.T())
+ (gradDU & gradDU.T())
)
)
)
+ (
0.5*lambdaf*
(
n*tr
(
fvc::interpolate
(
(gradU & gradDU.T())
+ (gradDU & gradU.T())
+ (gradDU & gradDU.T())
)
)
)
)
+ (
n & fvc::interpolate
(
(DSigma & gradU)
+ ((sigma + DSigma) & gradDU)
)
)
)
);
@ -48,8 +85,10 @@ if(divDSigmaNonLinExpMethod == "standard")
// + ( (sigma + DSigma) & gradDU )
// ) )
// );
}
else
{
FatalError << "divDSigmaExp method " << divDSigmaExpMethod << " not found!" << endl;
}
}
else
{
FatalErrorIn(args.executable())
<< "divDSigmaExp method " << divDSigmaExpMethod << " not found!"
<< abort(FatalError);
}

View file

@ -66,7 +66,7 @@ int main(int argc, char *argv[])
scalar initialResidual = 0;
lduMatrix::solverPerformance solverPerf;
scalar relativeResidual = GREAT;
lduMatrix::debug=0;
lduMatrix::debug = 0;
do
{
@ -88,7 +88,7 @@ int main(int argc, char *argv[])
//- fvc::div(2*mu*DEpsilonP, "div(sigma)")
- fvc::div
(
2*muf*( mesh.Sf() & fvc::interpolate(DEpsilonP))
2*muf*(mesh.Sf() & fvc::interpolate(DEpsilonP))
)
);
@ -108,7 +108,7 @@ int main(int argc, char *argv[])
(
(2*mu + lambda)*gradU, DU, "laplacian(DDU,DU)"
)
- fvc::div( (2*mu + lambda)*(gradU&gradDU), "div(sigma)");
- fvc::div((2*mu + lambda)*(gradU & gradDU), "div(sigma)");
}
if (nonLinearSemiImplicit)
@ -123,7 +123,7 @@ int main(int argc, char *argv[])
(
(2*mu + lambda)*gradDU, DU, "laplacian(DDU,DU)"
)
- fvc::div( (2*mu + lambda)*(gradDU&gradDU), "div(sigma)");
- fvc::div((2*mu + lambda)*(gradDU & gradDU), "div(sigma)");
}
solverPerf = DUEqn.solve();
@ -152,12 +152,13 @@ int main(int argc, char *argv[])
if (iCorr % infoFrequency == 0)
{
Info << "\tTime " << runTime.value()
Info<< "\tTime " << runTime.value()
<< ", Corrector " << iCorr
<< ", Solving for " << DU.name()
<< " using " << solverPerf.solverName()
<< ", res = " << solverPerf.initialResidual()
<< ", rel res = " << relativeResidual;
if (aitkenRelax)
{
Info << ", aitken = " << aitkenTheta;
@ -169,10 +170,11 @@ int main(int argc, char *argv[])
(
iCorr++ == 0
||
(//solverPerf.initialResidual() > convergenceTolerance
(
//solverPerf.initialResidual() > convergenceTolerance
relativeResidual > convergenceTolerance
&&
iCorr < nCorr)
&& iCorr < nCorr
)
);
Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()

View file

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

View file

@ -1,5 +1,5 @@
if (runTime.outputTime())
{
{
volScalarField epsilonEq
(
IOobject
@ -12,8 +12,7 @@ if (runTime.outputTime())
),
sqrt((2.0/3.0)*magSqr(dev(epsilon)))
);
Info<< "Max epsilonEq = " << max(epsilonEq).value()
<< endl;
Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
volScalarField epsilonPEq
(
@ -27,8 +26,7 @@ if (runTime.outputTime())
),
sqrt((2.0/3.0)*magSqr(dev(epsilonP)))
);
Info<< "Max epsilonPEq = " << max(epsilonPEq).value()
<< endl;
Info<< "Max epsilonPEq = " << max(epsilonPEq).value()<< endl;
volScalarField sigmaEq
(
@ -43,14 +41,13 @@ if (runTime.outputTime())
sqrt((3.0/2.0)*magSqr(dev(sigma)))
);
Info<< "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
// deformation gradient
// Deformation gradient
volTensorField F = I + gradU;
volScalarField J = det(F);
//- Calculate Cauchy stress
// Calculate Cauchy stress
volSymmTensorField sigmaCauchy
(
IOobject
@ -61,10 +58,10 @@ if (runTime.outputTime())
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
(1/J) * symm(F.T() & sigma & F)
(1/J)*symm(F.T() & sigma & F)
);
//- Cauchy von Mises stress
// Cauchy von Mises stress
volScalarField sigmaCauchyEq
(
IOobject
@ -78,11 +75,10 @@ if (runTime.outputTime())
sqrt((3.0/2.0)*magSqr(dev(sigmaCauchy)))
);
Info<< "Max sigmaCauchyEq = " << max(sigmaCauchyEq).value()
<< endl;
Info<< "Max sigmaCauchyEq = " << max(sigmaCauchyEq).value() << endl;
//volTensorField Finv = inv(F);
// volSymmTensorField epsilonAlmansi
// volTensorField Finv = inv(F);
// volSymmTensorField epsilonAlmansi
// (
// IOobject
// (
@ -95,6 +91,5 @@ if (runTime.outputTime())
// symm(Finv & epsilon & Finv.T())
// );
runTime.write();
}
}

View file

@ -0,0 +1,120 @@
if (runTime.outputTime())
{
// Displacement gradient
tetPointTensorField gradU = tetFec::grad(U);
// Stress tensor
tetPointSymmTensorField sigma =
rho*(2.0*mu*symm(gradU) + lambda*I*tr(gradU));
// Create pointMesh for field post-processing
const pointMesh& pMesh = pointMesh::New(mesh);
// U
pointVectorField Up
(
IOobject
(
"Up",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
U.dimensions()
);
Up.internalField() = vectorField::subField
(
U.internalField(),
pMesh.size()
);
// sigmaEq
pointScalarField sigmaEq
(
IOobject
(
"sigmaEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
sigma.dimensions()
);
sigmaEq.internalField() = scalarField::subField
(
sqrt((3.0/2.0)*magSqr(dev(sigma.internalField())))(),
pMesh.size()
);
// sigmaXX
pointScalarField sigmaXX
(
IOobject
(
"sigmaXX",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
sigma.dimensions()
);
sigmaXX.internalField() = scalarField::subField
(
sigma.component(symmTensor::XX)().internalField(),
pMesh.size()
);
// sigmaYY
pointScalarField sigmaYY
(
IOobject
(
"sigmaYY",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
sigma.dimensions()
);
sigmaYY.internalField() = scalarField::subField
(
sigma.component(symmTensor::YY)().internalField(),
pMesh.size()
);
// sigmaXY
pointScalarField sigmaXY
(
IOobject
(
"sigmaXY",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
sigma.dimensions()
);
sigmaXY.internalField() = scalarField::subField
(
sigma.component(symmTensor::XY)().internalField(),
pMesh.size()
);
runTime.write();
}

View file

@ -46,8 +46,9 @@ template
class PointPatch,
template<class> class MatrixType
>
void TractionPointPatchVectorField<PatchField, Mesh, PointPatch, MatrixType>
::checkFieldSize() const
void
TractionPointPatchVectorField<PatchField, Mesh, PointPatch, MatrixType>::
checkFieldSize() const
{
if
(
@ -229,7 +230,7 @@ addBoundarySourceDiag
const vectorField& pointNormals = this->patch().pointNormals();
const label nFaces = this->patch().nFaces();
vectorField& source = matrix.source();
vectorField& b = matrix.b();
for (label faceI = 0; faceI < nFaces; faceI++)
{
@ -240,7 +241,7 @@ addBoundarySourceDiag
{
const triFace& tri = faceTriangles[triI];
source[meshPoints[tri[0]]] +=
b[meshPoints[tri[0]]] +=
tri.mag(points)*
(
traction_[faceI]/6.0
@ -251,7 +252,7 @@ addBoundarySourceDiag
- pressure_[faceI]*pointNormals[tri[2]]/12.0
);
source[meshPoints[tri[1]]] +=
b[meshPoints[tri[1]]] +=
tri.mag(points)*
(
traction_[faceI]/6.0
@ -262,7 +263,7 @@ addBoundarySourceDiag
- pressure_[faceI]*pointNormals[tri[0]]/12.0
);
source[meshPoints[tri[2]]] +=
b[meshPoints[tri[2]]] +=
tri.mag(points)*
(
traction_[faceI]/6.0

View file

@ -36,6 +36,8 @@ Description
#include "fvCFD.H"
#include "pointMesh.H"
#include "pointFields.H"
#include "tetPolyMesh.H"
#include "tetPointFields.H"
#include "tetFem.H"
@ -69,7 +71,7 @@ int main(int argc, char *argv[])
+ tetFem::laplacianTrace(lambda, U)
);
solve(UEqn);
UEqn.solve();
# include "calculateStress.H"

View file

@ -19,28 +19,28 @@ convertToMeters 1;
vertices
(
(0.5 0 0)
(1 0 0)
(2 0 0)
(2 0.707107 0)
(0.707107 0.707107 0)
(0.353553 0.353553 0)
(2 2 0)
(0.707107 2 0)
(0 2 0)
(0 1 0)
(0 0.5 0)
(0.5 0 0.5)
(1 0 0.5)
(2 0 0.5)
(2 0.707107 0.5)
(0.707107 0.707107 0.5)
(0.353553 0.353553 0.5)
(2 2 0.5)
(0.707107 2 0.5)
(0 2 0.5)
(0 1 0.5)
(0 0.5 0.5)
(0.5 0 -0.1)
(1 0 -0.1)
(2 0 -0.1)
(2 0.707107 -0.1)
(0.707107 0.707107 -0.1)
(0.353553 0.353553 -0.1)
(2 2 -0.1)
(0.707107 2 -0.1)
(0 2 -0.1)
(0 1 -0.1)
(0 0.5 -0.1)
(0.5 0 0.1)
(1 0 0.1)
(2 0 0.1)
(2 0.707107 0.1)
(0.707107 0.707107 0.1)
(0.353553 0.353553 0.1)
(2 2 0.1)
(0.707107 2 0.1)
(0 2 0.1)
(0 1 0.1)
(0 0.5 0.1)
);
blocks
@ -54,14 +54,14 @@ blocks
edges
(
arc 0 5 (0.469846 0.17101 0)
arc 5 10 (0.17101 0.469846 0)
arc 1 4 (0.939693 0.34202 0)
arc 4 9 (0.34202 0.939693 0)
arc 11 16 (0.469846 0.17101 0.5)
arc 16 21 (0.17101 0.469846 0.5)
arc 12 15 (0.939693 0.34202 0.5)
arc 15 20 (0.34202 0.939693 0.5)
arc 0 5 (0.469846 0.17101 -0.1)
arc 5 10 (0.17101 0.469846 -0.1)
arc 1 4 (0.939693 0.34202 -0.1)
arc 4 9 (0.34202 0.939693 -0.1)
arc 11 16 (0.469846 0.17101 0.1)
arc 16 21 (0.17101 0.469846 0.1)
arc 12 15 (0.939693 0.34202 0.1)
arc 15 20 (0.34202 0.939693 0.1)
);
patches

View file

@ -23,13 +23,13 @@ startTime 0;
stopAt endTime;
endTime 100;
endTime 1;
deltaT 1;
writeControl timeStep;
writeInterval 20;
writeInterval 1;
purgeWrite 0;

View file

@ -19,11 +19,15 @@ solvers
{
U
{
solver PCG;
preconditioner DIC;
solver CG;
preconditioner Cholesky;
minIter 1;
maxIter 2000;
tolerance 1e-06;
relTol 0.1;
relTol 0;
}
}
// ************************************************************************* //