1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/modules/module-hunt-crossley/module-hunt-crossley.cc,v 1.8 2017/01/12 14:52:34 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31 /*
32 * Author: Pierangelo Masarati <masarati@aero.polimi.it>
33 */
34
35 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
36
37 #include <cmath>
38 #include <cfloat>
39
40 #include "dataman.h"
41 #include "constltp_impl.h"
42
43 class HuntCrossleyCL
44 : public ElasticConstitutiveLaw<doublereal, doublereal> {
45 protected:
46 doublereal m_dSign;
47 doublereal m_dAlpha;
48 doublereal m_dK;
49 doublereal m_dExp;
50 bool m_bActive; // is contact ongoing?
51 bool m_bToggling; // toggle m_bActive
52 doublereal m_InitialEpsPrime; // initial contact velocity
53
54 public:
HuntCrossleyCL(const TplDriveCaller<doublereal> * pTplDC,doublereal dSign,doublereal dAlpha,doublereal dK,doublereal dExp)55 HuntCrossleyCL(const TplDriveCaller<doublereal> *pTplDC,
56 doublereal dSign, doublereal dAlpha, doublereal dK, doublereal dExp)
57 : ElasticConstitutiveLaw<doublereal, doublereal>(pTplDC, 0.),
58 m_dSign(dSign), m_dAlpha(dAlpha), m_dK(dK), m_dExp(dExp),
59 m_bActive(false), m_bToggling(false), m_InitialEpsPrime(0.)
60 {
61 NO_OP;
62 };
63
~HuntCrossleyCL(void)64 virtual ~HuntCrossleyCL(void) {
65 NO_OP;
66 };
67
GetConstLawType(void) const68 virtual ConstLawType::Type GetConstLawType(void) const {
69 return ConstLawType::VISCOELASTIC;
70 };
71
pCopy(void) const72 virtual ConstitutiveLaw<doublereal, doublereal>* pCopy(void) const {
73 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
74
75 // pass parameters to copy constructor
76 SAFENEWWITHCONSTRUCTOR(pCL, HuntCrossleyCL,
77 HuntCrossleyCL(pGetDriveCaller()->pCopy(),
78 m_dSign, m_dAlpha, m_dK, m_dExp));
79 return pCL;
80 };
81
Restart(std::ostream & out) const82 virtual std::ostream& Restart(std::ostream& out) const {
83 out << "hunt crossley"
84 << ", sign, " << m_dSign
85 << ", alpha, " << m_dAlpha
86 << ", kappa, " << m_dK
87 << ", exp, " << m_dExp
88 << ", ", ElasticConstitutiveLaw<doublereal, doublereal>::Restart_int(out);
89 return out;
90 };
91
Update(const doublereal & Eps,const doublereal & EpsPrime)92 virtual void Update(const doublereal& Eps, const doublereal& EpsPrime) {
93 ConstitutiveLaw<doublereal, doublereal>::Epsilon = Eps - ElasticConstitutiveLaw<doublereal, doublereal>::Get();
94 ConstitutiveLaw<doublereal, doublereal>::EpsilonPrime = EpsPrime;
95
96 doublereal x = m_dSign*ConstitutiveLaw<doublereal, doublereal>::Epsilon;
97 if (x < 0.) {
98 if (m_bActive) {
99 if (!m_bToggling) {
100 m_bToggling = true;
101 }
102 }
103
104 ConstitutiveLaw<doublereal, doublereal>::F = 0.;
105 ConstitutiveLaw<doublereal, doublereal>::FDE = 0.;
106 ConstitutiveLaw<doublereal, doublereal>::FDEPrime = 0.;
107
108 } else {
109 doublereal xp = m_dSign*ConstitutiveLaw<doublereal, doublereal>::EpsilonPrime;
110
111 if (!m_bActive) {
112 if (!m_bToggling) {
113 m_bToggling = true;
114 m_InitialEpsPrime = xp;
115 }
116 }
117
118 doublereal xn = std::pow(x, m_dExp);
119 doublereal xnm1 = std::pow(x, m_dExp - 1.);
120
121 ConstitutiveLaw<doublereal, doublereal>::F = m_dSign*(m_dK*xn + 1.5*m_dAlpha*m_dK*xn*xp);
122 ConstitutiveLaw<doublereal, doublereal>::FDE = m_dSign*(m_dExp*m_dK*xnm1 + 1.5*m_dExp*m_dAlpha*m_dK*xnm1*xp);
123 ConstitutiveLaw<doublereal, doublereal>::FDEPrime = m_dSign*(1.5*m_dAlpha*m_dK*xn);
124 }
125 };
126
AfterConvergence(const doublereal & Eps,const doublereal & EpsPrime=0.)127 virtual void AfterConvergence(const doublereal& Eps, const doublereal& EpsPrime = 0.) {
128 if (m_bToggling) {
129 #if 0
130 silent_cout(">> HuntCrossleyCL::AfterConvergence() "
131 "m_bToggling=" << (m_bToggling ? "true" : "false") << " "
132 "m_bActive=" << (m_bActive ? "true" : "false") << " "
133 "m_InitialEpsPrime=" << m_InitialEpsPrime << std::endl);
134 #endif
135 if (m_bActive) {
136 m_bActive = false;
137 m_InitialEpsPrime = 0.;
138
139 } else {
140 m_bActive = true;
141 }
142 m_bToggling = false;
143 #if 0
144 silent_cout("<< HuntCrossleyCL::AfterConvergence() "
145 "m_bToggling=" << (m_bToggling ? "true" : "false") << " "
146 "m_bActive=" << (m_bActive ? "true" : "false") << " "
147 "m_InitialEpsPrime=" << m_InitialEpsPrime << std::endl);
148 #endif
149 }
150 };
151 };
152
153 /* specific functional object(s) */
154 struct HuntCrossleyCLR : public ConstitutiveLawRead<doublereal, doublereal> {
155 virtual ConstitutiveLaw<doublereal, doublereal> *
ReadHuntCrossleyCLR156 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
157 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
158
159 CLType = ConstLawType::VISCOELASTIC;
160
161 if (HP.IsKeyWord("help")) {
162 silent_cerr("HuntCrossleyCL:\n"
163 " hunt crossley,\n"
164 " [ , sign, { negative | positive | <sign> } , ]\n"
165 " alpha, <alpha>,\n"
166 " kappa, <kappa>,\n"
167 " exp, <exp>,\n"
168 " [ , prestress, <prestress> ]\n"
169 " [ , prestrain, (DriveCaller)<prestrain> ]\n"
170 << std::endl);
171
172 if (!HP.IsArg()) {
173 throw NoErr(MBDYN_EXCEPT_ARGS);
174 }
175 }
176
177 doublereal dSign = -1.;
178 if (HP.IsKeyWord("sign")) {
179 if (HP.IsKeyWord("positive")) {
180 dSign = 1.;
181 } else if (HP.IsKeyWord("negative")) {
182 dSign = -1.;
183 } else {
184 doublereal d = HP.GetReal();
185 if (d == 0.) {
186 silent_cerr("HuntCrossleyCLR: invalid sign " << d
187 << " at line " << HP.GetLineData() << std::endl);
188 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
189 }
190 dSign = copysign(1., d);
191 }
192 }
193
194 if (!HP.IsKeyWord("alpha")) {
195 silent_cerr("HuntCrossleyCLR: \"alpha\" expected at line " << HP.GetLineData() << std::endl);
196 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
197 }
198 doublereal dAlpha = HP.GetReal();
199 if (dAlpha < 0.) {
200 silent_cerr("HuntCrossleyCLR: invalid \"alpha\" at line " << HP.GetLineData() << std::endl);
201 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
202 }
203
204 if (!HP.IsKeyWord("kappa")) {
205 silent_cerr("HuntCrossleyCLR: \"kappa\" expected at line " << HP.GetLineData() << std::endl);
206 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
207 }
208 doublereal dK = HP.GetReal();
209 if (dK <= 0.) {
210 silent_cerr("HuntCrossleyCLR: invalid \"kappa\" at line " << HP.GetLineData() << std::endl);
211 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
212 }
213
214 if (!HP.IsKeyWord("exp")) {
215 silent_cerr("HuntCrossleyCLR: \"exp\" expected at line " << HP.GetLineData() << std::endl);
216 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
217 }
218 doublereal dExp = HP.GetReal();
219 if (dExp < 1.) {
220 silent_cerr("HuntCrossleyCLR: invalid \"exp\" at line " << HP.GetLineData() << std::endl);
221 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
222 }
223
224 /* Prestress and prestrain */
225 // doublereal PreStress(0.);
226 // GetPreStress(HP, PreStress);
227 TplDriveCaller<doublereal> *pTplDC = GetPreStrain<doublereal>(pDM, HP);
228
229 SAFENEWWITHCONSTRUCTOR(pCL, HuntCrossleyCL,
230 HuntCrossleyCL(pTplDC, dSign, dAlpha, dK, dExp));
231
232 return pCL;
233 };
234 };
235
236 extern "C" int
module_init(const char * module_name,void * pdm,void * php)237 module_init(const char *module_name, void *pdm, void *php)
238 {
239 #if 0
240 DataManager *pDM = (DataManager *)pdm;
241 MBDynParser *pHP = (MBDynParser *)php;
242 #endif
243
244 ConstitutiveLawRead<doublereal, doublereal> *rf1D = new HuntCrossleyCLR;
245 if (!SetCL1D("hunt" "crossley", rf1D)) {
246 delete rf1D;
247
248 silent_cerr("HuntCrossleyCL: "
249 "module_init(" << module_name << ") "
250 "failed" << std::endl);
251
252 return -1;
253 }
254
255 return 0;
256 }
257
258