yade-users team mailing list archive

[Idan <question706221@xxxxxxxxxxxxxxxxxxxxx>]

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.