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

Tutorial fix elasticThermalSolidFoam

Browse source
This commit is contained in:
Hrvoje Jasak 2014-04-10 20:54:24 +01:00
parent 45576c2dea
commit e732980a53
8 changed files with 181 additions and 190 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

57
applications/solvers/solidMechanics/elasticThermalSolidFoam/createFields.H
View file

@ -12,21 +12,20 @@
mesh mesh
); );
// volTensorField gradU = fvc::grad(U);
volTensorField gradU volTensorField gradU
( (
IOobject IOobject
( (
"grad(U)", "grad(U)",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedTensor("zero", dimless, tensor::zero) dimensionedTensor("zero", dimless, tensor::zero)
); );
//surfaceVectorField snGradU = fvc::snGrad(U);
surfaceVectorField snGradU surfaceVectorField snGradU
( (
IOobject IOobject
@ -38,7 +37,7 @@
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimless, vector::zero) dimensionedVector("zero", dimless, vector::zero)
); );
volSymmTensorField epsilon volSymmTensorField epsilon
@ -115,31 +114,33 @@
thermalModel thermal(T); thermalModel thermal(T);
volScalarField C = thermal.C(); volScalarField C = thermal.C();
volScalarField k = thermal.k(); volScalarField k
(
"DT",
thermal.k()
);
volScalarField threeKalpha = rheology.threeK()*rho*thermal.alpha(); volScalarField threeKalpha = rheology.threeK()*rho*thermal.alpha();
surfaceScalarField threeKalphaf = fvc::interpolate(threeKalpha, "threeKalpha"); surfaceScalarField threeKalphaf =
fvc::interpolate(threeKalpha, "threeKalpha");
volScalarField T0 = thermal.T0(); volScalarField T0 = thermal.T0();
volScalarField rhoC = rho*C; volScalarField rhoC = rho*C;
// for aitken relaxation // for aitken relaxation
volVectorField aitkenDelta volVectorField aitkenDelta
( (
IOobject IOobject
( (
"aitkenDelta", "aitkenDelta",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensionedVector("zero", dimLength, vector::zero) dimensionedVector("zero", dimLength, vector::zero)
); );
// aitken relaxation factor
scalar aitkenInitialRes = 1.0; // aitken relaxation factor
scalar aitkenTheta = 0.01; scalar aitkenInitialRes = 1.0;
// if(mesh.relax(U.name())) scalar aitkenTheta = 0.01;
// {
// aitkenTheta = mesh.relaxationFactor(U.name());
// }

171
applications/solvers/solidMechanics/elasticThermalSolidFoam/elasticThermalSolidFoam.C
View file

@ -46,100 +46,99 @@ Author
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "setRootCase.H" # include "setRootCase.H"
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
# include "createFields.H" # include "createFields.H"
# include "readDivSigmaExpMethod.H" # include "readDivSigmaExpMethod.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; int iCorr = 0;
scalar initialResidual = 1.0; scalar initialResidual = 1.0;
scalar relResT = 1.0; scalar relResT = 1.0;
scalar relResU = 1.0; scalar relResU = 1.0;
lduMatrix::solverPerformance solverPerfU; lduMatrix::solverPerformance solverPerfU;
lduMatrix::solverPerformance solverPerfT; lduMatrix::solverPerformance solverPerfT;
lduMatrix::debug = 0; lduMatrix::debug = 0;
// solve energy equation for temperature // solve energy equation for temperature
// the loop is for non-orthogonal corrections // the loop is for non-orthogonal corrections
Info<< "Solving for " << T.name() << nl; Info<< "Solving for " << T.name() << nl;
do do
{ {
T.storePrevIter(); T.storePrevIter();
fvScalarMatrix TEqn fvScalarMatrix TEqn
( (
rhoC*fvm::ddt(T) == fvm::laplacian(k, T, "laplacian(k,T)") rhoC*fvm::ddt(T) == fvm::laplacian(k, T, "laplacian(k,T)")
); );
solverPerfT = TEqn.solve(); solverPerfT = TEqn.solve();
T.relax(); T.relax();
# include "calculateRelResT.H" # include "calculateRelResT.H"
if (iCorr % infoFrequency == 0) if (iCorr % infoFrequency == 0)
{ {
Info<< "\tCorrector " << iCorr Info<< "\tCorrector " << iCorr
<< ", residual = " << solverPerfT.initialResidual() << ", residual = " << solverPerfT.initialResidual()
<< ", relative res = " << relResT << ", relative res = " << relResT
<< ", inner iters = " << solverPerfT.nIterations() << endl; << ", inner iters = " << solverPerfT.nIterations() << endl;
} }
} }
while while
( (
relResT > convergenceToleranceT relResT > convergenceToleranceT
&& && ++iCorr < nCorr
++iCorr < nCorr );
);
Info<< "Solved for " << T.name() Info<< "Solved for " << T.name()
<< " using " << solverPerfT.solverName() << " using " << solverPerfT.solverName()
<< " in " << iCorr << " iterations" << " in " << iCorr << " iterations"
<< ", residual = " << solverPerfT.initialResidual() << ", residual = " << solverPerfT.initialResidual()
<< ", relative res = " << relResT << nl << ", relative res = " << relResT << nl
<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< ", ClockTime = " << runTime.elapsedClockTime() << " s" << ", ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
// solve momentum equation for displacement // Solve momentum equation for displacement
iCorr = 0; iCorr = 0;
volVectorField gradThreeKalphaDeltaT = volVectorField gradThreeKalphaDeltaT =
fvc::grad(threeKalpha*(T-T0), "grad(threeKalphaDeltaT)"); fvc::grad(threeKalpha*(T-T0), "grad(threeKalphaDeltaT)");
surfaceVectorField threeKalphaDeltaTf = surfaceVectorField threeKalphaDeltaTf =
mesh.Sf()*threeKalphaf*fvc::interpolate(T-T0, "deltaT"); mesh.Sf()*threeKalphaf*fvc::interpolate(T-T0, "deltaT");
Info<< "Solving for " << U.name() << nl; Info<< "Solving for " << U.name() << nl;
do do
{ {
U.storePrevIter(); U.storePrevIter();
# include "calculateDivSigmaExp.H" # include "calculateDivSigmaExp.H"
// Linear momentum equaiton // Linear momentum equaiton
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
rho*fvm::d2dt2(U) rho*fvm::d2dt2(U)
== ==
fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)") fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+ divSigmaExp + divSigmaExp
); );
solverPerfU = UEqn.solve(); solverPerfU = UEqn.solve();
if (aitkenRelax) if (aitkenRelax)
{ {
# include "aitkenRelaxation.H" # include "aitkenRelaxation.H"
} }
else else
{ {
@ -148,7 +147,7 @@ int main(int argc, char *argv[])
gradU = fvc::grad(U); gradU = fvc::grad(U);
# include "calculateRelResU.H" # include "calculateRelResU.H"
if (iCorr == 0) if (iCorr == 0)
{ {
@ -160,6 +159,7 @@ int main(int argc, char *argv[])
Info<< "\tCorrector " << iCorr Info<< "\tCorrector " << iCorr
<< ", residual = " << solverPerfU.initialResidual() << ", residual = " << solverPerfU.initialResidual()
<< ", relative res = " << relResU; << ", relative res = " << relResU;
if (aitkenRelax) if (aitkenRelax)
{ {
Info << ", aitken = " << aitkenTheta; Info << ", aitken = " << aitkenTheta;
@ -167,25 +167,24 @@ int main(int argc, char *argv[])
Info<< ", inner iters = " << solverPerfU.nIterations() << endl; Info<< ", inner iters = " << solverPerfU.nIterations() << endl;
} }
} }
while while
( (
iCorr++ == 0 iCorr++ == 0
|| || (
(//solverPerfU.initialResidual() > convergenceTolerance relResU > convergenceToleranceU
relResU > convergenceToleranceU && iCorr < nCorr
&& )
iCorr < nCorr) );
);
Info<< "Solved for " << U.name() Info<< "Solved for " << U.name()
<< " using " << solverPerfU.solverName() << " using " << solverPerfU.solverName()
<< " in " << iCorr << " iterations" << " in " << iCorr << " iterations"
<< ", initial res = " << initialResidual << ", initial res = " << initialResidual
<< ", final res = " << solverPerfU.initialResidual() << ", final res = " << solverPerfU.initialResidual()
<< ", final rel res = " << relResU << nl << ", final rel res = " << relResU << nl
<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< ", ClockTime = " << runTime.elapsedClockTime() << " s" << ", ClockTime = " << runTime.elapsedClockTime() << " s"
<< endl; << endl;
# include "calculateEpsilonSigma.H" # include "calculateEpsilonSigma.H"
# include "writeFields.H" # include "writeFields.H"

18
tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/0/T
View file

@ -32,27 +32,23 @@ boundaryField
inside inside
{ {
type fixedValue; type fixedValue;
value uniform 100; value uniform 100;
} }
outside outside
{ {
type fixedValue; type fixedValue;
value uniform 100; value uniform 100;
} }
front front
{ {
type empty; type empty;
//type symmetryPlane;
} }
back back
{ {
type empty; type empty;
//type symmetryPlane;
} }
} }
// ************************************************************************* // // ************************************************************************* //

19
tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/0/U
View file

@ -29,32 +29,31 @@ boundaryField
{ {
type symmetryPlane; type symmetryPlane;
} }
inside inside
{ {
type solidTraction; type solidTraction;
traction uniform ( 0 0 0 ); traction uniform ( 0 0 0 );
pressure uniform 50e6; pressure uniform 50e6;
value uniform (0 0 0); DT k;
value uniform (0 0 0);
} }
outside outside
{ {
type solidTraction; type solidTraction;
traction uniform ( 0 0 0 ); traction uniform ( 0 0 0 );
pressure uniform 0.1e6; pressure uniform 0.1e6;
DT k;
value uniform (0 0 0); value uniform (0 0 0);
} }
front front
{ {
type empty; type empty;
} }
back back
{ {
type empty; type empty;
} }
} }
// ************************************************************************* // // ************************************************************************* //

