How to Iteratively changed C during Eulers Method ?

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.')

Answer 1

[James Browne]

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~