applications/solvers/solidMechanics/utilities/surfaceTractions/surfaceTractions.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 Application
  25     surfaceTractions
  26
  27 Description
  28     Calculates and writes the surface tractions as a volVectorField, using
  29     the sigma volSymmTensorField
  30
  31 Author
  32     Philip Cardiff UCD
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "fvCFD.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 int main(int argc, char *argv[])
  41 {
  42   Foam::argList::validOptions.insert("nonLinear", "");
  43
  44 # include "addTimeOptions.H"
  45 # include "setRootCase.H"
  46 # include "createTime.H"
  47
  48   bool nonLinear = args.optionFound("nonLinear");
  49
  50   // Get times list
  51   instantList Times = runTime.times();
  52
  53   // set startTime and endTime depending on -time and -latestTime options
  54 # include "checkTimeOptions.H"
  55
  56   runTime.setTime(Times[startTime], startTime);
  57
  58 # include "createMesh.H"
  59
  60   for (label i=startTime; i<endTime; i++)
  61   {
  62       runTime.setTime(Times[i], i);
  63
  64       Info<< "Time = " << runTime.timeName() << endl;
  65
  66       mesh.readUpdate();
  67
  68       IOobject sigmaheader
  69           (
  70               "sigma",
  71               runTime.timeName(),
  72               mesh,
  73               IOobject::MUST_READ
  74               );
  75
  76       // Check sigma exists
  77       if (sigmaheader.headerOk())
  78       {
  79           mesh.readUpdate();
  80
  81           Info<< "    Reading sigma" << endl;
  82           volSymmTensorField sigma(sigmaheader, mesh);
  83
  84           surfaceVectorField n = mesh.Sf()/mesh.magSf();
  85
  86           volVectorField totalTraction
  87               (
  88                   IOobject
  89                   (
  90                       "totalTraction",
  91                       runTime.timeName(),
  92                       mesh,
  93                       IOobject::NO_READ,
  94                       IOobject::AUTO_WRITE
  95                       ),
  96                   mesh,
  97                   dimensionedVector("zero", dimForce/dimArea, vector::zero)
  98                   );
  99           volScalarField normalTraction
 100               (
 101                   IOobject
 102                   (
 103                       "normalTraction",
 104                       runTime.timeName(),
 105                       mesh,
 106                       IOobject::NO_READ,
 107                       IOobject::AUTO_WRITE
 108                       ),
 109                   mesh,
 110                   dimensionedScalar("zero", dimForce/dimArea, 0.0)
 111                   );
 112           volVectorField shearTraction
 113               (
 114                   IOobject
 115                   (
 116                       "shearTraction",
 117                       runTime.timeName(),
 118                       mesh,
 119                       IOobject::NO_READ,
 120                       IOobject::AUTO_WRITE
 121                       ),
 122                   mesh,
 123                   dimensionedVector("zero", dimForce/dimArea, vector::zero)
 124                   );
 125
 126           volTensorField* gradUPtr = NULL;
 127           if (nonLinear)
 128           {
 129               gradUPtr = new volTensorField
 130                   (
 131                       IOobject
 132                       (
 133                           "grad(U)",
 134                           runTime.timeName(),
 135                           mesh,
 136                           IOobject::MUST_READ,
 137                           IOobject::NO_WRITE
 138                           ),
 139                       mesh
 140                       );
 141           }
 142
 143           forAll(totalTraction.boundaryField(), patchi)
 144           {
 145               const vectorField& nb = n.boundaryField()[patchi];
 146               const symmTensorField& sigmab = sigma.boundaryField()[patchi];
 147
 148               if (nonLinear)
 149               {
 150                   tensorField F = I + gradUPtr->boundaryField()[patchi];
 151                   totalTraction.boundaryField()[patchi] = nb & (sigmab & F);
 152               }
 153               else
 154               {
 155                   totalTraction.boundaryField()[patchi] = nb & sigmab;
 156               }
 157               normalTraction.boundaryField()[patchi] =
 158                   nb & totalTraction.boundaryField()[patchi];
 159               shearTraction.boundaryField()[patchi] =
 160                   (I -sqr(nb)) & totalTraction.boundaryField()[patchi];
 161           }
 162           totalTraction.write();
 163           normalTraction.write();
 164           shearTraction.write();
 165       }
 166       else
 167       {
 168           Info<< "    No sigma field" << endl;
 169       }
 170       Info<< endl;
 171   }
 172
 173   Info<< "End" << endl;
 174
 175   return(0);
 176 }
 177
 178
 179 // ************************************************************************* //