1/* Label Indexes 2Axle model example for the wheel4 module of MBDyn 3Louis Gagnon 2012, Louis.Gagnon.10@ulaval.ca 4Axle is moved according to a set function but in a vehicle model it would be attached the vehicle 5 6 7 8*************** CURRENT COORDS SYSTEM BELOW**************** 9 10 ^ x, points to front of car 11 | 12 | 7 (-z), points to center of the Earth 13 | / 14 | / 15 |/ points to right (UK driver side) 16 +-------------------> (-y) 17 origin located at GROUND 18 19*/ 20 21begin: data; 22 integrator: initial value; 23end: data; 24 25# FOR 60kph 26#set: real endtime = 8.00; 27#set: real Vxi = 16.667; #initial x velocity of center of gravity and all connected nodes 28# FOR 40kph 29#set: real endtime = 10.00; 30#set: real Vxi = 11.111; #initial x velocity of center of gravity and all connected nodes 31# FOR 20kph 32#set: real Vxi = 5.5556; #initial x velocity of center of gravity and all connected nodes 33#set: real endtime = 20.00; 34# FOR 2mph 35#set: real Vxi = 3.5;#0.88889; #initial x velocity of center of gravity and all connected nodes 36#set: real endtime = 125.00; 37# FOR iterator 38set: real Vxi = 16.6666666666; #iterZ #initial x velocity of center of gravity and all connected nodes 39set: real endtime = 45./Vxi; #iterZ 40 41set: real delaytime = 0.0000201; 42 43set: real timestep = delaytime; 44set: real maxstep = 0.01; 45set: real minstep = 0.00002; 46 47begin: initial value; 48 initial time: 0.0; 49 final time: endtime; 50 time step: timestep;#0.001/50.; 51 strategy: change, postponed, 99; ## defined later 52 max time step : maxstep;## this allows to cap overall timestep even if the wheel asks for a higher one 53 min time step : minstep; ## this caps the lowest timestep value ##ATTENTION: initial timestep must be between those two bounds 54 max iterations: 100; 55 derivatives max iterations: 1000; # // useful to help discontinuities convergence 56 tolerance: 1e-04; 57 derivatives tolerance : 1.;#2.0;#18 58 derivatives coefficient: 0.00001;#0.001; 59 method: ms, 0.6; 60# linear solver: naive, colamd, mt, 1; 61# linear solver: klu; 62# linear solver: umfpack; 63# output: all; 64# output: iterations; 65# output: residual; 66 threads: disable; # disable or num. of threads; 67 68end: initial value; 69 70set: integer EXTRA = 1; # number of extra driving wheels (how many more than 1) 71 72begin: control data; 73 default orientation: orientation matrix; 74 default output: accelerations; 75 output results : netcdf, notext; 76 output meter: meter, 0., forever, steps, 1; 77 structural nodes: 783 79 ; 80 joints: 815 82 ; 83 rigid bodies: 842 85 ; 86 forces: 87 88 ; 89 loadable elements: 902 91 ; 92# file drivers: 93#1 94#; 95 gravity; 96 print: dof stats; 97 print: dof description; 98end: control data; 99 100module load: "libmodule-wheel4.la"; 101#module load: "libmodule-wheel4.la", wheel4; 102 103set: integer GROUND = 0; #ground 104#set: integer MOVING_G = 4; #moving ground 105set: integer CG = 1; #center of gravity 106set: integer RW = 3; #rear wheel 107set: integer RWNR = 51; #rear wheel non rotating hub 108set: integer RWR = 50; #rear wheel ring 109set: integer RWA = 55; #real wheels axle (first) 110set: integer ROAD_LEFT = 7; #road profile driver (based on position) 111set: integer ROAD_RIGHT = 8; #road profile 112set: integer CE; 113 114 115### tire model data, 116 117 118 119set: real TRH = 1.E-9; # prevents division by zero at null x-velocity at the price of losing validity for velocities above -TRH*1.1 and below TRH*1.1 120set: real TRHA = 1.E-9; # buffer used to prevent division by zero 121set: real TRHT = 1.E-9; # prevents divison by zero for computation of the angle of the car (and wheels) 122set: real TRHTA = 1.E-9; # buffer used on angle zero division prevention 123set: real TRHC = 1.; # cap on kappa 124set: real TdLs = 0.001; # minimum value that the half-contact length will take 125set: real TdReDiv = 0.001; # minimum value that the wheel angular velocity will take in the calculation of the effective rolling radius 126set: real TdRe = 5.; # maximum (cap) ratio of abs((effective rolling radius)/(ring radius)) 127set: bool dtOn = 1; # bool to enable or disable adjustable timestep calculation (use only if you have few bumps on a very smooth road, otherwise it will only slow down the simulation) 128set: real TmaxH = 0.0001; # maximum height change wanted on the road profile for one step, this will govern the time step driver 129set: real dtRes = 0.001; # resolution of bump search, keep sufficiently smaller than ring radius (ie: look at every dt_Res meters for a bump in front of the tire) NOTE: will not look behind because it is assumed that the model is not made to move rearwards and that a small move rearward would not influence the timestep because it should already have been adjusted for the bump when driving towards it) Attention, try to avoid making it to low because the search loop will clog the resolution. 130##set: integer TnumA = 30; # number of steps to look ahead for a bump in the road, for timestep controller (the higher, the smoother and also the least chances of waiting until it's too late to reduce the timestep) 131set: real TmaxF = 1.43; # maximum divison of dt per timestep 132set: real TminF = 0.99; # minimum division of dt per timestep 133set: integer TminS = 12; # minimum number of steps wanted in a force cycle (approx., influences dt, see TdivF/3/4 below) 134set: real TdivF = 0.99; # factor by which to divide the timestep if the force oscillates more than wanted (as determined by TminS) (a number greater than one implies a reduction of the timestep and a number smaller than one implies an increase of the timestep only if the road profile also allows it, thus forcing a slower return to the max timestep if the forces oscillate alot) 135set: real TdivF3 = 1.2; # factor if force changes sign 3 times in last TminS steps 136set: real TdivF4 = 1.43; # factor if force changes sign 4 or more times in last TminS steps 137 138set: real R1 = 0.51; #front wheel radius 139set: real R1i = 0.49; #front wheel inner radius (tread+belt thickness taken out) 140set: real R1r = 0.29; # Contact Radius of the wheel (ie: inner tire radius) 141set: real dR_a1 = 0.4; # r_a1 contact length parameter from Besselink eq. 4.85 142set: real dR_a2 = 2.25;# dR_a2 143set: real rRatio = atan(0.09/R1)/pi; #rRatio, # length of the contact patch as ratio of full circle (to calculate how much mass contributes to the ring's centripetal forces) 144 145set: real pctMR = 0.95; #optimpctMR pct of mass to attribute to ring 146#set: real pct 147set: real MT = 50.; # mass of tire 148set: real MP = 2.; #optimMP mass of patch 149set: real MR = pctMR*MT; # #mass of ring 150set: real Jrx = MR*(R1^3.-R1i^3.)/(R1^2.-R1i^2.)*4./(3.*pi); #moment of inertia of ring about x axis 151set: real Jry = MR*(R1^3.-R1i^3.)/(R1^2.-R1i^2.)*2./3.; #moment of inertia of ring about y axis 152set: real Jrz = Jrx; #moment of inertia of ring about z axis 153set: real MTs = (1-pctMR)*MT; # mass of the tire sidewall portion considered to be attached to the rim 154set: real Jtsx = MTs*(R1i^3.-R1r^3.)/(R1i^2.-R1r^2.)*4./(3.*pi); #moment of inertia of tire sidewall about x axis 155set: real Jtsy = MTs*(R1i^3.-R1r^3.)/(R1i^2.-R1r^2.)*2./3.; #moment of inertia of tire sidewall about y axis 156#set: real Jtsz = Jrx; #moment of inertia of tire sidewall about z axis 157set: MR = MR - MP; # removing mass of patch from ring mass 158 159##temp mod: 160#set: MR=MR/3; 161 162set: real S_hf = -0.003; # residual aligning moment angle modifier 163set: real S_ht = -0.02; # pneumatic trail, for aligning moment, angle modifier 164set: real q_sy1 = 0.008; # tire rolling resistance linear vel coef 165set: real q_sy3 = 0.; # rolling resist other coef 166set: real dvao = 10.; # reference velocity probably as in ISO 28580 (only if q_sy3 is not null) 167 168set: real Pls = 1.1; #optimPls 169 170 171 172### tire data 173# vertical 174set: real KPAdz = 0.035; # vertical deflection used for stiffness calculation (m) 175set: real KPAlz = 30000.; # optimKPAlz 176set: real KPAfwrz = 0.15; # optimKPAfwrz 177set: real KPAfrpz = 1.-KPAfwrz; # fraction of z-deflection due to ring-patch 178# longitudinal 179set: real KPAdx = 0.03; # vertical deflection used for stiffness calculation (m) 180set: real KPAlx = 10000.;# 181set: real KPAfwrx = 0.9; # 182set: real KPAfwrxay = 0.9; # 183set: real KPAfrpx = 0.7; # 184 185 186set: real KPAfwrDen = KPAfwrx + KPAfwrxay + KPAfrpx; 187set: KPAfwrx = KPAfwrx/KPAfwrDen; # 188set: KPAfwrxay = KPAfwrxay/KPAfwrDen; # 189set: KPAfrpx = KPAfrpx/KPAfwrDen; # 190 191 192 193# lateral 194set: real KPAdy = 0.015; # vertical deflection used for stiffness calculation (m) 195set: real KPAly = 4500.; # load at specified vertical deflection (N) 196set: real KPAfwry = 0.15; # 197set: real KPAfwryax = 0.60; # 198set: real KPAfrpy = 0.50; # 199 200set: real PCTD = 0.06; # percent of critical damping for swift model 201 202#SWIFT coefficients for ring-patch connection 203set: real KPArpx = 1*KPAlx/(KPAfrpx*KPAdx); #; 204set: real KPArpy = 1*KPAly/(KPAfrpy*KPAdy); #; 205set: real KPArpz = 1*KPAlz/(KPAfrpz*KPAdz); #; 206set: real CPArpx = 0.2*1*PCTD*2*sqrt(MP*KPArpx);# 207set: real CPArpy = 0.9*1*PCTD*2*sqrt(MP*KPArpy);# 208set: real CPArpz = 1.5*1*PCTD*2*sqrt(MP*KPArpz);# 209#SWIFT coefficients for wheel-ring connection 210set: real KPAwrx = 1*KPAlx/(KPAfwrx*KPAdx); #; 211set: real KPAwry = 1*KPAly/(KPAfwry*KPAdy); #; 212set: real KPAwrz = 1*KPAlz/(KPAfwrz*KPAdz); #; 213set: real CPAwrx = 5.*1*PCTD*2*sqrt(MR*KPAwrx);# 214set: real CPAwry = 1.*1*PCTD*2*sqrt(MR*KPAwry);# 215set: real CPAwrz = 2.*1*PCTD*2*sqrt(MR*KPAwrz);# 216#SWIFT angular coefficients for wheel-ring connection 217set: real KPAwrax = 1*R1^2*KPAly/(KPAfwryax*KPAdy); # # angular stiffness about x-axis 218set: real KPAwraz = 1*KPAwrax; # forced to be equal by definition of the joints 219set: real KPAwray = 1*R1^2*KPAlx/(KPAfwrxay*KPAdx); # # angular elasticity about wheel spin axis 220set: real CPAwrax = 0.2*1*PCTD*2*sqrt(Jrx*KPAwrax);# 221set: real CPAwraz = CPAwrax; # forced to be equal by definition of the joints 222set: real CPAwray = 4.*1*PCTD*2*sqrt(Jry*KPAwray);# 223 224 225 226 227 228set: real MW = 40. + MTs; #mass of wheel (rim + mass of tire assumed on rim [mass tire sidewall]) 229set: real Jwx = 1.7 + Jtsx; # rim + tire sidewall contribution 230set: real Jwy = 2.4 + Jtsy; # rim + tire sidewall contribution 231set: real Jwz = Jwx; 232 233 234### ,tire model data 235 236 237 238 239set: real TF = 27440.; # target force at axle when tire is at rest 240set: real G = 9.81 ; # accel. of gravity 241set: real HG = R1 - ( ( TF + ( MR + MW ) * G ) / KPArpz * ( 1. + KPArpz / KPAwrz ) - MR * G / KPAwrz ) ; #axle height for lab test simulations 242 243 244set: real MI = 5.; #mass of intermediate node (non rotating) 245set: real RWx = -2.5; #x pos of R wheel (rear) 246set: real XPAz = 0.; #initial z pos of wheels... 247set: real CEx; #pos of wheel CE (Current Element) relative to reference 248set: real CEy; #pos of wheel CE (Current Element) relative to reference 249set: real CEz; #pos of wheel CE (Current Element) relative to reference 250set: real Vyi = 0.; #initial x velocity of center of gravity and all connected nodes 251set: real Vzi = 0.; #initial x velocity of center of gravity and all connected nodes 252set: real TORQUE = 100.; #torque value applied to wheels on demand 253 254 255 256 257####Note: could be usefull to put a cap on kappa when angular velocity is low in order to reduce oscillations of the wheel... 258 259### road conditions, 260set: real RDA = 5.; # road offset (null before that value) 261set: real RDBl = 1.; # road length for transition between initial contact patch height and road initial height 262set: real RDB = RDA+RDBl; # road offset (interpolated between RDA and that value) 263set: real RDL = 20.; # road loop condition (will loop after RDB+RDL over RDL) loop transfer zone must be done over a smooth zone (at least the size of half the contact patch length!) (ie: want, ideally, matching slopes on two first and two last points of the profile and matchin heights on first and alst points) 264### ,road conditions 265 266#######NOTE: (smoothens the results even with kappa being capped) this buffer needs to be tuned for best performance and has a significant impact on the stability of the tire, depending on the values chosen the tire will start to oscillate between positive and negative values or not.####### 267 268 269 270reference: GROUND, 271 null, 272 eye, 273 null, 274 null; 275reference: CG, 276 reference, GROUND, 0., 0., HG, 277 reference, GROUND, eye, 278 reference, GROUND, Vxi , Vyi , Vzi, 279 reference, GROUND, null; 280 281reference: RW, 282 reference, CG, RWx, 0., 0., 283 reference, CG, eye, 284 reference, CG, null, 285 reference, CG, 0., 0., 0.; # 50.943 for v =27. 286 287begin: nodes; 288 structural: RWA, static, # rear wheel axle (first) 289 reference, CG, RWx, 0., 0., 290 reference, CG, eye, 291 reference, CG, null, 292 reference, CG, null; 293 294set: CE = 200; 295set: CEx = 0.; #pos of wheel CE (Current Element) relative to reference 296set: CEy = 1.5; #pos of wheel CE (Current Element) relative to reference 297set: CEz = 0.; #pos of wheel CE (Current Element) relative to reference 298 structural: CE, dynamic, # Rear wheel 299 reference, RW, CEx, CEy, CEz, 300 reference, RW, eye, 301 reference, RW, null, 302 reference, RW, null; 303 structural: CE+1, dynamic, # Rear wheel ring 304 reference, RW, CEx, CEy, CEz, 305 reference, RW, eye, 306 reference, RW, null, 307 reference, RW, null; 308 309end: nodes; 310 311#begin: drivers; 312 313#end: drivers; 314 315 316### road definition, 317 318############################################################################### 319## the following command can be useful when feeding a large profile,########### 320#include: "road_cleat10x20.dat"; #iterZ 321############################################################################### 322 323### or for a flat road, 324#drive caller: ROAD_LEFT, scalar function, "road_left", const, 0; 325 326### example road, 327drive caller: ROAD_LEFT, scalar function, "road_left", 328multilinear, 3290.000000, 0, 3300.005000, 1e-03, 3310.010000, 3.5e-03, 3320.015000, 4.5e-3, 3330.020000, 7.5e-03, 3340.025000, 9e-03, 3350.030000, 0.01, 3360.035000, 0.01001; 337 338### ,road definition 339 340 341 342 343 344 345begin: elements; 346 347 joint: 998, total pin joint, # fixes the AXLE 348 RWA, 349 position constraint, 350 0,1,1, 351 0, 0, HG, 352 const, 1, 353 orientation constraint, 354 1,1,1, 355 null; 356 357 358 driven: 999, string, "Time>=delaytime", # this causes a delay in the ouput file for this joint 359 hint, "position-drive3{1., 0., 0., ramp, Vxi, delaytime, forever, model::xposition(RWA) }", 360 joint: 999, total pin joint, # fixes the AXLE velocity 361 RWA, 362 position constraint, 363 active, inactive, inactive, 364 null; # not used, will be replaced by hint 365 366 367 368 369set: integer ROAD_CE = ROAD_LEFT; 370 371 body: CE, #wheel 372 CE, 373 MW, 374 reference, node, 0.0, 0.0, 0.0, 375 diag, Jwx, Jwy, Jwz; 376 body: CE+1,#ring 377 CE+1, 378 MR, 379 reference, node, 0.0, 0.0, 0.0, 380 diag, Jrx, Jry, Jrz; 381 joint: CE+2, deformable hinge, 382 CE, 383 CE+1, 384 linear time variant viscoelastic generic, 385 diag, KPAwrax, KPAwray, KPAwraz, 386 string, "1", 387 diag, CPAwrax, CPAwray, CPAwraz, 388 string, "1"; 389 joint: CE+3, deformable displacement joint, 390 RWA, reference, node, CEx, CEy, XPAz+CEz, 391 CE+1, reference, node, null, 392 linear time variant viscoelastic generic, 393 diag, KPAwrx, KPAwry, KPAwrz, 394 string, "1", 395 diag, CPAwrx, CPAwry, CPAwrz, 396 string, "1"; 397 joint: CE+4, revolute hinge, 398 RWA, 399 position, reference, node, 0, CEy, XPAz+CEz, 400 orientation, 401 1, 0., 0., 1., 402 3, 0., 1., 0., 403 CE, 404 position, reference, node, null, 405 orientation, 406 1, 0., 0., 1., 407 3, 0., 1., 0.; 408 409 410 411set: integer EID = CE; # module element ID 412set: integer WNL = CE; # wheel structural node label 413set: integer WBL = CE; # wheel body label 414set: integer GSNL = RWA; # ground structural node label 415set: integer RNL = CE+1; # ring structural node label 416set: integer RBL = CE+1; # ring body label 417set: integer RPD = ROAD_CE; # driver for the profile of the road 418 419##include: "WUE.elm"; ## Wheel user elem 420 421user defined: EID, rigid ring tire, 422 WNL, # wheel structural node label 423 WBL, # wheel body label 424 0.,1.,0., # wheel axle direction (does sign matter??) 425## GSNL, # ground structural node label 426 R1, # wheel radius 427 swift, 428 RNL, # ring structural node label 429 RBL, # ring body label 430 KPArpx, KPArpy, KPArpz, # Kpa (stiffness) of contact patch (in the ring reference frame) 431 string, "1", # driver for time variant multiplier on the coefficient above 432 CPArpx, CPArpy, CPArpz, # Cpa # should be aboout 6 to 10% of crit damping? also need to distribute this damping over ring and patch and etc... 433 string, "1", # driver for time variant multiplier on the coefficient above 434 0., 0., 0., # // initial velocity of patch in absolute reference frame 435 MP, # // mass of patch 436 reference, RPD, #driver for the profile of the road 437 Pls, # patch contact-length to elleptical cam tandem base parameter (Schmeitz eq. 4.15) dPls = l_s/(2a) 438 dR_a1, # r_a1 contact length parameter from Besselink eq. 4.85 439 dR_a2, # r_a2 contact length parameter from Besselink eq. 4.85 440 rRatio, #length of the contact patch as ratio of full circle (to calculate how much mass contributes to the ring's centripetal forces) 441 KPAwrz, # vertical wheel-ring stiffness (to calculate applied centripetal forces) 442 loadedRadius, # this keyword enables the use of the (home-made) loaded radius as the brake lever arm (instead of the home-made effective rolling radius) comment out to use the effective effective rolling radius 443 slip, 444 ginac, "125*sin(0.01*atan(10.*Var-1.*(10.*Var-atan(10.*Var))))", # Fx/Fz 445 ginac, "0.7*sin(1.5*atan(7.*Var+0.8*(7.*Var-atan(7.*Var))))-0.02", # Fy/Fz without non-linear dependency on Fz (could be added later in module) 446 ginac, "(-0.05*cos(1.3*atan(8.*Var+1.5*(8.*Var-atan(8.*Var)))))/27440.", # Mz 447 ginac, "(-18.*cos(atan(43.*Var)))/27440.", # Mzr 448 S_ht, S_hf, q_sy1, q_sy3, dvao, 449 threshold, TRH, TRHA, TRHT, TRHTA, TRHC, TdLs, TdReDiv, TdRe, dtOn, TmaxH, dtRes, maxstep, minstep, TmaxF, TminF, TminS, TdivF, TdivF3, TdivF4, RDA, RDB, RDL 450 ; 451 452 453 454 455 456user defined: 204, timestep; 457drive caller: 99, element, 204, loadable, string, "1", direct; 458 459 460gravity: 461 0.0, 0.0, -1.0, G; 462end: elements; 463 464 465