How to put together different cases in one MATLAB code?

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

[Guillaume]

"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]

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.

[Dragana Skoko]

Exactly what I want :).

Thanks :))))))

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

[Dragana Skoko]

Actually , there is something wrong with this code.It works perfect when I run it for the first time.But next time when I try to run code in same window I get this message "Subscripted assignment dimension mismatch."I have to close MATLAB and run it again to get result...Do you have any advice how to solve it ?

[Rik]

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]

@Rik — Thank you!

I do not have any problems running it several times in succession.

What line throws the error?