New question #261724 on Yade:
https://answers.launchpad.net/yade/+question/261724
Hi guys,
I created a new class of material to present viscoelasticity by Burgers model.
This class of material has been derived from CohFrictMat and it's called
CohBurgersMat. I verified the response of CohBurgersMat-CohBurgersMat
interaction by comparing the load-displacement response of analytical
equations.
However, I need to have FrictMat-CohBurgersMat interaction as well. To do so, I
created a new Ip2 (Ip2_FrictMat_CohBurgersMat_CohBurgersPhys).
However, the model fails to create the required interaction physics. Can you
please have a look on my script? Thanks.
======== CohBurgers.cpp ==============================================
// 2014 © Behzad Majidi <behzad.maj...@gmail.com>
// 2014 © Ricardo Pieralisi <ricpieral...@gmail.com>
#include"CohBurgers.hpp"
#include<core/State.hpp>
#include<pkg/dem/ScGeom.hpp>
#include<core/Omega.hpp>
#include<core/Scene.hpp>
#include<pkg/dem/Shop.hpp>
#include<pkg/common/InteractionLoop.hpp>
#include<pkg/common/Sphere.hpp>
#include<pkg/common/Facet.hpp>
#include<pkg/common/Wall.hpp>
YADE_PLUGIN((CohBurgersMat)(CohBurgersPhys)(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys)(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys)(Law2_ScGeom_CohBurgersPhys_CohesiveBurgers));
//********************** Ip2_FrictMat_CohBurgersMat_CohBurgersPhys
****************************/
void Ip2_FrictMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>&
b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction)
{
// if(interaction->phys) return;
shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys());
Calculate_FrictMat_CohBurgersMat_CohBurgersPhys(b1, b2, interaction, phys);
interaction->phys = phys;
phys->cohesionBroken = false;
phys->initCohesion = true;
cohesionDefinitionIteration = -1;
setCohesionNow = true;
}
//********************** Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys
****************************/
void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::go(const
shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const
shared_ptr<Interaction>& interaction)
{
if(interaction->phys) return;
shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys());
Calculate_CohBurgersMat_CohBurgersMat_CohBurgersPhys(b1, b2, interaction,
phys);
interaction->phys = phys;
phys->cohesionBroken = false;
phys->initCohesion = true;
cohesionDefinitionIteration = -1;
setCohesionNow = true;
}
//********************** CohesiveBurgersContactLaw ****************************/
bool Law2_ScGeom_CohBurgersPhys_CohesiveBurgers::go(shared_ptr<IGeom>& _geom,
shared_ptr<IPhys>& _phys, Interaction* I)
{
Vector3r force = Vector3r::Zero();
Vector3r torque1 = Vector3r::Zero();
Vector3r torque2 = Vector3r::Zero();
CohBurgersPhys& phys=*static_cast<CohBurgersPhys*>(_phys.get());
if (computeForceTorqueCohBurgers(_geom, _phys, I, force, torque1, torque2)
and (I->isActive)){
const int id1 = I->getId1();
const int id2 = I->getId2();
addForce (id1,-force,scene);
addForce (id2, force,scene);
addTorque(id1, torque1,scene);
addTorque(id2, torque2,scene);
return true;
} else return false;
}
/*!Calculate_CohBurgersPhys (1)*/
void
Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::Calculate_CohBurgersMat_CohBurgersMat_CohBurgersPhys(const
shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const
shared_ptr<Interaction>& interaction, shared_ptr<CohBurgersPhys> phys) {
if(interaction->phys) return;
CohBurgersMat* mat1 = static_cast<CohBurgersMat*>(b1.get());
CohBurgersMat* mat2 = static_cast<CohBurgersMat*>(b2.get());
Real mass1 = 1.0;
Real mass2 = 1.0;
mass1=Body::byId(interaction->getId1())->state->mass;
mass2=Body::byId(interaction->getId2())->state->mass;
//Load the rheological parameters
Real kmn1 = (mat1->kmn)*2; Real kkn1 = (mat1->kkn)*2;
Real cmn1 = (mat1->cmn)*2; Real ckn1 = (mat1->ckn)*2;
Real kms1 = (mat1->kms)*2; Real kks1 = (mat1->kks)*2;
Real cms1 = (mat1->cms)*2; Real cks1 = (mat1->cks)*2;
Real kmn2 = (mat2->kmn)*2; Real kkn2 = (mat2->kkn)*2;
Real cmn2 = (mat2->cmn)*2; Real ckn2 = (mat2->ckn)*2;
Real kms2 = (mat2->kms)*2; Real kks2 = (mat2->kks)*2;
Real cms2 = (mat2->cms)*2; Real cks2 = (mat2->cks)*2;
phys->kmn = contactParameterCalculation(kmn1,kmn2);
phys->kms = contactParameterCalculation(kms1,kms2);
phys->cmn = contactParameterCalculation(cmn1,cmn2);
phys->cms = contactParameterCalculation(cms1,cms2);
phys->kkn = contactParameterCalculation(kkn1,kkn2);
phys->kks = contactParameterCalculation(kks1,kks2);
phys->ckn = contactParameterCalculation(ckn1,ckn2);
phys->cks = contactParameterCalculation(cks1,cks2);
phys->normalAdhesion= (mat1->normalCohesion);
interaction->phys = shared_ptr<CohBurgersPhys>(phys);
}
/*!Calculate_CohBurgersPhys (2)*/
void
Ip2_FrictMat_CohBurgersMat_CohBurgersPhys::Calculate_FrictMat_CohBurgersMat_CohBurgersPhys(const
shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const
shared_ptr<Interaction>& interaction, shared_ptr<CohBurgersPhys> phys) {
if(interaction->phys) return;
const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
const shared_ptr<CohBurgersMat>& mat2 = YADE_PTR_CAST<CohBurgersMat>(b2);
Real mass1 = 1.0;
Real mass2 = 1.0;
mass1=Body::byId(interaction->getId1())->state->mass;
mass2=Body::byId(interaction->getId2())->state->mass;
//Load the rheological parameters
Real kn1 = (mat1->young); Real ks1 = (mat1->young)*(mat1->poisson);
Real kmn2 = (mat2->kmn); Real kkn2 = (mat2->kkn);
Real cmn2 = (mat2->cmn); Real ckn2 = (mat2->ckn);
Real kms2 = (mat2->kms); Real kks2 = (mat2->kks);
Real cms2 = (mat2->cms); Real cks2 = (mat2->cks);
phys->kmn = contactParameterCalculation(kn1,kmn2);
phys->kms = contactParameterCalculation(ks1,kms2);
phys->cmn = cmn2;
phys->cms = cms2;
phys->kkn = kkn2;
phys->kks = kks2;
phys->ckn = ckn2;
phys->cks = cks2;
phys->normalAdhesion= ((mat2->normalCohesion));
interaction->phys = shared_ptr<CohBurgersPhys>(phys);
}
/*! computeForceTorqueBurgers */
bool computeForceTorqueCohBurgers(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>&
_phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r &
torque2) {
const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
CohBurgersPhys& phys=*static_cast<CohBurgersPhys*>(_phys.get());
Scene* scene=Omega::instance().getScene().get();
const int id1 = I->getId1();
const int id2 = I->getId2();
const BodyContainer& bodies = *scene->bodies;
const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
const Real& dt = scene->dt;
if (phys.normalForce.norm() > phys.normalAdhesion) {
return false;
}
if (I->isFresh(scene)) {
phys.ukn = Vector3r(0,0,0);
phys.uks = Vector3r(0,0,0);
phys.shearForce = Vector3r(0,0,0);
phys.normalForce = phys.kmn * geom.penetrationDepth * geom.normal;
} else {
Real An = 1.0 + phys.kkn * dt / (2.0 * phys.ckn);
Real Bn = 1.0 - phys.kkn * dt / (2.0 * phys.ckn);
Real Cn = dt / (2.0 * phys.ckn * An) + 1.0 / phys.kmn + dt / (2.0 *
phys.cmn);
Real Dn = dt / (2.0 * phys.ckn * An) - 1.0 / phys.kmn + dt / (2.0 *
phys.cmn);
Real As = 1.0 + phys.kks * dt / (2.0 * phys.cks);
Real Bs = 1.0 - phys.kks * dt / (2.0 * phys.cks);
Real Cs = dt / (2.0 * phys.cks * As) + 1.0 / phys.kms + dt / (2.0 *
phys.cms);
Real Ds = dt / (2.0 * phys.cks * As) - 1.0 / phys.kms + dt / (2.0 *
phys.cms);
const Vector3r shift2 = scene->isPeriodic ?
scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero();
const Vector3r shiftVel = scene->isPeriodic ?
scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero();
const Vector3r c1x = (geom.contactPoint - de1.pos);
const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) -
(de2.vel+de2.angVel.cross(c2x)) + shiftVel;
const Vector3r normalVelocity = geom.normal.dot(relativeVelocity) *
geom.normal;
const Vector3r shearVelocity = relativeVelocity-normalVelocity;
const Vector3r normalForceT = (normalVelocity * dt + phys.ukn * (1.0 - Bn
/ An) - phys.normalForce * Dn) / Cn;
phys.ukn = (phys.ukn * Bn + (normalForceT + phys.normalForce) * dt / (2.0
* phys.ckn)) / An;
phys.normalForce = normalForceT;
const Vector3r shearForceT = -1 * (shearVelocity * dt + phys.uks * (1.0 -
Bs / As) - phys.shearForce * Ds) / Cs;
phys.uks = (phys.uks * Bs + (shearForceT + phys.shearForce) * dt / (2.0 *
phys.cks)) / As;
phys.shearForce = shearForceT;
const Real maxFs = phys.normalForce.squaredNorm() *
std::pow(phys.tangensOfFrictionAngle,2);
if( phys.shearForce.squaredNorm() > maxFs ) {
const Real ratio = sqrt(maxFs) / phys.shearForce.norm();
phys.shearForce *= ratio;
}
force = phys.normalForce + phys.shearForce;
torque1 = -c1x.cross(force);
torque2 = c2x.cross(force);
return true;
}
}
==========================================================================
--
You received this question notification because you are a member of
yade-users, which 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
