← Back to branch summary

~ubuntu-branches/ubuntu/trusty/yade/trusty

« back to all changes in this revision

Viewing changes to pkg/dem/HertzMindlin.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky, cf3f8d9
  • Date: 2013-10-30 20:56:33 UTC
  • mfrom: (20.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20131030205633-1f01r7hjce17d723
Tags: 1.05.0-2
http://bugs.debian.org/726009
[cf3f8d9] Pass -ftrack-macro-expansion=0 only if gcc>=4.8. (Closes: #726009)

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/dem/HertzMindlin.cpp
31
31
 
32
32
void Ip2_FrictMat_FrictMat_MindlinPhys::go(const shared_ptr<Material>& b1,const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){
33
33
        if(interaction->phys) return; // no updates of an already existing contact necessary
34
 
        shared_ptr<MindlinPhys> mindlinPhys(new MindlinPhys());
35
 
        interaction->phys = mindlinPhys;
 
34
        shared_ptr<MindlinPhys> contactPhysics(new MindlinPhys());
 
35
        interaction->phys = contactPhysics;
36
36
        FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
37
37
        FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
38
38
        
64
64
        Real Rmean = (Da+Db)/2.; // mean radius
65
65
        Real Kno = 4./3.*E*sqrt(R); // coefficient for normal stiffness
66
66
        Real Kso = 2*sqrt(4*R)*G/(2-V); // coefficient for shear stiffness
67
 
        Real frictionAngle = std::min(fa,fb);
 
67
        Real frictionAngle = (!frictAngle) ? std::min(fa,fb) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle);
68
68
 
69
69
        Real Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory
70
70
 
71
71
        /* pass values calculated from above to MindlinPhys */
72
 
        mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle); 
73
 
        //mindlinPhys->prevNormal = scg->normal; // used to compute relative rotation
74
 
        mindlinPhys->kno = Kno; // this is just a coeff
75
 
        mindlinPhys->kso = Kso; // this is just a coeff
76
 
        mindlinPhys->adhesionForce = Adhesion;
 
72
        contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle); 
 
73
        //contactPhysics->prevNormal = scg->normal; // used to compute relative rotation
 
74
        contactPhysics->kno = Kno; // this is just a coeff
 
75
        contactPhysics->kso = Kso; // this is just a coeff
 
76
        contactPhysics->adhesionForce = Adhesion;
77
77
        
78
 
        mindlinPhys->kr = krot;
79
 
        mindlinPhys->ktw = ktwist;
80
 
        mindlinPhys->maxBendPl = eta*Rmean; // does this make sense? why do we take Rmean?
 
78
        contactPhysics->kr = krot;
 
79
        contactPhysics->ktw = ktwist;
 
80
        contactPhysics->maxBendPl = eta*Rmean; // does this make sense? why do we take Rmean?
81
81
 
82
82
        /* compute viscous coefficients */
83
83
        if(en && betan) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinPhys: only one of en, betan can be specified.");
86
86
        // en or es specified, just compute alpha, otherwise alpha remains 0
87
87
        if(en || es){
88
88
                Real logE = log((*en)(mat1->id,mat2->id));
89
 
                mindlinPhys->alpha = -sqrt(5/6.)*2*logE/sqrt(pow(logE,2)+pow(Mathr::PI,2))*sqrt(2*E*sqrt(R)); // (see Tsuji, 1992)
 
89
                contactPhysics->alpha = -sqrt(5/6.)*2*logE/sqrt(pow(logE,2)+pow(Mathr::PI,2))*sqrt(2*E*sqrt(R)); // (see Tsuji, 1992), also [Antypov2011] eq. 17
90
90
        }
91
91
        
92
92
        // betan specified, use that value directly; otherwise give zero
93
93
        else{   
94
 
                mindlinPhys->betan=betan ? (*betan)(mat1->id,mat2->id) : 0; 
95
 
                mindlinPhys->betas=betas ? (*betas)(mat1->id,mat2->id) : mindlinPhys->betan;
 
94
                contactPhysics->betan=betan ? (*betan)(mat1->id,mat2->id) : 0; 
 
95
                contactPhysics->betas=betas ? (*betas)(mat1->id,mat2->id) : contactPhysics->betan;
96
96
        }
97
97
}
98
98
 
319
319
        }
320
320
        else if (useDamping){ // (see Tsuji, 1992)
321
321
                Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body
322
 
                cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient
 
322
                cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient, see also [Antypov2011] eq. 10
323
323
                cs = cn; // same value for shear viscous coefficient
324
324
        }
325
325
 
346
346
        /**************************************/
347
347
        
348
348
        // normal force must be updated here before we apply the Mohr-Coulomb criterion
349
 
        if (useDamping){ // get normal viscous component 
 
349
        if (useDamping){ // get normal viscous component
350
350
                phys->normalViscous = cn*incidentVn;
351
 
                // add normal viscous component if damping is included
352
 
                phys->normalForce -= phys->normalViscous;
 
351
                Vector3r normTemp = phys->normalForce - phys->normalViscous; // temporary normal force
 
352
                // viscous force should not exceed the value of current normal force, i.e. no attraction force should be permitted if particles are non-adhesive
 
353
                // if particles are adhesive, then fixed the viscous force at maximum equal to the adhesion force
 
354
                // *** enforce normal force to zero if no adhesion is permitted ***
 
355
                if (phys->adhesionForce==0.0 || !includeAdhesion){
 
356
                                                if (normTemp.dot(scg->normal)<0.0){
 
357
                                                                                phys->normalForce = Vector3r::Zero();
 
358
                                                                                phys->normalViscous = phys->normalViscous + normTemp; // normal viscous force is such that the total applied force is null - it is necessary to compute energy correctly!
 
359
                                                }
 
360
                                                else{phys->normalForce -= phys->normalViscous;}
 
361
                }
 
362
                else if (includeAdhesion && phys->adhesionForce!=0.0){
 
363
                                                // *** limit viscous component to the max adhesive force ***
 
364
                                                if (normTemp.dot(scg->normal)<0.0 && (phys->normalViscous.norm() > phys->adhesionForce) ){
 
365
                                                                                Real normVisc = phys->normalViscous.norm(); Vector3r normViscVector = phys->normalViscous/normVisc;
 
366
                                                                                phys->normalViscous = phys->adhesionForce*normViscVector;
 
367
                                                                                phys->normalForce -= phys->normalViscous;
 
368
                                                }
 
369
                                                // *** apply viscous component - in the presence of adhesion ***
 
370
                                                else {phys->normalForce -= phys->normalViscous;}
 
371
                }
353
372
                if (calcEnergy) {normDampDissip += phys->normalViscous.dot(incidentVn*dt);} // calc dissipation of energy due to normal damping
354
373
        }
 
374
        
355
375
 
356
376
        /*************************************/
357
377
        /* SHEAR DISPLACEMENT (elastic only) */