How to get the propeller aerodynamic coefficient ?

[Dobid]

Hello !

As explained in #843, I'm trying to implement an autopilot for the x8 UAV based on a PID controller by following this book : Small Unmanned Aircrafts by Beard and McLain (here's the fdm).

In the book, in the longitudinal autopilot part, they give a set of formulas to compute PID gains for the airspeed controller :
$k_{iv} = {\omega_{nv}}^2/a_{v2}$ ; $k_{pv} = (2\xi_{v}\omega_{nv} - a_{v1})/a_{v2}$
Those formulas show the following coefficients : $a_{v1}, a_{v2}$.

Those are transfer function coefficients, their formulas are :

With $S_{prop}$ being the area swept by the propeller and $C_{prop}$ the aerodynamical coefficient of the propeller, for more context, they appear here, when finding the formula of the thrust produced by the propeller/engine :

So in order to achieve my PID, I need to find a way to get or compute this value of $C_{prop}$ and I have no idea of where to look for. I checked the fdm (propulsion tag) referencing the engine's xml file and the direct thrust file.

Can someone help me find a way to compute this $C_{prop}$ value ? The x8 UAV description table gives the dimensions of the propeller, but searching the internet about a corresponding $C_{prop}$ doesn't bring anything. This $C_{prop}$ coefficient notation seems to be specific to the book I'm following and I can't link it to information I find online (UAV specs, propeller specs, fdm files...)

Thanks in advance !

[seanmcleod]

I don't have a definite answer, but take a look at the following code:

jsbsim/src/models/propulsion/FGPropeller.cpp

Lines 207 to 298 in ebb770e

	double FGPropeller::Calculate(double EnginePower)
	{
	FGColumnVector3 vDXYZ = MassBalance->StructuralToBody(vXYZn);
	const FGMatrix33& mT = Transform();
	// Local air velocity is obtained from Stevens & Lewis' "Aircraft Control and
	// Simualtion (3rd edition)" eqn 8.2-1
	// Variables in.AeroUVW and in.AeroPQR include the wind and turbulence effects
	// as computed by FGAuxiliary.
	FGColumnVector3 localAeroVel = mT.Transposed() * (in.AeroUVW + in.AeroPQR*vDXYZ);
	double omega, PowerAvailable;
	double Vel = localAeroVel(eU);
	double rho = in.Density;
	double RPS = RPM/60.0;
	// Calculate helical tip Mach
	double Area = 0.25*Diameter*Diameter*M_PI;
	double Vtip = RPS * Diameter * M_PI;
	HelicalTipMach = sqrt(Vtip*Vtip + Vel*Vel) / in.Soundspeed;
	if (RPS > 0.01) J = Vel / (Diameter * RPS); // Calculate J normally
	else J = Vel / Diameter;
	PowerAvailable = EnginePower - GetPowerRequired();
	if (MaxPitch == MinPitch) { // Fixed pitch prop
	ThrustCoeff = cThrust->GetValue(J);
	} else { // Variable pitch prop
	ThrustCoeff = cThrust->GetValue(J, Pitch);
	}
	// Apply optional scaling factor to Ct (default value = 1)
	ThrustCoeff *= CtFactor;
	// Apply optional Mach effects from CT_MACH table
	if (CtMach) ThrustCoeff *= CtMach->GetValue(HelicalTipMach);
	Thrust = ThrustCoeff*RPS*RPS*D4*rho;
	// Induced velocity in the propeller disk area. This formula is obtained
	// from momentum theory - see B. W. McCormick, "Aerodynamics, Aeronautics,
	// and Flight Mechanics" 1st edition, eqn. 6.15 (propeller analysis chapter).
	// Since Thrust and Vel can both be negative we need to adjust this formula
	// To handle sign (direction) separately from magnitude.
	double Vel2sum = Vel*abs(Vel) + 2.0*Thrust/(rho*Area);
	if( Vel2sum > 0.0)
	Vinduced = 0.5 * (-Vel + sqrt(Vel2sum));
	else
	Vinduced = 0.5 * (-Vel - sqrt(-Vel2sum));
	// P-factor is simulated by a shift of the acting location of the thrust.
	// The shift is a multiple of the angle between the propeller shaft axis
	// and the relative wind that goes through the propeller disk.
	if (P_Factor > 0.0001) {
	double tangentialVel = localAeroVel.Magnitude(eV, eW);
	if (tangentialVel > 0.0001) {
	// The angle made locally by the air flow with respect to the propeller
	// axis is influenced by the induced velocity. This attenuates the
	// influence of a string cross wind and gives a more realistic behavior.
	double angle = atan2(tangentialVel, Vel+Vinduced);
	double factor = Sense * P_Factor * angle / tangentialVel;
	SetActingLocationY( GetLocationY() + factor * localAeroVel(eW));
	SetActingLocationZ( GetLocationZ() + factor * localAeroVel(eV));
	}
	}
	omega = RPS*2.0*M_PI;
	vFn(eX) = Thrust;
	vTorque(eX) = -Sense*EnginePower / max(0.01, omega);
	// The Ixx value and rotation speed given below are for rotation about the
	// natural axis of the engine. The transform takes place in the base class
	// FGForce::GetBodyForces() function.
	FGColumnVector3 vH(Ixx*omega*Sense*Sense_multiplier, 0.0, 0.0);
	if (omega > 0.01) ExcessTorque = PowerAvailable / omega;
	else ExcessTorque = PowerAvailable / 1.0;
	RPM = (RPS + ((ExcessTorque / Ixx) / (2.0 * M_PI)) * in.TotalDeltaT) * 60.0;
	if (RPM < 0.0) RPM = 0.0; // Engine won't turn backwards
	// Transform Torque and momentum first, as PQR is used in this
	// equation and cannot be transformed itself.
	vMn = in.PQRi*(mT*vH) + mT*vTorque;
	return Thrust; // return thrust in pounds
	}

I wonder if, without having a proper look whether $C_{prop}$ is the thrust_coefficient.

jsbsim/src/models/propulsion/FGPropeller.cpp

Lines 154 to 155 in ebb770e

	property_name = base_property_name + "/thrust-coefficient";
	PropertyManager->Tie( property_name.c_str(), this, &FGPropeller::GetThrustCoefficient );

jsbsim/src/models/propulsion/FGPropeller.cpp

Lines 118 to 120 in ebb770e

	if( (cPower == 0) || (cThrust == 0)){
	cerr << "Propeller configuration must contain C_THRUST and C_POWER tables!" << endl;
	}

jsbsim/src/models/propulsion/FGPropeller.cpp

Lines 103 to 104 in ebb770e

	if (name == "C_THRUST") {
	cThrust = new FGTable(PropertyManager, table_element);

[Dobid]

Thanks for the quick answer, I will go through the code and identify the laws used with the book mentionned in the comments to try and find some similarities with my $C_{prop}$ coefficient. I got assigned some other work before that, I'll keep this thread updated when I find something, or if I have another question.