src/lagrangian/dieselSpray/spraySubModels/wallModel/reflectParcel/reflectParcel.C2

   1 // The FOAM Project // File: reflectParcel.C
   2 /*
   3 -------------------------------------------------------------------------------
   4  =========         | Class Implementation
   5  \\      /         |
   6   \\    /          | Name:   reflectParcel
   7    \\  /           | Family: wallModel
   8     \\/            |
   9     F ield         | FOAM version: 2.3
  10     O peration     |
  11     A and          | Copyright (C) 1991-2008 OpenCFD Ltd.
  12     M anipulation  |          All Rights Reserved.
  13 -------------------------------------------------------------------------------
  14 DESCRIPTION
  15
  16 AUTHOR
  17     Henry G Weller.
  18
  19 -------------------------------------------------------------------------------
  20 */
  21
  22 #include "reflectParcel.H"
  23 #include "addToRunTimeSelectionTable.H"
  24 #include "wallPolyPatch.H"
  25
  26 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  27
  28 namespace Foam
  29 {
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 defineTypeNameAndDebug(reflectParcel, 0);
  34
  35 addToRunTimeSelectionTable
  36 (
  37     wallModel,
  38     reflectParcel,
  39     dictionary
  40 );
  41
  42
  43 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  44
  45 // Construct from components
  46 reflectParcel::reflectParcel
  47 (
  48     const dictionary& dict,
  49     const fvMesh& mesh,
  50     spray& sm
  51 )
  52 :
  53     wallModel(dict, mesh, sm),
  54     coeffsDict_(dict.subDict(typeName + "Coeffs")),
  55     elasticity_(readScalar(coeffsDict_.lookup("elasticity")))
  56 {}
  57
  58
  59 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  60
  61 reflectParcel::~reflectParcel()
  62 {}
  63
  64
  65 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  66
  67 // Return 'keepParcel'
  68 bool reflectParcel::wallTreatment
  69 (
  70     parcel& p
  71 ) const
  72 {
  73     label patchi = p.patch();
  74     label facei = p.patchFace(patchi);
  75
  76     const polyMesh& mesh = spray_.mesh();
  77
  78     if (typeid(mesh_.boundaryMesh()[patchi]) == typeid(wallPolyPatch))
  79     {
  80         // wallNormal defined to point outwards of domain
  81         vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
  82         Sf /= mag(Sf);
  83
  84         // adjust velocity when wall is moving
  85         if (!mesh.moving())
  86         {
  87             scalar Un = p.U() & Sf;
  88
  89             if (Un > 0)
  90             {
  91                 p.U() -= (1.0 + elasticity_)*Un*Sf;
  92             }
  93
  94         }
  95         else
  96         {
  97             label facep = p.face();
  98
  99             scalar dt = spray_.runTime().deltaT().value();
 100             scalar fraction = p.t0()/dt;
 101             const vectorField& oldPoints = mesh.oldPoints();
 102             const vector& Cf1 = mesh.faceCentres()[facep];
 103
 104             vector Cf0 = mesh.faces()[facep].centre(oldPoints);
 105             vector Cf = Cf0 + fraction*(Cf1 - Cf0);
 106             vector Sf0 = mesh.faces()[facep].normal(oldPoints);
 107             Sf0 /= mag(Sf0);
 108
 109             scalar magSfDiff = mag(Sf - Sf0);
 110
 111             if (magSfDiff > SMALL)
 112             {
 113                 bool pureRotation = false;
 114
 115                 if (pureRotation)
 116                 {
 117                     // rotation
 118
 119                     // find center of rotation
 120                     vector omega = Sf0 ^ Sf;
 121                     scalar magOmega = mag(omega);
 122                     omega /= magOmega+SMALL;
 123                     vector n0 = omega ^ Sf0;
 124                     scalar lam = ((Cf1 - Cf0) & Sf)/(n0 & Sf);
 125                     vector r0 = Cf0 + lam*n0;
 126
 127                     scalar phiVel = ::asin(magOmega)/dt;
 128
 129                     vector pos = p.position() - r0;
 130                     vector v = phiVel*(omega ^ pos);
 131
 132                     vector Sfp = Sf0 + fraction*(Sf - Sf0);
 133                     //vector Sfp = Sf;
 134
 135                     vector Ur = p.U() - v;
 136                     scalar Urn = Ur & Sfp;
 137                     /*
 138                     scalar dd = (p.position() - r0) & Sfp;
 139                     Info << "Urn = " << Urn
 140                         << ", dd = " << dd
 141                         << ", pos = " << p.position()
 142                         << ", Sfp = " << Sfp
 143                         << ", omega = " << omega
 144                         << ", r0 = " << r0
 145                         << endl;
 146                         */
 147                     if (Urn > 0.0)
 148                     {
 149                         p.U() -= (1.0 + elasticity_)*Urn*Sfp;
 150                     }
 151                 }
 152                 else
 153                 {
 154                     vector Sfp = Sf0 + fraction*(Sf - Sf0);
 155                     //vector Sfp = Sf;
 156                     vector omega = Sf0 ^ Sf;
 157                     scalar magOmega = mag(omega);
 158                     omega /= magOmega+SMALL;
 159
 160                     scalar phiVel = ::asin(magOmega)/dt;
 161
 162                     scalar dist = (p.position() - Cf) & Sfp;
 163                     vector pos = p.position() - dist*Sfp;
 164                     vector vrot = phiVel*(omega ^ (pos - Cf));
 165                     vector vtrans = (Cf1 - Cf0)/dt;
 166                     vector v = vtrans + vrot;
 167
 168                     scalar Un = ((p.U() - v) & Sfp);
 169
 170                     if (Un > 0.0)
 171                     {
 172                         p.U() -= (1.0 + elasticity_)*Un*Sfp;
 173                     }
 174                 }
 175             }
 176             else
 177             {
 178                 // translation
 179
 180                 vector v = (Cf1 - Cf0)/dt;
 181                 vector Ur = p.U() - v;
 182                 scalar Urn = Ur & Sf;
 183
 184                 if (Urn > 0.0)
 185                 {
 186                     Ur -= (1.0 + elasticity_)*Urn*Sf;
 187                     p.U() = Ur + v;
 188                 }
 189
 190
 191             }
 192         }
 193
 194     }
 195     else
 196     {
 197         FatalError
 198             << "bool reflectParcel::wallTreatment(parcel& parcel) const "
 199                 << " parcel has hit a boundary "
 200                 << mesh_.boundary()[patchi].type()
 201                 << " which not yet has been implemented."
 202             << abort(FatalError);
 203     }
 204     return true;
 205 }
 206
 207 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 208
 209 } // End namespace Foam
 210
 211 // ************************************************************************* //