1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
|
#include "Ip2_ElastMat.hpp"
#include <pkg/common/NormShearPhys.hpp>
#include <pkg/dem/DemXDofGeom.hpp>
namespace yade { // Cannot have #include directive inside.
YADE_PLUGIN((Ip2_ElastMat_ElastMat_NormPhys)(Ip2_ElastMat_ElastMat_NormShearPhys));
void Ip2_ElastMat_ElastMat_NormPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction)
{
if (interaction->phys) return;
const shared_ptr<ElastMat>& mat1 = YADE_PTR_CAST<ElastMat>(b1);
const shared_ptr<ElastMat>& mat2 = YADE_PTR_CAST<ElastMat>(b2);
Real Ea = mat1->young;
Real Eb = mat2->young;
interaction->phys = shared_ptr<NormPhys>(new NormPhys());
const shared_ptr<NormPhys>& phys = YADE_PTR_CAST<NormPhys>(interaction->phys);
Real Kn;
const GenericSpheresContact* geom = dynamic_cast<GenericSpheresContact*>(interaction->geom.get());
if (geom) {
Real Ra, Rb; //Vector3r normal;
Ra = geom->refR1 > 0 ? geom->refR1 : geom->refR2;
Rb = geom->refR2 > 0 ? geom->refR2 : geom->refR1;
//harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
Kn = 2 * Ea * Ra * Eb * Rb / (Ea * Ra + Eb * Rb);
} else {
Kn = 2 * Ea * Eb / (Ea + Eb);
}
phys->kn = Kn;
};
void Ip2_ElastMat_ElastMat_NormShearPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction)
{
if (interaction->phys) return;
const shared_ptr<ElastMat>& mat1 = YADE_PTR_CAST<ElastMat>(b1);
const shared_ptr<ElastMat>& mat2 = YADE_PTR_CAST<ElastMat>(b2);
Real Ea = mat1->young;
Real Eb = mat2->young;
Real Va = mat1->poisson;
Real Vb = mat2->poisson;
interaction->phys = shared_ptr<NormShearPhys>(new NormShearPhys());
const shared_ptr<NormShearPhys>& phys = YADE_PTR_CAST<NormShearPhys>(interaction->phys);
Real Kn = 0.0, Ks = 0.0;
GenericSpheresContact* geom = dynamic_cast<GenericSpheresContact*>(interaction->geom.get());
if (geom) {
Real Ra, Rb; //Vector3r normal;
Ra = geom->refR1 > 0 ? geom->refR1 : geom->refR2;
Rb = geom->refR2 > 0 ? geom->refR2 : geom->refR1;
//harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
Kn = 2 * Ea * Ra * Eb * Rb / (Ea * Ra + Eb * Rb);
Ks = 2 * Ea * Ra * Va * Eb * Rb * Vb / (Ea * Ra * Va + Eb * Rb * Vb);
} else {
Kn = 2 * Ea * Eb / (Ea + Eb);
Kn = 2 * Ea * Va * Eb * Vb / (Ea * Va + Eb * Vb);
}
phys->kn = Kn;
phys->ks = Ks;
};
} // namespace yade
|
