Algoritmo Matlab analisi modale2D

NOTE: è necessario sostituire il comando "assemblaMatrici_M_K_lumped" con

"assemblaMatrici_M_K_1piano"

ovvero l'altro algoritmo che potete scaricare sempre da me, che consente di

assemblare le

matrici di rigidezza e massa per un portale 1 piano. Quando caricate questo

algoritmo in matlab

dovrete salvarlo con il nome di "assemblaMatrici_M_K_1piano", cosi chè questo

algoritmo può

richiamarlo e fare l'analisi modale piana.

--------------------------------------------------------------------------------

------------------------------------------------

Algoritmo:

close all

clear all

clc

%% ANALISI MODALE

% rendo il sistema di ED delle equazioni di equilibrio, tutto disaccoppiato

assemblaMatrici_M_K_lumped %cambiare nel caso in cui si usa Mconsistent

%equazione del moto Ma+Ku=-M l ag(t)

%vettore di spostamento pseudostatico, il sisma non agisce sulle rotazioni

quindi zero

l=[1 1 1];

%fattori di eccitazione modale

for yy=1:length(freq)

eccit_mod(yy).L=modo(yy).n'*M_tt*l';

end

%Masse Modali

for yy=1:length(freq)

masse_mod(yy).M=modo(yy).n'* M_tt *modo(yy).n;

end

%fattori di partecipazione modale (quantizza il grado di partecipazione modo n

al totale)

for yy=1:length(freq)

gamma(yy).g=eccit_mod(yy).L/ masse_mod(yy).M;

end

%espansione modale del vettore di eccitazione --> s, peff(t)=s*p(t)=M l ag(t)

s=M_tt*l';

for yy=1:length(freq)

esp_s(yy).sn=gamma(yy).g*M_tt * modo(yy).n;

end

%% RISPOSTE STATICHE MODALI

for yy=1:length(freq)

mod=modo(yy).n; %ricavo l'autovettore corrispondente

Lh=0; %inizializzazione

Lteta=0;

for k=1:length(ut)

Lh= Lh + M_tt(k,k)*mod(k);

Lteta=Lteta + k*L*M_tt(k,k)*mod(k);

end Lh_n(yy)=Lh;

Lteta_n(yy)=Lteta;

end

for yy=1:length(freq)

h_star(yy)= Lteta_n(yy)/ Lh_n(yy); %altezza modale effettiva

M_star(yy)=gamma(yy).g*Lh_n(yy); %massa modale effettiva

%taglio alla base statico

Vb_st(yy)=M_star(yy); %Kg

%taglio ad ogni piano

V_piano_st(yy).piano(1)=Vb_st(yy);

for s=1:length(ut)-1

V_piano_st(yy).piano(s+1)=V_piano_st(yy).piano(s) - esp_s(yy).sn(s);

%kg end

%momento alla base statico

Mb_st(yy)=M_star(yy)*h_star(yy); %Kgm

%spostamenti degli impalcati

u_piani_st(yy).piano=(gamma(yy).g/(freq(yy)^2))*modo(yy).n;

%drift

for h=1:(length(ut)-1)

drift_st(yy).piani(h)=(gamma(yy).g/(freq(yy)^2))*(modo(yy).n(h+1)-

modo(yy).n(h));

end

end

%% RISPOSTE DINAMICHE - SISTEMA LINEARE - METODO NUMERICO DELL'INTERPOLAZIONE

DELLA FORZANTE

% ricavo le risposte di ogni SDOF

%accelerogramma

load Laquila_NS_AQA.txt

Laquila_NS_AQA(8001:length(Laquila_NS_AQA))=[];

%elcentro(:,1)=[]; %tolgo la colonna del tempo, qst

istruzione la posso togliere

acc=Laquila_NS_AQA/100; %perchè leggo dal file txt, prima riga che

l'unità di misura è cm/s^2, però g è in m/s^2

n=length(acc);

dt=0.005; %passo d'integrazione=passo di misura

t=0:dt:dt*(n-1);

