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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/////////////////////////////////////////////////////////////
// //
// Copyright (c) 2003-2014 by The University of Queensland //
// Centre for Geoscience Computing //
// http://earth.uq.edu.au/centre-geoscience-computing //
// //
// Primary Business: Brisbane, Queensland, Australia //
// Licensed under the Open Software License version 3.0 //
// http://www.opensource.org/licenses/osl-3.0.php //
// //
/////////////////////////////////////////////////////////////
//--- project includes ---
#include "Model/AdhesiveFriction.h"
#include "tml/message/packed_message_interface.h"
/*!
Default constructor for CAdhesiveFriction interaction
Zero all coefficients
*/
CAdhesiveFriction::CAdhesiveFriction()
{
m_k=0.0;
m_ks=0.0;
m_r0=0.0;
m_dt=0.0;
m_r_cut=0.0;
m_r_cut_h=0.0;
}
/*!
Construct a CAdhesiveFriction interaction from 2 particle pointers and the parameters
\param p1 pointer to the first particle
\param p2 pointer to the second particle
\param param the interaction parameters
*/
CAdhesiveFriction::CAdhesiveFriction(CParticle* p1,CParticle* p2,const CAdhesiveFrictionIGP& param):CFrictionInteraction(p1,p2)
{
m_k=param.k;
m_ks=param.k_s;
m_r0=p1->getRad()+p2->getRad();
m_dt=param.dt;
m_r_cut=param.r_cut;
m_r_cut_h=1.0+(param.r_cut-1.0)*0.5;
}
/*!
destruct a CAdehsiveFriction interaction, i.e.do nothing
*/
CAdhesiveFriction::~CAdhesiveFriction()
{}
void CAdhesiveFriction::calcForces()
{
Vec3 pos;
Vec3 force;
// calculate distance
Vec3 D=m_p1->getPos()-m_p2->getPos();
double dist=D*D;
// check if there is contact
if(dist<(m_r0*m_r0)){ // contact -> calculate forces as for normal frictional interaction
CFrictionInteraction::calcForces();
} else if (dist<(m_r0*m_r0*m_r_cut_h*m_r_cut_h)){ // between eq and half cut -> increase
//--- elastic force ---
dist=sqrt(dist);
force=D*(m_k*(dist-m_r0)/dist);
pos=m_p2->getPos()+(m_p2->getRad()/m_r0)*D;
m_Ffric=Vec3(0.0,0.0,0.0);
m_normal_force=Vec3(0.0,0.0,0.0);
// apply elastic force
m_p2->applyForce(force,pos);
m_p1->applyForce(-1.0*force,pos);
} else if (dist<(m_r0*m_r0*m_r_cut*m_r_cut)){ // between half cut and cut -> decrease
//--- elastic force ---
dist=sqrt(dist);
force=D*(m_k*((m_r0*m_r_cut)-dist)/dist);
pos=m_p2->getPos()+(m_p2->getRad()/m_r0)*D;
m_Ffric=Vec3(0.0,0.0,0.0);
m_normal_force=Vec3(0.0,0.0,0.0);
// apply elastic force
m_p2->applyForce(force,pos);
m_p1->applyForce(-1.0*force,pos);
}
}
/*!
Pack a CAdhesiveFriction into a TML packed message
\param I the interaction
*/
template<>
void TML_PackedMessageInterface::pack<CAdhesiveFriction>(const CAdhesiveFriction& I)
{
append(I.m_k);
append(I.m_r0);
append(I.m_mu);
append(I.m_ks);
append(I.m_dt);
append(I.m_r_cut);
append(I.m_id[0]);
append(I.m_id[1]);
}
/*!
Unpack a CAdhesiveFriction from a TML packed message
\param I the interaction
*/
template<>
void TML_PackedMessageInterface::unpack<CAdhesiveFriction>(CAdhesiveFriction& I)
{
I.m_k=pop_double();
I.m_r0=pop_double();
I.m_mu=pop_double();
I.m_ks=pop_double();
I.m_dt=pop_double();
I.m_r_cut=pop_double();
I.m_r_cut_h=1.0+(I.m_r_cut-1.0)*0.5;
I.m_id.erase(I.m_id.begin(),I.m_id.end());
I.m_id.push_back(pop_int());
I.m_id.push_back(pop_int());
}
|