help solving an equation using fzero

6 次查看（过去 30 天）
Terry Poole
Terry Poole 2021-2-9
回答： Star Strider 2021-2-9
I need to solve an equation
CT = (1/2)*sigma*r.^2*a.*(theta_75+(theta_tw*(r-0.75))-(lambda_i/r'))
by first using the trapz function to integrate this equation, and then using the fzero function to solve CT-CT_pres = 0, using theta_75 as an independent variable.
This is my code, I'm not sure if i'm using trapz correctly, and I can't get fzero to work at all, it's just one error after the other. Also i'm basically solving all this three times simultaniously by declaring theta_tw = [0,-10,-20] so I've tried to solve each theta_tw individually but still no luck.
Any help would be appreciated.
clear all; close all; clc;
Weight = 20000; % lbs
T = Weight; % Thrust = Weight (lbs)
rho = 0.00237; % Density (slug/ft^3)
R = 30; % Radius (ft)
A = 2827; % Area (ft^2)
v_tip = 650; % Tip Speed (ft/sec)
kappa = 1.15; % Induced Power Factor
Nb = 4; % Number of Blades
c = 2; % Chord (ft)
CD = 0.01; % Drag Coefficient
sigma = 0.085; % Solidity
a = 6; % Lift Curve Slope
%% Problem 2
% Rotor Disk Loading
DL = T/A; % lbs/ft^2
% Induced Velocity
v = sqrt(DL/(2*rho)); % ft/sec
% Ideal Induced Power
Pi_Ideal = (T*v)/550; % hp
% Ideal Power Loading
PL_Ideal = T/Pi_Ideal; % lbs/hp
% Non-Dimensional Coefficients at Hover
CT_pres = T/(rho*A*v_tip^2); % Thrust
CPi = kappa*(CT_pres*sqrt(CT_pres/2)); % Induced Power
CP0 = (1/8)*sigma*CD; % Rotor Profile Power
CQ0 = CP0; % Torque
CP = CPi+CP0; % Rotor Total Power
CQ = CP; % Torque
% Dimensional Coefficients at Hover
Pi = (CPi*rho*A*v_tip^3)/550; % Induced Power (hp)
P0 = (CP0*rho*A*v_tip^3)/550; % Profile Power (hp)
P = Pi+P0; % Rotor Total Power (hp)
% Figure of Merit
FM = Pi/P;
% Actual Power Loading
PL = T/P; % lbs/hp
fprintf(['----------- Problem 2 Answers -----------\n\n'...
'Rotor Disk Loading = %f lbs/ft^2 \n'...
'Ideal Induced Power = %f hp \nIdeal Power Loading = %f hp \n'...
'Thrust Coefficient = %f \nTorque Coefficient = %f \n'...
'Figure of Merit = %f \nProfile Power = %f hp \n'...
'Actual Power Loading = %f lbs/hp \n\n'],...
DL, Pi_Ideal,PL_Ideal, CT_pres, CQ, FM, P0, PL);
%% Problem 3
lambda_i = sqrt(CT_pres/2); % Uniform Inflow
theta_tw = [0;-10;-20]; % Twist Rate (Degrees)
r = 0:0.1:1; % Rotor section from root(0) to tip(1).
% Estimated Twist angle at 75% Raidus. (Degrees)
theta_75 = (((6*CT_pres)/(sigma*a))+((3/2)*sqrt(CT_pres/2)))
% Integrating for CT
dct = dCT(sigma,r,a,theta_75,theta_tw,lambda_i)'
ct = trapz(dct)
% Solving for theta_.75
fun = @(sigma,r,a,X,theta_tw,lambda_i)dCT;
X = 0;
fzero = (fun-Ct_pres,X);
%theta_r = (((6*CT)/(sigma*a))+((3/2)*sqrt(CT/2)))+(theta_tw*(r-0.75))
function lambda = lambda_r(sigma,a,theta_r)
lambda = (1/16)*(sqrt(1+(32/(sigma*a))*theta_r)-1);
end
function CT = dCT(sigma,r,a,theta_75,theta_tw,lambda_i)
CT = (1/2)*sigma*r.^2*a.*(theta_75+(theta_tw*(r-0.75))-(lambda_i/r'));
end

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2021-2-9
Since you specifically asked for help on fzero, this works:
X = 1;
for k = 1:numel(theta_tw)
T_75(k) = fzero(@(theta_75)dCT(sigma,R,a,theta_75,theta_tw(k),lambda_i)-CT_pres, X)
end
producing:
T_75 =
0.00201198547991055 292.50201198548 585.00201198548
Note that since ‘theta_tw’ is a vector, it’s necessary to iterate through its elements. The values for ‘T_75’ correspond to those eleements. Also, it is never appropriate to begin with an initial parameter estimate of 0, especially in a root-finding function such as fzero. I set ‘X’ to 1 for this reason.
The rest of your code appears to run without error with the fzero change I am posting here. You will need to determine if it produces the correct results. Note that it may be necessary to iterate through ‘T_75’ (name it whatever you want) if you are using non-vectorised code with its results.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by