src/thermophysicalModels/basicSolidThermo/directionalKSolidThermo/directionalKSolidThermo.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "directionalKSolidThermo.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "transform.H"
  29 #include "transformField.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35     defineTypeNameAndDebug(directionalKSolidThermo, 0);
  36     addToRunTimeSelectionTable
  37     (
  38         basicSolidThermo,
  39         directionalKSolidThermo,
  40         mesh
  41     );
  42
  43     addToRunTimeSelectionTable
  44     (
  45         basicSolidThermo,
  46         directionalKSolidThermo,
  47         dictionary
  48     );
  49
  50 }
  51
  52
  53 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  54
  55 Foam::directionalKSolidThermo::directionalKSolidThermo
  56 (
  57     const fvMesh& mesh,
  58     const dictionary& dict
  59 )
  60 :
  61     interpolatedSolidThermo(mesh, typeName + "Coeffs", dict),
  62     directionalK_
  63     (
  64         IOobject
  65         (
  66             "K",
  67             mesh_.time().timeName(),
  68             mesh_,
  69             IOobject::NO_READ,
  70             IOobject::NO_WRITE
  71         ),
  72         mesh_,
  73         dimEnergy/dimTime/(dimLength*dimTemperature)
  74     ),
  75     ccTransforms_
  76     (
  77         IOobject
  78         (
  79             "ccTransforms",
  80             mesh.time().timeName(),
  81             mesh,
  82             IOobject::NO_READ,
  83             IOobject::NO_WRITE
  84         ),
  85         mesh,
  86         dimLength
  87     )
  88 {
  89     init();
  90 }
  91
  92
  93 Foam::directionalKSolidThermo::directionalKSolidThermo(const fvMesh& mesh)
  94 :
  95     interpolatedSolidThermo(mesh, typeName + "Coeffs"),
  96     directionalK_
  97     (
  98         IOobject
  99         (
 100             "K",
 101             mesh_.time().timeName(),
 102             mesh_,
 103             IOobject::NO_READ,
 104             IOobject::NO_WRITE
 105         ),
 106         mesh_,
 107         dimEnergy/dimTime/(dimLength*dimTemperature)
 108     ),
 109     ccTransforms_
 110     (
 111         IOobject
 112         (
 113             "ccTransforms",
 114             mesh.time().timeName(),
 115             mesh,
 116             IOobject::NO_READ,
 117             IOobject::NO_WRITE
 118         ),
 119         mesh,
 120         dimLength
 121     )
 122 {
 123     init();
 124 }
 125
 126
 127 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 128
 129 Foam::directionalKSolidThermo::~directionalKSolidThermo()
 130 {}
 131
 132
 133 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 134
 135 void Foam::directionalKSolidThermo::init()
 136 {
 137     KValues_ = Field<vector>(subDict(typeName + "Coeffs").lookup("KValues"));
 138
 139     const fvMesh& mesh = directionalK_.mesh();
 140
 141     // Determine transforms for cell centres
 142     forAll(mesh.C(), cellI)
 143     {
 144         vector dir = mesh.C()[cellI] - coordSys_.origin();
 145         dir /= mag(dir);
 146
 147         // Define local coordinate system with
 148         // - e1 : axis from cc to centre
 149         // - e3 : rotation axis
 150         coordinateSystem cs
 151         (
 152             "cc",
 153             coordSys_.origin(),
 154             coordSys_.e3(),     //z',e3
 155             dir                 //x',e1
 156         );
 157
 158         ccTransforms_[cellI] = cs.R();
 159     }
 160
 161     forAll(mesh.C().boundaryField(), patchI)
 162     {
 163         const fvPatchVectorField& patchC = mesh.C().boundaryField()[patchI];
 164         fvPatchTensorField& patchT = ccTransforms_.boundaryField()[patchI];
 165
 166         tensorField tc(patchT.size());
 167         forAll(tc, i)
 168         {
 169             vector dir = patchC[i] - coordSys_.origin();
 170             dir /= mag(dir);
 171
 172             coordinateSystem cs
 173             (
 174                 "cc",
 175                 coordSys_.origin(),
 176                 coordSys_.e3(),     //z',e3
 177                 dir                 //x',e1
 178             );
 179
 180             tc[i] = cs.R();
 181         }
 182         patchT = tc;
 183     }
 184
 185     if (debug)
 186     {
 187         Info<< "directionalKSolidThermo : dumping converted Kxx, Kyy, Kzz"
 188             << endl;
 189         {
 190             volVectorField Kxx
 191             (
 192                 IOobject
 193                 (
 194                     "Kxx",
 195                     mesh.time().timeName(),
 196                     mesh,
 197                     IOobject::NO_READ,
 198                     IOobject::AUTO_WRITE,
 199                     false
 200                 ),
 201                 mesh,
 202                 dimless
 203             );
 204             Kxx.internalField() = transform
 205             (
 206                 ccTransforms_.internalField(),
 207                 vectorField
 208                 (
 209                     ccTransforms_.internalField().size(),
 210                     point(1, 0, 0)
 211                 )
 212             );
 213             forAll(Kxx.boundaryField(), patchI)
 214             {
 215                 Kxx.boundaryField()[patchI] = transform
 216                 (
 217                     ccTransforms_.boundaryField()[patchI],
 218                     vectorField
 219                     (
 220                         ccTransforms_.boundaryField()[patchI].size(),
 221                         point(1, 0, 0)
 222                     )
 223                 );
 224             }
 225             Kxx.write();
 226         }
 227         {
 228             volVectorField Kyy
 229             (
 230                 IOobject
 231                 (
 232                     "Kyy",
 233                     mesh.time().timeName(),
 234                     mesh,
 235                     IOobject::NO_READ,
 236                     IOobject::AUTO_WRITE,
 237                     false
 238                 ),
 239                 mesh,
 240                 dimless
 241             );
 242             Kyy.internalField() = transform
 243             (
 244                 ccTransforms_.internalField(),
 245                 vectorField
 246                 (
 247                     ccTransforms_.internalField().size(),
 248                     point(0, 1, 0)
 249                 )
 250             );
 251             forAll(Kyy.boundaryField(), patchI)
 252             {
 253                 Kyy.boundaryField()[patchI] = transform
 254                 (
 255                     ccTransforms_.boundaryField()[patchI],
 256                     vectorField
 257                     (
 258                         ccTransforms_.boundaryField()[patchI].size(),
 259                         point(0, 1, 0)
 260                     )
 261                 );
 262             }
 263             Kyy.write();
 264         }
 265         {
 266             volVectorField Kzz
 267             (
 268                 IOobject
 269                 (
 270                     "Kzz",
 271                     mesh.time().timeName(),
 272                     mesh,
 273                     IOobject::NO_READ,
 274                     IOobject::AUTO_WRITE,
 275                     false
 276                 ),
 277                 mesh,
 278                 dimless
 279             );
 280             Kzz.internalField() = transform
 281             (
 282                 ccTransforms_.internalField(),
 283                 vectorField
 284                 (
 285                     ccTransforms_.internalField().size(),
 286                     point(0, 0, 1)
 287                 )
 288             );
 289             forAll(Kzz.boundaryField(), patchI)
 290             {
 291                 Kzz.boundaryField()[patchI] = transform
 292                 (
 293                     ccTransforms_.boundaryField()[patchI],
 294                     vectorField
 295                     (
 296                         ccTransforms_.boundaryField()[patchI].size(),
 297                         point(0, 0, 1)
 298                     )
 299                 );
 300             }
 301             Kzz.write();
 302         }
 303     }
 304
 305     correct();
 306 }
 307
 308
 309 Foam::symmTensor Foam::directionalKSolidThermo::transformPrincipal
 310 (
 311     const tensor& tt,
 312     const vector& st
 313 ) const
 314 {
 315     return symmTensor
 316     (
 317         tt.xx()*st.x()*tt.xx()
 318       + tt.xy()*st.y()*tt.xy()
 319       + tt.xz()*st.z()*tt.xz(),
 320
 321         tt.xx()*st.x()*tt.yx()
 322       + tt.xy()*st.y()*tt.yy()
 323       + tt.xz()*st.z()*tt.yz(),
 324
 325         tt.xx()*st.x()*tt.zx()
 326       + tt.xy()*st.y()*tt.zy()
 327       + tt.xz()*st.z()*tt.zz(),
 328
 329         tt.yx()*st.x()*tt.yx()
 330       + tt.yy()*st.y()*tt.yy()
 331       + tt.yz()*st.z()*tt.yz(),
 332
 333         tt.yx()*st.x()*tt.zx()
 334       + tt.yy()*st.y()*tt.zy()
 335       + tt.yz()*st.z()*tt.zz(),
 336
 337         tt.zx()*st.x()*tt.zx()
 338       + tt.zy()*st.y()*tt.zy()
 339       + tt.zz()*st.z()*tt.zz()
 340     );
 341 }
 342
 343
 344 void Foam::directionalKSolidThermo::transformField
 345 (
 346     symmTensorField& fld,
 347     const tensorField& tt,
 348     const vectorField& st
 349 ) const
 350 {
 351     fld.setSize(tt.size());
 352     forAll(fld, i)
 353     {
 354         fld[i] = transformPrincipal(tt[i], st[i]);
 355     }
 356 }
 357
 358
 359 void Foam::directionalKSolidThermo::correct()
 360 {
 361     calculate();
 362     interpolatedSolidThermo::calculate();
 363 }
 364
 365
 366 const Foam::volSymmTensorField&
 367 Foam::directionalKSolidThermo::directionalK() const
 368 {
 369     return directionalK_;
 370 }
 371
 372
 373 void Foam::directionalKSolidThermo::calculate()
 374 {
 375     // Correct directionalK
 376     Field<vector> localK
 377     (
 378         interpolateXY
 379         (
 380             T_.internalField(),
 381             TValues_,
 382             KValues_
 383         )
 384     );
 385
 386     // Transform into global coordinate system
 387     transformField
 388     (
 389         directionalK_.internalField(),
 390         ccTransforms_.internalField(),
 391         localK
 392     );
 393
 394     forAll(directionalK_.boundaryField(), patchI)
 395     {
 396         directionalK_.boundaryField()[patchI] == this->directionalK(patchI)();
 397     }
 398 }
 399
 400
 401 const Foam::volScalarField& Foam::directionalKSolidThermo::K() const
 402 {
 403     forAll(KValues_, i)
 404     {
 405         const vector& v = KValues_[i];
 406         if
 407         (
 408             v.x() != v.y()
 409          || v.x() != v.z()
 410          || v.y() != v.z()
 411         )
 412         {
 413             FatalErrorIn("directionalKSolidThermo::K() const")
 414                 << "Supplied K values " << KValues_
 415                 << " are not isotropic." << exit(FatalError);
 416         }
 417     }
 418
 419     // Get temperature interpolated properties (principal directions)
 420     Field<vector> localK
 421     (
 422         interpolateXY
 423         (
 424             T_.internalField(),
 425             TValues_,
 426             KValues_
 427         )
 428     );
 429
 430     tmp<volScalarField> tK
 431     (
 432         new volScalarField
 433         (
 434             IOobject
 435             (
 436                 "K",
 437                 mesh_.time().timeName(),
 438                 mesh_,
 439                 IOobject::NO_READ,
 440                 IOobject::NO_WRITE
 441             ),
 442             mesh_,
 443             dimEnergy/dimTime/(dimLength*dimTemperature)
 444         )
 445     );
 446     volScalarField& K = tK();
 447
 448     K.internalField() = interpolateXY
 449     (
 450         T_.internalField(),
 451         TValues_,
 452         KValues_.component(0)()
 453     );
 454
 455     forAll(K.boundaryField(), patchI)
 456     {
 457         K.boundaryField()[patchI] == this->K(patchI)();
 458     }
 459
 460     return tK;
 461 }
 462
 463
 464 Foam::tmp<Foam::scalarField> Foam::directionalKSolidThermo::K
 465 (
 466     const label patchI
 467 ) const
 468 {
 469     forAll(KValues_, i)
 470     {
 471         const vector& v = KValues_[i];
 472         if
 473         (
 474             v.x() != v.y()
 475          || v.x() != v.z()
 476          || v.y() != v.z()
 477         )
 478         {
 479             FatalErrorIn("directionalKSolidThermo::K() const")
 480                 << "Supplied K values " << KValues_
 481                 << " are not isotropic." << exit(FatalError);
 482         }
 483     }
 484
 485     return tmp<scalarField>
 486     (
 487         new scalarField
 488         (
 489             interpolateXY
 490             (
 491                 T_.boundaryField()[patchI],
 492                 TValues_,
 493                 KValues_.component(0)()
 494             )
 495         )
 496     );
 497 }
 498
 499
 500 Foam::tmp<Foam::symmTensorField> Foam::directionalKSolidThermo::directionalK
 501 (
 502     const label patchI
 503 ) const
 504 {
 505     const fvPatchScalarField& patchT = T_.boundaryField()[patchI];
 506
 507     Field<vector> localK(interpolateXY(patchT, TValues_, KValues_));
 508
 509     tmp<symmTensorField> tglobalK(new symmTensorField(localK.size()));
 510     transformField(tglobalK(), ccTransforms_.boundaryField()[patchI], localK);
 511
 512     return tglobalK;
 513 }
 514
 515
 516 bool Foam::directionalKSolidThermo::read()
 517 {
 518     return read(subDict(typeName + "Coeffs"));
 519 }
 520
 521
 522 bool Foam::directionalKSolidThermo::read(const dictionary& dict)
 523 {
 524     coordSys_ = coordinateSystem(dict, mesh_);
 525     KValues_  = Field<vector>(subDict(typeName + "Coeffs").lookup("KValues"));
 526     return true;
 527 }
 528
 529
 530 bool Foam::directionalKSolidThermo::writeData(Ostream& os) const
 531 {
 532     bool ok = interpolatedSolidThermo::writeData(os);
 533     os.writeKeyword("KValues") << KValues_ << token::END_STATEMENT << nl;
 534     return ok && os.good();
 535 }
 536
 537
 538 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 539
 540 Foam::Ostream& Foam::operator<<(Ostream& os, const directionalKSolidThermo& s)
 541 {
 542     s.writeData(os);
 543     return os;
 544 }
 545
 546
 547 // ************************************************************************* //