temporal discretisation and time step value

Hello,

I know I need to add temporal discretisation defined by the value dt (delta t) to my code. I should have it as an input parameter. I know which value to add, just not how to do it in matlab. I haven't found a way to do it yet so a little help would be helpful ! (I'll add a plot in the code when I've figured out how to do the rest)

Thx a lot.

clc;clear;
X= input ('distance X=');
a=9*X^2;
b=-183*X;
c=1000;
F= (a+b+c); %Fpeak
a1=20;
b1=-0.12*X^2;
c1=4.2*X;

Td= (a1+b1+c1); %loading duration td

disp('z(t) equation');
%have delta t as input, should vary depending on k and m
dt=0.0003; %this is the delta t
tmax=100; %max time, arbitrary value for now
t=0:dt:tmax;% part i'm struggling with
M= input ('M kg=');
K=input('K in MN/m=');
w= sqrt(K/M);
for t= 1:length(t) % <- dunno if this is right
if t <= Td
    z(t)=(F/K)*(1-cos(w*t))+(F/(K*Td))*((sin(w*t)/w)-t);
else
    z(t)=(F/(K*w*Td))*(sin(w*t)-sin(w*(t-Td)))-(F/K)*cos(w*t);
end
end
T= 2*pi/w; %blast wall's natural period

disp ('zt')
disp (z(t))
disp ('T in seconds');
disp (T)
disp ('displacement must be >100mm')
%plot(z(t));

[Anthony Gurunian]

if you know dt as a function of M and K then you have to make a function which returns dt after the users inputs M and K.

Answer 1

[Mathieu NOE]

hello

some suggestions / corrections for your code

do not forget to pay attention to units , otherwise you may have wrong pulsation calculation (see my comments in the code)

I tested it with arbitrary inputs , as I don't know what are the typical values for X, M, K

clc;clear;
% X= input ('distance X=');  % units ?? 
X= 0.1;  % MN debug 

a=9*X^2;
b=-183*X;
c=1000;
F= (a+b+c); %Fpeak
a1=20;
b1=-0.12*X^2;
c1=4.2*X;
Td= (a1+b1+c1); %loading duration td % units ??
disp('z(t) equation');

% M= input ('M kg='); % unit must be kg 
% K=input('K in MN/m='); % unit must be kg mega N / m - see comment line below 

M= 1 ; % unit must be kg  % MN debug 
K= 0.0001 ; % unit must be megaN / m - see comment line below  % MN debug 

K = K*1e6; % MN : do not forget to convert to N/m to get the correct w pulsation
w= sqrt(K/M); % in rad/s 
period = (2*pi) / w; % in second
dt = period/5; % 5 samples / per
tmax = 50; %max time, arbitrary value for now
t=0:dt:tmax;% this is now fixed

for ci= 1:length(t) % <- dunno if this is right %% do not use t as loop index !! 
    ti = t(ci); % t at current time index (avoid doing this multiple times in the eqs below)
    
    if ti <= Td
        z(ci)=(F/K)*(1-cos(w*ti))+(F/(K*Td))*((sin(w*ti)/w)-ti);
    else
        z(ci)=(F/(K*w*Td))*(sin(w*ti)-sin(w*(ti-Td)))-(F/K)*cos(w*ti);
    end
end

T= 2*pi/w; %blast wall's natural period : we are in line ! 
disp ('zt')
disp (z(ci))
disp ('T in seconds');
disp (T)
disp ('displacement must be >100mm')
plot(t,z);

[Michael]

Hello,

Thanks for your help.

The default values are: k= 1/MN/m, M= 200kg and X= 0m. For tmax, how do I know what is the correct value to use? For dt I changed and I put period/1000 rather than period over 5 (I thought it would be more precise). I don't understand what you've put on lines 24 and 25? What is ci for?

Thx a lot,

M

