My integral results are not correct.

I'm writing a code for the propeller design, but the integral values ​​and results are not coming out as I want. I use this publicition. https://doi.org/10.2514/3.23779.

newzeta_ps needs to be equal to zeta. Also, Tc should be 0.0498.

There are 4 derivative equations in the code.

I looked many times but couldn't solve the problem.

clear

clc

V = 49.1744; % [m/s]

epsilon = 0.016789342857143;

zeta = 0.214700101676;

B = 2;

alpha = deg2rad(3.5007);

Cl = 0.7;

rho = 1.225;

D = 1.7526; % [m]

hub = 0.1524; % [m]

R = D/2; % [m]

r = linspace(hub, R, 7);

PropRPM = 2400;

i = 1;

PropRPS = PropRPM/60; %[dev/s]

omega = 2*pi*PropRPS;

x = omega/V*r(i);

xi = r(i)/R;

lambda = V/(omega*R);

J = V/(PropRPS*D);

phi_t = lambda*(1+(zeta/2)); rad2deg(phi_t);

phi = atan((tan(phi_t))/xi); rad2deg(phi);

f = (B/2)*(1-xi)/(sin(phi_t));

F = (2/pi)*acos(exp(-f));

a = (zeta/2)*((cos(phi))^2)*(1-(epsilon*tan(phi)));

aprime = (zeta/(2*x))*((cos(phi))*(sin(phi)))*(1+(epsilon/tan(phi)));

W = (V*(1+a))/sin(phi);

Mach = W/343;

G = F*x*cos(phi)*sin(phi);%dogru

gamma = (2*pi*(V^2)*zeta*G)/(B*omega);

c = (gamma*2)/(W*Cl);

c_ft = c/0.3048;

sigma = (B*c)/(2*pi*r(i));

xint = xi;

I1 = integral(@(xint) 4*xint*G*(1-(epsilon*tan(phi))),xint,1);

I2 = integral(@(xint) lambda*(4*xint*G*(1-epsilon*tan(phi)))/2.*xint*(1+epsilon/tan(phi))*sin(phi)*cos(phi),xint,1);

J1 = integral(@(xint) 4*xint*G*(1+epsilon/tan(phi)),xint,1);

J2 = integral(@(xint) ((4*xint*G*(1+epsilon/tan(phi)))/2)*(1-epsilon*tan(phi))*((cos(phi))^2),xint,1);

% Pc = (SHP/(2*rho*(PropRPS^3)*(D^5)))

Pc = 0.0402;

ccc = J1/2*J2;

newzeta_ps = (-1*(ccc))+sqrt((ccc^2)+Pc/J2);

Tc = (I1*zeta)-(I2*(zeta^2));

eta = (Tc/Pc)*J;

formatSpec = '%9.2f %9.3f %9.3f %9.4f %9.2f %9.4f %9.4f \n';

disp(' phi zeta newzeta Tc eta a aprime')

disp(' ______ ______ ______ ______ ______ ______ ______')

fprintf(formatSpec,rad2deg(phi),zeta, newzeta_ps, Tc, eta, a, aprime);

[Torsten]

xint as lower bound of integration for your integrals is wrong.

It mustn't be zeta, but zeta0 which is a fixed value.

[Emrah Can Ucar]

I think you should have said xi instead of zeta. because that's how it's shown in the article.

but anyway, my xi value never changes. always xi = 0.17.

I get the same values ​​when I change it to xi = 0.17.