src/lagrangian/dieselSpray/spraySubModels/wallModel/reflectParcel/reflectParcel.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 "reflectParcel.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "wallPolyPatch.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 defineTypeNameAndDebug(reflectParcel, 0);
  39
  40 addToRunTimeSelectionTable
  41 (
  42     wallModel,
  43     reflectParcel,
  44     dictionary
  45 );
  46
  47
  48 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  49
  50 // Construct from components
  51 reflectParcel::reflectParcel
  52 (
  53     const dictionary& dict,
  54     const volVectorField& U,
  55     spray& sm
  56 )
  57 :
  58     wallModel(dict, U, sm),
  59     U_(U),
  60     coeffsDict_(dict.subDict(typeName + "Coeffs")),
  61     elasticity_(readScalar(coeffsDict_.lookup("elasticity")))
  62 {}
  63
  64
  65 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  66
  67 reflectParcel::~reflectParcel()
  68 {}
  69
  70
  71 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  72
  73 // Return 'keepParcel'
  74 bool reflectParcel::wallTreatment
  75 (
  76     parcel& p,
  77     const label globalFacei
  78 ) const
  79 {
  80     label patchi = p.patch(globalFacei);
  81     label facei = p.patchFace(patchi, globalFacei);
  82
  83     const polyMesh& mesh = spray_.mesh();
  84
  85     if (isType<wallPolyPatch>(mesh_.boundaryMesh()[patchi]))
  86     {
  87         // wallNormal defined to point outwards of domain
  88         vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
  89         Sf /= mag(Sf);
  90
  91         if (!mesh.moving())
  92         {
  93             // static mesh
  94             scalar Un = p.U() & Sf;
  95
  96             if (Un > 0)
  97             {
  98                 p.U() -= (1.0 + elasticity_)*Un*Sf;
  99             }
 100
 101         }
 102         else
 103         {
 104             // moving mesh
 105             vector Ub1 = U_.boundaryField()[patchi][facei];
 106             vector Ub0 = U_.oldTime().boundaryField()[patchi][facei];
 107
 108             scalar dt = spray_.runTime().deltaT().value();
 109             const vectorField& oldPoints = mesh.oldPoints();
 110
 111             const vector& Cf1 = mesh.faceCentres()[globalFacei];
 112
 113             vector Cf0 = mesh.faces()[globalFacei].centre(oldPoints);
 114             vector Cf = Cf0 + p.stepFraction()*(Cf1 - Cf0);
 115             vector Sf0 = mesh.faces()[globalFacei].normal(oldPoints);
 116
 117             // for layer addition Sf0 = vector::zero and we use Sf
 118             if (mag(Sf0) > SMALL)
 119             {
 120                 Sf0 /= mag(Sf0);
 121             }
 122             else
 123             {
 124                 Sf0 = Sf;
 125             }
 126
 127             scalar magSfDiff = mag(Sf - Sf0);
 128
 129             vector Ub = Ub0 + p.stepFraction()*(Ub1 - Ub0);
 130
 131             if (magSfDiff > SMALL)
 132             {
 133                 // rotation + translation
 134                 vector Sfp = Sf0 + p.stepFraction()*(Sf - Sf0);
 135
 136                 vector omega = Sf0 ^ Sf;
 137                 scalar magOmega = mag(omega);
 138                 omega /= magOmega+SMALL;
 139
 140                 scalar phiVel = ::asin(magOmega)/dt;
 141
 142                 scalar dist = (p.position() - Cf) & Sfp;
 143                 vector pos = p.position() - dist*Sfp;
 144                 vector vrot = phiVel*(omega ^ (pos - Cf));
 145
 146                 vector v = Ub + vrot;
 147
 148                 scalar Un = ((p.U() - v) & Sfp);
 149
 150                 if (Un > 0.0)
 151                 {
 152                     p.U() -= (1.0 + elasticity_)*Un*Sfp;
 153                 }
 154             }
 155             else
 156             {
 157                 // translation
 158                 vector Ur = p.U() - Ub;
 159                 scalar Urn = Ur & Sf;
 160                 /*
 161                 if (mag(Ub-v) > SMALL)
 162                 {
 163                     Info << "reflectParcel:: v = " << v
 164                         << ", Ub = " << Ub
 165                         << ", facei = " << facei
 166                         << ", patchi = " << patchi
 167                         << ", globalFacei = " << globalFacei
 168                         << ", name = " << mesh_.boundaryMesh()[patchi].name()
 169                         << endl;
 170                 }
 171                     */
 172                 if (Urn > 0.0)
 173                 {
 174                     p.U() -= (1.0 + elasticity_)*Urn*Sf;
 175                 }
 176             }
 177         }
 178
 179     }
 180     else
 181     {
 182         FatalError
 183             << "bool reflectParcel::wallTreatment(parcel& parcel) const "
 184                 << " parcel has hit a boundary "
 185                 << mesh_.boundary()[patchi].type()
 186                 << " which not yet has been implemented."
 187             << abort(FatalError);
 188     }
 189     return true;
 190 }
 191
 192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 193
 194 } // End namespace Foam
 195
 196 // ************************************************************************* //