1 /***************************************************************************
2
3 file : aero.cpp
4 created : Sun Mar 19 00:04:50 CET 2000
5 copyright : (C) 2000 by Eric Espie
6 email : torcs@free.fr
7 version : $Id: aero.cpp,v 1.25.2.1 2008/12/31 03:53:56 berniw Exp $
8
9 ***************************************************************************/
10
11 /***************************************************************************
12 * *
13 * This program is free software; you can redistribute it and/or modify *
14 * it under the terms of the GNU General Public License as published by *
15 * the Free Software Foundation; either version 2 of the License, or *
16 * (at your option) any later version. *
17 * *
18 ***************************************************************************/
19
20
21
22 #include "sim.h"
23
24 void
SimAeroConfig(tCar * car)25 SimAeroConfig(tCar *car)
26 {
27 void *hdle = car->params;
28 tdble Cx, FrntArea;
29
30 Cx = GfParmGetNum(hdle, SECT_AERODYNAMICS, PRM_CX, (char*)NULL, 0.4);
31 FrntArea = GfParmGetNum(hdle, SECT_AERODYNAMICS, PRM_FRNTAREA, (char*)NULL, 2.5);
32 car->aero.Clift[0] = GfParmGetNum(hdle, SECT_AERODYNAMICS, PRM_FCL, (char*)NULL, 0.0);
33 car->aero.Clift[1] = GfParmGetNum(hdle, SECT_AERODYNAMICS, PRM_RCL, (char*)NULL, 0.0);
34 float aero_factor = car->options->aero_factor;
35
36 car->aero.SCx2 = 0.5f * AIR_DENSITY * Cx * FrntArea;
37 car->aero.Clift[0] *= aero_factor / 4.0f;
38 car->aero.Clift[1] *= aero_factor / 4.0f;
39 float max_lift = MaximumLiftGivenDrag (car->aero.SCx2, FrntArea);
40 float current_lift = 2.0f * (car->aero.Clift[0] + car->aero.Clift[1]);
41 if (current_lift > max_lift) {
42 fprintf (stderr, "Warning: car %s, driver %s: lift coefficients (%f, %f), generate a lift of %f, while maximum theoretical value is %f\n",
43 car->carElt->_carName,
44 car->carElt->_name,
45 car->aero.Clift[0],
46 car->aero.Clift[1],
47 current_lift,
48 max_lift);
49 }
50
51 GfParmSetNum(hdle, SECT_AERODYNAMICS, PRM_FCL, (char*)NULL, car->aero.Clift[0]);
52 GfParmSetNum(hdle, SECT_AERODYNAMICS, PRM_RCL, (char*)NULL, car->aero.Clift[1]);
53 //printf ("%f %f\n", GfParmGetNum(hdle, SECT_AERODYNAMICS, PRM_FCL, (char*)NULL, 0.0), GfParmGetNum(hdle, SECT_AERODYNAMICS, PRM_RCL, (char*)NULL, 0.0));
54 //printf ("cl: %f\n", car->aero.Clift[0]+car->aero.Clift[1]);
55 car->aero.Cd += car->aero.SCx2;
56 car->aero.rot_front[0] = 0.0;
57 car->aero.rot_front[1] = 0.0;
58 car->aero.rot_front[2] = 0.0;
59 car->aero.rot_lateral[0] = 0.0;
60 car->aero.rot_lateral[1] = 0.0;
61 car->aero.rot_lateral[2] = 0.0;
62 car->aero.rot_vertical[0] = 0.0;
63 car->aero.rot_vertical[1] = 0.0;
64 car->aero.rot_vertical[2] = 0.0;
65
66 }
67
68
69 void
SimAeroUpdate(tCar * car,tSituation * s)70 SimAeroUpdate(tCar *car, tSituation *s)
71 {
72 //tdble hm;
73 int i;
74 tdble airSpeed;
75 tdble dragK = 1.0;
76
77 airSpeed = car->DynGC.vel.x;
78
79 if (airSpeed > 10.0) {
80 tdble x = car->DynGC.pos.x;
81 tdble y = car->DynGC.pos.y;
82 // tdble x = car->DynGC.pos.x + cos(yaw)*wing->staticPos.x;
83 // tdble y = car->DynGC.pos.y + sin(yaw)*wing->staticPos.x;
84 tdble yaw = car->DynGC.pos.az;
85 tdble spdang = atan2(car->DynGCg.vel.y, car->DynGCg.vel.x);
86 for (i = 0; i < s->_ncars; i++) {
87 if (i == car->carElt->index) {
88 continue;
89 }
90
91
92 tdble tmpas = 1.00;
93
94 tCar* otherCar = &(SimCarTable[i]);
95 tdble otherYaw = otherCar->DynGC.pos.az;
96 tdble tmpsdpang = spdang - atan2(y - otherCar->DynGC.pos.y, x - otherCar->DynGC.pos.x);
97 NORM_PI_PI(tmpsdpang);
98 tdble dyaw = yaw - otherYaw;
99 NORM_PI_PI(dyaw);
100
101 if ((otherCar->DynGC.vel.x > 10.0) &&
102 (fabs(dyaw) < 0.1396)) {
103 if (fabs(tmpsdpang) > 2.9671) { /* 10 degrees */
104 /* behind another car - reduce overall airflow */
105 tdble factor = (fabs(tmpsdpang)-2.9671)/(M_PI-2.9671);
106
107 tmpas = 1.0 - factor * exp(- 2.0 * DIST(x, y, otherCar->DynGC.pos.x, otherCar->DynGC.pos.y)/(otherCar->aero.Cd * otherCar->DynGC.vel.x));
108 airSpeed = airSpeed * tmpas;
109 } else if (fabs(tmpsdpang) < 0.1396f) { /* 8 degrees */
110 tdble factor = 0.5f * (0.1396f-fabs(tmpsdpang))/(0.1396f);
111 /* before another car - breaks down rear eddies, reduces only drag*/
112 tmpas = 1.0f - factor * exp(- 8.0 * DIST(x, y, otherCar->DynGC.pos.x, otherCar->DynGC.pos.y) / (car->aero.Cd * car->DynGC.vel.x));
113 dragK = dragK * tmpas;
114 }
115 }
116 }
117 }
118
119 car->airSpeed2 = airSpeed * airSpeed;
120
121 tdble v2 = car->airSpeed2;
122 tdble dmg_coef = ((tdble)car->dammage / 10000.0);
123
124 car->aero.drag = -SIGN(car->DynGC.vel.x) * car->aero.SCx2 * v2 * (1.0 + dmg_coef) * dragK * dragK;
125
126
127 // Since we have the forces ready, we just multiply.
128 // Should insert constants here.
129 // Also, no torque is produced since the effect can be
130 // quite dramatic. Interesting idea to make all drags produce
131 // torque when the car is damaged.
132 car->aero.Mx = car->aero.drag * dmg_coef * car->aero.rot_front[0];
133 car->aero.My = car->aero.drag * dmg_coef * car->aero.rot_front[1];
134 car->aero.Mz = car->aero.drag * dmg_coef * car->aero.rot_front[2];
135
136
137 v2 = car->DynGC.vel.y;
138 car->aero.lateral_drag = -SIGN(v2)*v2*v2*0.7;
139 car->aero.Mx += car->aero.lateral_drag * dmg_coef * car->aero.rot_lateral[0];
140 car->aero.My += car->aero.lateral_drag * dmg_coef * car->aero.rot_lateral[1];
141 car->aero.Mz += car->aero.lateral_drag * dmg_coef * car->aero.rot_lateral[2];
142
143 v2 = car->DynGC.vel.z;
144 car->aero.vertical_drag = -SIGN(v2)*v2*v2*1.5;
145 car->aero.Mx += car->aero.vertical_drag * dmg_coef * car->aero.rot_vertical[0];
146 car->aero.My += car->aero.vertical_drag * dmg_coef * car->aero.rot_vertical[1];
147 car->aero.Mz += car->aero.vertical_drag * dmg_coef * car->aero.rot_vertical[2];
148
149
150
151
152
153 }
154
SimAeroDamage(tCar * car,sgVec3 poc,tdble F)155 void SimAeroDamage(tCar *car, sgVec3 poc, tdble F)
156 {
157 tAero* aero = &car->aero;
158 tdble dmg = F*0.0001;
159
160 aero->rot_front[0] += dmg*(urandom()-.5);
161 aero->rot_front[1] += dmg*(urandom()-.5);
162 aero->rot_front[2] += dmg*(urandom()-.5);
163 if (sgLengthVec3(car->aero.rot_front) > 1.0) {
164 sgNormaliseVec3 (car->aero.rot_front);
165 }
166 aero->rot_lateral[0] += dmg*(urandom()-.5);
167 aero->rot_lateral[1] += dmg*(urandom()-.5);
168 aero->rot_lateral[2] += dmg*(urandom()-.5);
169 if (sgLengthVec3(car->aero.rot_lateral) > 1.0) {
170 sgNormaliseVec3 (car->aero.rot_lateral);
171 }
172 aero->rot_vertical[0] += dmg*(urandom()-.5);
173 aero->rot_vertical[1] += dmg*(urandom()-.5);
174 aero->rot_vertical[2] += dmg*(urandom()-.5);
175 if (sgLengthVec3(car->aero.rot_vertical) > 1.0) {
176 sgNormaliseVec3 (car->aero.rot_vertical);
177 }
178
179 //printf ("aero damage:%f (->%f %f %f)\n", dmg, sgLengthVec3(car->aero.rot_front),
180 // sgLengthVec3(car->aero.rot_lateral), sgLengthVec3(car->aero.rot_vertical));
181
182
183
184 }
185
186 static const char *WingSect[2] = {SECT_FRNTWING, SECT_REARWING};
187
188 void
SimWingConfig(tCar * car,int index)189 SimWingConfig(tCar *car, int index)
190 {
191 void *hdle = car->params;
192 tWing *wing = &(car->wing[index]);
193 tdble area;
194
195 area = GfParmGetNum(hdle, WingSect[index], PRM_WINGAREA, (char*)NULL, 0);
196 // we need also the angle
197 wing->angle = GfParmGetNum(hdle, WingSect[index], PRM_WINGANGLE, (char*)NULL, 0);
198 wing->staticPos.x = GfParmGetNum(hdle, WingSect[index], PRM_XPOS, (char*)NULL, 0);
199 wing->staticPos.z = GfParmGetNum(hdle, WingSect[index], PRM_ZPOS, (char*)NULL, 0);
200
201 switch (car->options->aeroflow_model) {
202 case SIMPLE:
203 wing->Kx = -AIR_DENSITY * area; ///< \bug: there should be a 1/2 here.
204 //wing->Kz = 4.0 * wing->Kx;
205 wing->Kz = car->options->aero_factor * wing->Kx;
206 break;
207 case PLANAR:
208 wing->Kx = -AIR_DENSITY * area * 16.0f;
209 wing->Kz = wing->Kx;
210 break;
211 case OPTIMAL:
212 fprintf (stderr, "Using optimal wings\n");
213 wing->Kx = -AIR_DENSITY * area; ///< \bug: there should be a 1/2 here.
214 wing->Kz = car->options->aero_factor * wing->Kx;
215 break;
216 default:
217 fprintf (stderr, "Unimplemented option %d for aeroflow model\n", car->options->aeroflow_model);
218 }
219
220
221 if (index == 1) {
222 car->aero.Cd -= wing->Kx*sin(wing->angle);
223 }
224 }
225
226
227 void
SimWingUpdate(tCar * car,int index,tSituation * s)228 SimWingUpdate(tCar *car, int index, tSituation* s)
229 {
230 tWing *wing = &(car->wing[index]);
231 tdble vt2 = car->DynGC.vel.x;
232 tdble i_flow = 1.0;
233 // rear wing should not get any flow.
234
235 // compute angle of attack
236 // we don't add ay anymore since DynGC.vel.x,z are now in the correct frame of reference (see car.cpp)
237 tdble aoa = atan2(car->DynGC.vel.z, car->DynGC.vel.x); //+ car->DynGC.pos.ay;
238 // The flow to the rear wing can get cut off at large negative
239 // angles of attack. (so it won't produce lift because it will be
240 // completely shielded by the car's bottom)
241 // The value -0.4 should depend on the positioning of the wing.
242 // we also make this be like that.
243 if (index==1) {
244 i_flow = PartialFlowSmooth (-0.4, aoa);
245 }
246 // Flow to the wings gets cut off by other cars.
247 tdble airSpeed = car->DynGC.vel.x;
248
249 if (airSpeed > 10.0) {
250 tdble yaw = car->DynGC.pos.az;
251 tdble x = car->DynGC.pos.x + cos(yaw)*wing->staticPos.x;
252 tdble y = car->DynGC.pos.y + sin(yaw)*wing->staticPos.x;
253 tdble spdang = atan2(car->DynGCg.vel.y, car->DynGCg.vel.x);
254
255 int i;
256 for (i = 0; i < s->_ncars; i++) {
257 if (i == car->carElt->index) {
258 continue;
259 }
260 tdble tmpas = 1.00;
261 tCar* otherCar = &(SimCarTable[i]);
262 tdble otherYaw = otherCar->DynGC.pos.az;
263 tdble tmpsdpang = spdang - atan2(y - otherCar->DynGC.pos.y, x - otherCar->DynGC.pos.x);
264 NORM_PI_PI(tmpsdpang);
265 tdble dyaw = yaw - otherYaw;
266 NORM_PI_PI(dyaw);
267 if ((otherCar->DynGC.vel.x > 10.0) &&
268 (fabs(dyaw) < 0.1396)) {
269 if (fabs(tmpsdpang) > 2.9671) { /* 10 degrees */
270 /* behind another car - reduce overall airflow */
271 tdble factor = (fabs(tmpsdpang)-2.9671)/(M_PI-2.9671);
272 tmpas = 1.0 - factor*exp(- 2.0 * DIST(x, y, otherCar->DynGC.pos.x, otherCar->DynGC.pos.y) /
273 (otherCar->aero.Cd * otherCar->DynGC.vel.x));
274 i_flow = i_flow * tmpas;
275 }
276 }
277 }
278 }
279 //if (index==1) { -- thrown away so that we have different downforce
280 // reduction for front and rear parts.
281 if (1) {
282 // downforce due to body and ground effect.
283 tdble alpha = 0.0f;
284 tdble vt2b = vt2 * (alpha+(1-alpha)*i_flow);
285 vt2b = vt2b * vt2b;
286 tdble hm = 1.5 * (car->wheel[0].rideHeight + car->wheel[1].rideHeight + car->wheel[2].rideHeight + car->wheel[3].rideHeight);
287 hm = hm*hm;
288 hm = hm*hm;
289 hm = 1.0 + exp(-3.0*hm);
290 car->aero.lift[index] = - car->aero.Clift[index] * vt2b * hm;
291 //car->aero.lift[1] = - car->aero.Clift[1] * vt2b * hm;
292 //printf ("%f\n", car->aero.lift[0]+car->aero.lift[1]);
293 }
294
295
296 vt2=vt2*i_flow;
297 vt2=vt2*vt2;
298
299 aoa += wing->angle;
300
301 // the sinus of the angle of attack
302 tdble sinaoa = sin(aoa);
303 tdble cosaoa = cos(aoa);
304
305
306 if (car->DynGC.vel.x > 0.0f) {
307 switch (car->options->aeroflow_model) {
308 case SIMPLE:
309 wing->forces.x = wing->Kx * vt2 * (1.0f + (tdble)car->dammage / 10000.0f) * sinaoa;
310 wing->forces.z = wing->Kz * vt2 * sinaoa;
311 break;
312 case PLANAR:
313 wing->forces.x = wing->Kx * vt2 * (1.0f + (tdble)car->dammage / 10000.0f) * sinaoa * sinaoa * sinaoa;
314 wing->forces.z = wing->Kz * vt2 * sinaoa * sinaoa * cosaoa;
315 break;
316 case OPTIMAL:
317 wing->forces.x = wing->Kx * vt2 * (1.0f + (tdble)car->dammage / 10000.0f) * (1.0f - cosaoa);
318 wing->forces.x = wing->Kx * vt2 * (1.0f + (tdble)car->dammage / 10000.0f) * sinaoa;
319 break;
320 default:
321 fprintf (stderr, "Unimplemented option %d for aeroflow model\n", car->options->aeroflow_model);
322 }
323 } else {
324 wing->forces.x = wing->forces.z = 0.0f;
325 }
326 }
327
PartialFlowRectangle(tdble theta,tdble psi)328 tdble PartialFlowRectangle(tdble theta, tdble psi)
329 {
330 if (psi>0)
331 return 1.0;
332 if (fabs(psi)>fabs(2.0*theta))
333 return 0.0;
334 return (1.0-(1.0-sin(psi)/sin(2.0*theta)));
335 }
336
PartialFlowSmooth(tdble theta,tdble psi)337 tdble PartialFlowSmooth(tdble theta, tdble psi)
338 {
339 if (psi>0)
340 return 1.0;
341 if (fabs(psi)>fabs(2.0*theta))
342 return 0.0;
343 return (0.5*(1.0+tanh((theta-psi)/(fabs(1.0-(psi/theta))-1.0))));
344 }
345
PartialFlowSphere(tdble theta,tdble psi)346 tdble PartialFlowSphere(tdble theta, tdble psi)
347 {
348 if (psi>0)
349 return 1.0;
350 if (fabs(psi)>fabs(2.0*theta))
351 return 0.0;
352 return (0.5*(1.0+sin((theta-psi)*M_PI/(2.0*theta))));
353 }
354
Max_Cl_given_Cd(tdble Cd)355 tdble Max_Cl_given_Cd (tdble Cd)
356 {
357 // if Cd = 1, then all air hitting the surface is stopped.
358 // In any case, horizontal speed of air particles is given by
359 tdble ux = 1 - Cd;
360
361 // We assume no energy lost and thus can calculate the maximum
362 // possible vertical speed imparted to the paticles
363 tdble uy = sqrt(1 - ux*ux);
364
365 // So now Cl is just the imparted speed
366 return uy;
367 }
368
Max_SCl_given_Cd(tdble Cd,tdble A)369 tdble Max_SCl_given_Cd (tdble Cd, tdble A)
370 {
371 tdble Cl = Max_Cl_given_Cd (Cd);
372 return A * Cl * AIR_DENSITY / 2.0f;
373 }
374
MaximumLiftGivenDrag(tdble drag,tdble A)375 tdble MaximumLiftGivenDrag (tdble drag, tdble A)
376 {
377 // We know the drag, C/2 \rho A.
378 // We must calculate the drag coefficient
379 tdble Cd = (drag / A) * 2.0f / AIR_DENSITY;
380 return Max_SCl_given_Cd (Cd, A);
381 }
382