59
tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/polyMesh/blockMeshDict
View file

@ -19,15 +19,14 @@ convertToMeters 1;
vertices vertices
( (
(0 0.5 0) (0 0.5 0)
(0.5 0 0) (0.5 0 0)
(0.7 0 0) (0.7 0 0)
(0 0.7 0) (0 0.7 0)
(0 0.5 0.1) (0 0.5 0.1)
(0.5 0 0.1) (0.5 0 0.1)
(0.7 0 0.1) (0.7 0 0.1)
(0 0.7 0.1) (0 0.7 0.1)
); );
blocks blocks
@ -45,34 +44,30 @@ edges
patches patches
( (
symmetryPlane left
(
(4 7 3 0)
)
symmetryPlane bottom
(
(1 2 6 5)
)
patch inside
(
(0 1 5 4)
)
patch outside
(
(7 6 2 3)
)
empty back empty back
( (
(3 2 1 0) (3 2 1 0)
) )
empty front empty front
( (
(4 5 6 7) (4 5 6 7)
) )
symmetryPlane left
(
(4 7 3 0)
)
symmetryPlane bottom
(
(1 2 6 5)
)
patch inside
(
(0 1 5 4)
)
patch outside
(
(7 6 2 3)
)
); );
mergePatchPairs mergePatchPairs

