How to put together different cases in one MATLAB code?

Question

Dragana Skoko 2019-7-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/472619-how-to-put-together-different-cases-in-one-matlab-code

评论： 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 个评论
显示 -1更早的评论隐藏 -1更早的评论

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2019-7-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/472619-how-to-put-together-different-cases-in-one-matlab-code#answer_384057

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

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.