fly/flightgear
1
0
Fork 0
Code Activity

Helicopter update from Maik:

Browse source 
More realistic calculation of vortices at the blades and therefore
 real airfoil parameters can be used now (not to be mixed up with the
 vortex ring state which is still not simulated), ground effect is now
 continuous e. g. at buildings, calculating of the rotor in more than 4
 directions, better documentation of the airfoil parameters.
This commit is contained in:
andy 2006-09-14 18:18:33 +00:00
parent 2c193be0ca
commit 77ecc04676
8 changed files with 543 additions and 192 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
src/FDM/YASim/FGFDM.cpp
View file

@ -180,13 +180,17 @@ void FGFDM::startElement(const char* name, const XMLAttributes &atts)
Rotorgear* r = _airplane.getModel()->getRotorgear();
_currObj = r;
#define p(x) if (a->hasAttribute(#x)) r->setParameter((char *)#x,attrf(a,#x) );
#define p2(x,y) if (a->hasAttribute(#x)) {r->setParameter((char *)#y,attrf(a,#x) ); SG_LOG(SG_INPUT, SG_ALERT,"Warning: Aircraft defnition files uses outdated tag '" <<#x<<"', use '"<<#y<<"' instead (see README.YASim)" <<endl);}
p(max_power_engine)
p(engine_prop_factor)
p(yasimdragfactor)
p(yasimliftfactor)
p(max_power_rotor_brake)
p(engine_accell_limit)
p(rotorgear_friction)
p(engine_accel_limit)
p2(engine_accell_limit,engine_accel_limit)
#undef p
#undef p2
r->setInUse();
} else if(eq(name, "wing")) {
_airplane.setWing(parseWing(a, name));
@ -493,8 +497,10 @@ void FGFDM::setOutputProperties(float dt)
char b[256];
while((j = r->getValueforFGSet(j, b, &f)))
if(b[0]) fgSetFloat(b,f);
for(j=0; j < r->numRotorparts(); j++) {
j=0;
while((j = _airplane.getRotorgear()->getValueforFGSet(j, b, &f)))
if(b[0]) fgSetFloat(b,f);
for(j=0; j < r->numRotorparts(); j+=r->numRotorparts()>>2) {
Rotorpart* s = (Rotorpart*)r->getRotorpart(j);
char *b;
int k;
@ -692,6 +698,7 @@ Rotor* FGFDM::parseRotor(XMLAttributes* a, const char* type)
p(vortex_state_e2)
p(twist)
p(number_of_segments)
p(number_of_parts)
p(rel_len_where_incidence_is_measured)
p(chord)
p(taper)
@ -705,7 +712,8 @@ Rotor* FGFDM::parseRotor(XMLAttributes* a, const char* type)
p(airfoil_lift_coefficient)
p(airfoil_drag_coefficient0)
p(airfoil_drag_coefficient1)
p(cyclic_factor)
p(rotor_correction_factor)
#undef p
_currObj = w;
return w;

5
src/FDM/YASim/Math.cpp
View file

@ -80,6 +80,11 @@ float Math::pow(double base, double exp)
return (float)::pow(base, exp);
}
float Math::exp(float f)
{
return (float)::exp(f);
}
double Math::ceil(double f)
{
return ::ceil(f);

2
src/FDM/YASim/Math.hpp
View file

@ -20,6 +20,8 @@ public:
static float atan2(float y, float x);
static float asin(float f);
static float acos(float f);
static float exp(float f);
static float sqr(float f) {return f*f;}
// Takes two args and runs afoul of the Koenig rules.
static float pow(double base, double exp);

16
src/FDM/YASim/Model.cpp
View file

@ -312,21 +312,7 @@ void Model::updateGround(State* s)
}
for(i=0; i<_rotorgear.getRotors()->size(); i++) {
Rotor* r = (Rotor*)_rotorgear.getRotors()->get(i);
// Get the point of the rotor center
float pos[3];
r->getPosition(pos);
// Transform the local coordinates to
// global coordinates.
double pt[3];
s->posLocalToGlobal(pos, pt);
// Ask for the ground plane in the global coordinate system
double global_ground[4];
float global_vel[3];
_ground_cb->getGroundPlane(pt, global_ground, global_vel);
r->setGlobalGround(global_ground, global_vel);
r->findGroundEffectAltitude(_ground_cb,s);
}
// The arrester hook

509
src/FDM/YASim/Rotor.cpp
View file

@ -3,6 +3,7 @@
#include "Math.hpp"
#include "Surface.hpp"
#include "Rotorpart.hpp"
#include "Ground.hpp"
#include "Rotor.hpp"
#include STL_IOSTREAM
@ -10,7 +11,14 @@
SG_USING_STD(setprecision);
#ifdef TEST_DEBUG
#include <stdio.h>
#endif
#include <string.h>
#include <iostream>
#include <sstream>
namespace yasim {
@ -50,6 +58,8 @@ Rotor::Rotor()
_mincyclicail=-10./180*pi;
_mincyclicele=-10./180*pi;
_min_pitch=-.5/180*pi;
_cyclicele=0;
_cyclicail=0;
_name[0] = 0;
_normal[0] = _normal[1] = 0;
_normal[2] = 1;
@ -81,12 +91,14 @@ Rotor::Rotor()
_vortex_state_e1=1;
_vortex_state_e2=1;
_vortex_state_e3=1;
_vortex_state=0;
_lift_factor=1;
_liftcoef=0.1;
_dragcoef0=0.1;
_dragcoef1=0.1;
_twist=0;
_number_of_segments=1;
_number_of_parts=4;
_rel_len_where_incidence_is_measured=0.7;
_torque_of_inertia=1;
_torque=0;
@ -110,6 +122,13 @@ Rotor::Rotor()
_stall_v2sum=1;
_collective=0;
_yaw=_roll=0;
for (int k=0;k<4;k++)
for (i=0;i<3;i++)
_groundeffectpos[k][i]=0;
_ground_effect_altitude=1;
_cyclic_factor=1;
_lift_factor=_f_ge=_f_vs=_f_tl=1;
_rotor_correction_factor=.65;
}
Rotor::~Rotor()
@ -130,10 +149,13 @@ void Rotor::inititeration(float dt,float omegarel,float ddt_omegarel,float *rot)
_ddt_omega=_omegan*ddt_omegarel;
int i;
for(i=0; i<_rotorparts.size(); i++) {
float s = Math::sin(2*pi*i/_number_of_parts);
float c = Math::cos(2*pi*i/_number_of_parts);
Rotorpart* r = (Rotorpart*)_rotorparts.get(i);
r->setOmega(_omega);
r->setDdtOmega(_ddt_omega);
r->inititeration(dt,rot);
r->setCyclic(_cyclicail*c+_cyclicele*s);
}
//calculate the normal of the rotor disc, for calcualtion of the downwash
@ -169,11 +191,29 @@ float Rotor::calcStall(float incidence,float speed)
float Rotor::getLiftCoef(float incidence,float speed)
{
float stall=calcStall(incidence,speed);
/* the next shold look like this, but this is the inner loop of
the rotor simulation. For small angles (and we hav only small
angles) the first order approximation works well
float c1= Math::sin(incidence-_airfoil_incidence_no_lift)*_liftcoef;
for c2 we would need higher order, because in stall the angle can be large
*/
float i2;
if (incidence > (pi/2))
i2 = incidence-pi;
else if (incidence <-(pi/2))
i2 = (incidence+pi);
else
i2 = incidence;
float c1= (i2-_airfoil_incidence_no_lift)*_liftcoef;
if (stall > 0)
{
float c2= Math::sin(2*(incidence-_airfoil_incidence_no_lift))
*_liftcoef*_lift_factor_stall;
return (1-stall)*c1 + stall *c2;
}
else
return c1;
}
float Rotor::getDragCoef(float incidence,float speed)
{
@ -187,14 +227,14 @@ float Rotor::getDragCoef(float incidence,float speed)
int Rotor::getValueforFGSet(int j,char *text,float *f)
{
if (_name[0]==0) return 0;
if (4!=numRotorparts()) return 0; //compile first!
if (4>numRotorparts()) return 0; //compile first!
if (j==0)
{
sprintf(text,"/rotors/%s/cone", _name);
*f=( ((Rotorpart*)getRotorpart(0))->getrealAlpha()
+((Rotorpart*)getRotorpart(1))->getrealAlpha()
+((Rotorpart*)getRotorpart(2))->getrealAlpha()
+((Rotorpart*)getRotorpart(3))->getrealAlpha()
+((Rotorpart*)getRotorpart(1*(_number_of_parts>>2)))->getrealAlpha()
+((Rotorpart*)getRotorpart(2*(_number_of_parts>>2)))->getrealAlpha()
+((Rotorpart*)getRotorpart(3*(_number_of_parts>>2)))->getrealAlpha()
)/4*180/pi;
}
else
@ -202,7 +242,7 @@ int Rotor::getValueforFGSet(int j,char *text,float *f)
{
sprintf(text,"/rotors/%s/roll", _name);
_roll = ( ((Rotorpart*)getRotorpart(0))->getrealAlpha()
-((Rotorpart*)getRotorpart(2))->getrealAlpha()
-((Rotorpart*)getRotorpart(2*(_number_of_parts>>2)))->getrealAlpha()
)/2*(_ccw?-1:1);
*f=_roll *180/pi;
}
@ -210,8 +250,8 @@ int Rotor::getValueforFGSet(int j,char *text,float *f)
if (j==2)
{
sprintf(text,"/rotors/%s/yaw", _name);
_yaw=( ((Rotorpart*)getRotorpart(1))->getrealAlpha()
-((Rotorpart*)getRotorpart(3))->getrealAlpha()
_yaw=( ((Rotorpart*)getRotorpart(1*(_number_of_parts>>2)))->getrealAlpha()
-((Rotorpart*)getRotorpart(3*(_number_of_parts>>2)))->getrealAlpha()
)/2;
*f=_yaw*180/pi;
}
@ -271,18 +311,19 @@ int Rotor::getValueforFGSet(int j,char *text,float *f)
float rk,rl,p;
p=(*f/90);
k=int(p);
l=int(p+1);
l=k+1;
rk=l-p;
rk=Math::clamp(rk,0,1);//Delete this
rl=1-rk;
if(w==2) {k+=2;l+=2;}
else
if(w==1) {k+=1;l+=1;}
k%=4;
l%=4;
if (w==1) *f=rk*((Rotorpart*) getRotorpart(k))->getrealAlpha()*180/pi
+rl*((Rotorpart*) getRotorpart(l))->getrealAlpha()*180/pi;
else if(w==2) *f=rk*((Rotorpart*)getRotorpart(k))->getIncidence()*180/pi
+rl*((Rotorpart*)getRotorpart(l))->getIncidence()*180/pi;
if (w==1) *f=rk*((Rotorpart*) getRotorpart(k*(_number_of_parts>>2)))->getrealAlpha()*180/pi
+rl*((Rotorpart*) getRotorpart(l*(_number_of_parts>>2)))->getrealAlpha()*180/pi;
else if(w==2) *f=rk*((Rotorpart*)getRotorpart(k*(_number_of_parts>>2)))->getIncidence()*180/pi
+rl*((Rotorpart*)getRotorpart(l*(_number_of_parts>>2)))->getIncidence()*180/pi;
}
return j+1;
}
@ -557,6 +598,7 @@ void Rotor::setParameter(char *parametername, float value)
p(vortex_state_e3,1)
p(twist,pi/180.)
p(number_of_segments,1)
p(number_of_parts,1)
p(rel_len_where_incidence_is_measured,1)
p(chord,1)
p(taper,1)
@ -570,8 +612,11 @@ void Rotor::setParameter(char *parametername, float value)
p(airfoil_lift_coefficient,1)
p(airfoil_drag_coefficient0,1)
p(airfoil_drag_coefficient1,1)
cout << "internal error in parameter set up for rotor: '"
<< parametername <<"'" << endl;
p(cyclic_factor,1)
p(rotor_correction_factor,1)
SG_LOG(SG_INPUT, SG_ALERT,
"internal error in parameter set up for rotor: '" <<
parametername <<"'" << endl);
#undef p
}
@ -598,23 +643,15 @@ void Rotor::setCollective(float lval)
void Rotor::setCyclicele(float lval,float rval)
{
rval = Math::clamp(rval, -1, 1);
lval = Math::clamp(lval, -1, 1);
float col=_mincyclicele+(lval+1)/2*(_maxcyclicele-_mincyclicele);
if (_rotorparts.size()!=4) return;
((Rotorpart*)_rotorparts.get(1))->setCyclic(lval);
((Rotorpart*)_rotorparts.get(3))->setCyclic(-lval);
_cyclicele=_mincyclicele+(lval+1)/2*(_maxcyclicele-_mincyclicele);
}
void Rotor::setCyclicail(float lval,float rval)
{
lval = Math::clamp(lval, -1, 1);
rval = Math::clamp(rval, -1, 1);
float col=_mincyclicail+(lval+1)/2*(_maxcyclicail-_mincyclicail);
if (_rotorparts.size()!=4) return;
if (_ccw) lval *=-1;
((Rotorpart*)_rotorparts.get(0))->setCyclic(-lval);
((Rotorpart*)_rotorparts.get(2))->setCyclic( lval);
_cyclicail=-(_mincyclicail+(lval+1)/2*(_maxcyclicail-_mincyclicail));
}
void Rotor::getPosition(float* out)
@ -634,16 +671,8 @@ void Rotor::calcLiftFactor(float* v, float rho, State *s)
_f_tl=1;
_f_vs=1;
// find h, the distance to the ground
// The ground plane transformed to the local frame.
float ground[4];
s->planeGlobalToLocal(_global_ground, ground);
float h = ground[3] - Math::dot3(_base, ground);
// Now h is the distance from the rotor center to ground
// Calculate ground effect
_f_ge=1+_diameter/h*_ground_effect_constant;
_f_ge=1+_diameter/_ground_effect_altitude*_ground_effect_constant;
// Now calculate translational lift
float v_vert = Math::dot3(v,_normal);
@ -656,9 +685,104 @@ void Rotor::calcLiftFactor(float* v, float rho, State *s)
_lift_factor = _f_ge*_f_tl*_f_vs;
}
float Rotor::getGroundEffect(float* posOut)
void Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s)
{
return _diameter*0; //ground effect for rotor is calcualted not here
_ground_effect_altitude=findGroundEffectAltitude(ground_cb,s,
_groundeffectpos[0],_groundeffectpos[1],
_groundeffectpos[2],_groundeffectpos[3]);
}
float Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s,
float *pos0,float *pos1,float *pos2,float *pos3,
int iteration,float a0,float a1,float a2,float a3)
{
float a[5];
float *p[5],pos4[3];
a[0]=a0;
a[1]=a1;
a[2]=a2;
a[3]=a3;
a[4]=-1;
p[0]=pos0;
p[1]=pos1;
p[2]=pos2;
p[3]=pos3;
p[4]=pos4;
Math::add3(p[0],p[2],p[4]);
Math::mul3(0.5,p[4],p[4]);//the center
float mina=100*_diameter;
float suma=0;
for (int i=0;i<5;i++)
{
if (a[i]==-1)//in the first iteration,(iteration==0) no height is
//passed to this function, these missing values are
//marked by ==-1
{
double pt[3];
s->posLocalToGlobal(p[i], pt);
// Ask for the ground plane in the global coordinate system
double global_ground[4];
float global_vel[3];
ground_cb->getGroundPlane(pt, global_ground, global_vel);
// find h, the distance to the ground
// The ground plane transformed to the local frame.
float ground[4];
s->planeGlobalToLocal(global_ground, ground);
a[i] = ground[3] - Math::dot3(p[i], ground);
// Now a[i] is the distance from p[i] to ground
}
suma+=a[i];
if (a[i]<mina)
mina=a[i];
}
if (mina>=10*_diameter)
return mina; //the ground effect will be zero
//check if further recursion is neccessary
//if the height does not differ more than 20%, than
//we can return then mean height, if not split
//zhe square to four parts and calcualte the height
//for each part
//suma * 0.2 is the mean
//0.15 is the maximum allowed difference from the mean
//to the height at the center
if ((iteration>2)
||(Math::abs(suma*0.2-a[4])<(0.15*0.2*suma*(1<<iteration))))
return suma*0.2;
suma=0;
float pc[4][3],ac[4]; //pc[i]=center of pos[i] and pos[(i+1)&3]
for (int i=0;i<4;i++)
{
Math::add3(p[i],p[(i+1)&3],pc[i]);
Math::mul3(0.5,pc[i],pc[i]);
double pt[3];
s->posLocalToGlobal(pc[i], pt);
// Ask for the ground plane in the global coordinate system
double global_ground[4];
float global_vel[3];
ground_cb->getGroundPlane(pt, global_ground, global_vel);
// find h, the distance to the ground
// The ground plane transformed to the local frame.
float ground[4];
s->planeGlobalToLocal(global_ground, ground);
ac[i] = ground[3] - Math::dot3(p[i], ground);
// Now ac[i] is the distance from pc[i] to ground
}
return 0.25*
(findGroundEffectAltitude(ground_cb,s,p[0],pc[1],p[4],pc[3],
iteration+1,a[0],ac[0],a[4],ac[3])
+findGroundEffectAltitude(ground_cb,s,p[1],pc[0],p[4],pc[1],
iteration+1,a[1],ac[0],a[4],ac[1])
+findGroundEffectAltitude(ground_cb,s,p[2],pc[1],p[4],pc[2],
iteration+1,a[2],ac[1],a[4],ac[2])
+findGroundEffectAltitude(ground_cb,s,p[3],pc[2],p[4],pc[3],
iteration+1,a[3],ac[2],a[4],ac[3])
);
}
void Rotor::getDownWash(float *pos, float *v_heli, float *downwash)
@ -743,12 +867,12 @@ void Rotor::compile()
// Have we already been compiled?
if(_rotorparts.size() != 0) return;
//rotor is divided into 4 pointlike parts
//rotor is divided into _number_of_parts parts
//each part is calcualted at _number_of_segments points
SG_LOG(SG_FLIGHT, SG_DEBUG, "debug: e "
<< _mincyclicele << "..." <<_maxcyclicele << ' '
<< _mincyclicail << "..." << _maxcyclicail << ' '
<< _min_pitch << "..." << _max_pitch);
//clamp to 4..256
//and make it a factor of 4
_number_of_parts=(int(Math::clamp(_number_of_parts,4,256))>>2)<<2;
_dynamic=_dynamic*(1/ //inverse of the time
( (60/_rotor_rpm)/4 //for rotating 90 deg
@ -761,7 +885,6 @@ void Rotor::compile()
directions[0][1]=_forward[1];
directions[0][2]=_forward[2];
int i;
SG_LOG(SG_FLIGHT, SG_DEBUG, "Rotor rotating ccw? " << _ccw);
for (i=1;i<5;i++)
{
if (!_ccw)
@ -771,17 +894,23 @@ void Rotor::compile()
Math::unit3(directions[i],directions[i]);
}
Math::set3(directions[4],directions[0]);
float rotorpartmass = _weight_per_blade*_number_of_blades/4*.453;
for (i=0;i<4;i++)
{
Math::mul3(_diameter*0.7,directions[i],_groundeffectpos[i]);
Math::add3(_base,_groundeffectpos[i],_groundeffectpos[i]);
}
float rotorpartmass = _weight_per_blade*_number_of_blades/_number_of_parts*.453;
//was pounds -> now kg
_torque_of_inertia = 1/12. * ( 4 * rotorpartmass) * _diameter
_torque_of_inertia = 1/12. * ( _number_of_parts * rotorpartmass) * _diameter
* _diameter * _rel_blade_center * _rel_blade_center /(0.5*0.5);
float speed=_rotor_rpm/60*_diameter*_rel_blade_center*pi;
float lentocenter=_diameter*_rel_blade_center*0.5;
float lentoforceattac=_diameter*_rel_len_hinge*0.5;
float zentforce=rotorpartmass*speed*speed/lentocenter;
float pitchaforce=_force_at_pitch_a/4*.453*9.81;
//was pounds of force, now N, devided by 4 (so its now per rotorpart)
float pitchaforce=_force_at_pitch_a/_number_of_parts*.453*9.81;
// was pounds of force, now N, devided by _number_of_parts
//(so its now per rotorpart)
float torque0=0,torquemax=0,torqueb=0;
float omega=_rotor_rpm/60*2*pi;
@ -790,14 +919,14 @@ void Rotor::compile()
_delta*=pitchaforce/(_pitch_a*omega*lentocenter*2*rotorpartmass);
float phi=Math::atan2(2*omega*_delta,omega0*omega0-omega*omega);
float relamp=omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
+4*_delta*_delta*omega*omega));
float relamp=(omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
+4*_delta*_delta*omega*omega)))*_cyclic_factor;
if (!_no_torque)
{
torque0=_power_at_pitch_0/4*1000/omega;
torque0=_power_at_pitch_0/_number_of_parts*1000/omega;
// f*r=p/w ; p=f*s/t; r=s/t/w ; r*w*t = s
torqueb=_power_at_pitch_b/4*1000/omega;
torquemax=_power_at_pitch_b/4*1000/omega/_pitch_b*_max_pitch;
torqueb=_power_at_pitch_b/_number_of_parts*1000/omega;
torquemax=_power_at_pitch_b/_number_of_parts*1000/omega/_pitch_b*_max_pitch;
if(_ccw)
{
@ -807,36 +936,30 @@ void Rotor::compile()
}
}
SG_LOG(SG_FLIGHT, SG_DEBUG,
"spd: " << setprecision(8) << speed
<< " lentoc: " << lentocenter
<< " dia: " << _diameter
<< " rbl: " << _rel_blade_center
<< " hing: " << _rel_len_hinge
<< " lfa: " << lentoforceattac);
SG_LOG(SG_FLIGHT, SG_DEBUG,
"tq: " << setprecision(8) << torque0 << ".." << torquemax
<< " d3: " << _delta3);
SG_LOG(SG_FLIGHT, SG_DEBUG,
"o/o0: " << setprecision(8) << omega/omega0
<< " phi: " << phi*180/pi
<< " relamp: " << relamp
<< " delta: " <<_delta);
Rotorpart* rps[4];
for (i=0;i<4;i++)
Rotorpart* rps[256];
for (i=0;i<_number_of_parts;i++)
{
float lpos[3],lforceattac[3],lspeed[3],dirzentforce[3];
float s = Math::sin(2*pi*i/_number_of_parts);
float c = Math::cos(2*pi*i/_number_of_parts);
float direction[3],nextdirection[3],help[3];
Math::mul3(c ,directions[0],help);
Math::mul3(s ,directions[1],direction);
Math::add3(help,direction,direction);
Math::mul3(lentocenter,directions[i],lpos);
Math::mul3(c ,directions[1],help);
Math::mul3(s ,directions[2],nextdirection);
Math::add3(help,nextdirection,nextdirection);
Math::mul3(lentocenter,direction,lpos);
Math::add3(lpos,_base,lpos);
Math::mul3(lentoforceattac,directions[i+1],lforceattac);
//i+1: +90deg (gyro)!!!
Math::mul3(lentoforceattac,nextdirection,lforceattac);
//nextdirection: +90deg (gyro)!!!
Math::add3(lforceattac,_base,lforceattac);
Math::mul3(speed,directions[i+1],lspeed);
Math::mul3(1,directions[i+1],dirzentforce);
Math::mul3(speed,nextdirection,lspeed);
Math::mul3(1,nextdirection,dirzentforce);
float maxcyclic=(i&1)?_maxcyclicele:_maxcyclicail;
float mincyclic=(i&1)?_mincyclicele:_mincyclicail;
@ -845,19 +968,24 @@ void Rotor::compile()
lspeed,dirzentforce,zentforce,pitchaforce, _max_pitch,_min_pitch,
mincyclic,maxcyclic,_delta3,rotorpartmass,_translift,
_rel_len_hinge,lentocenter);
rp->setAlphaoutput(_alphaoutput[i&1?i:(_ccw?i^2:i)],0);
rp->setAlphaoutput(_alphaoutput[4+(i&1?i:(_ccw?i^2:i))],1+(i>1));
int k = i*4/_number_of_parts;
rp->setAlphaoutput(_alphaoutput[k&1?k:(_ccw?k^2:k)],0);
rp->setAlphaoutput(_alphaoutput[4+(k&1?k:(_ccw?k^2:k))],1+(k>1));
_rotorparts.add(rp);
rp->setTorque(torquemax,torque0);
rp->setRelamp(relamp);
rp->setDirectionofRotorPart(directions[i]);
rp->setTorqueOfInertia(_torque_of_inertia/4.);
rp->setDirectionofRotorPart(direction);
rp->setTorqueOfInertia(_torque_of_inertia/_number_of_parts);
}
for (i=0;i<4;i++)
for (i=0;i<_number_of_parts;i++)
{
rps[i]->setlastnextrp(rps[(i+3)%4],rps[(i+1)%4],rps[(i+2)%4]);
rps[i]->setlastnextrp(rps[(i-1+_number_of_parts)%_number_of_parts],
rps[(i+1)%_number_of_parts],
rps[(i+_number_of_parts/2)%_number_of_parts],
rps[(i-_number_of_parts/4+_number_of_parts)%_number_of_parts],
rps[(i+_number_of_parts/4)%_number_of_parts]);
}
for (i=0;i<4;i++)
for (i=0;i<_number_of_parts;i++)
{
rps[i]->setCompiled();
}
@ -868,25 +996,27 @@ void Rotor::compile()
if (_airfoil_lift_coefficient==0)
{
//calculate the lift and drag coefficients now
_liftcoef=0;
_dragcoef0=1;
_dragcoef1=1;
_liftcoef=1;
rps[0]->calculateAlpha(v_wind,rho_null,_pitch_a,0,0,
&(torque[0]),&(lift[0])); //max_pitch a
_liftcoef = pitchaforce/lift[0];
_dragcoef0=1;
_dragcoef1=0;
rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,&(torque[0]),&(lift[0]));
//0 degree, c0
_liftcoef=0;
_dragcoef0=0;
_dragcoef1=1;
rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,&(torque[1]),&(lift[1]));
//0 degree, c1
_liftcoef=0;
_dragcoef0=1;
_dragcoef1=0;
rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,&(torque[2]),&(lift[2]));
//picth b, c0
_liftcoef=0;
_dragcoef0=0;
_dragcoef1=1;
rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,&(torque[3]),&(lift[3]));
@ -903,17 +1033,12 @@ void Rotor::compile()
/(torque[1]/torque[0]-torque[3]/torque[2]);
_dragcoef0=(torqueb-_dragcoef1*torque[3])/torque[2];
}
_liftcoef=1;
rps[0]->calculateAlpha(v_wind,rho_null,_pitch_a,0,0,
&(torque[0]),&(lift[0])); //max_pitch a
_liftcoef = pitchaforce/lift[0];
}
else
{
_liftcoef=_airfoil_lift_coefficient/4*_number_of_blades;
_dragcoef0=_airfoil_drag_coefficient0/4*_number_of_blades;
_dragcoef1=_airfoil_drag_coefficient1/4*_number_of_blades;
_liftcoef=_airfoil_lift_coefficient/_number_of_parts*_number_of_blades;
_dragcoef0=_airfoil_drag_coefficient0/_number_of_parts*_number_of_blades*_c2;
_dragcoef1=_airfoil_drag_coefficient1/_number_of_parts*_number_of_blades*_c2;
}
//Check
@ -923,28 +1048,145 @@ void Rotor::compile()
&(torque[1]),&(lift[1])); //pitch b
rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,
&(torque[3]),&(lift[3])); //pitch 0
SG_LOG(SG_FLIGHT, SG_DEBUG,
SG_LOG(SG_GENERAL, SG_DEBUG,
"Rotor: coefficients for airfoil:" << endl << setprecision(6)
<< " drag0: " << _dragcoef0*4/_number_of_blades
<< " drag1: " << _dragcoef1*4/_number_of_blades
<< " lift: " << _liftcoef*4/_number_of_blades
<< " drag0: " << _dragcoef0*_number_of_parts/_number_of_blades/_c2
<< " drag1: " << _dragcoef1*_number_of_parts/_number_of_blades/_c2
<< " lift: " << _liftcoef*_number_of_parts/_number_of_blades
<< endl
<< "at 10 deg:" << endl
<< "drag: " << (Math::sin(10./180*pi)*_dragcoef1+_dragcoef0)
*4/_number_of_blades
<< "lift: " << Math::sin(10./180*pi)*_liftcoef*4/_number_of_blades
*_number_of_parts/_number_of_blades/_c2
<< " lift: " << Math::sin(10./180*pi)*_liftcoef*_number_of_parts/_number_of_blades
<< endl
<< "Some results (Pitch [degree], Power [kW], Lift [N]" << endl
<< 0.0f << "deg " << Math::abs(torque[3]*4*_omegan/1000) << "kW "
<< lift[3] << endl
<< _pitch_a*180/pi << "deg " << Math::abs(torque[0]*4*_omegan/1000)
<< "kW " << lift[0] << endl
<< _pitch_b*180/pi << "deg " << Math::abs(torque[1]*4*_omegan/1000)
<< "kW " << lift[1] << endl << endl);
<< "Some results (Pitch [degree], Power [kW], Lift [N])" << endl
<< 0.0f << "deg " << Math::abs(torque[3]*_number_of_parts*_omegan/1000) << "kW "
<< lift[3]*_number_of_parts << endl
<< _pitch_a*180/pi << "deg " << Math::abs(torque[0]*_number_of_parts*_omegan/1000)
<< "kW " << lift[0]*_number_of_parts << endl
<< _pitch_b*180/pi << "deg " << Math::abs(torque[1]*_number_of_parts*_omegan/1000)
<< "kW " << lift[1]*_number_of_parts << endl << endl );
rps[0]->setOmega(0);
//first calculation of relamp is wrong
//it used pitchaforce, but this was unknown and
//on the default value
_delta*=lift[0]/pitchaforce;
relamp=(omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
+4*_delta*_delta*omega*omega)))*_cyclic_factor;
for (i=0;i<_number_of_parts;i++)
{
rps[i]->setRelamp(relamp);
}
rps[0]->setOmega(0);
writeInfo();
}
std::ostream & operator<<(std::ostream & out, /*const*/ Rotor& r)
{
#define i(x) << #x << ":" << r.x << endl
#define iv(x) << #x << ":" << r.x[0] << ";" << r.x[1] << ";" <<r.x[2] << ";" << endl
out << "Writing Info on Rotor "
i(_name)
i(_torque)
i(_omega) i(_omegan) i(_omegarel) i(_ddt_omega) i(_omegarelneu)
i (_chord)
i( _taper)
i( _airfoil_incidence_no_lift)
i( _collective)
i( _airfoil_lift_coefficient)
i( _airfoil_drag_coefficient0)
i( _airfoil_drag_coefficient1)
i( _ccw)
i( _number_of_segments)
i( _number_of_parts)
iv( _base)
iv( _groundeffectpos[0])iv( _groundeffectpos[1])iv( _groundeffectpos[2])iv( _groundeffectpos[3])
i( _ground_effect_altitude)
iv( _normal)
iv( _normal_with_yaw_roll)
iv( _forward)
i( _diameter)
i( _number_of_blades)
i( _weight_per_blade)
i( _rel_blade_center)
i( _min_pitch)
i( _max_pitch)
i( _force_at_pitch_a)
i( _pitch_a)
i( _power_at_pitch_0)
i( _power_at_pitch_b)
i( _no_torque)
i( _sim_blades)
i( _pitch_b)
i( _rotor_rpm)
i( _rel_len_hinge)
i( _maxcyclicail)
i( _maxcyclicele)
i( _mincyclicail)
i( _mincyclicele)
i( _delta3)
i( _delta)
i( _dynamic)
i( _translift)
i( _c2)
i( _stepspersecond)
i( _engineon)
i( _alphamin) i(_alphamax) i(_alpha0) i(_alpha0factor)
i( _teeterdamp) i(_maxteeterdamp)
i( _rellenteeterhinge)
i( _translift_ve)
i( _translift_maxfactor)
i( _ground_effect_constant)
i( _vortex_state_lift_factor)
i( _vortex_state_c1)
i( _vortex_state_c2)
i( _vortex_state_c3)
i( _vortex_state_e1)
i( _vortex_state_e2)
i( _vortex_state_e3)
i( _lift_factor) i(_f_ge) i(_f_vs) i(_f_tl)
i( _vortex_state)
i( _liftcoef)
i( _dragcoef0)
i( _dragcoef1)
i( _twist) //outer incidence = inner inner incidence + _twist
i( _rel_len_where_incidence_is_measured)
i( _torque_of_inertia)
i( _rel_len_blade_start)
i( _incidence_stall_zero_speed)
i( _incidence_stall_half_sonic_speed)
i( _lift_factor_stall)
i( _stall_change_over)
i( _drag_factor_stall)
i( _stall_sum)
i( _stall_v2sum)
i( _yaw)
i( _roll)
i( _cyclicail)
i( _cyclicele)
i( _cyclic_factor) <<endl;
int j;
for(j=0; j<r._rotorparts.size(); j++) {
out << *((Rotorpart*)r._rotorparts.get(j));
}
out <<endl << endl;
#undef i
#undef iv
return out;
}
void Rotor:: writeInfo()
{
#ifdef TEST_DEBUG
std::ostringstream buffer;
buffer << *this;
FILE*f=fopen("c:\\fgmsvc\\bat\\log.txt","at");
if (!f) f=fopen("c:\\fgmsvc\\bat\\log.txt","wt");
if (f)
{
fprintf(f,"%s",(const char *)buffer.str().c_str());
fclose (f);
}
#endif
}
Rotorpart* Rotor::newRotorpart(float* pos, float *posforceattac, float *normal,
float* speed,float *dirzentforce, float zentforce,float maxpitchforce,
float maxpitch, float minpitch, float mincyclic,float maxcyclic,
@ -980,27 +1222,8 @@ Rotorpart* Rotor::newRotorpart(float* pos, float *posforceattac, float *normal,
p(number_of_segments)
p(rel_len_where_incidence_is_measured)
p(rel_len_blade_start)
p(rotor_correction_factor)
#undef p
SG_LOG(SG_FLIGHT, SG_DEBUG, setprecision(8)
<< "newrp: pos("
<< pos[0] << ' ' << pos[1] << ' ' << pos[2]
<< ") pfa ("
<< posforceattac[0] << ' ' << posforceattac[1] << ' '
<< posforceattac[2] << ')');
SG_LOG(SG_FLIGHT, SG_DEBUG, setprecision(8)
<< " nor("
<< normal[0] << ' ' << normal[1] << ' ' << normal[2]
<< ") spd ("
<< speed[0] << ' ' << speed[1] << ' ' << speed[2] << ')');
SG_LOG(SG_FLIGHT, SG_DEBUG, setprecision(8)
<< " dzf("
<< dirzentforce[0] << ' ' << dirzentforce[1] << dirzentforce[2]
<< ") zf (" << zentforce << ") mpf (" << maxpitchforce << ')');
SG_LOG(SG_FLIGHT, SG_DEBUG, setprecision(8)
<< " pit(" << minpitch << ".." << maxpitch
<< ") mcy (" << mincyclic << ".." << maxcyclic
<< ") d3 (" << delta3 << ')');
return r;
}
@ -1067,13 +1290,14 @@ void Rotorgear::calcForces(float* torqueOut)
else
rel_torque_engine=0;
//add the rotor brake
//add the rotor brake and the gear fritcion
float dt=0.1f;
if (r0->_rotorparts.size()) dt=((Rotorpart*)r0->_rotorparts.get(0))->getDt();
float rotor_brake_torque;
rotor_brake_torque=_rotorbrake*_max_power_rotor_brake;
rotor_brake_torque=_rotorbrake*_max_power_rotor_brake+_rotorgear_friction;
//clamp it to the value you need to stop the rotor
//to avod accelerate the rotor to neagtive rpm:
rotor_brake_torque=Math::clamp(rotor_brake_torque,0,
total_torque_of_inertia/dt*omegarel);
max_torque_of_engine-=rotor_brake_torque;
@ -1092,7 +1316,7 @@ void Rotorgear::calcForces(float* torqueOut)
float lim1=-total_torque/total_torque_of_inertia;
//accel. by autorotation
if (lim1<_engine_accell_limit) lim1=_engine_accell_limit;
if (lim1<_engine_accel_limit) lim1=_engine_accel_limit;
//if the accel by autorotation greater than the max. engine
//accel, then this is the limit, if not: the engine is the limit
if (_ddt_omegarel>lim1) _ddt_omegarel=lim1;
@ -1125,6 +1349,7 @@ void Rotorgear::calcForces(float* torqueOut)
}
}
}
_total_torque_on_engine=total_torque+_ddt_omegarel*total_torque_of_inertia;
}
}
@ -1173,24 +1398,36 @@ void Rotorgear::setParameter(char *parametername, float value)
p(yasimdragfactor,1)
p(yasimliftfactor,1)
p(max_power_rotor_brake,1000)
p(engine_accell_limit,0.01)
cout << "internal error in parameter set up for rotorgear: '"
<< parametername <<"'" << endl;
p(rotorgear_friction,1000)
p(engine_accel_limit,0.01)
SG_LOG(SG_INPUT, SG_ALERT,
"internal error in parameter set up for rotorgear: '"
<< parametername <<"'" << endl);
#undef p
}
int Rotorgear::getValueforFGSet(int j,char *text,float *f)
{
if (j==0)
{
sprintf(text,"/rotors/gear/total_torque");
*f=_total_torque_on_engine;
} else return 0;
return j+1;
}
Rotorgear::Rotorgear()
{
_in_use=0;
_engineon=0;
_rotorbrake=0;
_max_power_rotor_brake=1;
_rotorgear_friction=1;
_max_power_engine=1000*450;
_engine_prop_factor=0.05f;
_yasimdragfactor=1;
_yasimliftfactor=1;
_ddt_omegarel=0;
_engine_accell_limit=0.05f;
_engine_accel_limit=0.05f;
_total_torque_on_engine=0;
}
Rotorgear::~Rotorgear()

