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