help solving an equation using fzero

4 ビュー (過去 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 ExchangeNumerical Integration and Differentiation についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by