src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "CoEulerDdtScheme.H"
  28 #include "surfaceInterpolate.H"
  29 #include "fvcDiv.H"
  30 #include "fvMatrices.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace fv
  40 {
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 template<class Type>
  45 tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
  46 {
  47     surfaceScalarField cofrDeltaT = CofrDeltaT();
  48
  49     tmp<volScalarField> tcorDeltaT
  50     (
  51         new volScalarField
  52         (
  53             IOobject
  54             (
  55                 "CorDeltaT",
  56                 cofrDeltaT.instance(),
  57                 mesh()
  58             ),
  59             mesh(),
  60             dimensionedScalar("CorDeltaT", cofrDeltaT.dimensions(), 0.0),
  61             zeroGradientFvPatchScalarField::typeName
  62         )
  63     );
  64
  65     volScalarField& corDeltaT = tcorDeltaT();
  66
  67     const unallocLabelList& owner = mesh().owner();
  68     const unallocLabelList& neighbour = mesh().neighbour();
  69
  70     forAll(owner, faceI)
  71     {
  72         corDeltaT[owner[faceI]] =
  73             max(corDeltaT[owner[faceI]], cofrDeltaT[faceI]);
  74
  75         corDeltaT[neighbour[faceI]] =
  76             max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
  77     }
  78
  79     forAll(corDeltaT.boundaryField(), patchi)
  80     {
  81         const fvsPatchScalarField& pcofrDeltaT =
  82             cofrDeltaT.boundaryField()[patchi];
  83
  84         const fvPatch& p = pcofrDeltaT.patch();
  85         const unallocLabelList& faceCells = p.patch().faceCells();
  86
  87         forAll(pcofrDeltaT, patchFacei)
  88         {
  89             corDeltaT[faceCells[patchFacei]] = max
  90             (
  91                 corDeltaT[faceCells[patchFacei]],
  92                 pcofrDeltaT[patchFacei]
  93             );
  94         }
  95     }
  96
  97     corDeltaT.correctBoundaryConditions();
  98
  99     //corDeltaT = max(corDeltaT, max(corDeltaT)/100.0);
 100
 101     return tcorDeltaT;
 102 }
 103
 104
 105 template<class Type>
 106 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::CofrDeltaT() const
 107 {
 108     const dimensionedScalar& deltaT = mesh().time().deltaT();
 109
 110     const surfaceScalarField& phi =
 111         static_cast<const objectRegistry&>(mesh())
 112         .lookupObject<surfaceScalarField>(phiName_);
 113
 114     if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
 115     {
 116         surfaceScalarField Co
 117         (
 118             mesh().surfaceInterpolation::deltaCoeffs()
 119            *(mag(phi)/mesh().magSf())
 120            *deltaT
 121         );
 122
 123         return max(Co/maxCo_, scalar(1))/deltaT;
 124     }
 125     else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
 126     {
 127         const volScalarField& rho =
 128             static_cast<const objectRegistry&>(mesh())
 129            .lookupObject<volScalarField>(rhoName_).oldTime();
 130
 131         surfaceScalarField Co
 132         (
 133             mesh().surfaceInterpolation::deltaCoeffs()
 134            *(mag(phi)/(fvc::interpolate(rho)*mesh().magSf()))
 135            *deltaT
 136         );
 137
 138         return max(Co/maxCo_, scalar(1))/deltaT;
 139     }
 140     else
 141     {
 142         FatalErrorIn("CoEulerDdtScheme<Type>::CofrDeltaT() const")
 143             << "Incorrect dimensions of phi: " << phi.dimensions()
 144             << abort(FatalError);
 145
 146         return tmp<surfaceScalarField>(NULL);
 147     }
 148 }
 149
 150
 151 template<class Type>
 152 tmp<GeometricField<Type, fvPatchField, volMesh> >
 153 CoEulerDdtScheme<Type>::fvcDdt
 154 (
 155     const dimensioned<Type>& dt
 156 )
 157 {
 158     volScalarField rDeltaT = CorDeltaT();
 159
 160     IOobject ddtIOobject
 161     (
 162         "ddt("+dt.name()+')',
 163         mesh().time().timeName(),
 164         mesh()
 165     );
 166
 167     if (mesh().moving())
 168     {
 169         tmp<GeometricField<Type, fvPatchField, volMesh> > tdtdt
 170         (
 171             new GeometricField<Type, fvPatchField, volMesh>
 172             (
 173                 ddtIOobject,
 174                 mesh(),
 175                 dimensioned<Type>
 176                 (
 177                     "0",
 178                     dt.dimensions()/dimTime,
 179                     pTraits<Type>::zero
 180                 )
 181             )
 182         );
 183
 184         tdtdt().internalField() =
 185             rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
 186
 187         return tdtdt;
 188     }
 189     else
 190     {
 191         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 192         (
 193             new GeometricField<Type, fvPatchField, volMesh>
 194             (
 195                 ddtIOobject,
 196                 mesh(),
 197                 dimensioned<Type>
 198                 (
 199                     "0",
 200                     dt.dimensions()/dimTime,
 201                     pTraits<Type>::zero
 202                 ),
 203                 calculatedFvPatchField<Type>::typeName
 204             )
 205         );
 206     }
 207 }
 208
 209
 210 template<class Type>
 211 tmp<GeometricField<Type, fvPatchField, volMesh> >
 212 CoEulerDdtScheme<Type>::fvcDdt
 213 (
 214     const GeometricField<Type, fvPatchField, volMesh>& vf
 215 )
 216 {
 217     volScalarField rDeltaT = CorDeltaT();
 218
 219     IOobject ddtIOobject
 220     (
 221         "ddt("+vf.name()+')',
 222         mesh().time().timeName(),
 223         mesh()
 224     );
 225
 226     if (mesh().moving())
 227     {
 228         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 229         (
 230             new GeometricField<Type, fvPatchField, volMesh>
 231             (
 232                 ddtIOobject,
 233                 mesh(),
 234                 rDeltaT.dimensions()*vf.dimensions(),
 235                 rDeltaT.internalField()*
 236                 (
 237                     vf.internalField()
 238                   - vf.oldTime().internalField()*mesh().V0()/mesh().V()
 239                 ),
 240                 rDeltaT.boundaryField()*
 241                 (
 242                     vf.boundaryField() - vf.oldTime().boundaryField()
 243                 )
 244             )
 245         );
 246     }
 247     else
 248     {
 249         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 250         (
 251             new GeometricField<Type, fvPatchField, volMesh>
 252             (
 253                 ddtIOobject,
 254                 rDeltaT*(vf - vf.oldTime())
 255             )
 256         );
 257     }
 258 }
 259
 260
 261 template<class Type>
 262 tmp<GeometricField<Type, fvPatchField, volMesh> >
 263 CoEulerDdtScheme<Type>::fvcDdt
 264 (
 265     const dimensionedScalar& rho,
 266     const GeometricField<Type, fvPatchField, volMesh>& vf
 267 )
 268 {
 269     volScalarField rDeltaT = CorDeltaT();
 270
 271     IOobject ddtIOobject
 272     (
 273         "ddt("+rho.name()+','+vf.name()+')',
 274         mesh().time().timeName(),
 275         mesh()
 276     );
 277
 278     if (mesh().moving())
 279     {
 280         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 281         (
 282             new GeometricField<Type, fvPatchField, volMesh>
 283             (
 284                 ddtIOobject,
 285                 mesh(),
 286                 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
 287                 rDeltaT.internalField()*rho.value()*
 288                 (
 289                     vf.internalField()
 290                   - vf.oldTime().internalField()*mesh().V0()/mesh().V()
 291                 ),
 292                 rDeltaT.boundaryField()*rho.value()*
 293                 (
 294                     vf.boundaryField() - vf.oldTime().boundaryField()
 295                 )
 296             )
 297         );
 298     }
 299     else
 300     {
 301         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 302         (
 303             new GeometricField<Type, fvPatchField, volMesh>
 304             (
 305                 ddtIOobject,
 306                 rDeltaT*rho*(vf - vf.oldTime())
 307             )
 308         );
 309     }
 310 }
 311
 312
 313 template<class Type>
 314 tmp<GeometricField<Type, fvPatchField, volMesh> >
 315 CoEulerDdtScheme<Type>::fvcDdt
 316 (
 317     const volScalarField& rho,
 318     const GeometricField<Type, fvPatchField, volMesh>& vf
 319 )
 320 {
 321     volScalarField rDeltaT = CorDeltaT();
 322
 323     IOobject ddtIOobject
 324     (
 325         "ddt("+rho.name()+','+vf.name()+')',
 326         mesh().time().timeName(),
 327         mesh()
 328     );
 329
 330     if (mesh().moving())
 331     {
 332         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 333         (
 334             new GeometricField<Type, fvPatchField, volMesh>
 335             (
 336                 ddtIOobject,
 337                 mesh(),
 338                 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
 339                 rDeltaT.internalField()*
 340                 (
 341                     rho.internalField()*vf.internalField()
 342                   - rho.oldTime().internalField()
 343                    *vf.oldTime().internalField()*mesh().V0()/mesh().V()
 344                 ),
 345                 rDeltaT.boundaryField()*
 346                 (
 347                     rho.boundaryField()*vf.boundaryField()
 348                   - rho.oldTime().boundaryField()
 349                    *vf.oldTime().boundaryField()
 350                 )
 351             )
 352         );
 353     }
 354     else
 355     {
 356         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 357         (
 358             new GeometricField<Type, fvPatchField, volMesh>
 359             (
 360                 ddtIOobject,
 361                 rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
 362             )
 363         );
 364     }
 365 }
 366
 367
 368 template<class Type>
 369 tmp<fvMatrix<Type> >
 370 CoEulerDdtScheme<Type>::fvmDdt
 371 (
 372     GeometricField<Type, fvPatchField, volMesh>& vf
 373 )
 374 {
 375     tmp<fvMatrix<Type> > tfvm
 376     (
 377         new fvMatrix<Type>
 378         (
 379             vf,
 380             vf.dimensions()*dimVol/dimTime
 381         )
 382     );
 383
 384     fvMatrix<Type>& fvm = tfvm();
 385
 386     scalarField rDeltaT = CorDeltaT()().internalField();
 387
 388     fvm.diag() = rDeltaT*mesh().V();
 389
 390     if (mesh().moving())
 391     {
 392         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
 393     }
 394     else
 395     {
 396         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
 397     }
 398
 399     return tfvm;
 400 }
 401
 402
 403 template<class Type>
 404 tmp<fvMatrix<Type> >
 405 CoEulerDdtScheme<Type>::fvmDdt
 406 (
 407     const dimensionedScalar& rho,
 408     GeometricField<Type, fvPatchField, volMesh>& vf
 409 )
 410 {
 411     tmp<fvMatrix<Type> > tfvm
 412     (
 413         new fvMatrix<Type>
 414         (
 415             vf,
 416             rho.dimensions()*vf.dimensions()*dimVol/dimTime
 417         )
 418     );
 419     fvMatrix<Type>& fvm = tfvm();
 420
 421     scalarField rDeltaT = CorDeltaT()().internalField();
 422
 423     fvm.diag() = rDeltaT*rho.value()*mesh().V();
 424
 425     if (mesh().moving())
 426     {
 427         fvm.source() = rDeltaT
 428             *rho.value()*vf.oldTime().internalField()*mesh().V0();
 429     }
 430     else
 431     {
 432         fvm.source() = rDeltaT
 433             *rho.value()*vf.oldTime().internalField()*mesh().V();
 434     }
 435
 436     return tfvm;
 437 }
 438
 439
 440 template<class Type>
 441 tmp<fvMatrix<Type> >
 442 CoEulerDdtScheme<Type>::fvmDdt
 443 (
 444     const volScalarField& rho,
 445     GeometricField<Type, fvPatchField, volMesh>& vf
 446 )
 447 {
 448     tmp<fvMatrix<Type> > tfvm
 449     (
 450         new fvMatrix<Type>
 451         (
 452             vf,
 453             rho.dimensions()*vf.dimensions()*dimVol/dimTime
 454         )
 455     );
 456     fvMatrix<Type>& fvm = tfvm();
 457
 458     scalarField rDeltaT = CorDeltaT()().internalField();
 459
 460     fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
 461
 462     if (mesh().moving())
 463     {
 464         fvm.source() = rDeltaT
 465             *rho.oldTime().internalField()
 466             *vf.oldTime().internalField()*mesh().V0();
 467     }
 468     else
 469     {
 470         fvm.source() = rDeltaT
 471             *rho.oldTime().internalField()
 472             *vf.oldTime().internalField()*mesh().V();
 473     }
 474
 475     return tfvm;
 476 }
 477
 478
 479 template<class Type>
 480 tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
 481 CoEulerDdtScheme<Type>::fvcDdtPhiCorr
 482 (
 483     const volScalarField& rA,
 484     const GeometricField<Type, fvPatchField, volMesh>& U,
 485     const fluxFieldType& phi
 486 )
 487 {
 488     IOobject ddtIOobject
 489     (
 490         "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
 491         mesh().time().timeName(),
 492         mesh()
 493     );
 494
 495     if (mesh().moving())
 496     {
 497         return tmp<fluxFieldType>
 498         (
 499             new fluxFieldType
 500             (
 501                 ddtIOobject,
 502                 mesh(),
 503                 dimensioned<typename flux<Type>::type>
 504                 (
 505                     "0",
 506                     rA.dimensions()*phi.dimensions()/dimTime,
 507                     pTraits<typename flux<Type>::type>::zero
 508                 )
 509             )
 510         );
 511     }
 512     else
 513     {
 514         volScalarField rDeltaT = CorDeltaT();
 515
 516         return tmp<fluxFieldType>
 517         (
 518             new fluxFieldType
 519             (
 520                 ddtIOobject,
 521                 fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
 522                 (
 523                     fvc::interpolate(rDeltaT*rA)*phi.oldTime()
 524                   - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
 525                 )
 526             )
 527         );
 528     }
 529 }
 530
 531
 532 template<class Type>
 533 tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
 534 CoEulerDdtScheme<Type>::fvcDdtPhiCorr
 535 (
 536     const volScalarField& rA,
 537     const volScalarField& rho,
 538     const GeometricField<Type, fvPatchField, volMesh>& U,
 539     const fluxFieldType& phi
 540 )
 541 {
 542     IOobject ddtIOobject
 543     (
 544         "ddtPhiCorr("
 545       + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
 546         mesh().time().timeName(),
 547         mesh()
 548     );
 549
 550     if (mesh().moving())
 551     {
 552         return tmp<fluxFieldType>
 553         (
 554             new fluxFieldType
 555             (
 556                 ddtIOobject,
 557                 mesh(),
 558                 dimensioned<typename flux<Type>::type>
 559                 (
 560                     "0",
 561                     rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
 562                     pTraits<typename flux<Type>::type>::zero
 563                 )
 564             )
 565         );
 566     }
 567     else
 568     {
 569         volScalarField rDeltaT = CorDeltaT();
 570
 571         if
 572         (
 573             U.dimensions() == dimVelocity
 574          && phi.dimensions() == dimVelocity*dimArea
 575         )
 576         {
 577             return tmp<fluxFieldType>
 578             (
 579                 new fluxFieldType
 580                 (
 581                     ddtIOobject,
 582                     fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
 583                    *(
 584                         fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
 585                       - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
 586                       & mesh().Sf())
 587                     )
 588                 )
 589             );
 590         }
 591         else if
 592         (
 593             U.dimensions() == dimVelocity
 594          && phi.dimensions() == dimDensity*dimVelocity*dimArea
 595         )
 596         {
 597             return tmp<fluxFieldType>
 598             (
 599                 new fluxFieldType
 600                 (
 601                     ddtIOobject,
 602                     fvcDdtPhiCoeff
 603                     (
 604                         U.oldTime(),
 605                         phi.oldTime()/fvc::interpolate(rho.oldTime())
 606                     )
 607                    *(
 608                         fvc::interpolate(rDeltaT*rA*rho.oldTime())
 609                        *phi.oldTime()/fvc::interpolate(rho.oldTime())
 610                       - (
 611                             fvc::interpolate
 612                             (
 613                                 rDeltaT*rA*rho.oldTime()*U.oldTime()
 614                             ) & mesh().Sf()
 615                         )
 616                     )
 617                 )
 618             );
 619         }
 620         else if
 621         (
 622             U.dimensions() == dimDensity*dimVelocity
 623          && phi.dimensions() == dimDensity*dimVelocity*dimArea
 624         )
 625         {
 626             return tmp<fluxFieldType>
 627             (
 628                 new fluxFieldType
 629                 (
 630                     ddtIOobject,
 631                     fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
 632                    *(
 633                         fvc::interpolate(rDeltaT*rA)*phi.oldTime()
 634                       - (
 635                             fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
 636                         )
 637                     )
 638                 )
 639             );
 640         }
 641         else
 642         {
 643             FatalErrorIn
 644             (
 645                 "CoEulerDdtScheme<Type>::fvcDdtPhiCorr"
 646             )   << "dimensions of phi are not correct"
 647                 << abort(FatalError);
 648
 649             return fluxFieldType::null();
 650         }
 651     }
 652 }
 653
 654
 655 template<class Type>
 656 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::meshPhi
 657 (
 658     const GeometricField<Type, fvPatchField, volMesh>&
 659 )
 660 {
 661     return tmp<surfaceScalarField>
 662     (
 663         new surfaceScalarField
 664         (
 665             IOobject
 666             (
 667                 "meshPhi",
 668                 mesh().time().timeName(),
 669                 mesh()
 670             ),
 671             mesh(),
 672             dimensionedScalar("0", dimVolume/dimTime, 0.0)
 673         )
 674     );
 675 }
 676
 677
 678 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 679
 680 } // End namespace fv
 681
 682 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 683
 684 } // End namespace Foam
 685
 686 // ************************************************************************* //