32
src/FDM/YASim/Rotor.hpp
View file

@ -11,9 +11,11 @@ namespace yasim {
class Surface;
class Rotorpart;
class Ground;
const float rho_null=1.184f; //25DegC, 101325Pa
class Rotor {
friend std::ostream & operator<<(std::ostream & out, /*const*/ Rotor& r);
private:
float _torque;
float _omega,_omegan,_omegarel,_ddt_omega,_omegarelneu;
@ -25,6 +27,9 @@ private:
float _airfoil_drag_coefficient0;
float _airfoil_drag_coefficient1;
int _ccw;
int _number_of_blades;
int _number_of_segments;
int _number_of_parts;
public:
Rotor();
@ -78,7 +83,7 @@ public:
void getTip(float* tip);
void calcLiftFactor(float* v, float rho, State *s);
void getDownWash(float *pos, float * v_heli, float *downwash);
float getGroundEffect(float* posOut);
int getNumberOfBlades(){return _number_of_blades;}
// Query the list of Rotorpart objects
int numRotorparts();
@ -102,24 +107,36 @@ public:
float getOmegan() {return _omegan;}
float getTaper() { return _taper;}
float getChord() { return _chord;}
int getNumberOfParts() { return _number_of_parts;}
float getOverallStall()
{if (_stall_v2sum !=0 ) return _stall_sum/_stall_v2sum; else return 0;}
float getAirfoilIncidenceNoLift() {return _airfoil_incidence_no_lift;}
Vector _rotorparts;
void findGroundEffectAltitude(Ground * ground_cb,State *s);
void writeInfo();
private:
void strncpy(char *dest,const char *src,int maxlen);
void interp(float* v1, float* v2, float frac, float* out);
float calcStall(float incidence,float speed);
float findGroundEffectAltitude(Ground * ground_cb,State *s,
float *pos0,float *pos1,float *pos2,float *pos3,
int iteration=0,float a0=-1,float a1=-1,float a2=-1,float a3=-1);
Rotorpart* newRotorpart(float* pos, float *posforceattac, float *normal,
float* speed,float *dirzentforce, float zentforce,float maxpitchforce,
float maxpitch, float minpitch, float mincyclic,float maxcyclic,
float delta3,float mass,float translift,float rellenhinge,float len);
float _base[3];
float _groundeffectpos[4][3];
float _ground_effect_altitude;
//some postions, where to calcualte the ground effect
float _normal[3];//the normal vector (direction of rotormast, pointing up)
float _normal_with_yaw_roll[3];//the normal vector (perpendicular to rotordisc)
float _forward[3];
float _diameter;
int _number_of_blades;
float _weight_per_blade;
float _rel_blade_center;
float _min_pitch;
@ -166,7 +183,6 @@ private:
float _dragcoef0;
float _dragcoef1;
float _twist; //outer incidence = inner inner incidence + _twist
int _number_of_segments;
float _rel_len_where_incidence_is_measured;
float _torque_of_inertia;
float _rel_len_blade_start;
@ -179,7 +195,12 @@ private:
float _stall_v2sum;
float _yaw;
float _roll;
float _cyclicail;
float _cyclicele;
float _cyclic_factor;
float _rotor_correction_factor;
};
std::ostream & operator<<(std::ostream & out, /*const*/ Rotor& r);
class Rotorgear {
private:
@ -191,8 +212,10 @@ private:
float _yasimliftfactor;
float _rotorbrake;
float _max_power_rotor_brake;
float _rotorgear_friction;
float _ddt_omegarel;
float _engine_accell_limit;
float _engine_accel_limit;
float _total_torque_on_engine;
Vector _rotors;
public:
@ -218,6 +241,7 @@ public:
Vector* getRotors() { return &_rotors;}
void initRotorIteration(float *lrot,float dt);
void getDownWash(float *pos, float * v_heli, float *downwash);
int getValueforFGSet(int j,char *b,float *f);
};
}; // namespace yasim

