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