Spring + friction contribution (force density ∝ volj).
98 {
99 const auto yji = p.yj - p.yi;
100 const auto Rji = yji.length();
101 if (!(Rji > 0.) || !
util::isLess(Rji, p.deck.d_contactR))
102 return {};
103
104 auto en = yji / Rji;
105 auto scalar_f = p.deck.d_Kn * (Rji - p.deck.d_contactR) * p.volj;
106 if (scalar_f > 0.)
107 scalar_f = 0.;
108
110 if (!p.deck.d_frictionOn)
111 return f;
112
113 const double fn_mag = -scalar_f;
114 if (!(fn_mag > 0.) || !(p.dt > 0.))
115 return f;
116
117 const auto vji = p.vj - p.vi;
118 auto vt = vji - (vji * en) * en;
119
121 {
122 std::lock_guard<std::mutex> lock(
d_mutex);
125 hist.delta_t += vt * p.dt;
126
127 hist.delta_t -= (hist.delta_t * en) * en;
128
129
130 const double Kt = p.deck.d_Kn;
131 const double kt_vol = Kt * p.volj;
132 ft_trial = hist.delta_t * (-kt_vol);
133 const double ft_mag = ft_trial.length();
134 const double ft_max = p.deck.d_mu * fn_mag;
135
137 ft_trial *= (ft_max / ft_mag);
138 if (kt_vol > 0.)
139 hist.delta_t = ft_trial * (-1.0 / kt_vol);
140 }
141 }
142
143 f += ft_trial;
144 return f;
145}
bool isGreater(const double &a, const double &b)
Returns true if a > b.
bool isLess(const double &a, const double &b)
Returns true if a < b.
A structure to represent 3d vectors.
Public Member Functions inherited from contact::PairForce
1.9.8