Question #706221 on Yade changed:
https://answers.launchpad.net/yade/+question/706221
Idan gave more information on the question:
Hi,
I asked it to continue implementing the model with exponential
temperature dependence (it was lacking the nature of temperature
dependence as it does not exist in the original burgers formulation).
Here is the result:
Header:
#pragma once
#include<yade/pkg/dem/CohesiveFrictionalContactLaw.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/pkg/dem/FrictPhys.hpp>
#include<yade/pkg/dem/ThermalState.hpp>
class TempDependentBurgersViscoElasticModel: public CohFrictMat {
public:
TempDependentBurgersViscoElasticModel(Real young, Real poisson, Real
frictionAngle, Real density, Real cohesion, Real viscousDamping, Real
viscosity, Real tempRef, Real tempExpFactor);
virtual ~TempDependentBurgersViscoElasticModel();
// Temperature reference for the model
Real tempRef;
// Temperature exponential factor
Real tempExpFactor;
// Viscosity for the viscoelastic component of the model
Real viscosity;
YADE_CLASS_BASE_DOC_ATTRS_CTOR(TempDependentBurgersViscoElasticModel,
CohFrictMat, "Temperature-dependent Burger's viscoelastic model inheriting from
CohFrictMat",
((Real, viscousDamping, 0, ,"Viscous damping for the viscoelastic
component of the model"))
((Real, tempRef, 298.15, ,"Reference temperature for the model"))
((Real, viscosity, 0, ,"Viscosity for the viscoelastic component of the
model"))
((Real, tempExpFactor, 0, ,"Temperature exponential factor for the
temperature-dependent tangential force")),
/* ctor */ TempDependentBurgersViscoElasticModel(young, poisson,
frictionAngle, density, cohesion, viscousDamping, viscosity, tempRef,
tempExpFactor); );
};
REGISTER_SERIALIZABLE(TempDependentBurgersViscoElasticModel);
class TempDependentBurgersViscoElasticContactLaw: public
Law2_ScGeom_FrictPhys_CundallStrack {
public:
virtual void go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction*
contact);
YADE_CLASS_BASE_DOC(TempDependentBurgersViscoElasticContactLaw,
Law2_ScGeom_FrictPhys_CundallStrack, "Contact law for the temperature-dependent
Burger's viscoelastic model");
};
REGISTER_SERIALIZABLE(TempDependentBurgersViscoElasticContactLaw);
CCP file:
#include "TempDependentBurgersViscoElastic.hpp"
TempDependentBurgersViscoElasticModel::TempDependentBurgersViscoElasticModel(Real
young, Real poisson, Real frictionAngle, Real density, Real cohesion, Real
viscousDamping, Real viscosity, Real tempRef, Real tempExpFactor)
: CohFrictMat(young, poisson, frictionAngle, density, cohesion)
, viscousDamping(viscousDamping)
, viscosity(viscosity)
, tempRef(tempRef)
, tempExpFactor(tempExpFactor) {
}
TempDependentBurgersViscoElasticModel::~TempDependentBurgersViscoElasticModel()
{
}
YADE_PLUGIN((TempDependentBurgersViscoElasticModel));
void TempDependentBurgersViscoElasticContactLaw::go(shared_ptr<IGeom>& ig,
shared_ptr<IPhys>& ip, Interaction* contact) {
// Cast the geometry and physical parameters to the correct types
ScGeom* geom = static_cast<ScGeom*>(ig.get());
FrictPhys* phys = static_cast<FrictPhys*>(ip.get());
// Get the associated bodies and materials
Body* b1 = Body::byId(contact->getId1(), scene);
Body* b2 = Body::byId(contact->getId2(), scene);
TempDependentBurgersViscoElasticModel* mat1 =
static_cast<TempDependentBurgersViscoElasticModel*>(b1->material.get());
TempDependentBurgersViscoElasticModel* mat2 =
static_cast<TempDependentBurgersViscoElasticModel*>(b2->material.get());
// Compute the normal force based on Hooke's law
Real kn = 2 * phys->kn * geom->radius / (geom->radius + geom->radius);
Real Fn = kn * geom->penetrationDepth;
// Compute the tangential force based on a simplified version of the Burger's
model
Real kt = 2 * phys->kt * geom->radius / (geom->radius + geom->radius);
Real viscosity = 0.5 * (mat1->viscosity + mat2->viscosity);
Real viscousDamping = 0.5 * (mat1->viscousDamping + mat2->viscousDamping);
Vector3r relVel = b1->state->vel - b2->state->vel;
Vector3r relAngVel = b1->state->angVel - b2->state->angVel;
Vector3r relTangentialVel = relVel - geom->normal.dot(relVel) * geom->normal;
Vector3r relTangentialAngVel = geom->radius * 0.5 *
(relAngVel.cross(geom->normal));
Vector3r relTangentialTotalVel = relTangentialVel + relTangentialAngVel;
Vector3r tangentialForce = -(kt * geom->uT + viscosity * relTangentialTotalVel);
// Get the temperatures from the thermal state
shared_ptr<ThermalState> ts1 = YADE_PTR_CAST<ThermalState>(b1->state);
shared_ptr<ThermalState> ts2 = YADE_PTR_CAST<ThermalState>(b2->state);
Real avgTemp = 0.5 * (ts1->temp + ts2->temp);
// Apply temperature-dependent exponential factor
Real tempExpFactor = 0.5 * (mat1->tempExpFactor + mat2->tempExpFactor);
Real tempRef = 0.5 * (mat1->tempRef + mat2->tempRef);
Real tempFactor = std::exp(tempExpFactor * (avgTemp - tempRef));
tangentialForce *= tempFactor;
// Apply Coulomb friction law
Real maxTangentialForce = phys->tangensOfFrictionAngle * Fn;
if (tangentialForce.norm() > maxTangentialForce) {
tangentialForce = tangentialForce.normalized() * maxTangentialForce;
}
Vector3r normalForce = geom->normal * Fn;
// Apply the forces to the bodies
scene->forces.addForce(contact->getId1(), normalForce + tangentialForce);
scene->forces.addForce(contact->getId2(), -(normalForce + tangentialForce));
scene->forces.addTorque(contact->getId1(), geom->radius * 0.5 *
tangentialForce.cross(geom->normal));
scene->forces.addTorque(contact->getId2(), -geom->radius * 0.5 *
tangentialForce.cross(geom->normal));
geom->uT = geom->uT - geom->rotate(tangentialForce / kt);
}
YADE_PLUGIN((TempDependentBurgersViscoElasticContactLaw));
Hope it is more helpful now.
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
