← Back to team overview

yade-users team mailing list archive

Re: [Question #251218]: Modifying the Cohesive Frictional Contact Law

 

Question #251218 on Yade changed:
https://answers.launchpad.net/yade/+question/251218

    Status: Answered => Open

behzad is still having a problem:
Good!

and now my Ip2_CohBurgersMat_CohBurgersMat_CohbergersPhys.hpp:

================


#pragma once

#include<yade/pkg/common/Dispatching.hpp>
#include<yade/pkg/dem/CohBurgersMat.hpp>

class Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys : public IPhysFunctor
{
	public :
		virtual void go(	const shared_ptr<Material>& b1,
					const shared_ptr<Material>& b2,
					const shared_ptr<Interaction>& interaction);
		int cohesionDefinitionIteration;

		YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor,
		((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts in current time-step. The flag is turned false automatically, so that assignment is done in the current timestep only."))
		((bool,setCohesionOnNewContacts,true))	
		,
		cohesionDefinitionIteration = -1;
		);
	FUNCTOR2D(CohBurgersMat,CohBurgersMat);
};

REGISTER_SERIALIZABLE(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys);

=======================


and my Ip2_CohBurgersMat_CohBurgersMat_CohbergersPhys.cpp:

================================

#include"Ip2_CohBurgersMat_CohBurgersMat_CohbergersPhys.hpp"
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/pkg/dem/CohBurgersPhys.hpp>
#include<yade/pkg/dem/CohBurgersMat.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>

void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1    // CohBurgersMat
                                        , const shared_ptr<Material>& b2 // CohBurgersMat
                                        , const shared_ptr<Interaction>& interaction)
{
	CohBurgersMat* sdec1 = static_cast<CohBurgersMat*>(b1.get());
	CohBurgersMat* sdec2 = static_cast<CohBurgersMat*>(b2.get());
	ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());

	//Create cohesive interractions only once
	if (setCohesionNow && cohesionDefinitionIteration==-1) cohesionDefinitionIteration=scene->iter;
	if (setCohesionNow && cohesionDefinitionIteration!=-1 && cohesionDefinitionIteration!=scene->iter) {
		cohesionDefinitionIteration = -1;
		setCohesionNow = 0;}

	if (geom) {
		if (!interaction->phys) {
			interaction->phys = shared_ptr<CohBurgersPhys>(new CohBurgersPhys());
			CohBurgersPhys* contactPhysics = YADE_CAST<CohBurgersPhys*>(interaction->phys.get());
			Real Ea 	= sdec1->young;
			Real Eb 	= sdec2->young;
			Real Va 	= sdec1->poisson;
			Real Vb 	= sdec2->poisson;
			Real Da 	= geom->radius1;
			Real Db 	= geom->radius2;
			Real fa 	= sdec1->frictionAngle;
			Real fb 	= sdec2->frictionAngle;
			Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses

			// harmonic average of alphas parameters
			Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
			Real AlphaKtw;
			if (sdec1->alphaKtw && sdec2->alphaKtw) AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
			else AlphaKtw=0;

			Real Ks;
			if (Va && Vb) Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
			else Ks=0;

			// Jean-Patrick Plassiard, Noura Belhaine, Frederic
			// Victor Donze, "Calibration procedure for spherical
			// discrete elements using a local moemnt law".
			contactPhysics->kr = Da*Db*Ks*AlphaKr;
			contactPhysics->ktw = Da*Db*Ks*AlphaKtw;
			contactPhysics->tangensOfFrictionAngle		= std::tan(std::min(fa,fb));

			if ((setCohesionOnNewContacts || setCohesionNow) && sdec1->isCohesive && sdec2->isCohesive)
			{
				contactPhysics->cohesionBroken = false;
				contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(Db, Da),2);
				contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(Db, Da),2);
				geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
			}
			contactPhysics->kn = Kn;
			contactPhysics->ks = Ks;

			contactPhysics->maxRollPl = min(sdec1->etaRoll*Da,sdec2->etaRoll*Db);
			contactPhysics->momentRotationLaw=(sdec1->momentRotationLaw && sdec2->momentRotationLaw);
			//contactPhysics->elasticRollingLimit = elasticRollingLimit;
		}
		else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
			CohBurgersPhys* contactPhysics = YADE_CAST<CohBurgersPhys*>(interaction->phys.get());
			if ((setCohesionNow && sdec1->isCohesive && sdec2->isCohesive) || contactPhysics->initCohesion)
			{
				contactPhysics->cohesionBroken = false;
				contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
				contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(geom->radius2, geom->radius1),2);

				geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
				contactPhysics->initCohesion=false;
			}
		}
	}
};
YADE_PLUGIN((Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys));


==============================

Thanks for having a look on it! I appreciate it

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.