src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.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 "EulerD2dt2Scheme.H"
  28 #include "fvcDiv.H"
  29 #include "fvMatrices.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace fv
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 template<class Type>
  44 tmp<GeometricField<Type, fvPatchField, volMesh> >
  45 EulerD2dt2Scheme<Type>::fvcD2dt2
  46 (
  47     const GeometricField<Type, fvPatchField, volMesh>& vf
  48 )
  49 {
  50     dimensionedScalar rDeltaT2 =
  51         4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
  52
  53     IOobject d2dt2IOobject
  54     (
  55         "d2dt2("+vf.name()+')',
  56         mesh().time().timeName(),
  57         mesh(),
  58         IOobject::NO_READ,
  59         IOobject::NO_WRITE
  60     );
  61
  62     scalar deltaT = mesh().time().deltaT().value();
  63     scalar deltaT0 = mesh().time().deltaT0().value();
  64
  65     scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
  66     scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
  67     scalar coefft0  = coefft + coefft00;
  68
  69     if (mesh().moving())
  70     {
  71         scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
  72
  73         scalarField VV0 = mesh().V() + mesh().V0();
  74         scalarField V0V00 = mesh().V0() + mesh().V00();
  75
  76         return tmp<GeometricField<Type, fvPatchField, volMesh> >
  77         (
  78             new GeometricField<Type, fvPatchField, volMesh>
  79             (
  80                 d2dt2IOobject,
  81                 mesh(),
  82                 rDeltaT2.dimensions()*vf.dimensions(),
  83                 halfRdeltaT2*
  84                 (
  85                     coefft*VV0*vf.internalField()
  86
  87                   - (coefft*VV0 + coefft00*V0V00)
  88                    *vf.oldTime().internalField()
  89
  90                   + (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
  91                 )/mesh().V(),
  92                 rDeltaT2.value()*
  93                 (
  94                     coefft*vf.boundaryField()
  95                   - coefft0*vf.oldTime().boundaryField()
  96                   + coefft00*vf.oldTime().oldTime().boundaryField()
  97                 )
  98             )
  99         );
 100     }
 101     else
 102     {
 103         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 104         (
 105             new GeometricField<Type, fvPatchField, volMesh>
 106             (
 107                 d2dt2IOobject,
 108                 rDeltaT2*
 109                 (
 110                     coefft*vf
 111                   - coefft0*vf.oldTime()
 112                   + coefft00*vf.oldTime().oldTime()
 113                 )
 114             )
 115         );
 116     }
 117 }
 118
 119
 120 template<class Type>
 121 tmp<GeometricField<Type, fvPatchField, volMesh> >
 122 EulerD2dt2Scheme<Type>::fvcD2dt2
 123 (
 124     const volScalarField& rho,
 125     const GeometricField<Type, fvPatchField, volMesh>& vf
 126 )
 127 {
 128     dimensionedScalar rDeltaT2 =
 129         4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
 130
 131     IOobject d2dt2IOobject
 132     (
 133         "d2dt2("+rho.name()+','+vf.name()+')',
 134         mesh().time().timeName(),
 135         mesh(),
 136         IOobject::NO_READ,
 137         IOobject::NO_WRITE
 138     );
 139
 140     scalar deltaT = mesh().time().deltaT().value();
 141     scalar deltaT0 = mesh().time().deltaT0().value();
 142
 143     scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
 144     scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 145
 146     if (mesh().moving())
 147     {
 148         scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
 149         scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
 150
 151         scalarField VV0rhoRho0 =
 152             (mesh().V() + mesh().V0())
 153            *(rho.internalField() + rho.oldTime().internalField());
 154
 155         scalarField V0V00rho0Rho00 =
 156             (mesh().V0() + mesh().V00())
 157            *(
 158                rho.oldTime().internalField()
 159              + rho.oldTime().oldTime().internalField()
 160             );
 161
 162         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 163         (
 164             new GeometricField<Type, fvPatchField, volMesh>
 165             (
 166                 d2dt2IOobject,
 167                 mesh(),
 168                 rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
 169                 quarterRdeltaT2*
 170                 (
 171                     coefft*VV0rhoRho0*vf.internalField()
 172
 173                   - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
 174                    *vf.oldTime().internalField()
 175
 176                   + (coefft00*V0V00rho0Rho00)
 177                    *vf.oldTime().oldTime().internalField()
 178                 )/mesh().V(),
 179                 halfRdeltaT2*
 180                 (
 181                     coefft
 182                    *(rho.boundaryField() + rho.oldTime().boundaryField())
 183                    *vf.boundaryField()
 184
 185                   - (
 186                         coefft
 187                        *(
 188                            rho.boundaryField()
 189                          + rho.oldTime().boundaryField()
 190                         )
 191                       + coefft00
 192                        *(
 193                            rho.oldTime().boundaryField()
 194                          + rho.oldTime().oldTime().boundaryField()
 195                         )
 196                     )*vf.oldTime().boundaryField()
 197
 198                   + coefft00
 199                    *(
 200                        rho.oldTime().boundaryField()
 201                      + rho.oldTime().oldTime().boundaryField()
 202                     )*vf.oldTime().oldTime().boundaryField()
 203                 )
 204             )
 205         );
 206     }
 207     else
 208     {
 209         dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
 210
 211         volScalarField rhoRho0 = rho + rho.oldTime();
 212         volScalarField rho0Rho00 = rho.oldTime() +rho.oldTime().oldTime();
 213
 214         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 215         (
 216             new GeometricField<Type, fvPatchField, volMesh>
 217             (
 218                 d2dt2IOobject,
 219                 halfRdeltaT2*
 220                 (
 221                     coefft*rhoRho0*vf
 222                   - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
 223                   + coefft00*rho0Rho00*vf.oldTime().oldTime()
 224                 )
 225             )
 226         );
 227     }
 228 }
 229
 230
 231 template<class Type>
 232 tmp<fvMatrix<Type> >
 233 EulerD2dt2Scheme<Type>::fvmD2dt2
 234 (
 235     GeometricField<Type, fvPatchField, volMesh>& vf
 236 )
 237 {
 238     tmp<fvMatrix<Type> > tfvm
 239     (
 240         new fvMatrix<Type>
 241         (
 242             vf,
 243             vf.dimensions()*dimVol/dimTime/dimTime
 244         )
 245     );
 246
 247     fvMatrix<Type>& fvm = tfvm();
 248
 249     scalar deltaT = mesh().time().deltaT().value();
 250     scalar deltaT0 = mesh().time().deltaT0().value();
 251
 252     scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
 253     scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 254     scalar coefft0  = coefft + coefft00;
 255
 256     scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
 257
 258     if (mesh().moving())
 259     {
 260         scalar halfRdeltaT2 = rDeltaT2/2.0;
 261
 262         scalarField VV0 = mesh().V() + mesh().V0();
 263         scalarField V0V00 = mesh().V0() + mesh().V00();
 264
 265         fvm.diag() = (coefft*halfRdeltaT2)*VV0;
 266
 267         fvm.source() = halfRdeltaT2*
 268         (
 269             (coefft*VV0 + coefft00*V0V00)
 270            *vf.oldTime().internalField()
 271
 272           - (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
 273         );
 274     }
 275     else
 276     {
 277         fvm.diag() = (coefft*rDeltaT2)*mesh().V();
 278
 279         fvm.source() = rDeltaT2*mesh().V()*
 280         (
 281             coefft0*vf.oldTime().internalField()
 282           - coefft00*vf.oldTime().oldTime().internalField()
 283         );
 284     }
 285
 286     return tfvm;
 287 }
 288
 289
 290 template<class Type>
 291 tmp<fvMatrix<Type> >
 292 EulerD2dt2Scheme<Type>::fvmD2dt2
 293 (
 294     const dimensionedScalar& rho,
 295     GeometricField<Type, fvPatchField, volMesh>& vf
 296 )
 297 {
 298     tmp<fvMatrix<Type> > tfvm
 299     (
 300         new fvMatrix<Type>
 301         (
 302             vf,
 303             rho.dimensions()*vf.dimensions()*dimVol
 304             /dimTime/dimTime
 305         )
 306     );
 307
 308     fvMatrix<Type>& fvm = tfvm();
 309
 310     scalar deltaT = mesh().time().deltaT().value();
 311     scalar deltaT0 = mesh().time().deltaT0().value();
 312
 313     scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
 314     scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 315
 316     scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
 317
 318     if (mesh().moving())
 319     {
 320         scalar halfRdeltaT2 = 0.5*rDeltaT2;
 321
 322         scalarField VV0 = mesh().V() + mesh().V0();
 323
 324         scalarField V0V00 = mesh().V0() + mesh().V00();
 325
 326         fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;
 327
 328         fvm.source() = halfRdeltaT2*rho.value()*
 329         (
 330             (coefft*VV0 + coefft00*V0V00)
 331            *vf.oldTime().internalField()
 332
 333           - (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
 334         );
 335     }
 336     else
 337     {
 338         fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();
 339
 340         fvm.source() = rDeltaT2*mesh().V()*rho.value()*
 341         (
 342             (coefft + coefft00)*vf.oldTime().internalField()
 343           - coefft00*vf.oldTime().oldTime().internalField()
 344         );
 345     }
 346
 347     return tfvm;
 348 }
 349
 350
 351 template<class Type>
 352 tmp<fvMatrix<Type> >
 353 EulerD2dt2Scheme<Type>::fvmD2dt2
 354 (
 355     const volScalarField& rho,
 356     GeometricField<Type, fvPatchField, volMesh>& vf
 357 )
 358 {
 359     tmp<fvMatrix<Type> > tfvm
 360     (
 361         new fvMatrix<Type>
 362         (
 363             vf,
 364             rho.dimensions()*vf.dimensions()*dimVol
 365             /dimTime/dimTime
 366         )
 367     );
 368
 369     fvMatrix<Type>& fvm = tfvm();
 370
 371     scalar deltaT = mesh().time().deltaT().value();
 372     scalar deltaT0 = mesh().time().deltaT0().value();
 373
 374     scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
 375     scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 376
 377     scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
 378
 379     if (mesh().moving())
 380     {
 381         scalar quarterRdeltaT2 = 0.25*rDeltaT2;
 382
 383         scalarField VV0rhoRho0 =
 384             (mesh().V() + mesh().V0())
 385            *(rho.internalField() + rho.oldTime().internalField());
 386
 387         scalarField V0V00rho0Rho00 =
 388             (mesh().V0() + mesh().V00())
 389            *(
 390                rho.oldTime().internalField()
 391              + rho.oldTime().oldTime().internalField()
 392             );
 393
 394         fvm.diag() = (coefft*quarterRdeltaT2)*VV0rhoRho0;
 395
 396         fvm.source() = quarterRdeltaT2*
 397         (
 398             (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
 399            *vf.oldTime().internalField()
 400
 401           - (coefft00*V0V00rho0Rho00)
 402            *vf.oldTime().oldTime().internalField()
 403         );
 404     }
 405     else
 406     {
 407         scalar halfRdeltaT2 = 0.5*rDeltaT2;
 408
 409         scalarField rhoRho0 =
 410             (rho.internalField() + rho.oldTime().internalField());
 411
 412         scalarField rho0Rho00 =
 413         (
 414             rho.oldTime().internalField()
 415           + rho.oldTime().oldTime().internalField()
 416         );
 417
 418         fvm.diag() = (coefft*halfRdeltaT2)*mesh().V()*rhoRho0;
 419
 420         fvm.source() = halfRdeltaT2*mesh().V()*
 421         (
 422             (coefft*rhoRho0 + coefft00*rho0Rho00)
 423            *vf.oldTime().internalField()
 424
 425           - (coefft00*rho0Rho00)
 426            *vf.oldTime().oldTime().internalField()
 427         );
 428     }
 429
 430     return tfvm;
 431 }
 432
 433
 434 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 435
 436 } // End namespace fv
 437
 438 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 439
 440 } // End namespace Foam
 441
 442 // ************************************************************************* //