g=9.807; %controllare se l'accelerogramma fornito è in

m/s2 in caso cambiare

figure(100)

plot(t,acc/g)

title('Accelerogramma - L Aquila, comp.N-S, staz.AQA');

xlabel('t (s)')

ylabel('a(t), g')

grid

% . . . . .

% D(t) + 2zitawnD(t)+wn^2 D(t)= - ug(t)=peff(t)

for yy=1:length(freq)

Tn(yy)=2*pi/freq(yy); %periodo naturale corrispondente ad ogni modo

zita=0.05;

wn=freq(yy); %frequenza di ogni modo

wd=wn*sqrt(1-zita*zita);

%M=m;

k=(wn^2); %poichè è fissata la frequenza trovo la rigidezza

peff=-acc; %sarebbe la p(t)

% Coefficienti della formula

%coefficienti dello spostamento

A=(exp(-zita*wn*dt))*((zita/sqrt(1-zita*zita))*sin (wd*dt) +

cos(wd*dt)); B=(exp(-zita*wn*dt))*((1/wd)*sin (wd*dt));

C=(1/k)*((2*zita/(wn*dt))+(exp(-zita*wn*dt))*((((1-2*zita*zita)/

(wd*dt))-(zita/sqrt(1-zita*zita)))*sin(wd*dt) - (1+((2*zita)/(wn*dt)))*cos

(wd*dt))); D=(1/k)*(1-(2*zita/(wn*dt))+(exp(-zita*wn*dt))*(((2*zita*zita - 1)/

(wd*dt))*sin (wd*dt)+(2*zita/(wn*dt))*cos(wd*dt)));

%coefficienti della velocità

A_=-(exp(-zita*wn*dt))*((wn/(sqrt(1-zita*zita)))*sin(wd*dt));

B_=(exp(-zita*wn*dt))*(cos(wd*dt)-(zita/(sqrt(1-

zita*zita)))*sin(wd*dt));

C_=(1/k)*(-1/dt +(exp(-zita*wn*dt))*((wn/(sqrt(1-zita*zita)) + zita/

(dt*sqrt(1-zita*zita)))*sin(wd*dt) +(1/dt)*cos(wd*dt)));

D_=(1/(k*dt))*(1-(exp(-zita*wn*dt))*((zita/(sqrt(1-

zita*zita)))*sin(wd*dt)+cos(wd*dt)));

%Risposte-uso la variabile struct cosi memorizzo la risposta A per ogni modo

spost(yy).D(1)=0;

velocita(yy).v(1)=0;

for i=1:n-1

spost(yy).D(i+1)=A*spost(yy).D(i)+B*velocita(yy).v(i)+C*peff(i)

+ D*peff(i+1); % Dn(t) che proviene dall'equazione modale

velocita(yy).v(i+1)=A_*spost(yy).D(i)+B_*velocita(yy).v(i)

+C_*peff(i) + D_*peff(i+1);

Pseud_Acc(yy).A(i+1)=(wn^2)*spost(yy).D(i+1);

%An(t) end

end

%rappresentazione delle risposte SDOF di ogni modo

for yy=1:length(freq)

figure(yy+10)

plot(t, spost(yy).D)

title(['Spostamento SDOF modo: ',num2str(yy)]);

xlabel('t (s)')

ylabel('D(t), m')

grid

figure (yy+16)

plot(t,Pseud_Acc(yy).A)

title(['PseudoAccelerazione spettrale SDOF modo: ',num2str(yy)]);

xlabel('t (s)')

ylabel('A(t), m/s^2')

grid

end

%% RISPOSTE MODALI rn(t)=rn_st*An(t)

for yy=1:length(freq)

%taglio fra i piani

for ii=1:length(ut) %ogni piano deve essere moltiplicato per la

risposta dinamica

V_piano(yy).piano(ii).piano=V_piano_st(yy).piano(ii)*Pseud_Acc(yy).A;

end

%taglio alla base

V_base(yy).Taglio=Vb_st(yy)*Pseud_Acc(yy).A;