From 995e56c8a703711723eb660bac3b83237564c7d3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 2 Dec 2013 11:17:49 +0000 Subject: [PATCH] Rewrite and clean-up --- .../analyticalHotCylinder.C | 564 ++++++++---------- 1 file changed, 244 insertions(+), 320 deletions(-) diff --git a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C index 39a3dd7a4..c89c7c5f3 100644 --- a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C +++ b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C @@ -26,7 +26,7 @@ Description Generate analytical solution for a thick-walled cylinder with a temperature gradient. Temperature field T and stress field sigma and generated. - Based on solution outlined in Timoshenko, Theory of Elasticity. + Based on solution outlined in Timoshenko, Theory of Elasticity. Author philip.cardiff@ucd.ie @@ -43,351 +43,275 @@ int main(int argc, char *argv[]) # include "createTime.H" # include "createMesh.H" - runTime++; + runTime++; - Info << "Writing analytical solution for a plain strain cylinder with concentric hole,\nwhere" - << "\n\tinner radius = 0.5" - << "\n\touter radius = 0.7" - << "\n\tinner temperature = 100" - << "\n\touter temperature = 0" - << "\n\tinner pressure = 0" - << "\n\touter pressure = 0" - << "\n\tE = 200e9" - << "\n\tu = 0.3" - << "\n\talpha = 1e-5" - << nl << endl; - - //- inner and outer radii and temperatures - scalar a = 0.5; - scalar b = 0.7; - scalar Ti = 100; - scalar To = 0; + Info<< "Writing analytical solution for a plain strain cylinder " + << "with concentric hole,\nwhere" + << "\n\tinner radius = 0.5" + << "\n\touter radius = 0.7" + << "\n\tinner temperature = 100" + << "\n\touter temperature = 0" + << "\n\tinner pressure = 0" + << "\n\touter pressure = 0" + << "\n\tE = 200e9" + << "\n\tu = 0.3" + << "\n\talpha = 1e-5" + << nl << endl; - //- mechanical and thermal properties - scalar E = 200e9; - scalar nu = 0.3; - scalar alpha = 1e-5; + //- inner and outer radii and temperatures + scalar a = 0.5; + scalar b = 0.7; + scalar Ti = 100; + scalar To = 0; - //- create T field - volScalarField T + //- mechanical and thermal properties + scalar E = 200e9; + scalar nu = 0.3; + scalar alpha = 1e-5; + + const volVectorField& C = mesh.C(); + + //- radial coordinate + volScalarField radii ( - IOobject - ( - "analyticalT", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimTemperature, 0.0) - ); - - const volVectorField& C = mesh.C(); + sqrt + ( + sqr(C.component(vector::X)) + + sqr(C.component(vector::Y)) + )/dimensionedScalar("one", dimLength, 1) + ); - //- radial coordinate - volScalarField radii = - C.component(vector::X)*C.component(vector::X) + C.component(vector::Y)*C.component(vector::Y); - forAll(radii.internalField(), celli) - { - radii.internalField()[celli] = ::sqrt(radii.internalField()[celli]); - } - forAll(radii.boundaryField(), patchi) - { - forAll(radii.boundaryField()[patchi], facei) - { - radii.boundaryField()[patchi][facei] = ::sqrt(radii.boundaryField()[patchi][facei]); - } - } + const scalarField& rIn = radii.internalField(); - forAll(T.internalField(), celli) - { - const scalar& r = radii[celli]; - - T.internalField()[celli] = - ( (Ti-To)/Foam::log(b/a) ) * Foam::log(b/r); - } - - forAll(T.boundaryField(), patchi) - { - forAll(T.boundaryField()[patchi], facei) - { - const scalar& r = radii.boundaryField()[patchi][facei]; - - T.boundaryField()[patchi][facei] = - ( (Ti-To)/Foam::log(b/a) ) * Foam::log(b/r); - } - } - - //- write temperature file - Info << "Writing analytical termpature field" << endl; - T.write(); - - - //- create sigma field - volScalarField sigmaR + Info << "Writing analytical termpature field" << endl; + //- create T field + volScalarField T ( - IOobject - ( - "sigmaR", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimForce/dimArea, 0.0) - ); + IOobject + ( + "analyticalT", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + ((Ti - To)/Foam::log(b/a))*Foam::log(b/radii) + ); + T.write(); - forAll(sigmaR.internalField(), celli) - { - const scalar& r = radii.internalField()[celli]; - - sigmaR.internalField()[celli] = - ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) * - (-Foam::log(b/r) -( a*a/(b*b - a*a))*(1 - (b*b)/(r*r))*Foam::log(b/a)); - } - - forAll(sigmaR.boundaryField(), patchi) - { - forAll(sigmaR.boundaryField()[patchi], facei) - { - const scalar& r = radii.boundaryField()[patchi][facei]; - - sigmaR.boundaryField()[patchi][facei] = - ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) * - ( -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 - (b*b)/(r*r))*Foam::log(b/a) ); - } - } - - //- write temperature file - Info << "\nWriting analytical sigmaR field" << endl; - sigmaR.write(); - - - volScalarField sigmaTheta + //- create sigma field + Info << "\nWriting analytical sigmaR field" << endl; + volScalarField sigmaR ( - IOobject - ( - "sigmaTheta", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimForce/dimArea, 0.0) - ); - - forAll(sigmaTheta.internalField(), celli) - { - const scalar& r = radii.internalField()[celli]; - - sigmaTheta.internalField()[celli] = - ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) * - (1 -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 + (b*b)/(r*r))*Foam::log(b/a) ); - } - - forAll(sigmaTheta.boundaryField(), patchi) - { - forAll(sigmaTheta.boundaryField()[patchi], facei) - { - const scalar& r = radii.boundaryField()[patchi][facei]; - - sigmaTheta.boundaryField()[patchi][facei] = - ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) * - (1 -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 + (b*b)/(r*r))*Foam::log(b/a) ); - } - } + IOobject + ( + "sigmaR", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))* + ( + -Foam::log(b/radii) + - (sqr(a)/(sqr(b) - sqr(a)))*(1 - sqr(b)/sqr(radii))*Foam::log(b/a) + ) + ); + sigmaR.write(); - //- write temperature file - Info << "\nWriting analytical sigmaTheta field" << endl; - sigmaTheta.write(); - - volScalarField sigmaZ + Info << "\nWriting analytical sigmaTheta field" << endl; + volScalarField sigmaTheta ( - IOobject - ( - "sigmaZ", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("zero", dimForce/dimArea, 0.0) - ); + IOobject + ( + "sigmaTheta", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))* + ( + 1 - Foam::log(b/radii) + - (sqr(a)/(sqr(b) - sqr(a)))*(1 + sqr(b)/sqr(radii))*Foam::log(b/a) + ) + ); + sigmaTheta.write(); - forAll(sigmaZ.internalField(), celli) - { - //- Timoshenko says this but I am not sure I am not sure the BCs in - //- the z direction - // sigmaZ.internalField()[celli] = - // ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) * - // (1 - 2*Foam::log(b/r) - ( 2*a*a/(b*b - a*a))*Foam::log(b/a)); - - sigmaZ.internalField()[celli] = - 0.3*(sigmaR.internalField()[celli] + sigmaTheta.internalField()[celli]) - - E*alpha*(T.internalField()[celli]); - } - - forAll(sigmaZ.boundaryField(), patchi) - { - forAll(sigmaZ.boundaryField()[patchi], facei) - { - //- Timoshenko says this but I am not sure I am not sure the BCs in - //- the z direction - //sigmaZ.boundaryField()[patchi][facei] = - //( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) * - //(1 - 2*Foam::log(b/r) - ( 2*a*a/(b*b - a*a))*Foam::log(b/a)); - - //-for general 2-D plain strain problems, the axial stress is given by this: - sigmaZ.boundaryField()[patchi][facei] = - nu*(sigmaR.boundaryField()[patchi][facei] + sigmaTheta.boundaryField()[patchi][facei]) - - E*alpha*(T.boundaryField()[patchi][facei]); - } - } - - - //- write temperature file - Info << "\nWriting analytical sigmaZ field" << endl; - sigmaZ.write(); - - - //- create analytical sigma tensor - - //- create theta field - volScalarField theta + Info << "\nWriting analytical sigmaZ field" << endl; + volScalarField sigmaZ ( - IOobject - ( - "theta", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimless, 0.0) - ); - forAll(theta.internalField(), celli) - { - const scalar& x = mesh.C().internalField()[celli][vector::X]; - const scalar& y = mesh.C().internalField()[celli][vector::Y]; + IOobject + ( + "sigmaZ", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + // Timoshenko says this but I am not sure I am not sure the BCs in + // the z direction + // ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))* + // (1 - 2*Foam::log(b/radii) - ( 2*sqr(a)/(sqr(b) - sqr(a)))*Foam::log(b/a)); + 0.3*(sigmaR + sigmaTheta) - E*alpha*(T) + ); + sigmaZ.write(); - theta.internalField()[celli] = Foam::atan(y/x); - } - - forAll(theta.boundaryField(), patchi) - { - forAll(theta.boundaryField()[patchi], facei) - { - const scalar& x = mesh.C().boundaryField()[patchi][facei][vector::X]; - const scalar& y = mesh.C().boundaryField()[patchi][facei][vector::Y]; - - theta.boundaryField()[patchi][facei] = Foam::atan(y/x); - } - } - - //- rotation matrix to convert polar stresses to cartesian - volTensorField rotMat + //- create theta field + volScalarField yOverX ( - IOobject - ( - "rotMat", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedTensor("zero", dimless, tensor::zero) - ); + "yOverX", + Foam::max + ( + scalar(-1), + Foam::min + ( + scalar(1), + mesh.C().component(vector::Y)/ + stabilise + ( + mesh.C().component(vector::X), + dimensionedScalar("small", dimLength, SMALL) + ) + ) + ) + ); - forAll(rotMat.internalField(), celli) - { - const scalar& t = theta.internalField()[celli]; - - rotMat.internalField()[celli] = tensor(::cos(t), ::sin(t), 0, - -::sin(t), ::cos(t), 0, - 0, 0, 1); - } - - forAll(rotMat.boundaryField(), patchi) - { - forAll(rotMat.boundaryField()[patchi], facei) - { - const scalar& t = theta.boundaryField()[patchi][facei]; - - rotMat.boundaryField()[patchi][facei] = tensor(::cos(t), ::sin(t), 0, - -::sin(t), ::cos(t), 0, - 0, 0, 1); - } - } - - volSymmTensorField sigma + volScalarField theta ( - IOobject - ( - "analyticalSigma", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) - ); + IOobject + ( + "theta", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + Foam::atan(yOverX) + ); - forAll(sigma.internalField(), celli) + //- rotation matrix to convert polar stresses to cartesian + volTensorField rotMat + ( + IOobject + ( + "rotMat", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedTensor("zero", dimless, tensor::zero) + ); + + tensorField& rotMatIn = rotMat.internalField(); + const scalarField tIn = theta.internalField(); + + forAll (rotMatIn, celli) { - const scalar& r = sigmaR.internalField()[celli]; - const scalar& t = sigmaTheta.internalField()[celli]; - const scalar& z = sigmaZ.internalField()[celli]; + const scalar& t = tIn[celli]; - const tensor& rot = rotMat.internalField()[celli]; - - symmTensor sigmaCart(r, 0, 0, - t, 0, - z); - - sigma.internalField()[celli] = - symm(rot.T() & sigmaCart & rot); - - //-for general 2-D plain strain problems, the axial stress is given by this: - //- (which is not equal to the solution by Timoshenko... hmmmnn) -// sigma.internalField()[celli][symmTensor::ZZ] = -// 0.3*(sigma.internalField()[celli][symmTensor::XX] + sigma.internalField()[celli][symmTensor::YY]) -// - E*alpha*(T.internalField()[celli]); - } - - forAll(sigma.boundaryField(), patchi) - { - forAll(sigma.boundaryField()[patchi], facei) - { - const scalar& r = sigmaR.boundaryField()[patchi][facei]; - const scalar& t = sigmaTheta.boundaryField()[patchi][facei]; - const scalar& z = sigmaZ.boundaryField()[patchi][facei]; - - const tensor& rot = rotMat.boundaryField()[patchi][facei]; - - symmTensor sigmaCart(r, 0, 0, - t, 0, - z); - sigma.boundaryField()[patchi][facei] = - symm(rot.T() & sigmaCart & rot); - } + rotMatIn[celli] = + tensor + ( + Foam::cos(t), Foam::sin(t), 0, + -Foam::sin(t), Foam::cos(t), 0, + 0, 0, 1 + ); } + forAll (rotMat.boundaryField(), patchi) + { + forAll (rotMat.boundaryField()[patchi], facei) + { + const scalar& t = theta.boundaryField()[patchi][facei]; - Info << "\nWriting analytical sigma tensor" << endl; - sigma.write(); + rotMat.boundaryField()[patchi][facei] = + tensor + ( + Foam::cos(t), Foam::sin(t), 0, + -Foam::sin(t), Foam::cos(t), 0, + 0, 0, 1 + ); + } + } - Info << nl << "End" << endl; - - return 0; + volSymmTensorField sigma + ( + IOobject + ( + "analyticalSigma", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero) + ); + + { + symmTensorField& sigmaIn = sigma.internalField(); + + const scalarField& rIn = sigmaR.internalField(); + const scalarField& tIn = sigmaTheta.internalField(); + const scalarField& zIn = sigmaZ.internalField(); + + forAll (sigmaIn, celli) + { + symmTensor sigmaCart + ( + rIn[celli], 0, 0, + tIn[celli], 0, + zIn[celli] + ); + + const tensor& rot = rotMatIn[celli]; + + sigmaIn[celli] = symm(rot.T() & sigmaCart & rot); + + // for general 2-D plain strain problems, the axial stress is: + // (which is not equal to the solution by Timoshenko... hmmmnn) + // sigmaIn[celli][symmTensor::ZZ] = + // 0.3*(sigmaIn[celli][symmTensor::XX] + // + sigmaIn[celli][symmTensor::YY]) + // - E*alpha*(T.internalField()[celli]); + } + } + + forAll (sigma.boundaryField(), patchi) + { + symmTensorField& pSigma = sigma.boundaryField()[patchi]; + const scalarField& pR = sigmaR.boundaryField()[patchi]; + const scalarField& pT = sigmaTheta.boundaryField()[patchi]; + const scalarField& pZ = sigmaZ.boundaryField()[patchi]; + + const tensorField pRot = rotMat.boundaryField()[patchi]; + + forAll (pSigma, facei) + { + const tensor& rot = pRot[facei]; + + symmTensor sigmaCart + ( + pR[facei], 0, 0, + pT[facei], 0, + pZ[facei] + ); + + pSigma[facei] = symm(rot.T() & sigmaCart & rot); + } + } + + Info << "\nWriting analytical sigma tensor" << endl; + sigma.write(); + + Info << nl << "End" << endl; + + return 0; }