applications/solvers/solidMechanics/utilities/forceDisp/forceDisp.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     Calculates the force and average displacement of the specified patch
  26     and writes it to a file
  27
  28 Author
  29     philip.cardiff@ucd.ie
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "fvCFD.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 int main(int argc, char *argv[])
  38 {
  39   Foam::argList::validOptions.insert("noMeshUpdate", "");
  40   Foam::argList::validOptions.insert("nonLinear", "");
  41   argList::validArgs.append("patchName");
  42
  43 # include "addTimeOptions.H"
  44
  45 # include "setRootCase.H"
  46
  47 # include "createTime.H"
  48
  49   word patchName(args.additionalArgs()[0]);
  50
  51   // Get times list
  52   instantList Times = runTime.times();
  53
  54   // set startTime and endTime depending on -time and -latestTime options
  55 # include "checkTimeOptions.H"
  56
  57   runTime.setTime(Times[startTime], startTime);
  58
  59 # include "createMesh.H"
  60
  61   bool noMeshUpdate = args.optionFound("noMeshUpdate");
  62   bool nonLinear = args.optionFound("nonLinear");
  63
  64   // check patch exists
  65   label patchID = mesh.boundaryMesh().findPatchID(patchName);
  66   if (patchID == -1)
  67     {
  68       FatalError
  69           << "Cannot find patch " << patchName << exit(FatalError);
  70     }
  71
  72   // open file
  73   OFstream forceDispFile("forceDisp_"+patchName+".dat");
  74   label width = 20;
  75   forceDispFile << "average disp";
  76   forceDispFile.width(width);
  77   forceDispFile << "totalForce";
  78   forceDispFile << endl;
  79
  80   for (label i=startTime; i<endTime; i++)
  81   {
  82       runTime.setTime(Times[i], i);
  83
  84       Info<< "Time = " << runTime.timeName() << endl;
  85
  86       if (!noMeshUpdate)
  87       {
  88           mesh.readUpdate();
  89       }
  90
  91       IOobject sigmaheader
  92           (
  93               "sigma",
  94               runTime.timeName(),
  95               mesh,
  96               IOobject::MUST_READ
  97               );
  98
  99       IOobject Uheader
 100           (
 101               "U",
 102               runTime.timeName(),
 103               mesh,
 104               IOobject::MUST_READ
 105               );
 106
 107       // Check sigma exists
 108       if (sigmaheader.headerOk() && Uheader.headerOk())
 109       {
 110           Info<< "    Reading sigma" << endl;
 111           volSymmTensorField sigma(sigmaheader, mesh);
 112           Info<< "    Reading U" << endl;
 113           volVectorField U(Uheader, mesh);
 114
 115           Info << nl;
 116
 117           // calculate patch force
 118           const vectorField n = mesh.boundary()[patchID].nf();
 119           const vectorField& Sf = mesh.boundary()[patchID].Sf();
 120           const symmTensorField& sigmaPatch = sigma.boundaryField()[patchID];
 121           vector totalForce = vector::zero;
 122           scalar normalForce = 0.0;
 123           scalar shearForce = 0.0;
 124           if (nonLinear)
 125           {
 126               // assuming sigma is the second Piola-Kirchhoff tensor
 127
 128               // get deformation gradient
 129               volTensorField gradU
 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               tensorField F = I + gradU.boundaryField()[patchID];
 142               vectorField force = Sf & (sigmaPatch & F);
 143               tensorField Finv = inv(F);
 144               vectorField nCur = Finv & n;
 145               nCur /= mag(nCur);
 146               normalForce = sum(nCur & force);
 147               shearForce = sum( mag((I - sqr(nCur)) & force) );
 148               totalForce = sum(force);
 149           }
 150           else
 151           {
 152               // small strain or UL large strain
 153               // (as the mesh is moved and sigma is Cauchy)
 154               //Info << "\tSmall Strain Total Lagrangian"<<nl<<endl;
 155               vectorField force = Sf & sigmaPatch;
 156               normalForce = sum(n & force);
 157               shearForce = sum( mag((I - sqr(n)) & force) );
 158               totalForce = sum(force);
 159           }
 160
 161           // calculate average displacement
 162           // these should be a weighted average but OK for most models
 163           vector disp = average(U.boundaryField()[patchID]);
 164
 165           // write to file
 166           forceDispFile
 167               << disp.x() << " " << disp.y() << " " << disp.z();
 168           forceDispFile.width(width);
 169           forceDispFile
 170               << totalForce.x() << " " << totalForce.y()
 171               << " " << totalForce.z();
 172           forceDispFile
 173               << " " << normalForce << " " << shearForce
 174               << endl;
 175       }
 176   }
 177
 178   Info << nl << "End" << endl;
 179
 180   return 0;
 181 }
 182
 183
 184 // ************************************************************************* //