yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #02939
Re: UniaxialTest
2010/4/30 Václav Šmilauer <eudoxos@xxxxxxxx>
> there is some problem with your contact law and I cannot guess without
> seeing the
> code.
I attached my files, please have a look at them.
Read this 10 times aloud:
>
> PARTICLE MATERIAL'S MODULUS IS NOT EQUAL TO MACROSCOPIC MODULUS AND
> THAT IS OK.
>
> (You _will_ have to do calibration, but don't worry about that now)
>
> v
>
> I use the Kn and Ks of this paper (FCC,relation (26));
Macroscopic elastic properties of regular lattices (Yucang Wang , Peter
Mora)
http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6TXB-4TF2J92-1&_user=1399990&_coverDate=12%2F31%2F2008&_rdoc=1&_fmt=high&_orig=search&_sort=d&_docanchor=&view=c&_acct=C000052576&_version=1&_urlVersion=0&_userid=1399990&md5=6bdea54ae28d243cac06403f6a3d4748
In these Kn and Ks, E and v are macroscopic properties.
/*************************************************************************
* Nasibeh Moradi 2010 *
*************************************************************************/
#include"NasibehMoradi.hpp"
#include<yade/pkg-dem/FrictPhys.hpp>
#include<yade/pkg-dem/CohFrictPhys.hpp>
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
#include<yade/pkg-common/ElastMat.hpp>
YADE_PLUGIN((Law2_Dem3DofGeom_CohFrictPhys_Continuum) (Ip2_FrictMat_FrictMat_CohFrictPhys_Continuum));
// same as Law2_Dem3DofGeom_FrictPhys_Basic, but using tension and cohision
void Law2_Dem3DofGeom_CohFrictPhys_Continuum::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene*){
Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
CohFrictPhys* phys=static_cast<CohFrictPhys*>(ip.get());
/* Calsculation of stiffneses when break occure*/
//Real Da=geom->refR1>0?geom->refR1:geom->refR2;
//Real Ea=74.1e9;
//Real Va=0.2;
//Real Kn = Ea*Da;
//Real Ks = Ea*Da*Va;
/* crossSection */
Real crossSection=sqrt(2)*pow(geom->refR1,2);
/* Normal Force */
Real displN=geom->displacementN();
Real Fn=phys->kn*displN;
Real SigmaC=350e6;
Real maxFnCompresion=SigmaC*crossSection;
//Real SigmaT=34.6e6;
//Real maxdisplN=crossSection*SigmaT/phys->kn;
if(displN>0 && (Fn>phys->normalAdhesion || phys->normalAdhesion==0))
{
/*Break due to tension*/
phys->cohesionBroken = true;
phys->normalAdhesion = 0;
phys->shearAdhesion = 0;
scene->interactions->requestErase(contact->getId1(),contact->getId2());
//if (Fn>phys->normalAdhesion)
//{
//phys->kn=Kn;
//phys->ks=Ks;
//}
return;
}
if (displN<0 && abs(Fn)>maxFnCompresion)
{
Fn=-maxFnCompresion;
}
phys->normalForce=Fn*geom->normal;
/* Shear Force */
/* Mohr-Coulomb criterion */
Real maxFs=std::max((Real) 0,phys->shearAdhesion-Fn*phys->tangensOfFrictionAngle);
Vector3r trialFs=phys->ks*geom->displacementT();
//std::cout<<"displacementT()"<<geom->displacementT()<<" "<<"trialFs="<<trialFs<<" "<<"maxFs="<<maxFs<<endl;
if(trialFs.squaredNorm()>(maxFs*maxFs))
{
/*Break due to Shear*/
//geom->slipToDisplacementTMax((sqrt(phys->normalForce.squaredNorm()*pow(phys->tangensOfFrictionAngle,2)))/phys->ks);
trialFs*=sqrt((phys->normalForce.squaredNorm()*pow(phys->tangensOfFrictionAngle,2))/(trialFs.squaredNorm()));
//std::cout<<"normalForce="<<phys->normalForce<<" "<<"trialFs="<<trialFs<<endl;
if (displN>0)
{
trialFs=Vector3r::Zero();
phys->normalForce=Vector3r::Zero();
}
phys->cohesionBroken = true;
phys->normalAdhesion = 0;
phys->shearAdhesion = 0;
//phys->kn=Kn;
//phys->ks=Ks;
}
phys->shearForce=trialFs;
/* Apply forces*/
applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,scene);
}
/*------------------ New Class -------------------*/
void Ip2_FrictMat_FrictMat_CohFrictPhys_Continuum::go( const shared_ptr<Material>& b1
, const shared_ptr<Material>& b2
, const shared_ptr<Interaction>& interaction)
{
if(interaction->interactionPhysics) return;
const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
if (!interaction->interactionPhysics)
interaction->interactionPhysics = shared_ptr<CohFrictPhys>(new CohFrictPhys());
const shared_ptr<CohFrictPhys>& contactPhysics = YADE_PTR_CAST<CohFrictPhys>(interaction->interactionPhysics);
Real Ea = mat1->young;
Real Eb = mat2->young;
Real Va = mat1->poisson;
Real Vb = mat2->poisson;
Real Da,Db; Vector3r normal;
//FIXME : dynamic casts here???!!!
Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
if(scg){Da=scg->radius1; Db=scg->radius2; normal=scg->normal;}
else if(d3dg){Da=d3dg->refR1>0?d3dg->refR1:d3dg->refR2; Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; normal=d3dg->normal;}
else throw runtime_error("geometry is neither ScGeom nor Dem3DofGeom");
//if ((scene->currentIteration<cohesiveTresholdIteration) && (normalAdhesion>0))
//if (scene->currentIteration<200)
//{
contactPhysics->cohesionBroken=false;
Real crossSection=sqrt(2)*pow(std::min(Da,Db),2);
contactPhysics->normalAdhesion=crossSection*tensileStrength;
contactPhysics->shearAdhesion=crossSection*cohesion;
Real Kna = sqrt(2)*Ea*Da/(2*(1-2*Va));
Real Knb = sqrt(2)*Eb*Db/(2*(1-2*Vb));
//harmonic average of the two stiffnesses when (Di.Ei) is the stiffness of sphere "i"
Real Kn = 2*Kna*Knb/(Kna+Knb);
//same for shear stiffness
Real Ksa = (1-3*Va)*Kna/(1+Va);
Real Ksb = (1-3*Vb)*Knb/(1+Vb);
Real Ks = 2*Ksa*Ksb/(Ksa+Ksb);
//}
//else
if (scene->currentIteration>cohesiveTresholdIteration)
{
contactPhysics->cohesionBroken=true;
contactPhysics->normalAdhesion = 0;
contactPhysics->shearAdhesion = 0;
//harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of sphere "i"
//Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);
//same for shear stiffness
//Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);
}
contactPhysics->frictionAngle = std::min(mat1->frictionAngle,mat2->frictionAngle);
contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
contactPhysics->prevNormal = normal;
contactPhysics->kn = Kn;
contactPhysics->ks = Ks;
return;
};
/*************************************************************************
* Nasibeh Moradi 2010 *
*************************************************************************/
#pragma once
#include<yade/core/GlobalEngine.hpp>
// only to see whether SCG_SHEAR is defined, may be removed in the future
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/ElastMat.hpp>
#include<yade/pkg-common/LawFunctor.hpp>
#include<yade/pkg-common/InteractionPhysicsFunctor.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>
/* Constitutive Law for linear compression and tension for glass with cohision*/
class Law2_Dem3DofGeom_CohFrictPhys_Continuum: public LawFunctor{
public:
virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
FUNCTOR2D(Dem3DofGeom,CohFrictPhys);
YADE_CLASS_BASE_DOC(Law2_Dem3DofGeom_CohFrictPhys_Continuum,LawFunctor,"Constitutive law for linear compression,tension,shear");
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CohFrictPhys_Continuum);
class Ip2_FrictMat_FrictMat_CohFrictPhys_Continuum: public InteractionPhysicsFunctor{
public:
virtual void go(const shared_ptr<Material>& b1,
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
FUNCTOR2D(FrictMat, FrictMat);
DECLARE_LOGGER;
YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_CohFrictPhys_Continuum,InteractionPhysicsFunctor,"Converts 2 Frictmat instances to CohesiveFrictionalContactInteraction with corresponding parameters.",
((int,cohesiveTresholdIteration,1,"Should new contacts be cohesive? They will before this iter, they won't afterward."))
((Real,tensileStrength,0,"Defines the maximum admissible normal force in traction FnMax=tensileStrength*crossSection. [Pa]"))
((Real,cohesion,0,"Defines the maximum admissible tangential force in shear FsMax=cohesion*crossSection. [Pa]")));
};
REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_CohFrictPhys_Continuum);
Follow ups
References
