temporal discretisation and time step value

12 views (last 30 days)
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));
  1 Comment
Anthony Gurunian
Anthony Gurunian on 7 Jan 2021
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.

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 8 Jan 2021
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);
  2 Comments
Michael
Michael on 12 Jan 2021
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
Mathieu NOE on 13 Jan 2021
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');

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!