How to Iteratively changed C during Eulers Method ?

2 views (last 30 days)
I am trying to match a force curve from excel to a pedulum curve by manipulating a damping constant c I want to run through the array Total Force which is a 921x1 and see if the next value is larger or smaller and then manipulate c by either adding or subratcing from it and then use that value during the euler solution. Right now the code changes the value but then it is constant and runs it at that number I need the number changing through out the entire iterative solution. Any help would be greatly appreciated, Thank you in advance.
%%%%%%%%% Force Output Data %%%%%%%%%
CFDoutput = readmatrix('Sphere_5mRadius_10%_2mAmp.csv');
time = CFDoutput(:,1);
xForce = smooth(CFDoutput(:,8));
yForce = smooth(CFDoutput(:,9));
zForce = smooth(CFDoutput(:,10));
TotalForce = smooth(sqrt(xForce.^2 + yForce.^2 + zForce.^2));
Average_Force = mean(TotalForce);
% CurveFit = polyfit(time,TotalForce,17);
% CurvePlot = polyval(CurveFit,time);
%%
%%%%%%%%% 2D Pendulum Model %%%%%%%%
% System Variables
mass=10000; % [kg] pendulum bob mass
L=0.65; % [m] pendulum rod length
g=9.81; % [m/s^2] acceleration due to gravity
t=6; % [s] simulation time
c = 0at tha;
dt=0.005; % Time step
i= t/dt; % Number of iterations
%%%%%%%%% Variable Damping Coefficient %%%%%%%%%
for j = 1:length(TotalForce)-1
if j <= length(TotalForce)
if TotalForce(j) >= TotalForce(j+1)
c = c + 0.01 ;
elseif TotalForce(j) <= TotalForce(j+1)
c = c - 0.01 ;
end
end
end
% Initial Conditons for Eulers Method
t=0:i;
theta=zeros(i,1);
omega=zeros(i,1);
alpha=zeros(i,1);
Fr=zeros(i,1);
theta(1,:)=deg2rad(90); % initial angular position
omega(1,:)=0; % initial angular velocity
alpha(1,:)=0; % initial angular acceleration
% Iteratively solve equations of motion using Euler's Method
for n=1:length(TotalForce)
theta(n+1,:)=theta(n,:)+omega(n,:)*dt; % new angular position
omega(n+1,:)=omega(n,:)+alpha(n,:)*dt; % new angular velocity
alpha(n+1,:)=(-g*sin(theta(n+1,:)))/L-c*omega(n+1,:); % new angular acceleration
Fr(n+1,:)=mass*g*cos(theta(n+1,:))+mass*L*(omega(n+1,:)).^2; % Reaction Force
end
%%%%%%%%% Force Plots %%%%%%%%%
hold on
%plot(time,TotalForce)%,time,CurvePlot)%time,xForce,time,yForce,time,zForce)
title('Pressure Forces')
xlabel('Time (s)')
ylabel('Force (N)')
%plot(t*dt,mass*g,'b.');
plot(time,TotalForce,t*dt,Fr','k.')

Answers (1)

James Browne
James Browne on 15 Jun 2019
Greetings,
Have you tried storing the values of c in a vector, in each iteration of the first for loop and the indexing said vector in the second for loop? If I understand your question correctly, I believe this should solve your issue. Try something like:
%%%%%%%%% Force Output Data %%%%%%%%%
CFDoutput = readmatrix('Sphere_5mRadius_10%_2mAmp.csv');
time = CFDoutput(:,1);
xForce = smooth(CFDoutput(:,8));
yForce = smooth(CFDoutput(:,9));
zForce = smooth(CFDoutput(:,10));
TotalForce = smooth(sqrt(xForce.^2 + yForce.^2 + zForce.^2));
Average_Force = mean(TotalForce);
% CurveFit = polyfit(time,TotalForce,17);
% CurvePlot = polyval(CurveFit,time);
%%
%%%%%%%%% 2D Pendulum Model %%%%%%%%
% System Variables
mass=10000; % [kg] pendulum bob mass
L=0.65; % [m] pendulum rod length
g=9.81; % [m/s^2] acceleration due to gravity
t=6; % [s] simulation time
c = 0at tha;
dt=0.005; % Time step
i= t/dt; % Number of iterations
%%%%%%%%% Variable Damping Coefficient %%%%%%%%%
for j = 1:length(TotalForce)-1
if j <= length(TotalForce)
if TotalForce(j) >= TotalForce(j+1)
c(j+1) = c(j) + 0.01 ;
elseif TotalForce(j) <= TotalForce(j+1)
c(j+1) = c(j) - 0.01 ;
end
end
end
% Initial Conditons for Eulers Method
t=0:i;
theta=zeros(i,1);
omega=zeros(i,1);
alpha=zeros(i,1);
Fr=zeros(i,1);
theta(1,:)=deg2rad(90); % initial angular position
omega(1,:)=0; % initial angular velocity
alpha(1,:)=0; % initial angular acceleration
% Iteratively solve equations of motion using Euler's Method
for n=1:length(TotalForce)
theta(n+1,:)=theta(n,:)+omega(n,:)*dt; % new angular position
omega(n+1,:)=omega(n,:)+alpha(n,:)*dt; % new angular velocity
alpha(n+1,:)=(-g*sin(theta(n+1,:)))/L-c(n)*omega(n+1,:); % new angular acceleration
Fr(n+1,:)=mass*g*cos(theta(n+1,:))+mass*L*(omega(n+1,:)).^2; % Reaction Force
end
%%%%%%%%% Force Plots %%%%%%%%%
hold on
%plot(time,TotalForce)%,time,CurvePlot)%time,xForce,time,yForce,time,zForce)
title('Pressure Forces')
xlabel('Time (s)')
ylabel('Force (N)')
%plot(t*dt,mass*g,'b.');
plot(time,TotalForce,t*dt,Fr','k.')
Please note: I cannot debug the modifications to your code as I do not have access to the file:
Sphere_5mRadius_10%_2mAmp.csv
so you may have some additional debugging to do, but I think this might be close to what you need to accomplish your goals~

Community Treasure Hunt

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

Start Hunting!

Translated by