src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-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 "WallSpringSliderDashpot.H"
  27
  28 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  29
  30 template<class CloudType>
  31 void Foam::WallSpringSliderDashpot<CloudType>::findMinMaxProperties
  32 (
  33     scalar& rMin,
  34     scalar& rhoMax,
  35     scalar& UMagMax
  36 ) const
  37 {
  38     rMin = VGREAT;
  39     rhoMax = -VGREAT;
  40     UMagMax = -VGREAT;
  41
  42     forAllConstIter(typename CloudType, this->owner(), iter)
  43     {
  44         const typename CloudType::parcelType& p = iter();
  45
  46         // Finding minimum diameter to avoid excessive arithmetic
  47
  48         scalar dEff = p.d();
  49
  50         if (useEquivalentSize_)
  51         {
  52             dEff *= cbrt(p.nParticle()*volumeFactor_);
  53         }
  54
  55         rMin = min(dEff, rMin);
  56
  57         rhoMax = max(p.rho(), rhoMax);
  58
  59         UMagMax = max
  60         (
  61             mag(p.U()) + mag(p.omega())*dEff/2,
  62             UMagMax
  63         );
  64     }
  65
  66     // Transform the minimum diameter into minimum radius
  67     //     rMin = dMin/2
  68
  69     rMin /= 2.0;
  70 }
  71
  72
  73 template<class CloudType>
  74 void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
  75 (
  76     typename CloudType::parcelType& p,
  77     const point& site,
  78     const WallSiteData<vector>& data,
  79     scalar pREff,
  80     scalar kN
  81 ) const
  82 {
  83     vector r_PW = p.position() - site;
  84
  85     vector U_PW = p.U() - data.wallData();
  86
  87     scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
  88
  89     vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
  90
  91     scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
  92
  93     vector fN_PW =
  94         rHat_PW
  95        *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
  96
  97     p.f() += fN_PW;
  98
  99     vector USlip_PW =
 100         U_PW - (U_PW & rHat_PW)*rHat_PW
 101       + (p.omega() ^ (pREff*-rHat_PW));
 102
 103     scalar deltaT = this->owner().mesh().time().deltaTValue();
 104
 105     vector& tangentialOverlap_PW =
 106         p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();
 107
 108     tangentialOverlap_PW += USlip_PW*deltaT;
 109
 110     scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
 111
 112     if (tangentialOverlapMag > VSMALL)
 113     {
 114         scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar_;
 115
 116         scalar etaT = etaN;
 117
 118         // Tangential force
 119         vector fT_PW;
 120
 121         if (kT*tangentialOverlapMag > mu_*mag(fN_PW))
 122         {
 123             // Tangential force greater than sliding friction,
 124             // particle slips
 125
 126             fT_PW = -mu_*mag(fN_PW)*USlip_PW/mag(USlip_PW);
 127
 128             tangentialOverlap_PW = vector::zero;
 129         }
 130         else
 131         {
 132             fT_PW =
 133                 -kT*tangentialOverlapMag
 134                *tangentialOverlap_PW/tangentialOverlapMag
 135               - etaT*USlip_PW;
 136         }
 137
 138         p.f() += fT_PW;
 139
 140         p.torque() += (pREff*-rHat_PW) ^ fT_PW;
 141     }
 142 }
 143
 144
 145 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 146
 147 template<class CloudType>
 148 Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
 149 (
 150     const dictionary& dict,
 151     CloudType& cloud
 152 )
 153 :
 154     WallModel<CloudType>(dict, cloud, typeName),
 155     Estar_(),
 156     Gstar_(),
 157     alpha_(readScalar(this->coeffDict().lookup("alpha"))),
 158     b_(readScalar(this->coeffDict().lookup("b"))),
 159     mu_(readScalar(this->coeffDict().lookup("mu"))),
 160     collisionResolutionSteps_
 161     (
 162         readScalar
 163         (
 164             this->coeffDict().lookup("collisionResolutionSteps")
 165         )
 166     ),
 167     volumeFactor_(1.0),
 168     useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
 169 {
 170     if (useEquivalentSize_)
 171     {
 172         volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
 173     }
 174
 175     scalar nu = readScalar(this->coeffDict().lookup("poissonsRatio"));
 176
 177     scalar E = readScalar(this->coeffDict().lookup("youngsModulus"));
 178
 179     scalar pNu = this->owner().constProps().poissonsRatio();
 180
 181     scalar pE = this->owner().constProps().youngsModulus();
 182
 183     Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
 184
 185     Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
 186 }
 187
 188
 189 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 190
 191 template<class CloudType>
 192 Foam::WallSpringSliderDashpot<CloudType>::~WallSpringSliderDashpot()
 193 {}
 194
 195
 196 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 197
 198 template<class CloudType>
 199 Foam::scalar Foam::WallSpringSliderDashpot<CloudType>::pREff
 200 (
 201     const typename CloudType::parcelType& p
 202 ) const
 203 {
 204     if (useEquivalentSize_)
 205     {
 206         return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
 207     }
 208     else
 209     {
 210         return p.d()/2;
 211     }
 212 }
 213
 214
 215 template<class CloudType>
 216 bool Foam::WallSpringSliderDashpot<CloudType>::controlsTimestep() const
 217 {
 218     return true;
 219 }
 220
 221
 222 template<class CloudType>
 223 Foam::label Foam::WallSpringSliderDashpot<CloudType>::nSubCycles() const
 224 {
 225     if (!(this->owner().size()))
 226     {
 227         return 1;
 228     }
 229
 230     scalar rMin;
 231     scalar rhoMax;
 232     scalar UMagMax;
 233
 234     findMinMaxProperties(rMin, rhoMax, UMagMax);
 235
 236     // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
 237     scalar minCollisionDeltaT =
 238         5.429675
 239        *rMin
 240        *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
 241        /collisionResolutionSteps_;
 242
 243     return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
 244 }
 245
 246
 247 template<class CloudType>
 248 void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
 249 (
 250     typename CloudType::parcelType& p,
 251     const List<point>& flatSitePoints,
 252     const List<WallSiteData<vector> >& flatSiteData,
 253     const List<point>& sharpSitePoints,
 254     const List<WallSiteData<vector> >& sharpSiteData
 255 ) const
 256 {
 257     scalar pREff = this->pREff(p);
 258
 259     scalar kN = (4.0/3.0)*sqrt(pREff)*Estar_;
 260
 261     forAll(flatSitePoints, siteI)
 262     {
 263         evaluateWall
 264         (
 265             p,
 266             flatSitePoints[siteI],
 267             flatSiteData[siteI],
 268             pREff,
 269             kN
 270         );
 271     }
 272
 273     forAll(sharpSitePoints, siteI)
 274     {
 275         // Treating sharp sites like flat sites
 276
 277         evaluateWall
 278         (
 279             p,
 280             sharpSitePoints[siteI],
 281             sharpSiteData[siteI],
 282             pREff,
 283             kN
 284         );
 285     }
 286 }
 287
 288
 289 // ************************************************************************* //