How to put together different cases in one MATLAB code?

2 ビュー (過去 30 日間)
Dragana Skoko
Dragana Skoko 2019 年 7 月 20 日
コメント済み: Star Strider 2019 年 7 月 21 日
Hello everyone, here is exact matlab code for solving dymanic response using HHT direct integration method .But is written for just one case when alpha is -0.3. Yet, I have to solve same problem, but for more different values of parametar alpha (alpha =-0.16 , alpha =0)..Can anyone explain me how to put it together in one code (all threee cases)?Which funktion to use?
close all
clear all
clc
%----------------------------------------------------------------
% Hilbe-Hughes-Taylor (alpha) metod
%----------------------------------------------------------------
%----------------------------------------------------------------
% Definisanje parametara sistema
%----------------------------------------------------------------
M=2500 %masa [kg]
K=2250000 %krutost [N/m']
C=7500 %konstanta prigušenja [Ns/m']
dt=0.002 %vremenski korak integracije
Td=0.2 %ukupno trajanje impulsa [sec]
%--------------------------------------
% Početni uslovi
%--------------------------------------
X0=0 %početno pomjeranje [m]
X0d=0 %početna brzina [m/s]
F0=20000 %sila u trenutku t=0
X0dd=inv(M)*(F0-C*X0d-K*X0) %početno ubrzanje [m/s2]
%-----------------------------------------
% Definisanje parametara integracije
%-----------------------------------------
alpha=-0.3 %koeficijent prigušenja
beta=((1-alpha)^2)/4
gamma=(1-2*alpha)/2
%-------------------------------------------
% Definisanje integracionih konstanti (ai)
%-------------------------------------------
a0=1/(beta*(dt)^2);a1=gamma/(beta*dt);a2=1/(beta*dt);
a3=(1/(2*beta))-1; a4=(gamma/beta)-1; a5=(dt/2)*((gamma/beta)-2);
a6=dt*(1-gamma); a7=gamma*dt;
%-------------------------------------------------
% Formiranje efektivne matrice krutosti
%-------------------------------------------------
Keff=(1+alpha)*K+a0*M+a1*(1+alpha)*C
%---------------------------------------------------------
% Definisanje funkcije opterećenja *trougaoni impuls
%---------------------------------------------------------
i=1
v(1)=X0d
a(1)=X0dd
U(1)=X0
F(1)=F0
t=0
for t=dt:dt:Td
i=i+1
if t>=0.0 F(i)=F0
else if t<0.0 F(i)=0
end
end
%---------------------------------------------------------
% Proračun za svaki naredni vremenski korak t+delta t
%---------------------------------------------------------
Reff=(1+alpha)*F(i)-alpha*F(i-1)+M*(a0*U(i-1)+a2*v(i-1)+a3*a(i-1))+C*(1+alpha)*(a1*U(i-1)+a4*v(i-1)+a5*a(i-1))+alpha*(C*v(i-1)+K*U(i-1))
U(i)=inv(Keff)*Reff
a(i)=a0*(U(i)-U(i-1))-a2*v(i-1)-a3*a(i-1)
v(i)=v(i-1)+a6*a(i-1)+a7*a(i)
end
fprintf('\nvrijeme\t\tpomjeranje\tbrzina\tubrzanje\n')
i=1;
for t=0:dt:Td
fprintf('%f\t%f\t%f\t%f\n',t, U(i), v(i), a(i));
i=i+1
end
Umax=max(U)
vmax=max(v)
amax=max(a)
figure (1)
t=[0:dt:Td]
plot(t,U,'-')
xlabel('vrijeme t[s]')
ylabel('pomjeranje U[m]')
grid on
t=[0:dt:Td]
figure (2)
plot(t,v,'-')
xlabel('vrijeme t[s]')
ylabel('brzina V[m/s]')
grid on
t=[0:dt:Td]
figure (3)
plot(t,a,'-')
xlabel('vrijeme t[s]')
ylabel('ubrzanje a[m/s^2]')
grid on
  1 件のコメント
Guillaume
Guillaume 2019 年 7 月 20 日
"But is written for just one case when alpha is -1/3"
It looks to me like it's written for alpha = -0.3, which is a far cry from -1/3 (10% error).

サインインしてコメントする。

回答 (1 件)

Star Strider
Star Strider 2019 年 7 月 20 日
With this change:
alpha= [-1/3, -1/6, 0]; %koeficijent prigušenja <— CREATE VECTOR FOR ‘alpha’
and others to accommodate its use as a vector in the rest of your code), see if this does what you want:
%----------------------------------------------------------------
% Hilbe-Hughes-Taylor (alpha) metod
%----------------------------------------------------------------
%----------------------------------------------------------------
% Definisanje parametara sistema
%----------------------------------------------------------------
M=2500; %masa [kg]
K=2250000; %krutost [N/m']
C=7500; %konstanta prigušenja [Ns/m']
dt=0.002; %vremenski korak integracije
Td=0.2; %ukupno trajanje impulsa [sec]
%--------------------------------------
% Početni uslovi
%--------------------------------------
X0=0; %početno pomjeranje [m]
X0d=0; %početna brzina [m/s]
F0=20000; %sila u trenutku t=0
X0dd=inv(M)*(F0-C*X0d-K*X0); %početno ubrzanje [m/s2]
%-----------------------------------------
% Definisanje parametara integracije
%-----------------------------------------
alpha= [-1/3, -1/6, 0]; %koeficijent prigušenja <— CREATE VECTOR FOR ‘alpha’
beta=((1-alpha).^2)/4;
gamma=(1-2*alpha)/2;
%-------------------------------------------
% Definisanje integracionih konstanti (ai)
%-------------------------------------------
a0=1./(beta*(dt).^2);
a1=gamma./(beta*dt);
a2=1./(beta*dt);
a3=(1./(2*beta))-1;
a4=(gamma./beta)-1;
a5=(dt/2)*((gamma./beta)-2);
a6=dt*(1-gamma);
a7=gamma*dt;
%-------------------------------------------------
% Formiranje efektivne matrice krutosti
%-------------------------------------------------
Keff=(1+alpha).*K+a0*M+a1.*(1+alpha)*C;
%---------------------------------------------------------
% Definisanje funkcije opterećenja *trougaoni impuls
%---------------------------------------------------------
i=1;
v(1,:)=X0d*[1 1 1];
a(1,:)=X0dd*[1 1 1];
U(1,:)=X0*[1 1 1];
F(1,:)=F0*[1 1 1];
t=0;
for t=dt:dt:Td
i=i+1;
if t>=0.0
F(i)=F0
else if t<0.0
F(i)=0
end
end
%---------------------------------------------------------
% Proračun za svaki naredni vremenski korak t+delta t
%---------------------------------------------------------
Reff=(1+alpha).*F(i)-alpha.*F(i-1)+M*(a0*U(i-1)+a2*v(i-1)+a3*a(i-1))+C*(1+alpha).*(a1*U(i-1)+a4*v(i-1)+a5*a(i-1))+alpha*(C*v(i-1)+K*U(i-1));
% U(i)=inv(Keff)*Reff;
Q1 = (Keff).\Reff;
U(i,:)=(Keff).\Reff;
a(i,:)=a0*(U(i)-U(i-1))-a2*v(i-1)-a3*a(i-1);
v(i,:)=v(i-1)+a6*a(i-1)+a7*a(i);
end
fprintf('\nvrijeme\t\tpomjeranje\tbrzina\tubrzanje\n')
t=0:dt:Td;
for i = 1:numel(t)
for k = 1:numel(alpha)
fprintf('%f\t%f\t%f\t%f\n',t(i), U(i,k), v(i,k), a(i,k));
end
fprintf('\n')
end
Umax=max(U);
vmax=max(v);
amax=max(a);
figure (1)
t=[0:dt:Td];
plot(t,U,'-')
xlabel('vrijeme t[s]')
ylabel('pomjeranje U[m]')
grid on
t=[0:dt:Td];
figure (2)
plot(t,v,'-')
xlabel('vrijeme t[s]')
ylabel('brzina V[m/s]')
grid on
t=[0:dt:Td];
figure (3)
plot(t,a,'-')
xlabel('vrijeme t[s]')
ylabel('ubrzanje a[m/s^2]')
grid on
This runs without error. Make appropriate changes to get the result you want.
  5 件のコメント
Rik
Rik 2019 年 7 月 21 日
Run this as a function. That will ensure a clean workspace every time you run this code. Scripts should only be used for really basic things and debugging.
Star Strider
Star Strider 2019 年 7 月 21 日
@Rik — Thank you!
I do not have any problems running it several times in succession.
What line throws the error?

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by