140
src/FDM/YASim/Rotorpart.cpp
View file

@ -45,6 +45,8 @@ Rotorpart::Rotorpart()
_lastrp=0;
_nextrp=0;
_oppositerp=0;
_last90rp=0;
_next90rp=0;
_translift=0;
_dynamic=100;
_c2=0;
@ -62,6 +64,8 @@ Rotorpart::Rotorpart()
_diameter=10;
_torque_of_inertia=0;
_rel_len_blade_start=0;
_torque=0;
_rotor_correction_factor=0.6;
}
void Rotorpart::inititeration(float dt,float *rot)
@ -73,10 +77,10 @@ void Rotorpart::inititeration(float dt,float *rot)
float a=Math::dot3(rot,_normal);
if(a>0)
_alphaalt=_alpha*Math::cos(a)
+_nextrp->getrealAlpha()*Math::sin(a);
+_next90rp->getrealAlpha()*Math::sin(a);
else
_alphaalt=_alpha*Math::cos(a)
+_lastrp->getrealAlpha()*Math::sin(-a);
+_last90rp->getrealAlpha()*Math::sin(-a);
//calculate the rotation of the fuselage, determine
//the part in the same direction as alpha
//and add it ro _alphaalt
@ -102,8 +106,10 @@ void Rotorpart::setParameter(char *parametername, float value)
p(number_of_segments)
p(rel_len_where_incidence_is_measured)
p(rel_len_blade_start)
cout << "internal error in parameter set up for rotorpart: '"
<< parametername <<"'" << endl;
p(rotor_correction_factor)
SG_LOG(SG_INPUT, SG_ALERT,
"internal error in parameter set up for rotorpart: '"
<< parametername <<"'" << endl);
#undef p
}
@ -281,14 +287,14 @@ float Rotorpart::getAlpha(int i)
i=i&1;
if (i==0)
return _alpha*180/3.14;//in Grad = 1
return _alpha*180/pi;//in Grad = 1
else
{
if (_alpha2type==1) //yaw or roll
return (getAlpha(0)-_oppositerp->getAlpha(0))/2;
else //collective
return (getAlpha(0)+_oppositerp->getAlpha(0)+
_nextrp->getAlpha(0)+_lastrp->getAlpha(0))/4;
_next90rp->getAlpha(0)+_last90rp->getAlpha(0))/4;
}
}
float Rotorpart::getrealAlpha(void)
@ -298,11 +304,7 @@ float Rotorpart::getrealAlpha(void)
void Rotorpart::setAlphaoutput(char *text,int i)
{
SG_LOG(SG_FLIGHT, SG_DEBUG, "setAlphaoutput rotorpart ["
<< text << "] typ" << i);
strncpy(_alphaoutputbuf[i>0],text,255);
if (i>0) _alpha2type=i;
}
@ -339,11 +341,13 @@ void Rotorpart::setCyclic(float pos)
}
void Rotorpart::setlastnextrp(Rotorpart*lastrp,Rotorpart*nextrp,
Rotorpart *oppositerp)
Rotorpart *oppositerp,Rotorpart*last90rp,Rotorpart*next90rp)
{
_lastrp=lastrp;
_nextrp=nextrp;
_oppositerp=oppositerp;
_last90rp=last90rp;
_next90rp=next90rp;
}
void Rotorpart::strncpy(char *dest,const char *src,int maxlen)
@ -373,7 +377,10 @@ float Rotorpart::calculateAlpha(float* v_rel_air, float rho,
if((_nextrp==NULL)||(_lastrp==NULL)||(_rotor==NULL))
return 0.0;//not initialized. Can happen during startupt of flightgear
if (returnlift!=NULL) *returnlift=0;
float flap_omega=(_nextrp->getrealAlpha()-_lastrp->getrealAlpha())
/*float flap_omega=(_nextrp->getrealAlpha()-_lastrp->getrealAlpha())
*_omega / pi*_rotor->getNumberOfParts()/4; hier mal die alte version probieren?
*/
float flap_omega=(_next90rp->getrealAlpha()-_last90rp->getrealAlpha())
*_omega / pi;
float local_width=_diameter*(1-_rel_len_blade_start)/2.
/(float (_number_of_segments));
@ -382,8 +389,8 @@ float Rotorpart::calculateAlpha(float* v_rel_air, float rho,
float rel = (n+.5)/(float (_number_of_segments));
float r= _diameter *0.5 *(rel*(1-_rel_len_blade_start)
+_rel_len_blade_start);
float local_incidence=incidence+_twist *rel - _twist
*_rel_len_where_incidence_is_measured;
float local_incidence=incidence+_twist *rel -
_twist *_rel_len_where_incidence_is_measured;
float local_chord = _rotor->getChord()*rel+_rotor->getChord()
*_rotor->getTaper()*(1-rel);
float A = local_chord * local_width;
@ -408,12 +415,23 @@ float Rotorpart::calculateAlpha(float* v_rel_air, float rho,
//split into direction and magnitude
v_local_scalar=Math::mag3(v_local);
if (v_local_scalar!=0)
Math::unit3(v_local,v_local);
//Math::unit3(v_local,v_help);
Math::mul3(1/v_local_scalar,v_local,v_help);
float incidence_of_airspeed = Math::asin(Math::clamp(
Math::dot3(v_local,_normal),-1,1)) + local_incidence;
Math::dot3(v_help,_normal),-1,1)) + local_incidence;
ias = incidence_of_airspeed;
float lift_wo_cyc =
_rotor->getLiftCoef(incidence_of_airspeed-cyc,v_local_scalar)
//reduce the ias (Prantl factor)
float prantl_factor=2/pi*Math::acos(Math::exp(
-_rotor->getNumberOfBlades()/2.*(1-rel)
*Math::sqrt(1+1/Math::sqr(Math::tan(
pi/2-Math::abs(incidence_of_airspeed-local_incidence))))));
incidence_of_airspeed = (incidence_of_airspeed+
_rotor->getAirfoilIncidenceNoLift())*prantl_factor
*_rotor_correction_factor-_rotor->getAirfoilIncidenceNoLift();
ias = incidence_of_airspeed;
float lift_wo_cyc = _rotor->getLiftCoef(incidence_of_airspeed
-cyc*_rotor_correction_factor,v_local_scalar)
* v_local_scalar * v_local_scalar * A *rho *0.5;
float lift_with_cyc =
_rotor->getLiftCoef(incidence_of_airspeed,v_local_scalar)
@ -427,15 +445,28 @@ float Rotorpart::calculateAlpha(float* v_rel_air, float rho,
//angle between blade movement caused by rotor-rotation and the
//total movement of the blade
/* the next shold look like this, but this is the inner loop of
the rotor simulation. For small angles (and we hav only small
angles) the first order approximation works well
lift_moment += r*(lift * Math::cos(angle)
- drag * Math::sin(angle));
*torque += r*(drag * Math::cos(angle)
+ lift * Math::sin(angle));
*/
lift_moment += r*(lift * (1-angle*angle)
- drag * angle);
*torque += r*(drag * (1-angle*angle)
+ lift * angle);
if (returnlift!=NULL) *returnlift+=lift;
}
float alpha=Math::atan2(lift_moment,_zentipetalforce * _len);
//as above, use 1st order approximation
//float alpha=Math::atan2(lift_moment,_zentipetalforce * _len);
float alpha;
if ((_zentipetalforce >1e-8) || (_zentipetalforce <-1e-8))
alpha=lift_moment/(_zentipetalforce * _len);
else
alpha=0;
return (alpha);
}
@ -461,12 +492,13 @@ void Rotorpart::calcForce(float* v, float rho, float* out, float* torque,
//Angle of blade which would produce no vertical force (where the
//effective incidence is zero)
float cyc=_mincyclic+(_cyclic+1)/2*(_maxcyclic-_mincyclic);
//float cyc=_mincyclic+(_cyclic+1)/2*(_maxcyclic-_mincyclic);
float cyc=_cyclic;
float col=_minpitch+(_collective+1)/2*(_maxpitch-_minpitch);
_incidence=(col+cyc)-_delta3*_alphaalt;
//the incidence of the rotorblade due to control input reduced by the
//delta3 effect, see README.YASIM
float beta=_relamp*cyc+col;
//float beta=_relamp*cyc+col;
//the incidence of the rotorblade which is used for the calculation
float alpha,factor; //alpha is the flapping angle
@ -490,7 +522,7 @@ void Rotorpart::calcForce(float* v, float rho, float* out, float* torque,
&scalar_torque);
//calculate drag etc., e. g. for deccelrating the rotor if engine
//is off and omega <10%
alpha = Math::clamp(alpha,_alphamin,_alphamax);
float rel =_omega*10 / _omegan;
alpha=rel * alpha + (1-rel)* _alpha0;
factor=_dt*_dynamic/10;
@ -505,11 +537,11 @@ void Rotorpart::calcForce(float* v, float rho, float* out, float* torque,
alpha=_alphaalt+(alpha-_alphaalt)*factor;
_alpha=alpha;
float meancosalpha=(1*Math::cos(_lastrp->getrealAlpha())
+1*Math::cos(_nextrp->getrealAlpha())
float meancosalpha=(1*Math::cos(_last90rp->getrealAlpha())
+1*Math::cos(_next90rp->getrealAlpha())
+1*Math::cos(_oppositerp->getrealAlpha())
+1*Math::cos(alpha))/4;
float schwenkfactor=1-(Math::cos(_lastrp->getrealAlpha())-meancosalpha);
float schwenkfactor=1-(Math::cos(_lastrp->getrealAlpha())-meancosalpha)*_rotor->getNumberOfParts()/4;
//missing: consideration of rellenhinge
float xforce = Math::cos(alpha)*_zentipetalforce;
@ -531,8 +563,62 @@ void Rotorpart::getAccelTorque(float relaccel,float *t)
int i;
float f=_rotor->getCcw()?1:-1;
for(i=0; i<3; i++) {
t[i]=f*-1* _normal[i]*relaccel*_omegan* _torque_of_inertia;
t[i]=f*-1* _normal[i]*relaccel*_omegan* _torque_of_inertia;// *_omeagan ?
_rotor->addTorque(-relaccel*_omegan* _torque_of_inertia);
}
}
std::ostream & operator<<(std::ostream & out, const Rotorpart& rp)
{
#define i(x) << #x << ":" << rp.x << endl
#define iv(x) << #x << ":" << rp.x[0] << ";" << rp.x[1] << ";" <<rp.x[2] << ";" << endl
out << "Writing Info on Rotorpart " << endl
i( _dt)
iv( _last_torque)
i( _compiled)
iv( _pos) // position in local coords
iv( _posforceattac) // position in local coords
iv( _normal) //direcetion of the rotation axis
i( _torque_max_force)
i( _torque_no_force)
iv( _speed)
iv( _direction_of_movement)
iv( _directionofzentipetalforce)
iv( _directionofrotorpart)
i( _zentipetalforce)
i( _maxpitch)
i( _minpitch)
i( _maxcyclic)
i( _mincyclic)
i( _cyclic)
i( _collective)
i( _delta3)
i( _dynamic)
i( _translift)
i( _c2)
i( _mass)
i( _alpha)
i( _alphaalt)
i( _alphamin) i(_alphamax) i(_alpha0) i(_alpha0factor)
i( _rellenhinge)
i( _relamp)
i( _omega) i(_omegan) i(_ddt_omega)
i( _phi)
i( _len)
i( _incidence)
i( _twist) //outer incidence = inner inner incidence + _twist
i( _number_of_segments)
i( _rel_len_where_incidence_is_measured)
i( _rel_len_blade_start)
i( _diameter)
i( _torque_of_inertia)
i( _torque)
i (_alphaoutputbuf[0])
i (_alphaoutputbuf[1])
i( _alpha2type)
i(_rotor_correction_factor)
<< endl;
#undef i
#undef iv
return out;
}
}; // namespace yasim

13
src/FDM/YASim/Rotorpart.hpp
View file

@ -1,10 +1,12 @@
#ifndef _ROTORPART_HPP
#define _ROTORPART_HPP
#include <sstream>
#include <iostream>
namespace yasim {
class Rotor;
class Rotorpart
{
friend std::ostream & operator<<(std::ostream & out, const Rotorpart& rp);
private:
float _dt;
float _last_torque[3];
@ -54,7 +56,7 @@ namespace yasim {
float calculateAlpha(float* v, float rho, float incidence, float cyc,
float alphaalt, float *torque,float *returnlift=0);
void setlastnextrp(Rotorpart*lastrp,Rotorpart*nextrp,
Rotorpart *oppositerp);
Rotorpart *oppositerp,Rotorpart*last90rp,Rotorpart*next90rp);
void setTorque(float torque_max_force,float torque_no_force);
void setOmega(float value);
void setOmegaN(float value);
@ -69,10 +71,11 @@ namespace yasim {
void setParameter(char *parametername, float value);
void setRotor(Rotor *rotor);
void setTorqueOfInertia(float toi);
void writeInfo(std::ostringstream &buffer);
private:
void strncpy(char *dest,const char *src,int maxlen);
Rotorpart *_lastrp,*_nextrp,*_oppositerp;
Rotorpart *_lastrp,*_nextrp,*_oppositerp,*_last90rp,*_next90rp;
Rotor *_rotor;
float _pos[3]; // position in local coords
@ -109,14 +112,14 @@ namespace yasim {
int _number_of_segments;
float _rel_len_where_incidence_is_measured;
float _rel_len_blade_start;
float _rel_len_blade_measured;
float _diameter;
float _torque_of_inertia;
float _torque;
// total torque of rotor (scalar) for calculating new rotor rpm
char _alphaoutputbuf[2][256];
int _alpha2type;
float _rotor_correction_factor;
};
std::ostream & operator<<(std::ostream & out, const Rotorpart& rp);
}; // namespace yasim
#endif // _ROTORPART_HPP