How to change initial conditions and run loop again until R^2 is above 0.8
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
I am trying to create a combined pendulum model to model the force from a CFD solution and instead of guessing and checking I want to create a script to change conditions such as lengths (L and L2), initial angle(A1 and A2) and damping coefficient (c and c2) and re run the (2 pendulum forces) Fr and Fr2 in my script and compare the FR (combined force pendulum) to the Total Force Variable.
I currently an running a while loop with a nested if statement at the end to chnage the varaible but its not doing anything and the constants remain zero any help would be greatly appreciated.
Thanks, Willie
%%
CFDoutput = readmatrix('Cyl_3T_4_15_Fill_90sec.csv');
% Cut first 5 seconds from CFD
T = CFDoutput(:,1);
loc = find(T<5);
start = length(loc);
loc2 = find(T==80.0130);
endend = length(loc2);
time = T(start:end);
XForce = CFDoutput(:,8);
xForce = XForce(start:end);
YForce = CFDoutput(:,9);
yForce = YForce(start:end);
ZForce = CFDoutput(:,10);
zForce = ZForce(start:end);
TotalForce = sqrt(xForce.^2 + yForce.^2 + zForce.^2);
%% 2D pendulum motion
R2 = 0;
while R2 < 0.8
% System 1 Variables
mass=67500; % [kg] pendulum bob mass
L=5.9; % [m] pendulum rod length
g=9.81; % [m/s^2] acceleration due to gravity
t=90; % [s] simulation time
c =0;% 0.05; % Damping coefficient
dt= 2/91; % Time step
i= t/dt; % Number of iterations
A1 = 0; %initial angle deg
% System 2 Variables
mass2=67500; % [kg] pendulum bob mass
L2=2; % [m] pendulum rod length
c2 =0;% 0.1; % Damping coefficient
A2 = 0; % initial angle deg
% Allocation of Variables
theta=zeros(round(i),1);
omega=zeros(round(i),1);
alpha=zeros(round(i),1);
Fr=zeros(round(i),1);
theta(1,:)=deg2rad(A1); % initial angular position
omega(1,:)=0; % initial angular velocity
alpha(1,:)=0; % initial angular acceleration
% Solve Equations of Motion
for n=1:i
theta(n+1,:)= theta(n,:)+omega(n,:)*dt;
omega(n+1,:)= omega(n,:)+alpha(n,:)*dt;
alpha(n+1,:)=(-g*sin(theta(n+1,:)))/L-c*omega(n+1,:);
Fr(n+1,:)= mass*g*cos(theta(n+1,:))+ mass*L*(omega(n+1,:)).^2;
end
% % 2D pendulum motion 2
% Allocation of variables 2
t=0:i;
theta2=zeros(round(i),1);
omega2=zeros(round(i),1);
alpha2=zeros(round(i),1);
Fr2=zeros(round(i),1);
theta2(1,:)=deg2rad(A2); % initial angular position
omega2(1,:)=0; % initial angular velocity
alpha2(1,:)=0; % initial angular acceleration
% Solve Equations of Motion
for n=1:i
theta2(n+1,:)= theta2(n,:)+omega2(n,:)*dt;
omega2(n+1,:)= omega2(n,:)+alpha2(n,:)*dt;
alpha2(n+1,:)=(-g*sin(theta2(n+1,:)))/L2-c2*omega2(n+1,:);
Fr2(n+1,:)= mass2*g*cos(theta2(n+1,:))+ mass2*L2*(omega2(n+1,:)).^2;
end
% Combined Pendulum Force
FR = (Fr + Fr2)/2;
%R squared Value curve fitting
Fit = fitlm(TotalForce,FR);
R2 = Fit.Rsquared.Ordinary;
if R2 < 0.8
for j=1:i
c = c + 0.001;
c2 = c2 + 0.001;
L = L + 0.001;
L2 = L2 + 0.001;
break
end
end
2 Kommentare
Walter Roberson
am 3 Jul. 2019
Put everything except the initial conditions into a function. Have the function return R2. Now loop around providing new initial conditions and calling the function until R2 meets your criteria.
Antworten (0)
Siehe auch
Kategorien
Mehr zu Classical Mechanics finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!