_cx = _cy = _cz = 1;
_cz0 = 0;
_peaks[0] = _peaks[1] = 1;
- for(int i=0; i<4; i++)
- _stalls[i] = _widths[i] = 0;
+ int i;
+ for(i=0; i<4; i++) {
+ _stalls[i] = 0;
+ _widths[i] = 0.01; // half a degree
+ }
_orient[0] = 1; _orient[1] = 0; _orient[2] = 0;
_orient[3] = 0; _orient[4] = 1; _orient[5] = 0;
_orient[6] = 0; _orient[7] = 0; _orient[8] = 1;
+ _chord = 0;
_incidence = 0;
+ _twist = 0;
_slatPos = _spoilerPos = _flapPos = 0;
_slatDrag = _spoilerDrag = _flapDrag = 1;
+
_flapLift = 0;
_slatAlpha = 0;
_spoilerLift = 1;
+ _inducedDrag = 1;
}
void Surface::setPosition(float* p)
{
- for(int i=0; i<3; i++) _pos[i] = p[i];
+ int i;
+ for(i=0; i<3; i++) _pos[i] = p[i];
}
void Surface::getPosition(float* out)
{
- for(int i=0; i<3; i++) out[i] = _pos[i];
+ int i;
+ for(i=0; i<3; i++) out[i] = _pos[i];
}
void Surface::setChord(float chord)
void Surface::setOrientation(float* o)
{
- for(int i=0; i<9; i++)
+ int i;
+ for(i=0; i<9; i++)
_orient[i] = o[i];
}
_incidence = angle;
}
+void Surface::setTwist(float angle)
+{
+ _twist = angle;
+}
+
void Surface::setSlatParams(float stallDelta, float dragPenalty)
{
_slatAlpha = stallDelta;
// Handle the blowup condition. Zero velocity means zero force by
// definition.
if(vel == 0) {
- for(int i=0; i<3; i++) out[i] = torque[i] = 0;
+ int i;
+ for(i=0; i<3; i++) out[i] = torque[i] = 0;
return;
}
// "Rotate" by the incidence angle. Assume small angles, so we
// need to diddle only the Z component, X is relatively unchanged
// by small rotations.
- out[2] += _incidence * out[0]; // z' = z + incidence * x
+ float incidence = _incidence + _twist;
+ out[2] += incidence * out[0]; // z' = z + incidence * x
+ // Hold onto the local wind vector so we can multiply the induced
+ // drag at the end.
+ float lwind[3];
+ Math::set3(out, lwind);
+
// Diddle the Z force according to our configuration
float stallMul = stallFunc(out);
stallMul *= 1 + _spoilerPos * (_spoilerLift - 1);
float stallLift = (stallMul - 1) * _cz * out[2];
- float flapLift = _cz * _flapPos * (_flapLift-1);
+ float flaplift = flapLift(out[2]);
out[2] *= _cz; // scaling factor
out[2] += _cz*_cz0; // zero-alpha lift
out[2] += stallLift;
- out[2] += flapLift;
+ out[2] += flaplift;
// Airfoil lift (pre-stall and zero-alpha) torques "up" (negative
// torque) around the Y axis, while flap lift pushes down. Both
// forces are considered to act at one third chord from the
- // center. Convert to local (i.e. airplane) coordiantes and store
+ // edge. Convert to local (i.e. airplane) coordiantes and store
// into "torque".
torque[0] = 0;
- torque[1] = 0.33 * _chord * (flapLift - (_cz*_cz0 + stallLift));
+ torque[1] = 0.1667f * _chord * (flaplift - (_cz*_cz0 + stallLift));
torque[2] = 0;
Math::tmul33(_orient, torque, torque);
- // Diddle X (drag) and Y (side force) in the same manner
- out[0] *= _cx * controlDrag();
+ // The X (drag) force gets diddled for control deflection
+ out[0] = controlDrag(out[2], _cx * out[0]);
+
+ // Add in any specific Y (side force) coefficient.
out[1] *= _cy;
+ // Diddle the induced drag
+ Math::mul3(-1*_inducedDrag*out[2]*lwind[2], lwind, lwind);
+ Math::add3(lwind, out, out);
+
// Reverse the incidence rotation to get back to surface
// coordinates.
- out[2] -= _incidence * out[0];
+ out[2] -= incidence * out[0];
// Convert back to external coordinates
Math::tmul33(_orient, out, out);
// Add in the units to make a real force:
- float scale = 0.5*rho*vel*vel*_c0;
+ float scale = 0.5f*rho*vel*vel*_c0;
Math::mul3(scale, out, out);
Math::mul3(scale, torque, torque);
}
return 1;
// (note mask: we want to use the "positive" stall angle here)
- float scale = 0.5*_peaks[fwdBak]/_stalls[i&2];
+ float scale = 0.5f*_peaks[fwdBak]/_stalls[i&2];
// Before the stall
if(alpha <= stallAlpha)
return scale*(1-frac) + frac;
}
-float Surface::controlDrag()
+// Similar to the above -- interpolates out the flap lift past the
+// stall alpha
+float Surface::flapLift(float alpha)
{
- float d = 1;
- d *= 1 + _spoilerPos * (_spoilerDrag - 1);
- d *= 1 + _slatPos * (_slatDrag - 1);
+ float flapLift = _cz * _flapPos * (_flapLift-1);
+ if(alpha < 0) alpha = -alpha;
+ if(alpha < _stalls[0])
+ return flapLift;
+ else if(alpha > _stalls[0] + _widths[0])
+ return 1;
+
+ float frac = (alpha - _stalls[0]) / _widths[0];
+ frac = frac*frac*(3-2*frac);
+ return flapLift * (1-frac) + frac;
+}
+
+float Surface::controlDrag(float lift, float drag)
+{
// Negative flap deflections don't affect drag until their lift
- // multiplier exceeds the "camber" (cz0) of the surface.
+ // multiplier exceeds the "camber" (cz0) of the surface. Use a
+ // synthesized "fp" number instead of the actual flap position.
float fp = _flapPos;
if(fp < 0) {
fp = -fp;
fp -= _cz0/(_flapLift-1);
if(fp < 0) fp = 0;
}
-
- d *= 1 + fp * (_flapDrag - 1);
-
- return d;
+
+ // Calculate an "effective" drag -- this is the drag that would
+ // have been produced by an unflapped surface at the same lift.
+ float flapDragAoA = (_flapLift - 1 - _cz0) * _stalls[0];
+ float fd = Math::abs(lift * flapDragAoA * fp);
+ if(drag < 0) fd = -fd;
+ drag += fd;
+
+ // Now multiply by the various control factors
+ drag *= 1 + fp * (_flapDrag - 1);
+ drag *= 1 + _spoilerPos * (_spoilerDrag - 1);
+ drag *= 1 + _slatPos * (_slatDrag - 1);
+
+ return drag;
}
}; // namespace yasim
projects / flightgear.git / blobdiff