src/solidModels/constitutiveModel/rheologyLaws/orthotropicLinearElastic/orthotropicLinearElastic.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "orthotropicLinearElastic.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "zeroGradientFvPatchFields.H"
  29 #include "transformField.H"
  30 #include "transformGeometricField.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTypeNameAndDebug(orthotropicLinearElastic, 0);
  37     addToRunTimeSelectionTable
  38     (rheologyLaw, orthotropicLinearElastic, dictionary);
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 // Construct from dictionary
  48 Foam::orthotropicLinearElastic::orthotropicLinearElastic
  49 (
  50     const word& name,
  51     const volSymmTensorField& sigma,
  52     const dictionary& dict
  53 )
  54 :
  55     rheologyLaw(name, sigma, dict),
  56     rho_(dict.lookup("rho")),
  57     Ex_(dict.lookup("Ex")),
  58     Ey_(dict.lookup("Ey")),
  59     Ez_(dict.lookup("Ez")),
  60     E_(Ex_),
  61     nuxy_(dict.lookup("nuxy")),
  62     nuyz_(dict.lookup("nuyz")),
  63     nuzx_(dict.lookup("nuzx")),
  64     nuyx_(nuxy_*Ey_/Ex_),
  65     nuzy_(nuyz_*Ez_/Ey_),
  66     nuxz_(nuzx_*Ex_/Ez_),
  67     nu_(nuxy_),
  68     Gxy_(dict.lookup("Gxy")),
  69     Gyz_(dict.lookup("Gyz")),
  70     Gzx_(dict.lookup("Gzx")),
  71     J_(
  72         (1 - nuxy_*nuyx_ - nuyz_*nuzy_ - nuzx_*nuxz_ - 2*nuyx_*nuzy_*nuxz_)
  73         / (Ex_*Ey_*Ez_)
  74         ),
  75     A11_( (1 - nuyz_*nuzy_)/(J_*Ey_*Ez_) ),
  76     A22_( (1 - nuxz_*nuzx_)/(J_*Ex_*Ez_) ),
  77     A33_( (1 - nuyx_*nuxy_)/(J_*Ey_*Ex_) ),
  78     A12_( (nuxy_ + nuzy_*nuxz_)/(J_*Ex_*Ez_) ),
  79     A31_( (nuzx_ + nuyx_*nuzy_)/(J_*Ey_*Ez_) ),
  80     A23_( (nuyz_ + nuyx_*nuxz_)/(J_*Ex_*Ey_) ),
  81     A44_( 2*Gxy_ ),
  82     A55_( 2*Gyz_ ),
  83     A66_( 2*Gzx_ ),
  84     C_(
  85        IOobject
  86        (
  87            "rheologyLawStoredC",
  88            mesh().time().timeName(),
  89            mesh(),
  90            IOobject::NO_READ,
  91            IOobject::NO_WRITE
  92            ),
  93        mesh(),
  94        dimensionedSymmTensor4thOrder
  95        (
  96            "zero",
  97            dimForce/dimArea,
  98            symmTensor4thOrder(A11_.value(), A12_.value(), A31_.value(),
  99                               A22_.value(), A23_.value(),
 100                               A33_.value(),
 101                               A44_.value(),
 102                               A55_.value(),
 103                               A66_.value())
 104            ),
 105        zeroGradientFvPatchSymmTensor4thOrderField::typeName
 106         ),
 107     matDir_(
 108         IOobject
 109         (
 110             "materialDirections",
 111          mesh().time().timeName(),
 112          mesh(),
 113          IOobject::MUST_READ,
 114          IOobject::NO_WRITE
 115          ),
 116         mesh()
 117         )
 118 {
 119   //- check material properties lie within physical constraints
 120   //- ref Abaqus analysis user's manual orthotropic material
 121     Info<< "\tChecking physical constraints on the orthotropic"
 122         << " material properties" << endl;
 123
 124   //- E and G should be greater than zero
 125   if (Ex_.value() < 0.0 || Ey_.value() < 0.0 || Ez_.value() < 0.0
 126      || Gxy_.value() < 0.0 || Gyz_.value() < 0.0 || Gzx_.value() < 0.0)
 127     {
 128       FatalError << "Ex, Ey, Ez, Gxy, Gyz, Gzx should all be greater than zero!"
 129          << exit(FatalError);
 130     }
 131
 132   //- restriction on Poisson's ratio
 133   if (mag(nuxy_.value()) >= sqrt(Ex_.value()/Ey_.value())
 134      || mag(nuyz_.value()) >= sqrt(Ey_.value()/Ez_.value())
 135      || mag(nuzx_.value()) >= sqrt(Ez_.value()/Ex_.value()))
 136     {
 137       FatalError << "mag(nuij) should be less sqrt(Ei/Ej)"
 138          << exit(FatalError);
 139     }
 140
 141   if (
 142       dimensionedScalar
 143       (1 - nuxy_*nuyx_ - nuyz_*nuzy_
 144        - nuzx_*nuxz_ - 2*nuyx_*nuzy_*nuxz_).value()
 145       <= 0 )
 146     {
 147       FatalError
 148           << "(1 - nuxy*nuyx - nuyz*nuzy "
 149           << "- nuzx*nuxz - 2*nuyx*nuzy*nuxz) should be greater than zero!"
 150          << exit(FatalError);
 151     }
 152
 153
 154   Info<< "\tRotating local material properties to global coordinate system"
 155       << endl;
 156   //- rotate tensors
 157   volTensorField R
 158     (
 159      IOobject
 160      (
 161       "R",
 162       mesh().time().timeName(),
 163       mesh(),
 164       IOobject::NO_READ,
 165       IOobject::NO_WRITE
 166       ),
 167      mesh(),
 168      dimensionedTensor("zero", dimless, tensor::zero),
 169      zeroGradientFvPatchTensorField::typeName
 170      );
 171
 172   //- make sure matDir_ are unit directions
 173   forAll(matDir_, celli)
 174     {
 175       {
 176     scalar magVec =
 177       mag(vector(matDir_[celli][0], matDir_[celli][1], matDir_[celli][2]));
 178     matDir_[celli][0] /= magVec;
 179     matDir_[celli][1] /= magVec;
 180     matDir_[celli][2] /= magVec;
 181       }
 182       {
 183     scalar magVec =
 184       mag(vector(matDir_[celli][3], matDir_[celli][4], matDir_[celli][5]));
 185     matDir_[celli][3] /= magVec;
 186     matDir_[celli][4] /= magVec;
 187     matDir_[celli][5] /= magVec;
 188       }
 189       {
 190     scalar magVec =
 191       mag(vector(matDir_[celli][6], matDir_[celli][7], matDir_[celli][8]));
 192     matDir_[celli][6] /= magVec;
 193     matDir_[celli][7] /= magVec;
 194     matDir_[celli][8] /= magVec;
 195       }
 196     }
 197
 198   //- global axes
 199   vector e0(1,0,0);
 200   vector e1(0,1,0);
 201   vector e2(0,0,1);
 202
 203   forAll(R, celli)
 204     {
 205       // R_ij = xold_i & xnew_i;
 206       {
 207     vector mD(matDir_[celli][0], matDir_[celli][1], matDir_[celli][2]);
 208     R[celli][0] = e0 & mD;
 209     R[celli][1] = e1 & mD;
 210     R[celli][2] = e2 & mD;
 211       }
 212       {
 213     vector mD(matDir_[celli][3], matDir_[celli][4], matDir_[celli][5]);
 214     R[celli][3] = e0 & mD;
 215     R[celli][4] = e1 & mD;
 216     R[celli][5] = e2 & mD;
 217       }
 218       {
 219     vector mD(matDir_[celli][6], matDir_[celli][7], matDir_[celli][8]);
 220     R[celli][6] = e0 & mD;
 221     R[celli][7] = e1 & mD;
 222     R[celli][8] = e2 & mD;
 223       }
 224     }
 225
 226   //Info << "R is " << R.internalField() << endl;
 227   R.correctBoundaryConditions();
 228   //R.write();
 229
 230   //- rotate C to global corrdinate system
 231   //C_.correctBoundaryConditions();
 232   C_ = transform(R, C_);
 233   C_.correctBoundaryConditions();
 234 }
 235
 236
 237 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 238
 239 Foam::orthotropicLinearElastic::~orthotropicLinearElastic()
 240 {}
 241
 242
 243 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 244
 245 Foam::tmp<Foam::volScalarField> Foam::orthotropicLinearElastic::rho() const
 246 {
 247     tmp<volScalarField> tresult
 248     (
 249         new volScalarField
 250         (
 251             IOobject
 252             (
 253                 "rho",
 254                 mesh().time().timeName(),
 255                 mesh(),
 256                 IOobject::NO_READ,
 257                 IOobject::NO_WRITE
 258             ),
 259             mesh(),
 260             rho_,
 261             zeroGradientFvPatchScalarField::typeName
 262         )
 263     );
 264
 265     tresult().correctBoundaryConditions();
 266
 267     return tresult;
 268 }
 269
 270
 271 Foam::tmp<Foam::volScalarField> Foam::orthotropicLinearElastic::E() const
 272 {
 273     tmp<volScalarField> tresult
 274     (
 275         new volScalarField
 276         (
 277             IOobject
 278             (
 279                 "E",
 280                 mesh().time().timeName(),
 281                 mesh(),
 282                 IOobject::NO_READ,
 283                 IOobject::NO_WRITE
 284             ),
 285             mesh(),
 286             E_,
 287             zeroGradientFvPatchScalarField::typeName
 288         )
 289     );
 290
 291     tresult().correctBoundaryConditions();
 292
 293     return tresult;
 294 }
 295
 296
 297 Foam::tmp<Foam::volScalarField> Foam::orthotropicLinearElastic::nu() const
 298 {
 299     tmp<volScalarField> tresult
 300     (
 301         new volScalarField
 302         (
 303             IOobject
 304             (
 305                 "nu",
 306                 mesh().time().timeName(),
 307                 mesh(),
 308                 IOobject::NO_READ,
 309                 IOobject::NO_WRITE
 310             ),
 311             mesh(),
 312             nu_,
 313             zeroGradientFvPatchScalarField::typeName
 314         )
 315     );
 316
 317     tresult().correctBoundaryConditions();
 318
 319     return tresult;
 320 }
 321
 322 Foam::tmp<Foam::volScalarField> Foam::orthotropicLinearElastic::Ep() const
 323 {
 324 //     notImplemented(type() + "::Ep()");
 325
 326     return tmp<volScalarField>
 327     (
 328         new volScalarField
 329         (
 330             IOobject
 331             (
 332                 "Ep",
 333                 mesh().time().timeName(),
 334                 mesh(),
 335                 IOobject::NO_READ,
 336                 IOobject::NO_WRITE
 337             ),
 338             mesh(),
 339             dimensionedScalar("zeroEp", dimForce/dimArea, 0.0),
 340             zeroGradientFvPatchScalarField::typeName
 341         )
 342     );
 343 }
 344
 345 Foam::tmp<Foam::volScalarField> Foam::orthotropicLinearElastic::sigmaY() const
 346 {
 347 //     notImplemented(type() + "::sigmaY()");
 348
 349     return tmp<volScalarField>
 350     (
 351         new volScalarField
 352         (
 353             IOobject
 354             (
 355                 "sigmaY",
 356                 mesh().time().timeName(),
 357                 mesh(),
 358                 IOobject::NO_READ,
 359                 IOobject::NO_WRITE
 360             ),
 361             mesh(),
 362             dimensionedScalar("zeroSigmaY", dimForce/dimArea, GREAT),
 363             zeroGradientFvPatchScalarField::typeName
 364         )
 365     );
 366 }
 367
 368 Foam::scalar Foam::orthotropicLinearElastic::sigmaY
 369 (const scalar epsilonPEq, const label cellID) const
 370 {
 371   return GREAT;
 372 }
 373
 374 //Foam::tmp<Foam::volTensorField> Foam::orthotropicLinearElastic::K() const
 375 Foam::tmp<Foam::volDiagTensorField> Foam::orthotropicLinearElastic::K() const
 376 {
 377   tmp<volDiagTensorField> tresult
 378     (
 379         new volDiagTensorField
 380         (
 381             IOobject
 382             (
 383                 "K",
 384                 mesh().time().timeName(),
 385                 mesh(),
 386                 IOobject::NO_READ,
 387                 IOobject::NO_WRITE
 388             ),
 389             mesh(),
 390             dimensionedDiagTensor("zero", A11_.dimensions(), diagTensor::zero),
 391             zeroGradientFvPatchScalarField::typeName
 392         )
 393     );
 394
 395   volDiagTensorField& K = tresult();
 396
 397   forAll(K, celli)
 398     {
 399       K[celli].xx() = C_[celli].xxxx();
 400       K[celli].yy() = C_[celli].yyyy();
 401       K[celli].zz() = C_[celli].zzzz();
 402     }
 403
 404   K.correctBoundaryConditions();
 405
 406   return tresult;
 407 }
 408
 409 Foam::tmp<Foam::volSymmTensor4thOrderField>
 410 Foam::orthotropicLinearElastic::C() const
 411 {
 412   tmp<volSymmTensor4thOrderField> tresult
 413     (
 414         new volSymmTensor4thOrderField
 415         (
 416      C_
 417      )
 418     );
 419
 420   volSymmTensor4thOrderField& result = tresult();
 421
 422   result.correctBoundaryConditions();
 423
 424   return tresult;
 425 }
 426 // ************************************************************************* //