1
#include "Ip2_ElastMat.hpp"
3
#include <yade/pkg/common/NormShearPhys.hpp>
4
#include <yade/pkg/common/NormShearPhys.hpp>
5
#include <yade/pkg/dem/DemXDofGeom.hpp>
7
YADE_PLUGIN((Ip2_ElastMat_ElastMat_NormPhys)(Ip2_ElastMat_ElastMat_NormShearPhys));
10
void Ip2_ElastMat_ElastMat_NormPhys::go( const shared_ptr<Material>& b1
11
, const shared_ptr<Material>& b2
12
, const shared_ptr<Interaction>& interaction)
14
if(interaction->phys) return;
15
const shared_ptr<ElastMat>& mat1 = YADE_PTR_CAST<ElastMat>(b1);
16
const shared_ptr<ElastMat>& mat2 = YADE_PTR_CAST<ElastMat>(b2);
17
Real Ea = mat1->young;
18
Real Eb = mat2->young;
19
interaction->phys = shared_ptr<NormPhys>(new NormPhys());
20
const shared_ptr<NormPhys>& phys = YADE_PTR_CAST<NormPhys>(interaction->phys);
22
const GenericSpheresContact* geom=dynamic_cast<GenericSpheresContact*>(interaction->geom.get());
24
Real Ra,Rb;//Vector3r normal;
25
Ra=geom->refR1>0?geom->refR1:geom->refR2;
26
Rb=geom->refR2>0?geom->refR2:geom->refR1;
27
//harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
28
Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
36
void Ip2_ElastMat_ElastMat_NormShearPhys::go( const shared_ptr<Material>& b1
37
, const shared_ptr<Material>& b2
38
, const shared_ptr<Interaction>& interaction)
40
if(interaction->phys) return;
41
const shared_ptr<ElastMat>& mat1 = YADE_PTR_CAST<ElastMat>(b1);
42
const shared_ptr<ElastMat>& mat2 = YADE_PTR_CAST<ElastMat>(b2);
43
Real Ea = mat1->young;
44
Real Eb = mat2->young;
45
Real Va = mat1->poisson;
46
Real Vb = mat2->poisson;
47
interaction->phys = shared_ptr<NormShearPhys>(new NormShearPhys());
48
const shared_ptr<NormShearPhys>& phys = YADE_PTR_CAST<NormShearPhys>(interaction->phys);
50
GenericSpheresContact* geom=dynamic_cast<GenericSpheresContact*>(interaction->geom.get());
52
Real Ra,Rb;//Vector3r normal;
53
Ra=geom->refR1>0?geom->refR1:geom->refR2;
54
Rb=geom->refR2>0?geom->refR2:geom->refR1;
55
//harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
56
Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
57
Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);
60
Kn = 2*Ea*Va*Eb*Vb/(Ea*Va+Eb*Vb);