OpenLB 1.7
Loading...
Searching...
No Matches
Public Member Functions | List of all members
olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR > Class Template Reference

#include <linearContactForce3D.h>

+ Inheritance diagram for olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >:
+ Collaboration diagram for olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >:

Public Member Functions

 LinearContactForce3D (T G1, T G2, T v1, T v2, T scale1=T(1.), T scale2=T(1.), bool validationKruggelEmden=false)
 
 ~LinearContactForce3D () override
 
void applyForce (typename std::deque< PARTICLETYPE< T > >::iterator p, int pInt, ParticleSystem3D< T, PARTICLETYPE > &pSys) override
 
void computeForce (typename std::deque< PARTICLETYPE< T > >::iterator p, int pInt, ParticleSystem3D< T, PARTICLETYPE > &pSys, T force[3])
 
- Public Member Functions inherited from olb::Force3D< T, PARTICLETYPE >
 Force3D ()
 
 Force3D (Force3D< T, PARTICLETYPE > &)
 
 Force3D (const Force3D< T, PARTICLETYPE > &)
 
virtual ~Force3D ()
 

Additional Inherited Members

- Protected Attributes inherited from olb::Force3D< T, PARTICLETYPE >
OstreamManager clout
 

Detailed Description

template<typename T, template< typename U > class PARTICLETYPE, template< typename W > class DESCRIPTOR>
class olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >

Definition at line 49 of file linearContactForce3D.h.

Constructor & Destructor Documentation

◆ LinearContactForce3D()

template<typename T , template< typename U > class PARTICLETYPE, template< typename W > class DESCRIPTOR>
olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >::LinearContactForce3D ( T G1,
T G2,
T v1,
T v2,
T scale1 = T(1.),
T scale2 = T(1.),
bool validationKruggelEmden = false )

Definition at line 44 of file linearContactForce3D.hh.

45 :
46 Force3D<T, PARTICLETYPE>(), _G1(G1), _G2(G2), _v1(v1), _v2(v2), _scale1(
47 scale1), _scale2(scale2), _validationKruggelEmden(validationKruggelEmden)
48{
49 // E-Modul Particle
50 E1 = 2 * (1 + _v1) * _G1;
51 E2 = 2 * (1 + _v2) * _G2;
52
53 // equivalent combined E-Modul
54 eE = (1 - util::pow(_v1, 2)) / E1 + (1 - util::pow(_v2, 2)) / E2;
55 eE = 1 / eE;
56
57 // equivalent combined E-Modul
58 eG = (2.0 - _v1) / _G1 + (2 - _v2) / _G2;
59 eG = 1. / eG;
60}
olb::util::pow
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112

◆ ~LinearContactForce3D()

template<typename T , template< typename U > class PARTICLETYPE, template< typename W > class DESCRIPTOR>
olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >::~LinearContactForce3D ( )
inlineoverride

Definition at line 54 of file linearContactForce3D.h.

54{};

Member Function Documentation

◆ applyForce()

template<typename T , template< typename U > class PARTICLETYPE, template< typename W > class DESCRIPTOR>
void olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >::applyForce ( typename std::deque< PARTICLETYPE< T > >::iterator p,
int pInt,
ParticleSystem3D< T, PARTICLETYPE > & pSys )
overridevirtual

Implements olb::Force3D< T, PARTICLETYPE >.

Definition at line 64 of file linearContactForce3D.hh.

67{
68 T force[3] = {T(), T(), T()};
69
70 computeForce(p, pInt, pSys, force);
71}
olb::LinearContactForce3D::computeForce
void computeForce(typename std::deque< PARTICLETYPE< T > >::iterator p, int pInt, ParticleSystem3D< T, PARTICLETYPE > &pSys, T force[3])
Definition linearContactForce3D.hh:75

◆ computeForce()

template<typename T , template< typename U > class PARTICLETYPE, template< typename W > class DESCRIPTOR>
void olb::LinearContactForce3D< T, PARTICLETYPE, DESCRIPTOR >::computeForce ( typename std::deque< PARTICLETYPE< T > >::iterator p,
int pInt,
ParticleSystem3D< T, PARTICLETYPE > & pSys,
T force[3] )

Definition at line 75 of file linearContactForce3D.hh.

