tutorials/solidMechanics/elasticSolidFoam/plateHole/analyticalPlateHole/analyticalPlateHole.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Generate analytical solution for a infinite plaste with a circular
  26     hole.
  27     Stress field sigma is generated.
  28     Based on solution outlined in Timoshenko, Theory of Elasticity.
  29
  30 Author
  31     plateHoleSolution function by Z. Tukovic
  32     utility assembled by P. Cardiff
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "fvCFD.H"
  37 #include "volFields.H"
  38 #include "fvc.H"
  39 #include "fixedValueFvPatchFields.H"
  40 #include "coordinateSystem.H"
  41
  42 symmTensor plateHoleSolution(const vector& C);
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 int main(int argc, char *argv[])
  47 {
  48 #   include "setRootCase.H"
  49 #   include "createTime.H"
  50 #   include "createMesh.H"
  51
  52     runTime++;
  53
  54     Info << "Writing analytical solution for an infinite plate with a circular hole,\nwhere"
  55         << "\n\tradius = 0.5"
  56         << "\n\tdistant traction = (10,000 0 0 )"
  57         << nl << endl;
  58
  59     volSymmTensorField sigma
  60     (
  61         IOobject
  62         (
  63             "sigmaAnalyticalCylin",
  64            runTime.timeName(),
  65            mesh,
  66            IOobject::NO_READ,
  67            IOobject::AUTO_WRITE
  68         ),
  69         mesh,
  70         dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
  71     );
  72
  73     const volVectorField& C = mesh.C();
  74
  75     forAll(sigma.internalField(), celli)
  76     {
  77         vector curR = vector(C[celli].x(), C[celli].y(), 0);
  78
  79         sigma.internalField()[celli] = plateHoleSolution(curR);
  80     }
  81
  82     forAll(sigma.boundaryField(), patchi)
  83     {
  84         forAll(sigma.boundaryField()[patchi], facei)
  85         {
  86             vector curR = vector(C.boundaryField()[patchi][facei].x(), C.boundaryField()[patchi][facei].y(), 0);
  87
  88             sigma.boundaryField()[patchi][facei] = plateHoleSolution(curR);
  89         }
  90     }
  91
  92     Info << "Writing analytical sigma tensor" << endl;
  93     sigma.write();
  94
  95     Info << nl << "End" << endl;
  96
  97     return 0;
  98 }
  99
 100 // ************************************************************************* //
 101
 102
 103 symmTensor plateHoleSolution(const vector& C)
 104 {
 105     tensor sigma = tensor::zero;
 106
 107     scalar T = 10000;
 108     scalar a = 0.5;
 109
 110     scalar r = ::sqrt(sqr(C.x()) + sqr(C.y()));
 111     scalar theta = Foam::atan2(C.y(), C.x());
 112
 113     coordinateSystem cs("polarCS", C, vector(0, 0, 1), C/mag(C));
 114
 115     sigma.xx() =
 116         T*(1 - sqr(a)/sqr(r))/2
 117       + T*(1 + 3*pow(a,4)/pow(r,4) - 4*sqr(a)/sqr(r))*::cos(2*theta)/2;
 118
 119     sigma.xy() =
 120       - T*(1 - 3*pow(a,4)/pow(r,4) + 2*sqr(a)/sqr(r))*::sin(2*theta)/2;
 121
 122     sigma.yx() = sigma.xy();
 123
 124     sigma.yy() =
 125         T*(1 + sqr(a)/sqr(r))/2
 126       - T*(1 + 3*pow(a,4)/pow(r,4))*::cos(2*theta)/2;
 127
 128
 129     // Transformation to global coordinate system
 130     sigma = ((cs.R()&sigma)&cs.R().T());
 131
 132     symmTensor S = symmTensor::zero;
 133
 134     S.xx() = sigma.xx();
 135     S.xy() = sigma.xy();
 136     S.yy() = sigma.yy();
 137
 138     return S;
 139 }