src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.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 "EulerDdtScheme.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<GeometricField<Type, fvPatchField, volMesh> >
  46 EulerDdtScheme<Type>::fvcDdt
  47 (
  48     const dimensioned<Type>& dt
  49 )
  50 {
  51     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
  52
  53     IOobject ddtIOobject
  54     (
  55         "ddt("+dt.name()+')',
  56         mesh().time().timeName(),
  57         mesh()
  58     );
  59
  60     if (mesh().moving())
  61     {
  62         tmp<GeometricField<Type, fvPatchField, volMesh> > tdtdt
  63         (
  64             new GeometricField<Type, fvPatchField, volMesh>
  65             (
  66                 ddtIOobject,
  67                 mesh(),
  68                 dimensioned<Type>
  69                 (
  70                     "0",
  71                     dt.dimensions()/dimTime,
  72                     pTraits<Type>::zero
  73                 )
  74             )
  75         );
  76
  77         tdtdt().internalField() =
  78             rDeltaT.value()*dt.value()*(1.0 - mesh().V0()/mesh().V());
  79
  80         return tdtdt;
  81     }
  82     else
  83     {
  84         return tmp<GeometricField<Type, fvPatchField, volMesh> >
  85         (
  86             new GeometricField<Type, fvPatchField, volMesh>
  87             (
  88                 ddtIOobject,
  89                 mesh(),
  90                 dimensioned<Type>
  91                 (
  92                     "0",
  93                     dt.dimensions()/dimTime,
  94                     pTraits<Type>::zero
  95                 ),
  96                 calculatedFvPatchField<Type>::typeName
  97             )
  98         );
  99     }
 100 }
 101
 102
 103 template<class Type>
 104 tmp<GeometricField<Type, fvPatchField, volMesh> >
 105 EulerDdtScheme<Type>::fvcDdt
 106 (
 107     const GeometricField<Type, fvPatchField, volMesh>& vf
 108 )
 109 {
 110     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 111
 112     IOobject ddtIOobject
 113     (
 114         "ddt("+vf.name()+')',
 115         mesh().time().timeName(),
 116         mesh()
 117     );
 118
 119     if (mesh().moving())
 120     {
 121         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 122         (
 123             new GeometricField<Type, fvPatchField, volMesh>
 124             (
 125                 ddtIOobject,
 126                 mesh(),
 127                 rDeltaT.dimensions()*vf.dimensions(),
 128                 rDeltaT.value()*
 129                 (
 130                     vf.internalField()
 131                   - vf.oldTime().internalField()*mesh().V0()/mesh().V()
 132                 ),
 133                 rDeltaT.value()*
 134                 (
 135                     vf.boundaryField() - vf.oldTime().boundaryField()
 136                 )
 137             )
 138         );
 139     }
 140     else
 141     {
 142         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 143         (
 144             new GeometricField<Type, fvPatchField, volMesh>
 145             (
 146                 ddtIOobject,
 147                 rDeltaT*(vf - vf.oldTime())
 148             )
 149         );
 150     }
 151 }
 152
 153
 154 template<class Type>
 155 tmp<GeometricField<Type, fvPatchField, volMesh> >
 156 EulerDdtScheme<Type>::fvcDdt
 157 (
 158     const dimensionedScalar& rho,
 159     const GeometricField<Type, fvPatchField, volMesh>& vf
 160 )
 161 {
 162     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 163
 164     IOobject ddtIOobject
 165     (
 166         "ddt("+rho.name()+','+vf.name()+')',
 167         mesh().time().timeName(),
 168         mesh()
 169     );
 170
 171     if (mesh().moving())
 172     {
 173         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 174         (
 175             new GeometricField<Type, fvPatchField, volMesh>
 176             (
 177                 ddtIOobject,
 178                 mesh(),
 179                 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
 180                 rDeltaT.value()*rho.value()*
 181                 (
 182                     vf.internalField()
 183                   - vf.oldTime().internalField()*mesh().V0()/mesh().V()
 184                 ),
 185                 rDeltaT.value()*rho.value()*
 186                 (
 187                     vf.boundaryField() - vf.oldTime().boundaryField()
 188                 )
 189             )
 190         );
 191     }
 192     else
 193     {
 194         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 195         (
 196             new GeometricField<Type, fvPatchField, volMesh>
 197             (
 198                 ddtIOobject,
 199                 rDeltaT*rho*(vf - vf.oldTime())
 200             )
 201         );
 202     }
 203 }
 204
 205
 206 template<class Type>
 207 tmp<GeometricField<Type, fvPatchField, volMesh> >
 208 EulerDdtScheme<Type>::fvcDdt
 209 (
 210     const volScalarField& rho,
 211     const GeometricField<Type, fvPatchField, volMesh>& vf
 212 )
 213 {
 214     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 215
 216     IOobject ddtIOobject
 217     (
 218         "ddt("+rho.name()+','+vf.name()+')',
 219         mesh().time().timeName(),
 220         mesh()
 221     );
 222
 223     if (mesh().moving())
 224     {
 225         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 226         (
 227             new GeometricField<Type, fvPatchField, volMesh>
 228             (
 229                 ddtIOobject,
 230                 mesh(),
 231                 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
 232                 rDeltaT.value()*
 233                 (
 234                     rho.internalField()*vf.internalField()
 235                   - rho.oldTime().internalField()
 236                    *vf.oldTime().internalField()*mesh().V0()/mesh().V()
 237                 ),
 238                 rDeltaT.value()*
 239                 (
 240                     rho.boundaryField()*vf.boundaryField()
 241                   - rho.oldTime().boundaryField()
 242                    *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*(rho*vf - rho.oldTime()*vf.oldTime())
 255             )
 256         );
 257     }
 258 }
 259
 260
 261 template<class Type>
 262 tmp<fvMatrix<Type> >
 263 EulerDdtScheme<Type>::fvmDdt
 264 (
 265     GeometricField<Type, fvPatchField, volMesh>& vf
 266 )
 267 {
 268     tmp<fvMatrix<Type> > tfvm
 269     (
 270         new fvMatrix<Type>
 271         (
 272             vf,
 273             vf.dimensions()*dimVol/dimTime
 274         )
 275     );
 276
 277     fvMatrix<Type>& fvm = tfvm();
 278
 279     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 280
 281     fvm.diag() = rDeltaT*mesh().V();
 282
 283     if (mesh().moving())
 284     {
 285         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
 286     }
 287     else
 288     {
 289         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
 290     }
 291
 292     return tfvm;
 293 }
 294
 295
 296 template<class Type>
 297 tmp<fvMatrix<Type> >
 298 EulerDdtScheme<Type>::fvmDdt
 299 (
 300     const dimensionedScalar& rho,
 301     GeometricField<Type, fvPatchField, volMesh>& vf
 302 )
 303 {
 304     tmp<fvMatrix<Type> > tfvm
 305     (
 306         new fvMatrix<Type>
 307         (
 308             vf,
 309             rho.dimensions()*vf.dimensions()*dimVol/dimTime
 310         )
 311     );
 312     fvMatrix<Type>& fvm = tfvm();
 313
 314     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 315
 316     fvm.diag() = rDeltaT*rho.value()*mesh().V();
 317
 318     if (mesh().moving())
 319     {
 320         fvm.source() = rDeltaT
 321             *rho.value()*vf.oldTime().internalField()*mesh().V0();
 322     }
 323     else
 324     {
 325         fvm.source() = rDeltaT
 326             *rho.value()*vf.oldTime().internalField()*mesh().V();
 327     }
 328
 329     return tfvm;
 330 }
 331
 332
 333 template<class Type>
 334 tmp<fvMatrix<Type> >
 335 EulerDdtScheme<Type>::fvmDdt
 336 (
 337     const volScalarField& rho,
 338     GeometricField<Type, fvPatchField, volMesh>& vf
 339 )
 340 {
 341     tmp<fvMatrix<Type> > tfvm
 342     (
 343         new fvMatrix<Type>
 344         (
 345             vf,
 346             rho.dimensions()*vf.dimensions()*dimVol/dimTime
 347         )
 348     );
 349     fvMatrix<Type>& fvm = tfvm();
 350
 351     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 352
 353     fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
 354
 355     if (mesh().moving())
 356     {
 357         fvm.source() = rDeltaT
 358             *rho.oldTime().internalField()
 359             *vf.oldTime().internalField()*mesh().V0();
 360     }
 361     else
 362     {
 363         fvm.source() = rDeltaT
 364             *rho.oldTime().internalField()
 365             *vf.oldTime().internalField()*mesh().V();
 366     }
 367
 368     return tfvm;
 369 }
 370
 371
 372 template<class Type>
 373 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
 374 EulerDdtScheme<Type>::fvcDdtPhiCorr
 375 (
 376     const volScalarField& rA,
 377     const GeometricField<Type, fvPatchField, volMesh>& U,
 378     const fluxFieldType& phiAbs
 379 )
 380 {
 381     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 382
 383     IOobject ddtIOobject
 384     (
 385         "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phiAbs.name() + ')',
 386         mesh().time().timeName(),
 387         mesh()
 388     );
 389
 390     tmp<fluxFieldType> phiCorr =
 391         phiAbs.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
 392
 393     return tmp<fluxFieldType>
 394     (
 395         new fluxFieldType
 396         (
 397             ddtIOobject,
 398             fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime(), phiCorr())
 399            *fvc::interpolate(rDeltaT*rA)*phiCorr
 400         )
 401     );
 402 }
 403
 404
 405 template<class Type>
 406 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
 407 EulerDdtScheme<Type>::fvcDdtPhiCorr
 408 (
 409     const volScalarField& rA,
 410     const volScalarField& rho,
 411     const GeometricField<Type, fvPatchField, volMesh>& U,
 412     const fluxFieldType& phiAbs
 413 )
 414 {
 415     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 416
 417     IOobject ddtIOobject
 418     (
 419         "ddtPhiCorr("
 420       + rA.name() + ','
 421       + rho.name() + ','
 422       + U.name() + ','
 423       + phiAbs.name() + ')',
 424         mesh().time().timeName(),
 425         mesh()
 426     );
 427
 428     if
 429     (
 430         U.dimensions() == dimVelocity
 431      && phiAbs.dimensions() == dimVelocity*dimArea
 432     )
 433     {
 434         return tmp<fluxFieldType>
 435         (
 436             new fluxFieldType
 437             (
 438                 ddtIOobject,
 439                 rDeltaT
 440                *fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
 441                *(
 442                    fvc::interpolate(rA*rho.oldTime())*phiAbs.oldTime()
 443                  - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
 444                   & mesh().Sf())
 445                 )
 446             )
 447         );
 448     }
 449     else if
 450     (
 451         U.dimensions() == dimVelocity
 452      && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
 453     )
 454     {
 455         return tmp<fluxFieldType>
 456         (
 457             new fluxFieldType
 458             (
 459                 ddtIOobject,
 460                 rDeltaT
 461                *fvcDdtPhiCoeff
 462                 (
 463                     U.oldTime(),
 464                     phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
 465                 )
 466                *(
 467                    fvc::interpolate(rA*rho.oldTime())
 468                   *phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
 469                  - (
 470                        fvc::interpolate
 471                        (
 472                            rA*rho.oldTime()*U.oldTime()
 473                        ) & mesh().Sf()
 474                    )
 475                 )
 476             )
 477         );
 478     }
 479     else if
 480     (
 481         U.dimensions() == rho.dimensions()*dimVelocity
 482      && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
 483     )
 484     {
 485         return tmp<fluxFieldType>
 486         (
 487             new fluxFieldType
 488             (
 489                 ddtIOobject,
 490                 rDeltaT
 491                *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
 492                *(
 493                    fvc::interpolate(rA)*phiAbs.oldTime()
 494                  - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
 495                 )
 496             )
 497         );
 498     }
 499     else
 500     {
 501         FatalErrorIn
 502         (
 503             "EulerDdtScheme<Type>::fvcDdtPhiCorr"
 504         )   << "dimensions of phiAbs are not correct"
 505             << abort(FatalError);
 506
 507         return fluxFieldType::null();
 508     }
 509 }
 510
 511
 512 template<class Type>
 513 tmp<surfaceScalarField> EulerDdtScheme<Type>::meshPhi
 514 (
 515     const GeometricField<Type, fvPatchField, volMesh>&
 516 )
 517 {
 518     return mesh().phi();
 519 }
 520
 521
 522 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 523
 524 } // End namespace fv
 525
 526 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 527
 528 } // End namespace Foam
 529
 530 // ************************************************************************* //