78{
79
80 std::vector<std::pair<size_t, T>> ret_matches;
81 // kind of contactDetection has to be chosen in application
82 pSys.getContactDetection()->getMatches(pInt, ret_matches);
83
84 PARTICLETYPE<T>* p2 = NULL;
85
86 // iterator walks through number of neighbored particles = ret_matches
87 for (const auto& it : ret_matches) {
88
89 if (!util::nearZero(it.second)) {
90
91 p2 = &pSys[it.first];
92
93 // overlap
94 T delta = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
95
96 // equivalent mass
97 T M = p->getMass() * p2->getMass() / (p->getMass() + p2->getMass());
98 // equivalent radius
99 T R = p->getRad() * p2->getRad() / (p->getRad() + p2->getRad());
100 // relative velocity
101 std::vector < T > _velR(3, T());
102 _velR[0] = -(p2->getVel()[0] - p->getVel()[0]);
103 _velR[1] = -(p2->getVel()[1] - p->getVel()[1]);
104 _velR[2] = -(p2->getVel()[2] - p->getVel()[2]);
105
106 std::vector < T > _d(3, T());
107 std::vector < T > _normal(3, T());
108
109 //_d: vector from particle1 to particle2
110 _d[0] = p2->getPos()[0] - p->getPos()[0];
111 _d[1] = p2->getPos()[1] - p->getPos()[1];
112 _d[2] = p2->getPos()[2] - p->getPos()[2];
113
114 if ( !util::nearZero(util::norm(_d)) ) {
115 _normal = util::normalize(_d);
116 }
117 else {
118 return;
119 }
120
121
122 Vector<T, 3> d_(_d);
123 Vector<T, 3> velR_(_velR);
124 T dot = velR_[0] * _normal[0] + velR_[1] * _normal[1] + velR_[2] * _normal[2];
125
126 // normal part of relative velocity
127 // normal relative to surface of particles at contact point
128 std::vector < T > _velN(3, T());
129 _velN[0] = dot * _normal[0];
130 _velN[1] = dot * _normal[1];
131 _velN[2] = dot * _normal[2];
132
133 // tangential part of relative velocity
134 // tangential relative to surface of particles at contact point
135 std::vector < T > _velT(3, T());
136 _velT[0] = _velR[0] - _velN[0];
137 _velT[1] = _velR[1] - _velN[1];
138 _velT[2] = _velR[2] - _velN[2];
139
140 if (delta > 0.) {
141
142
143 // Force normal
144 // spring constant in normal direction
145 T A = M_PI * util::pow(util::min(p->getRad(),p2->getRad()), 2.);
146 T kn = util::pow(p->getRad()/(A*E1) + p2->getRad()/(A*E2), -1.);
147
148 // part of mechanical force of spring in normal direction
149 // Hertz Contact (P. A. Langston, Powder Technology 85 (1995))
150 std::vector < T > Fs_n(3, T());
151 Fs_n[0] = -kn * delta * _normal[0];
152 Fs_n[1] = -kn * delta * _normal[1];
153 Fs_n[2] = -kn * delta * _normal[2];
154
155 // part of mechanical force of damper in normal direction
156 // damped linear spring (Cundall, Strack 1979)
157 // (K.W. Chu, A.B. Yu, Powder Technology 179 (2008) 104 – 114)
158 // damper constant in normal direction
159 // constant eta_n from H. Kruggel-Endem
160 T eta_n = 0.3 * util::sqrt(4.5 * M * util::sqrt(delta) * kn);
161 if (_validationKruggelEmden) {
162 eta_n = 1.96e5; // to compare to Kruggel-Emden
163 }
164
165 std::vector < T > Fd_n(3, T());
166 Fd_n[0] = -eta_n * _velN[0] * util::sqrt(delta);
167 Fd_n[1] = -eta_n * _velN[1] * util::sqrt(delta);
168 Fd_n[2] = -eta_n * _velN[2] * util::sqrt(delta);
169
170 std::vector < T > F_n(3, T());
171 F_n[0] = Fs_n[0] + Fd_n[0];
172 F_n[1] = Fs_n[1] + Fd_n[1];
173 F_n[2] = Fs_n[2] + Fd_n[2];
174
175 // Force tangential
176 // spring constant in tangential direction
177 // (N.G. Deen, Chemical Engineering Science 62 (2007) 28 - 44)
178 T kt = 2 * util::sqrt(2 * R) * _G1 / (2 - _v1) * util::pow(delta, 0.5);
179
180 // damper constant in normal direction
181 T eta_t = 2 * util::sqrt(2. / 7. * M * kt);
182
183 // part of mechanical force of damper in tangential direction
184 std::vector < T > F_t(3, T());
185 F_t[0] = -eta_t * _velT[0];
186 F_t[1] = -eta_t * _velT[1];
187 F_t[2] = -eta_t * _velT[2];
188
189 // entire force
190 // factor _scale to prevent instability
191 force[0] = _scale1 * F_n[0] + _scale2 * F_t[0];
192 force[1] = _scale1 * F_n[1] + _scale2 * F_t[1];
193 force[2] = _scale1 * F_n[2] + _scale2 * F_t[2];
194
195 p->getForce()[0] += force[0] * 0.5 ;
196 p->getForce()[1] += force[1] * 0.5 ;
197 p->getForce()[2] += force[2] * 0.5 ;
198 p2->getForce()[0] -= force[0] * 0.5 ;
199 p2->getForce()[1] -= force[1] * 0.5 ;
200 p2->getForce()[2] -= force[2] * 0.5 ;
201 }
202 }
203 }
204}
M_PI
#define M_PI
Definition analyticalF.hh:43
olb::descriptors::mrt_data::M
platform_constant Fraction M[Q][Q]
Definition mrtLatticeDescriptors.h:65
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb::util::min
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
olb::util::norm
T norm(const std::vector< T > &a)
l2 norm of a vector of arbitrary length
Definition vectorHelpers.h:146
olb::util::normalize
Vector< T, D > normalize(const Vector< T, D > &a)
Definition vectorHelpers.h:193
olb::util::nearZero
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087

References M_PI.

The documentation for this class was generated from the following files:
Generated on Tue Mar 5 2024 09:58:02 for OpenLB by doxygen 1.10.0