|
HG-MD
1
|
|
HG-MD
1
|
00001 #ifndef CSPECIES_H 00002 #define CSPECIES_H 00003 00004 //This is a header with some extra standard maths function and constants that are required by the code 00005 #include "ExtendedMath.h" 00006 00007 #include "Vector.h" 00008 #include <vector> 00009 #include <iostream> 00010 #include <iomanip> 00011 #include <stdio.h> 00012 #include <stdlib.h> 00013 #include <string.h> 00014 00028 class CSpecies { 00029 public: 00030 CSpecies() { 00031 k = 0; 00032 kt = 0; 00033 krolling = 0; 00034 ktorsion = 0; 00035 disp = 0; 00036 dispt = 0; 00037 disprolling = 0; 00038 disptorsion = 0; 00039 mu = 0; 00040 mus = 0; 00041 mu = 0; 00042 murolling = 0; 00043 musrolling = 0; 00044 mutorsion = 0; 00045 mustorsion = 0; 00046 k2max = 0; 00047 kc = 0; 00048 depth = 0; 00049 rho = 0; 00050 dim_particle = 0; 00051 } 00052 00053 CSpecies(const CSpecies& S) { 00054 k = S.get_k(); 00055 kt = S.get_kt(); 00056 krolling = S.get_krolling(); 00057 ktorsion = S.get_ktorsion(); 00058 disp = S.get_dissipation(); 00059 dispt = S.get_dispt(); 00060 disprolling = S.get_disprolling(); 00061 disptorsion = S.get_disptorsion(); 00062 mu = S.get_mu(); 00063 mus = S.get_mus(); 00064 murolling = S.get_murolling(); 00065 musrolling = S.get_musrolling(); 00066 mutorsion = S.get_mutorsion(); 00067 mustorsion = S.get_mustorsion(); 00068 k2max = S.get_k2max(); 00069 kc = S.get_kc(); 00070 depth = S.get_depth(); 00071 rho = S.get_rho(); 00072 dim_particle = S.get_dim_particle(); 00073 } 00074 00076 void set_k(Mdouble new_k){if(new_k>=0) k=new_k; else { cerr << "Error in set_k" << endl; exit(-1); }} 00078 Mdouble get_k() const {return k;} 00079 00081 void set_kt(Mdouble new_kt){if(new_kt>=0) {kt=new_kt;} else { cerr << "Error in set_kt" << endl; exit(-1); }} 00083 Mdouble get_kt() const {return kt;} 00085 void set_krolling(Mdouble new_k){if(new_k>=0) {krolling=new_k;} else { cerr << "Error in set_krolling" << endl; exit(-1); }} 00087 Mdouble get_krolling() const {return krolling;} 00089 void set_ktorsion(Mdouble new_k){if(new_k>=0) {ktorsion=new_k;} else { cerr << "Error in set_ktorsion" << endl; exit(-1); }} 00091 Mdouble get_ktorsion() const {return ktorsion;} 00092 00095 void set_rho(Mdouble new_rho){if(new_rho>=0) rho=new_rho; else { cerr << "Error in set_rho" << endl; exit(-1); }} 00097 Mdouble get_rho() const {return rho;} 00098 00100 void set_dispt(Mdouble new_dispt){if (new_dispt>=0) dispt = new_dispt; else { cerr << "Error in set_dispt" << endl; exit(-1); }} 00102 Mdouble get_dispt() const {return dispt;} 00104 void set_disprolling(Mdouble new_disprolling){if (new_disprolling>=0) disprolling = new_disprolling; else { cerr << "Error in set_disprolling" << endl; exit(-1); }} 00106 Mdouble get_disprolling() const {return disprolling;} 00108 void set_disptorsion(Mdouble new_disptorsion){if (new_disptorsion>=0) disptorsion = new_disptorsion; else { cerr << "Error in set_disptorsion" << endl; exit(-1); }} 00110 Mdouble get_disptorsion() const {return disptorsion;} 00111 00113 void set_dissipation(Mdouble new_disp){if(new_disp>=0) disp=new_disp; else { cerr << "Error in set_dissipation(" << new_disp << ")" << endl; exit(-1); }} 00115 Mdouble get_dissipation() const {return disp;} 00116 00118 //mu has to be set to allow tangential forces (sets dispt=disp as default) 00119 void set_mu(Mdouble new_mu){if (new_mu>=0) {mu = new_mu; mus = mu;} else { cerr << "Error in set_mu" << endl; exit(-1); }} 00121 Mdouble get_mu() const {return mu;} 00123 void set_mus(Mdouble new_mu){if (new_mu>=0) {mus = new_mu;} else { cerr << "Error in set_mus" << endl; exit(-1); }} 00125 Mdouble get_mus() const {return mus;} 00126 00128 //murolling has to be set to allow tangential forces (sets dispt=disp as default) 00129 void set_murolling(Mdouble new_murolling){if (new_murolling>=0) {murolling = new_murolling; musrolling = murolling;} else { cerr << "Error in set_murolling" << endl; exit(-1); }} 00131 Mdouble get_murolling() const {return murolling;} 00133 void set_musrolling(Mdouble new_mu){if (new_mu>=0) {musrolling = new_mu;} else { cerr << "Error in set_musrolling" << endl; exit(-1); }} 00135 Mdouble get_musrolling() const {return musrolling;} 00136 00138 //mu has to be set to allow tangential forces (sets dispt=disp as default) 00139 void set_mutorsion(Mdouble new_mu){if (new_mu>=0) {mutorsion = new_mu; mustorsion = mutorsion;} else { cerr << "Error in set_mutorsion" << endl; exit(-1); }} 00141 Mdouble get_mutorsion() const {return mutorsion;} 00143 void set_mustorsion(Mdouble new_mu){if (new_mu>=0) {mustorsion = new_mu;} else { cerr << "Error in set_mustorsion" << endl; exit(-1); }} 00145 Mdouble get_mustorsion() const {return mustorsion;} 00146 00148 void set_dim_particle(int new_dim){if (new_dim>=1 && new_dim<=3) dim_particle = new_dim; else { cerr << "Error in set_dim_particle" << endl; exit(-1); }} 00150 int get_dim_particle() const {return dim_particle;} 00151 00153 friend inline std::ostream& operator<<(std::ostream& os, const CSpecies &s) { 00154 os << "k " << s.k 00155 << " disp " << s.disp 00156 << " kt " << s.kt; 00157 if (s.krolling) os << " krolling " << s.krolling; 00158 if (s.ktorsion) os << " ktorsion " << s.ktorsion; 00159 os << " dispt " << s.dispt; 00160 if (s.disprolling) os << " disprolling " << s.disprolling; 00161 if (s.disptorsion) os << " disptorsion " << s.disptorsion; 00162 os << " mu " << s.mu; 00163 //optional output 00164 if (s.mu!=s.mus) os << " mus " << s.mus; 00165 if (s.murolling) os << " murolling " << s.murolling; 00166 if (s.murolling!=s.musrolling) os << " musrolling " << s.musrolling; 00167 if (s.mutorsion) os << " mutorsion " << s.mutorsion; 00168 if (s.mutorsion!=s.mustorsion) os << " mustorsion " << s.mustorsion; 00169 if (s.depth) { 00170 os << " k2max " << s.k2max 00171 << " kc " << s.kc 00172 << " depth " << s.depth; 00173 } 00174 os << " rho " << s.rho 00175 << " dim_particle " << s.dim_particle; 00176 return os; 00177 } 00178 00180 friend inline std::istream& operator>>(std::istream& is, CSpecies &s) { 00181 is >> s.k >> s.disp >> s.kt >> s.krolling >> s.ktorsion >> s.dispt >> s.disprolling >> s.disptorsion >> s.mu >> s.mus >> s.murolling >> s.musrolling >> s.mutorsion >> s.mustorsion >> s.k2max >> s.kc >> s.depth >> s.rho >> s.dim_particle; 00182 //s.mus=s.mu; 00183 return is; 00184 } 00185 00186 void read(std::istream& is) { 00187 string dummy; 00188 is >> dummy >> k 00189 >> dummy >> disp 00190 >> dummy >> kt >> dummy; 00191 if (!strcmp(dummy.c_str(),"krolling")) is >> krolling >> dummy; else krolling = 0; 00192 if (!strcmp(dummy.c_str(),"ktorsion")) is >> ktorsion >> dummy; else ktorsion = 0; 00193 is >> dispt >> dummy; 00194 if (!strcmp(dummy.c_str(),"disprolling")) is >> disprolling >> dummy; else disprolling = 0; 00195 if (!strcmp(dummy.c_str(),"disptorsion")) is >> disptorsion >> dummy; else disptorsion = 0; 00196 is >> mu >> dummy; 00197 if (!strcmp(dummy.c_str(),"mus")) is >> mus >> dummy; else mus=mu; 00198 if (!strcmp(dummy.c_str(),"murolling")) is >> murolling >> dummy; else murolling = 0; 00199 if (!strcmp(dummy.c_str(),"musrolling")) is >> musrolling >> dummy; else musrolling = murolling; 00200 if (!strcmp(dummy.c_str(),"mutorsion")) is >> mutorsion >> dummy; else mutorsion = 0; 00201 if (!strcmp(dummy.c_str(),"mustorsion")) is >> mustorsion >> dummy; else mustorsion = mutorsion; 00202 00203 // checks if plastic model parameters are included (for backward compability) 00204 if (!strcmp(dummy.c_str(),"k2max")) { 00205 is >> k2max 00206 >> dummy >> kc 00207 >> dummy >> depth 00208 >> dummy; 00209 } 00210 is >> rho 00211 >> dummy >> dim_particle; 00212 //read rest of line 00213 getline(is,dummy); 00214 } 00215 00217 void print(std::ostream& os) { 00218 os << "k:" << k 00219 << ", disp:" << disp 00220 << ", kt:" << kt; 00221 if (krolling) os << ", krolling: " << krolling; 00222 if (ktorsion) os << ", ktorsion: " << ktorsion; 00223 os << ", dispt: " << dispt; 00224 if (disprolling) os << ", disprolling: " << disprolling; 00225 if (disptorsion) os << ", disptorsion: " << disptorsion; 00226 os << ", mu: " << mu; 00227 if (mu!=mus) os << ", mus: " << mus; 00228 if (murolling) os << ", murolling: " << murolling; 00229 if (murolling!=musrolling) os << ", musrolling: " << musrolling; 00230 if (mutorsion) os << ", mutorsion: " << mutorsion; 00231 if (mutorsion!=mustorsion) os << ", mustorsion: " << mustorsion; 00232 os << ", k2max:" << k2max 00233 << ", kc:" << kc 00234 << ", depth:" << depth; 00235 if (dim_particle) { 00236 os << ", rho:" << rho 00237 << ", dim_particle:" << dim_particle; 00238 } else { 00239 os << " (mixed)"; 00240 } 00241 } 00242 00244 Mdouble get_collision_time(Mdouble mass){ 00245 if(mass<=0) {cerr << "Error in get_collision_time(Mdouble mass) mass is not set or has an unexpected value, (get_collision_time("<<mass<<"))"<< endl; exit(-1);} 00246 if(k<=0) {cerr << "Error in get_collision_time(Mdouble mass) stiffness is not set or has an unexpected value, (get_collision_time("<<mass<<"), with stiffness="<<k<<")"<< endl; exit(-1);} 00247 if(disp<0) {cerr << "Error in get_collision_time(Mdouble mass) dissipation is not set or has an unexpected value, (get_collision_time("<<mass<<"), with dissipation="<<disp<<")"<< endl; exit(-1);} 00248 Mdouble tosqrt=k/(.5*mass) - sqr(disp/mass); 00249 //~ cout<<"tosqrt "<<tosqrt<<" 1e "<<k/(.5*mass)<<" 2e "<<sqr(disp/mass)<<endl; 00250 if (tosqrt<=0) {cerr << "Error in get_collision_time(Mdouble mass) values for mass, stiffness and dissipation would leads to an over damped system, (get_collision_time("<<mass<<"), with stiffness="<<k<<" and dissipation="<<disp<<")"<< endl; exit(-1);} 00251 return constants::pi / sqrt( tosqrt ); 00252 } 00253 00255 Mdouble get_restitution_coefficient(Mdouble mass){return exp(-disp/mass*get_collision_time(mass));} 00256 00258 Mdouble get_maximum_velocity(Mdouble radius, Mdouble mass){return radius * sqrt(k/(.5*mass));} 00259 00261 void set_k_and_restitution_coefficient(Mdouble k_, Mdouble eps, Mdouble mass){ 00262 k = k_; 00263 disp = - sqrt(2.0*mass*k/(constants::sqrt_pi+sqr(log(eps)))) * log(eps); 00264 } 00265 00267 void set_collision_time_and_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble mass){ 00268 disp = - mass / tc * log(eps); 00269 k = .5 * mass * (sqr(constants::pi/tc) + sqr(disp/mass)); 00270 } 00271 00275 void set_collision_time_and_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble mass1, Mdouble mass2) 00276 { 00277 Mdouble reduced_mass = mass1*mass2/(mass1+mass2); 00278 set_collision_time_and_restitution_coefficient(tc,eps,2.0*reduced_mass); 00279 } 00280 00282 void set_collision_time_and_normal_and_tangential_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass){ 00283 set_collision_time_and_restitution_coefficient(tc,eps,mass); 00284 //from Deen...Kuipers2006, eq. 43 and 30 00285 set_kt(2.0/7.0*get_k()*(sqr(constants::pi)+sqr(log(beta)))/(sqr(constants::pi)+sqr(log(eps)))); 00286 if (beta) set_dispt(-2*log(beta)*sqrt(1.0/7.0*mass*get_kt()/(sqr(constants::pi)+sqr(log(beta))))); 00287 else set_dispt(2.*sqrt(1.0/7.0*mass*get_kt())); 00288 } 00289 00291 void set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt(Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass){ 00292 set_collision_time_and_restitution_coefficient(tc,eps,mass); 00293 //from BeckerSchwagerPoeschel2008, eq. 56 00294 set_kt(2.0/7.0*get_k()*sqr(acos(-beta)/constants::pi)); 00295 set_dispt(0); 00296 } 00297 00299 void set_collision_time_and_normal_and_tangential_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2) 00300 { 00301 Mdouble reduced_mass = mass1*mass2/(mass1+mass2); 00302 set_collision_time_and_normal_and_tangential_restitution_coefficient(tc,eps,beta,2.0*reduced_mass); 00303 } 00304 00306 void set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt(Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2) 00307 { 00308 Mdouble reduced_mass = mass1*mass2/(mass1+mass2); 00309 set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt(tc,eps,beta,2.0*reduced_mass); 00310 } 00311 00312 Mdouble get_average(Mdouble a, Mdouble b) { 00313 return (a+b) ? (2.*(a*b)/(a+b)) : 0; 00314 } 00315 00317 void mix(CSpecies& S0, CSpecies& S1){ 00318 k = get_average(S0.get_k(),S1.get_k()); 00319 kt = get_average(S0.get_kt(),S1.get_kt()); 00320 krolling = get_average(S0.get_krolling(),S1.get_krolling()); 00321 ktorsion = get_average(S0.get_ktorsion(),S1.get_ktorsion()); 00322 disp = get_average(S0.get_dissipation(),S1.get_dissipation()); 00323 dispt = get_average(S0.get_dispt(),S1.get_dispt()); 00324 disprolling = get_average(S0.get_disprolling(),S1.get_disprolling()); 00325 disptorsion = get_average(S0.get_disptorsion(),S1.get_disptorsion()); 00326 mu = get_average(S0.get_mu(),S1.get_mu()); 00327 mus = get_average(S0.get_mus(),S1.get_mus()); 00328 murolling = get_average(S0.get_murolling(),S1.get_murolling()); 00329 musrolling = get_average(S0.get_musrolling(),S1.get_musrolling()); 00330 mutorsion = get_average(S0.get_mutorsion(),S1.get_mutorsion()); 00331 mustorsion = get_average(S0.get_mustorsion(),S1.get_mustorsion()); 00332 k2max = get_average(S0.get_k2max(),S1.get_k2max()); 00333 kc = get_average(S0.get_kc(),S1.get_kc()); 00334 depth = get_average(S0.get_depth(),S1.get_depth()); 00335 rho = 0; 00336 dim_particle = 0; //this will be used to distinguish a mixed species 00337 } 00338 00339 00340 Mdouble k; //(normal) spring constant 00341 Mdouble kt; 00342 Mdouble krolling; 00343 Mdouble ktorsion; //tangential spring constant 00344 Mdouble disp; //(normal) viscosity 00345 Mdouble dispt; //tangential viscosity: should satisfy 4*dispt*dt<mass, dispt \approx disp 00346 Mdouble disprolling; //tangential viscosity: should satisfy 4*dispt*dt<mass, dispt \approx disp 00347 Mdouble disptorsion; //tangential viscosity: should satisfy 4*dispt*dt<mass, dispt \approx disp 00348 Mdouble mu;//< (dynamic) Coulomb friction coefficient 00349 Mdouble mus; //static Coulomb friction coefficient (by default set equal to mu) 00350 Mdouble murolling;//< (dynamic) Coulomb friction coefficient 00351 Mdouble musrolling; //static Coulomb friction coefficient 00352 Mdouble mutorsion;//< (dynamic) Coulomb friction coefficient 00353 Mdouble mustorsion; //static Coulomb friction coefficient 00354 Mdouble rho; //density 00355 int dim_particle; //determines if 2D or 3D volume is used for mass calculations 00356 vector<CSpecies> MixedSpecies; 00357 00358 public: 00359 00361 void set_plastic_k1_k2max_kc_depth(Mdouble k1_, Mdouble k2max_, Mdouble kc_, Mdouble depth_){ 00362 set_k1(k1_); 00363 set_k2max(k2max_); 00364 set_kc(kc_); 00365 set_depth(depth_); 00366 } 00367 Mdouble get_k1() const {return k;} 00368 Mdouble get_k2max() const {return k2max;} 00369 Mdouble get_kc() const {return kc;} 00370 Mdouble get_depth() const {return depth;} 00371 void set_k1(Mdouble new_) {k = new_;} 00372 void set_k2max(Mdouble new_) {k2max = new_;} 00373 void set_kc(Mdouble new_) {kc = new_;} 00374 void set_depth(Mdouble new_) {depth = new_;} 00376 Mdouble get_plastic_dt(Mdouble mass){ 00377 return 0.02 * constants::pi / sqrt( k2max/(.5*mass) - sqr(disp/mass) ); 00378 } 00379 00380 Mdouble k2max; //<for plastic deformations; the maximum elastic constant 00381 Mdouble kc; //<for plastic deformations; the adhesive spring constant 00382 Mdouble depth; //<for plastic deformations; the depth (relative to the normalized radius) at which kpmax is used 00383 00384 }; 00385 00386 #endif
1.7.6.1
