# How to Iteratively changed C during Eulers Method ?

1 view (last 30 days)
Wilfredo Huaman on 14 Jun 2019
Answered: James Browne on 15 Jun 2019
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 %%%%%%%%%
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);
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.')

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 %%%%%%%%%
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);
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.')
so you may have some additional debugging to do, but I think this might be close to what you need to accomplish your goals~