MATLAB: Newton Raphson Method. Error: In an assignment A(I) = B, the number of elements in B and I must be the same. Error in Rotorcraft2 (line 61) theta(r)=x1;

??? in an assignment a(i) = berrorfuselage attitudehelicopternewton raphsonrotorcraftsea kingtrim angleswestland

I am writing a code to calculate the required trim angles for a sea king helicopter at velocities 0-70 m/s in steps of 5m/s. Data has been provided for me including an initial guess for fuselage attitude (theta).
I have completed my code and I am attempting to debug it. The first problem I'm having is getting the fuselage attitude angle iteration to run the entire way through. The code runs to line 61 suggesting that the code is working as it should, however it isn't as in line 43 (x1=theta_i(r);) a 15×1 column vector is produced when in fact x1 should be a single value assigned from the theta_i vector and correspond to the appropriate velocity that the iteration is being run for. Line 61 where the code fails should then store the final iterated value into the vector initialized before the loop, line 38 (theta=zeros(size(V)); %Empty Fuselage Attitude Matrix).
Any help to solving why x1 is not a single value would be great as this should solve the problem…
Thanks, the code is below.
Oliver
(The line calls may be incorrect since I spaced the lines so that reading is easier)
clear all
AUM=10000; %All up mass [kg]
AUW=AUM*9.81; %All up weight [N]
XG=0; ZG=0; %CoG coordinates [m]
XR=0; ZR=2; %Rotor coordinates [m]
XA=0; ZA=1; %Aerodynamic origin coordinates [m]
L_100=0; D_100=12000; M_100=0; %Lift/Drag/Moment at 100 m/s
N=5; %Number of blades
Lamba_beta=1; %Flapping Freqency
I_beta=2400; %Flapping inertia of blade [kg/m^3]
K_beta=0; %Spring stiffness constant
Ms=0; %Main rotor hub moment per unit disc tilt (a1)
gamma=11; %Lock number
V_T=200; %Rotor Tip Speed [m/s]
R=10; %Rotor Radius [m]
c=0.5; %Blade Chord [m]
alpha_s=0*(pi/180); %Shaft Tilt angle [rad]
sigma=1; %Density ratio
omega=20; %Rotor speed [rad/s]
a=5.28; %liftslope [1/rad]
s=0.079577472; %Solidity
Vf=0:5:70; %Forward Velocity Increments [m/s]
V=Vf';
%Lift, Drag and Pitching moment generated by the fuselage at each velocity increment
LF=zeros(size(V));
Df=[0 30 120 270 480 750 1080 1470 1920 2430 3000 3630 4320 5070 5880];
DF=Df';
Mf=zeros(size(V));
A=-DF; %Lateral Cyclic Pitch angle
B=2*AUW; %Longitudinal cyclic pitch angle
tf=atan(DF./(LF-AUW));%Thrust tilt angle
theta=zeros(size(V)); %Empty Fuselage Attitude Matrix
theta_initial=[0.0001 0 -0.000152905 -0.00611621 -0.001376146 -0.002446478 -0.003822611 -0.005504532 -0.007492215 -0.00978562 -0.012384688 -0.015289328 -0.018499418 -0.022014791 -0.025835229]; %Initial guess for theta
theta_i=theta_initial'; %initial guess for theta
for r=1:1:15
x1=theta_i(r);
itermax=100;
iter=0;
tol=0.0001;
converge=1;
while converge>tol
F=A(r)*cos(x1)+B*sin(x1)-Ms*x1+((Ms*tf(r))-(Ms*alpha_s)+Mf);
dF=-A(r)*sin(x1)+B*cos(x1)-Ms;
upd=-F./dF;
x1=x1+upd;
converge=abs(upd);
iter=iter+1;
if iter<itermax;
continue
else
break
end
end
theta(r)=x1;
end

Best Answer

  • That was the problem.
    Thankyou