34
tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/polyMesh/boundary
View file

@ -1,7 +1,7 @@
/*--------------------------------*- C++ -*----------------------------------*\ /*--------------------------------*- C++ -*----------------------------------*\
| ========= | | | ========= | |
| \\ / F ield | foam-extend: Open Source CFD | | \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.0 | | \\ / O peration | Version: 3.0 |
| \\ / A nd | Web: http://www.extend-project.de | | \\ / A nd | Web: http://www.extend-project.de |
| \\/ M anipulation | | | \\/ M anipulation | |
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -17,41 +17,41 @@ FoamFile
6 6
( (
back
{
type empty;
nFaces 600;
startFace 1130;
}
front
{
type empty;
nFaces 600;
startFace 1730;
}
left left
{ {
type symmetryPlane; type symmetryPlane;
nFaces 10; nFaces 10;
startFace 2330; startFace 1130;
} }
bottom bottom
{ {
type symmetryPlane; type symmetryPlane;
nFaces 10; nFaces 10;
startFace 2340; startFace 1140;
} }
inside inside
{ {
type patch; type patch;
nFaces 60; nFaces 60;
startFace 2350; startFace 1150;
} }
outside outside
{ {
type patch; type patch;
nFaces 60; nFaces 60;
startFace 2410; startFace 1210;
}
back
{
type empty;
nFaces 600;
startFace 1270;
}
front
{
type empty;
nFaces 600;
startFace 1870;
} }
) )

11
tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/thermalProperties
View file

@ -16,11 +16,12 @@ FoamFile
thermal thermal
{ {
type constant; type constant;
C C [0 2 -2 -1 0 0 0] 434; C C [0 2 -2 -1 0 0 0] 434;
k k [1 1 -3 -1 0 0 0] 250; k k [1 1 -3 -1 0 0 0] 250;
alpha alpha [0 0 0 -1 0 0 0] 2.3e-05; alpha alpha [0 0 0 -1 0 0 0] 2.3e-05;
T0 T0 [0 0 0 1 0 0 0] 0; T0 T0 [0 0 0 1 0 0 0] 0;
} }
// ************************************************************************* // // ************************************************************************* //