temporal discretisation and time step value

조회 수: 13(최근 30일)
Michael
Michael 2021년 1월 7일
편집: Michael 2021년 1월 15일
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
Anthony Gurunian
Anthony Gurunian 2021년 1월 7일
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.

댓글을 달려면 로그인하십시오.

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 1월 8일
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
Mathieu NOE
Mathieu NOE 2021년 1월 13일
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');

댓글을 달려면 로그인하십시오.

추가 답변(0개)

Community Treasure Hunt

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

Start Hunting!

Translated by