29#ifndef HERTZMINDLINDERESIEWICZ3D_HH
30#define HERTZMINDLINDERESIEWICZ3D_HH
37template<
typename T,
template<
typename U>
class PARTICLETYPE,
template<
38 typename W>
class DESCRIPTOR>
39HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::HertzMindlinDeresiewicz3D(
40 T G1, T G2, T v1, T v2, T scale1, T scale2, bool validationKruggelEmden) :
41 Force3D<T, PARTICLETYPE>(), _G1(G1), _G2(G2), _v1(v1), _v2(v2), _scale1(
42 scale1), _scale2(scale2), _validationKruggelEmden(validationKruggelEmden)
45 E1 = 2 * (1 + _v1) * _G1;
46 E2 = 2 * (1 + _v2) * _G2;
49 eE = (1 - util::pow(_v1, 2)) / E1 + (1 - util::pow(_v2, 2)) / E2;
53 eG = (2.0 - _v1) / _G1 + (2 - _v2) / _G2;
57template<
typename T,
template<
typename U>
class PARTICLETYPE,
template<
58 typename W>
class DESCRIPTOR>
59void HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::applyForce(
60 typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt,
61 ParticleSystem3D<T, PARTICLETYPE>& pSys)
63 T force[3] = {T(), T(), T()};
65 computeForce(p, pInt, pSys, force);
69template<
typename T,
template<
typename U>
class PARTICLETYPE,
template<
70 typename W>
class DESCRIPTOR>
71void HertzMindlinDeresiewicz3D<T, PARTICLETYPE, DESCRIPTOR>::computeForce(
72 typename std::deque<PARTICLETYPE<T> >::iterator p, int pInt,
73 ParticleSystem3D<T, PARTICLETYPE>& pSys, T force[3])
75 if (p->getSActivity() > 1) {
77 std::vector<std::pair<size_t, T>> ret_matches;
79 pSys.getContactDetection()->getMatches(pInt, ret_matches);
81 PARTICLETYPE<T>* p2 = NULL;
84 for (
const auto& it : ret_matches) {
86 if (!util::nearZero(it.second)) {
91 T delta = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
109 T
M = p->getMass() * p2->getMass() / (p->getMass() + p2->getMass());
111 T R = p->getRad() * p2->getRad() / (p->getRad() + p2->getRad());
113 std::vector < T > _velR(3, T());
114 _velR[0] = -(p2->getVel()[0] - p->getVel()[0]);
115 _velR[1] = -(p2->getVel()[1] - p->getVel()[1]);
116 _velR[2] = -(p2->getVel()[2] - p->getVel()[2]);
118 std::vector < T > _d(3, T());
119 std::vector < T > _normal(3, T());
122 _d[0] = p2->getPos()[0] - p->getPos()[0];
123 _d[1] = p2->getPos()[1] - p->getPos()[1];
124 _d[2] = p2->getPos()[2] - p->getPos()[2];
126 if ( !util::nearZero(util::norm(_d)) ) {
127 _normal = util::normalize(_d);
134 Vector<T, 3> velR_(_velR);
135 T dot = velR_[0] * _normal[0] + velR_[1] * _normal[1] + velR_[2] * _normal[2];
139 std::vector < T > _velN(3, T());
140 _velN[0] = dot * _normal[0];
141 _velN[1] = dot * _normal[1];
142 _velN[2] = dot * _normal[2];
146 std::vector < T > _velT(3, T());
147 _velT[0] = _velR[0] - _velN[0];
148 _velT[1] = _velR[1] - _velN[1];
149 _velT[2] = _velR[2] - _velN[2];
158 T kn = 4 / 3. * util::sqrt(R) * eE;
159 if (_validationKruggelEmden) {
165 std::vector < T > Fs_n(3, T());
166 Fs_n[0] = -kn * util::pow(delta, 1.5) * _normal[0];
167 Fs_n[1] = -kn * util::pow(delta, 1.5) * _normal[1];
168 Fs_n[2] = -kn * util::pow(delta, 1.5) * _normal[2];
175 T eta_n = 0.3 * util::sqrt(4.5 * M * util::sqrt(delta) * kn);
176 if (_validationKruggelEmden) {
180 std::vector < T > Fd_n(3, T());
181 Fd_n[0] = -eta_n * _velN[0] * util::sqrt(delta);
182 Fd_n[1] = -eta_n * _velN[1] * util::sqrt(delta);
183 Fd_n[2] = -eta_n * _velN[2] * util::sqrt(delta);
185 std::vector < T > F_n(3, T());
186 F_n[0] = Fs_n[0] + Fd_n[0];
187 F_n[1] = Fs_n[1] + Fd_n[1];
188 F_n[2] = Fs_n[2] + Fd_n[2];
193 T kt = 2 * util::sqrt(2 * R) * _G1 / (2 - _v1) * util::pow(delta, 0.5);
196 T eta_t = 2 * util::sqrt(2. / 7. * M * kt);
199 std::vector < T > F_t(3, T());
200 F_t[0] = -eta_t * _velT[0];
201 F_t[1] = -eta_t * _velT[1];
202 F_t[2] = -eta_t * _velT[2];
206 force[0] = _scale1 * F_n[0] + _scale2 * F_t[0];
207 force[1] = _scale1 * F_n[1] + _scale2 * F_t[1];
208 force[2] = _scale1 * F_n[2] + _scale2 * F_t[2];
210 p->getForce()[0] += force[0] * 0.5 ;
211 p->getForce()[1] += force[1] * 0.5 ;
212 p->getForce()[2] += force[2] * 0.5 ;
213 p2->getForce()[0] -= force[0] * 0.5 ;
214 p2->getForce()[1] -= force[1] * 0.5 ;
215 p2->getForce()[2] -= force[2] * 0.5 ;
217 if ((p->getSActivity() || p2->getSActivity()) == 3) {
221 if ((p->getSActivity() == 4) && (p2->getSActivity() != (4 || 3))) {
224 if ((p2->getSActivity() == 4) && (p->getSActivity() != (4 || 3))) {
platform_constant Fraction M[Q][Q]
Top level namespace for all of OpenLB.