src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-2011 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 "PairSpringSliderDashpot.H"
  27
  28 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  29
  30 template<class CloudType>
  31 void Foam::PairSpringSliderDashpot<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     // then rMin into minimum R,
  69     //     1/RMin = 1/rMin + 1/rMin,
  70     //     RMin = rMin/2 = dMin/4
  71
  72     RMin /= 4.0;
  73
  74     // Multiply by two to create the worst-case relative velocity
  75
  76     UMagMax = 2*UMagMax;
  77 }
  78
  79
  80 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  81
  82 template<class CloudType>
  83 Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
  84 (
  85     const dictionary& dict,
  86     CloudType& cloud
  87 )
  88 :
  89     PairModel<CloudType>(dict, cloud, typeName),
  90     Estar_(),
  91     Gstar_(),
  92     alpha_(readScalar(this->coeffDict().lookup("alpha"))),
  93     b_(readScalar(this->coeffDict().lookup("b"))),
  94     mu_(readScalar(this->coeffDict().lookup("mu"))),
  95     cohesionEnergyDensity_
  96     (
  97         readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
  98     ),
  99     cohesion_(false),
 100     collisionResolutionSteps_
 101     (
 102         readScalar(this->coeffDict().lookup("collisionResolutionSteps"))
 103     ),
 104     volumeFactor_(1.0),
 105     useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
 106 {
 107     if (useEquivalentSize_)
 108     {
 109         volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
 110     }
 111
 112     scalar nu = this->owner().constProps().poissonsRatio();
 113
 114     scalar E = this->owner().constProps().youngsModulus();
 115
 116     Estar_ = E/(2.0*(1.0 - sqr(nu)));
 117
 118     scalar G = E/(2.0*(1.0 + nu));
 119
 120     Gstar_ = G/(2.0*(2.0 - nu));
 121
 122     cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
 123 }
 124
 125
 126 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 127
 128 template<class CloudType>
 129 Foam::PairSpringSliderDashpot<CloudType>::~PairSpringSliderDashpot()
 130 {}
 131
 132
 133 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 134
 135 template<class CloudType>
 136 bool Foam::PairSpringSliderDashpot<CloudType>::controlsTimestep() const
 137 {
 138     return true;
 139 }
 140
 141
 142 template<class CloudType>
 143 Foam::label Foam::PairSpringSliderDashpot<CloudType>::nSubCycles() const
 144 {
 145     if (!(this->owner().size()))
 146     {
 147         return 1;
 148     }
 149
 150     scalar RMin;
 151     scalar rhoMax;
 152     scalar UMagMax;
 153
 154     findMinMaxProperties(RMin, rhoMax, UMagMax);
 155
 156     // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
 157     scalar minCollisionDeltaT =
 158         5.429675
 159        *RMin
 160        *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
 161        /collisionResolutionSteps_;
 162
 163     return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
 164 }
 165
 166
 167 template<class CloudType>
 168 void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
 169 (
 170     typename CloudType::parcelType& pA,
 171     typename CloudType::parcelType& pB
 172 ) const
 173 {
 174     vector r_AB = (pA.position() - pB.position());
 175
 176     scalar dAEff = pA.d();
 177
 178     if (useEquivalentSize_)
 179     {
 180         dAEff *= cbrt(pA.nParticle()*volumeFactor_);
 181     }
 182
 183     scalar dBEff = pB.d();
 184
 185     if (useEquivalentSize_)
 186     {
 187         dBEff *= cbrt(pB.nParticle()*volumeFactor_);
 188     }
 189
 190     scalar r_AB_mag = mag(r_AB);
 191
 192     scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
 193
 194     if (normalOverlapMag > 0)
 195     {
 196         //Particles in collision
 197
 198         vector rHat_AB = r_AB/(r_AB_mag + VSMALL);
 199
 200         vector U_AB = pA.U() - pB.U();
 201
 202         // Effective radius
 203         scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);
 204
 205         // Effective mass
 206         scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
 207
 208         scalar kN = (4.0/3.0)*sqrt(R)*Estar_;
 209
 210         scalar etaN = alpha_*sqrt(M*kN)*pow025(normalOverlapMag);
 211
 212         // Normal force
 213         vector fN_AB =
 214             rHat_AB
 215            *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
 216
 217         // Cohesion force
 218         if (cohesion_)
 219         {
 220             fN_AB +=
 221                 -cohesionEnergyDensity_
 222                 *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
 223                 *rHat_AB;
 224         }
 225
 226         pA.f() += fN_AB;
 227         pB.f() += -fN_AB;
 228
 229         vector USlip_AB =
 230             U_AB - (U_AB & rHat_AB)*rHat_AB
 231           + (pA.omega() ^ (dAEff/2*-rHat_AB))
 232           - (pB.omega() ^ (dBEff/2*rHat_AB));
 233
 234         scalar deltaT = this->owner().mesh().time().deltaTValue();
 235
 236         vector& tangentialOverlap_AB =
 237             pA.collisionRecords().matchPairRecord
 238             (
 239                 pB.origProc(),
 240                 pB.origId()
 241             ).collisionData();
 242
 243         vector& tangentialOverlap_BA =
 244             pB.collisionRecords().matchPairRecord
 245             (
 246                 pA.origProc(),
 247                 pA.origId()
 248             ).collisionData();
 249
 250         vector deltaTangentialOverlap_AB = USlip_AB*deltaT;
 251
 252         tangentialOverlap_AB += deltaTangentialOverlap_AB;
 253         tangentialOverlap_BA += -deltaTangentialOverlap_AB;
 254
 255         scalar tangentialOverlapMag = mag(tangentialOverlap_AB);
 256
 257         if (tangentialOverlapMag > VSMALL)
 258         {
 259             scalar kT = 8.0*sqrt(R*normalOverlapMag)*Gstar_;
 260
 261             scalar etaT = etaN;
 262
 263             // Tangential force
 264             vector fT_AB;
 265
 266             if (kT*tangentialOverlapMag > mu_*mag(fN_AB))
 267             {
 268                 // Tangential force greater than sliding friction,
 269                 // particle slips
 270
 271                 fT_AB = -mu_*mag(fN_AB)*USlip_AB/mag(USlip_AB);
 272
 273                 tangentialOverlap_AB = vector::zero;
 274                 tangentialOverlap_BA = vector::zero;
 275             }
 276             else
 277             {
 278                 fT_AB =
 279                     -kT*tangentialOverlapMag
 280                    *tangentialOverlap_AB/tangentialOverlapMag
 281                   - etaT*USlip_AB;
 282             }
 283
 284             pA.f() += fT_AB;
 285             pB.f() += -fT_AB;
 286
 287             pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
 288             pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
 289         }
 290     }
 291 }
 292
 293
 294 // ************************************************************************* //