clc;clear;
% X= input ('distance X in m=');
X= 0;  % MN debug 
a=9*X^2;
b=-183*X;
c=1000;
F= (a+b+c); %Fpeak
a1=20;
b1=-0.12*X^2;
c1=4.2*X;
Td= (a1+b1+c1); %loading duration td in ms
Td=Td/1000
disp('z(t) equation');
% M= input ('M in kg='); % unit must be kg 
% K=input('K in MN/m='); % unit must be kg mega N / m - see comment line below 
M= 200 ; % in kg
K= 1; % MN/m
K = K*1e6; % convert to N/m
w= sqrt(K/M); % in rad/s 
period = (2*pi) / w; % in s
dt = period/1000; % 5 samples / per % I changed it to 1000 to make it more precise
tmax = 50; %max time, arbitrary value for now
t=0:dt:tmax;
for ci= 1:length(t) % <- dunno if this is right %% do not use t as loop index !! 
    ti = t(ci); % t at current time index (avoid doing this multiple times in the eqs below)
    
    if ti <= Td
        z(ci)=(F/K)*(1-cos(w*ti))+(F/(K*Td))*((sin(w*ti)/w)-ti);
    else
        z(ci)=(F/(K*w*Td))*(sin(w*ti)-sin(w*(ti-Td)))-(F/K)*cos(w*ti);
    end
end
T= 2*pi/w; %blast wall's natural period in s : we are in line ! 
disp ('zt')
disp (z(ci))
disp ('T in seconds');
disp (T)
disp ('displacement must be >100mm')
%plot(t,z);

[Mathieu NOE]

hello

lines 24 and 25 : i changes your code because cannot be a time vector and a loop index (scalar) as the same time.

I guess you wanted to do a kind of vectorized code, which btw avoids to do a time consuming for loop. see the new code below.

Also the simulation time can be greatly reduced , as I don't see the need to go beyong t max = 1 period. Anyway you just looking at a repetitive sine wave, so plotting thousands of periods of it is useless - unless there is a good reason to do so.

I plotted in two colors the z values for before and after Td time - just for fun - maybe this can visually explain why no need to do this simulation on very long time vector.

all the best

clc;clear;
% X= input ('distance X in m=');
X= 0;  % MN debug 
a=9*X^2;
b=-183*X;
c=1000;
F= (a+b+c); %Fpeak
a1=20;
b1=-0.12*X^2;
c1=4.2*X;
Td= (a1+b1+c1); %loading duration td in ms
Td=Td/1000; % convertion ms to s 
disp('z(t) equation');
% M= input ('M in kg='); % unit must be kg 
% K=input('K in MN/m='); % unit must be kg mega N / m - see comment line below 
M= 200 ; % in kg
K= 1; % MN/m
K = K*1e6; % convert to N/m
w= sqrt(K/M); % in rad/s 
period = (2*pi) / w; % in s
dt = period/1000; % 5 samples / per % I changed it to 1000 to make it more precise
% tmax = 50; %max time, arbitrary value for now
tmax = period; %max time, arbitrary value for now

t=0:dt:tmax;
% for ci= 1:length(t) % <- dunno if this is right %% do not use t as loop index !! 
%     ti = t(ci); % t at current time index (avoid doing this multiple times in the eqs below)
%     
%     if ti <= Td
%         z(ci)=(F/K)*(1-cos(w*ti))+(F/(K*Td))*((sin(w*ti)/w)-ti);
%     else
%         z(ci)=(F/(K*w*Td))*(sin(w*ti)-sin(w*(ti-Td)))-(F/K)*cos(w*ti);
%     end
% end
%% instead of a for loop , this code can be easily vectorized 
% period 1 : t <= Td
t1 = t(t<=Td);
z1=(F/K)*(1-cos(w*t1))+(F/(K*Td))*((sin(w*t1)/w)-t1);
% period 2 : t > Td
t2 = t(t>Td);
z2=(F/(K*w*Td))*(sin(w*t2)-sin(w*(t2-Td)))-(F/K)*cos(w*t2);


T= 2*pi/w; %blast wall's natural period in s : we are in line ! 
disp ('zt')
disp (z2(end))
disp ('T in seconds');
disp (T)
disp ('displacement must be >100mm')
plot(t1,z1,'b',t2,z2,'r');