면진 구조물 해석에 관하여

#해당 구조물 해석을 진행할 시 계산을 간소화시키기 위하여 레일라이 감쇠를 고려한 비고전적 감쇠행렬을 고전적 감쇠행렬로 치환하여 가정하고 진행하였음을 미리 고지한다. 해당 가정은 지반을 강체로 고려할 수 있는(지상구조물의 주기가 충분히 큰) 경우에 수치적으로 유효한 가정이다. 만약 내진설계를 구조물의 강성을 높여서 진행하는 경우, 구조물이 매우 작은 고유주기를 가지게 된다면 지반-구조물 상호작용을 고려한 레일라이 감쇠 행렬을 그대로 적용하여야만 한다.

 

 

 

해당 내용은 Chopra의 'Dynamics of Structures Theory and Applications to Earthquake Engineering' 4판의 21장 'Earthquake Dynamics of Base-Isolated Buildings' 내용에 대해 공부하면서 과거에 정리 해 둔 것을 옮기는 내용이다.

 

 누구에게 설명하는 용도가 아닌, 내 개인 노트에 정리된 내용을 옮기는 과정이므로 관련 전문 지식이 없으면 이해하기 어렵다. 설명에서 이해하기 어려운 부분이나 누락 된 점 있을 수 있는것을 양해바란다. 내용 하단에 있는 코드도 책 내용의 계산 검토를 위해 작성한 코드라 책이 없으면 무엇을 의미하는지 알기가 힘들다. (해당 내용을 공부한지 시간이 많이 지나서... 최근 어떤 구조물의 해석 과정에서 발생한 이해가 안되는 부분 때문에 출퇴근길에 관련 내용을 되짚어 보려고 올린다.)

 

 

우선 주기 T=0.4에 해당하는 각속도 wn 값은 15.7080 [단위길이/sec] 에 해당한다. 또한 각 층 슬래브 중심마다 층의 질량 100kips/g을 집중시킨 질량 행렬은

이 그림과 같다. 또한 각 층의 층강성을 k라 하면 건물의 강성행렬은

여기에 k 값을 곱한 것과 같다.

이제 wn^2 * [m] = [k] 일때 건물의 고유치가 존재한다는 것을 이용하여, m,k 행렬에 대한 고유치 해석을 matlab 함수 eig를 사용하여 진행하거나 행렬식(determinant)이 0이 되게 만드는 k 값을 구할 수 있다. 이렇게 구해진 k값은 3.0456e+5 이 나온다. 이때 주의할 점은 wn값이 n차의 고유주기가 될 수 있는 만큼 k값도 그에 대응하여 n개의 값이 나오므로 wn이 1차 모드의 각속도라는 것을 wn 행렬에서 확인할 필요가 있다. 여기서 k값이 여러개라는건 위의 식을 만족하는 wn 배열을 달리하면 그렇게 보이는 것이지 실제 k행렬은 구조물이 결정될때 하나로 결정되므로 k행렬 자체가 여러개라는 뜻이 아니다. (이것을 위해서 eig 함수를 사용하는 것이 더 직관적이다.)

 

 건축물의 모드해석시에 고유벡터는 각 모드의 형상에 대응하며 [더 정확히는 여기서 구한 고유벡터는 modal matrix이며 질량 정규화한 [pi]'*[m]*[pi]가 단위행렬이 되는 [pi]가 mode shape 이 된다.] 고유값은 각 모드 형상과 연관되는 고유진동수에 대응된다.

 그리고 건물의 해석시 일반적으로 행렬 [A] 하나로 바로 구해지는 것이 아니라 질량 행렬[m], 강성 행렬[k]을 통해 k/m=wn^2 공식을 사용하여 해석을 진행하므로 앞에서 나온 표준식 (1)이 아닌 [k]*[pi]=wn^2 * [m]*[pi] 로 정리된 일반식을 사용하여 고유치 해석을 진행하게 된다.

 

-------------------------------------------------------[인용 자료]--------------------------------------------------------------

***참고로 고유치 해석은 고유값과 고유벡터를 구하는 것인데 

이 부분은 다른 블로거가 정리한 자료를 바탕으로 해설을 남긴다. 해당 내용 인용을 허락해 준 다크님에게 감사드린다.

해당 블로그 주소를 인용자료 하단에 같이 남기니 해당 내용이 궁금하신 분들이 있다면 찾아가서 읽어보길 추천 드린다. 고유치해석 이외에도 상당히 유용한 선형대수학 설명자료와 프로그래밍 내용들이 많이 포스팅 된 곳이다.

 

고유값(eigenvalue), 고유벡터(eigenvector)에 대한 수학적 정의는 비교적 간단하다.

 

행렬 A를 선형변환으로 봤을 때, 선형변환 A에 의한 변환 결과가 자기 자신의 상수배가 되는 0이 아닌 벡터를 고유벡터(eigenvector)라 하고 이 상수배 값을 고유값(eigenvalue)라 한다.

 

즉, n x n 정방행렬(고유값, 고유벡터는 정방행렬에 대해서만 정의된다) A에 대해 Av = λv를 만족하는 0이 아닌 열벡터 v를 고유벡터, 상수 λ를 고유값이라 정의한다.

 

 ---(1)

 

 --- (2)

 

좀더 정확한 용어로는 λ는 '행렬 A의 고유값', v는 '행렬 A의 λ에 대한 고유벡터'이다.

 

즉, 고유값과 고유벡터는 행렬에 따라 정의되는 값으로서 어떤 행렬은 이러한 고유값-고유벡터가 아에 존재하지 않을수도 있고 어떤 행렬은 하나만 존재하거나 또는 최대 n개까지 존재할 수 있다.

 

기하학적으로 봤을 때, 행렬(선형변환) A의 고유벡터는 선형변환 A에 의해 방향은 보존되고 스케일(scale)만 변화되는 방향 벡터를 나타내고 고유값은 그 고유벡터의 변화되는 스케일 정도를 나타내는 값이다.

 

 

예를 들어 지구의 자전운동과 같이 3차원 회전변환을 생각했을 때, 이 회전변환에 의해 변하지 않는 고유벡터는 회전축 벡터이고 그 고유값은 1이 될 것이다.

 

 

출처: https://darkpgmr.tistory.com/105

[선형대수학 #3] 고유값과 고유벡터 (eigenvalue & eigenvector)

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

 

 

 

가정된 감쇠계수의 행렬은(쵸플라 책 예제에서) k 행렬에 비례한 [c]=a1*[k] 이므로 레일리 감쇠에서 a0=0인 형태와 같으며, 이때 감쇠비의 공식은 (감쇠비)= a1*wn/2 가 된다. 따라서 a1=2*0.02/15.7080=0.0025가 된다.

면진장치의 감쇠정수와 강성값은 주어진 자료의 식 21.4.1a 와 21.4.1b를 참고하여 구하면, kb=5.9218e+03, cb=376.9911 값이 나온다.

이제 앞서 구한 k, a1, kb, cb 값을 이용하여 기존 구조물과 면진장치가 설치된 구조물의 강성행렬 및 감쇠행렬을 구하면 다음과 같다.

%% problem1
clear; clc; close all;
zeta=0.02;                            
T=0.4;
 
wn1=2*pi/T;            
m=diag(ones(1,5))*100; 
kshear=zeros(5);             
for i=1:6               
   if i~=6
       kshear(i,i)=2;
      if i>1
        kshear(i,i-1)=-1;
        kshear(i-1,i)=-1;      
      end
   else
       kshearb=kshear;          
       kshear(5,5)=1;            
       kshearb(1,1)=1; kshearb(6,6)=1;   
       kshearb(i,i-1)=-1;
       kshearb(i-1,i)=-1;
       
   end    
end
 
k0=wn1^2 * 10000;                      
k=fzero(@kk,k0,[],kshear,m,wn1)        
                                       
 
 
 
a1=2*zeta/wn1                         
Tb=2; zetab=0.1;                       
wb=2*pi/Tb;
mb=zeros(6,6); mb(6,6)=100;                  
mb(1:5, 1:5)=m;                                    
kb=wb^2 * sum(mb,'all')       
cb=zetab*  2* sum(mb,'all') * wb
 
kshearb=k*kshearb;
kshearb(1,1)=kshearb(1,1)+kb;
 
syms wn0;
eqn = kk(k,kshear,m,wn0)==0;                 
Wn=double(solve(eqn,wn0)); j=0;              
for i=1:length(Wn)                          
    if Wn(i)>0
        j=j+1;
        wn0(j)=Wn(i); 
    end
    if i==length(Wn)
        wn0=double(wn0);
        wn0=sort(wn0);
        for i=1:length(wn0)
            wn(i,i)=wn0(i);             
        end
    end
    
end
 
syms wn0;
eqn = kk(1,kshearb,mb,wn0)==0;               
Wn= double(solve(eqn,wn0)); j=0;            
for i=1:length(Wn)                          
    if Wn(i)>0
        j=j+1;
        wn0(j)=Wn(i); 
    end
    if i==length(Wn)
        wn0=double(wn0);
        wn0=sort(wn0);
        for i=1:length(wn0)
            wnb(i,i)=wn0(i);
        end
    end
  
end
 
c_matrix=a1*k*kshear;                             
cb_matrix=[zeros(1,5);c_matrix]; cb_matrix=[zeros(6,1), cb_matrix];
cb_matrix(1,1)=cb+c_matrix(5,5); cb_matrix(1,2)=c_matrix(1,2); cb_matrix(2,1)=c_matrix(2,1); 
 
 
 
Tn=diag(2*pi./diag(wn));                                 
Tnb=diag(2*pi./diag(wnb));                               
 
zetas= diag(wn)*a1/2                                     
 
[Phib,Omega2b]=eig(kshearb,mb);                         
Mb=Phib'*mb*Phib;
Cb=Phib'*cb_matrix*Phib;                                 
Kb=Phib'*kshearb*Phib;
zetasb= 0.5*diag(Cb)./ (diag(Mb).*diag(Kb)).^0.5   
 
tablea=table(round(diag(Tn),4),round(zetas,4))
tableb=table(round(diag(Tnb),4),round(zetasb,4))
 
 
%% problem 2 _1번 실행이후 실행
clc;  close all;          %clear안함 problem1의 변수 가져옴
 
[Phi,Omega2]=eig(k*kshear,m);
[Phib,Omega2b]=eig(kshearb,mb);
 
for i=1:6
    text(i)= sprintf("%d shape",i);
end    
Phi=[text(1:5);round(Phi,4)]
Phib=[text(1:6);round(Phib,4)]
 
%% problem 3 _1,2번 실행이후 실행
clc;  close all;          %clear안함 problem1,2의 변수 가져옴
[Phi,Omega2]=eig(k*kshear,m);
[Phib,Omega2b]=eig(kshearb,mb);
%modal participation factor/ effective modal mass ratio
for i=1:5
    PF(i)= (diag(m)'*Phi(:,i))/ (diag(m)'*(Phi(:,i).*Phi(:,i)));        
    Mn_(i)= (diag(m)'*Phi(:,i))^2/ (diag(m)'*(Phi(:,i).*Phi(:,i)));      
    Mratio(i)=Mn_(i)/sum(m,'all');                                      
end
for i=1:6
    PFb(i)= (diag(mb)'*Phib(:,i))/ (diag(mb)'*(Phib(:,i).*Phib(:,i)));
    Mn_b(i)= (diag(mb)'*Phib(:,i))^2/ (diag(mb)'*(Phib(:,i).*Phib(:,i)));
    Mratiob(i)=Mn_b(i)/sum(mb,'all');
end
 
 
for i=1:5
    sn(:,i)= Phi(:,i)*PF(i)  ;               %모드 참여율*모드세잎
    snm(:,i)=sn(:,i).*diag(m);
 end
for i=1:6
    snb(:,i)= Phib(:,i)*PFb(i);
    snmb(:,i)=snb(:,i).*diag(mb);
end
 
                       
                                  
snm                                
snmb                               
 
%% Problem 4 문제 1~3까지 실행후 실행
clc; close all;         %clear안함 위에 변수 가져옴
% Input
g=386; 
for i=1:6   
    if i~=6
        Sa1(i)=Design_Spectrum_KDS_damping('S3',0.176,Tn(i,i),zetas(i)) ;       
        Sa2(i)=Design_Spectrum_KDS_damping('S3',0.176,Tnb(i,i),zetasb(i)) ;    
    else
        Sa2(i)=Design_Spectrum_KDS_damping('S3',0.176,Tnb(i,i),zetasb(i))  ;   
    end   
end
 
Vbn= Mn_.*Sa1*g ;
Vbnb= Mn_b.*Sa2*g ;
Vst_bn= Mn_;
Vst_bnb= Mn_b;
SRSS= sqrt(sum((Vbn/500/g).*(Vbn/500/g),'all'));
SRSSb= sqrt(sum((Vbnb/600/g).*(Vbnb/600/g),'all'));
 
fprintf("An/g 기본 건물: %.4f, %.4f, %.4f, %.4f, %.4f\n",Sa1)
fprintf("An/g 면진 건물: %.4f, %.4f, %.4f, %.4f, %.4f, %.4f\n",Sa2)
 
fprintf("V^{st}_{bn}/m 기본 건물: %.4f, %.4f, %.4f, %.4f, %.4f\n",Vst_bn/100)
fprintf("V^{st}_{bn}/m 면진 건물: %.4f, %.4f, %.4f, %.4f, %.4f, %.4f \n \n",Vst_bnb/100)
fprintf("Mn*의 총합과 m의 총합이 같아야 하는걸 생각하면.\n") 
fprintf("Vbn/W 기본 건물: %.4f, %.4f, %.4f, %.4f, %.4f\n",Vbn/500/g)
fprintf("Vbn/W 면진 건물: %.4f, %.4f, %.4f, %.4f, %.4f, %.4f\n",Vbnb/600/g)
 
fprintf("SRSS 기본 건물 : %.4f \n",SRSS)
fprintf("SRSS 면진건물: %.4f \n",SRSSb)
 
%% problem 5 (1~4번 실행 이후) Cb 행렬:non-classical damping) [앞에서 푼 문제와 같음], diag(diag(Cb))행렬:calssical 사용시 응답비교 
 
 
 
% non classical
Matrix= readmatrix('Elcentro.xlsx');      
p=Matrix(:,2)';
tt=Matrix(:,1)'; 
 
 
 
u0=[0 0 0 0 0 0]; v0=[0 0 0 0 0 0];
l=[1 1 1 1 1 1]'; 
 
A=[zeros(6,6) eye(6,6); -inv(mb)*kshearb -inv(mb)*cb_matrix];
B=[zeros(6,1) ; inv(mb)*l];                                                              
                                                                     
C=[-inv(mb)*kshearb -inv(mb)*cb_matrix];
D=[ inv(mb)*l]; 
ync_lsim=lsim(A,B,C,D,p,tt,[u0 v0]');
 
 
 
cb_matrix2=diag(diag(cb_matrix));
A=[zeros(6,6) eye(6,6); -inv(mb)*kshearb -inv(mb)*cb_matrix2];
B=[zeros(6,1) ; inv(mb)*l];                                                              
 
C=[-inv(mb)*kshearb -inv(mb)*cb_matrix2];
D=[ inv(mb)*l]; 
 
yc_lsim=lsim(A,B,C,D,p,tt,[u0 v0]');
 
subplot(6,1,1)
pp1=plot(tt,ync_lsim(:,1));
hold on
pp2=plot(tt,yc_lsim(:,1));
hold on
legend('non classical damping','classical damping');
xlabel('sec')
ylabel('가속도[단위?]')
title('1번째 자유도에 대한 그래프')
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
subplot(6,1,2)
pp1=plot(tt,ync_lsim(:,2));
hold on
pp2=plot(tt,yc_lsim(:,2));
hold on
legend('non classical damping','classical damping');
xlabel('sec')
ylabel('가속도[단위?]')
title('2번째 자유도에 대한 그래프')
 
subplot(6,1,3)
pp1=plot(tt,ync_lsim(:,3));
hold on
pp2=plot(tt,yc_lsim(:,3));
hold on
legend('non classical damping','classical damping');
xlabel('sec')
ylabel('가속도[단위?]')
title('3번째 자유도에 대한 그래프')
 
subplot(6,1,4)
pp1=plot(tt,ync_lsim(:,4));
hold on
pp2=plot(tt,yc_lsim(:,4));
hold on
legend('non classical damping','classical damping');
xlabel('sec')
ylabel('가속도[단위?]')
title('4번째 자유도에 대한 그래프')
 
subplot(6,1,5)
pp1=plot(tt,ync_lsim(:,5));
hold on
pp2=plot(tt,yc_lsim(:,5));
hold on
legend('non classical damping','classical damping');
xlabel('sec')
ylabel('가속도[단위?]')
title('5번째 자유도에 대한 그래프')
 
subplot(6,1,6)
pp1=plot(tt,ync_lsim(:,6));
hold on
pp2=plot(tt,yc_lsim(:,6));
hold on
legend('non classical damping','classical damping');
xlabel('sec')
ylabel('가속도[단위?]')
title('6번째 자유도에 대한 그래프')
 
 
 
%% 사용 함수
 
%행렬식 함수
function k= kk(k,kshear,m,wn1)
k=det(k*kshear-(wn1^2)*m); 
end
 
function Sa=Design_Spectrum_KDS_damping(site,S,T,z_eff)
 
% 지반증폭계수
switch site
    case {'S1'}
        Fas=[1.12 1.12 1.12];
        Fvs=[0.84 0.84 0.84];
    case {'S2'}
        Fas=[1.4 1.4 1.3];
        Fvs=[1.5 1.4 1.3];
    case {'S3'}
        Fas=[1.7 1.5 1.3];
        Fvs=[1.7 1.6 1.5];
    case {'S4'}
        Fas=[1.6 1.4 1.2];
        Fvs=[2.2 2.0 1.8];
    case {'S5'}
        Fas=[1.8 1.3 1.3];
        Fvs=[3.0 2.7 2.4];
end
 
if S<=0.1
    Fa=Fas(1);
    Fv=Fvs(1);
elseif S<=0.2
    Fa=interp1([0.1 0.2],[Fas(1) Fas(2)],S);
    Fv=interp1([0.1 0.2],[Fvs(1) Fvs(2)],S);
elseif S<=0.3
    Fa=interp1([0.2 0.3],[Fas(2) Fas(3)],S);
    Fv=interp1([0.2 0.3],[Fvs(2) Fvs(3)],S);
else
    Fa=Fas(3);
    Fv=Fvs(3);
end
 
S_DS=2/3*S*2.5*Fa;   % 단주기 설계응답스펙트럼 가속도(g)
S_D1=2/3*S*Fv;       % 1초주기 설계응답스펙트럼 가속도(g)
 
% Ts=S_D1/S_DS; % 가속도 일정구간과 속도지배구간의 전이주기
% To=0.2*Ts; % 가속도 일정구간의 시작점
% Tl=5;
% 
% for ii=1:length(T)
%     if T(ii)<=To
%         Sa(ii)=0.6*S_DS*T(ii)/To+0.4*S_DS;
%     elseif T(ii)<=Ts
%         Sa(ii)=S_DS;
%     elseif T(ii)<=Tl
%         Sa(ii)=S_D1/T(ii);
%     else
%         Sa(ii)=S_D1*Tl/T(ii)^2;
%     end
% end
 
 
% 감쇠비 고려
z_eff=z_eff*100;    % 감쇠비(%)
if z_eff <=2
    Bs=0.8;B1=0.8;
elseif z_eff <=5
    Bs=interp1([2 5],[0.8 1.0],z_eff);
    B1=interp1([2 5],[0.8 1.0],z_eff);
elseif z_eff <=10
    Bs=interp1([5 10],[1.0 1.3],z_eff);
    B1=interp1([5 10],[1.0 1.2],z_eff);
elseif z_eff <=20
    Bs=interp1([10 20],[1.3 1.8],z_eff);
    B1=interp1([10 20],[1.2 1.5],z_eff);
elseif z_eff <=30
    Bs=interp1([20 30],[1.8 2.3],z_eff);
    B1=interp1([20 30],[1.5 1.7],z_eff);
elseif z_eff <=40
    Bs=interp1([30 40],[2.3 2.7],z_eff);
    B1=interp1([30 40],[1.7 1.9],z_eff);
elseif z_eff <=50
    Bs=interp1([40 50],[2.7 3.0],z_eff);
    B1=interp1([40 50],[1.9 2.0],z_eff);
else
    Bs=3.0;B1=2.0;
end
 
Ts=S_D1*Bs/(S_DS*B1);
To=0.2*Ts;
Tl=5;
 
for ii=1:length(T)
    if T(ii)<=To
        Sa(ii)=S_DS*((5/Bs-2)*T(ii)/Ts+0.4);
    elseif T(ii)<=Ts
        Sa(ii)=S_DS/Bs;
    elseif T(ii)<=Tl
        Sa(ii)=S_D1/(B1*T(ii));
    else
        Sa(ii)=S_D1*Tl/(B1*T(ii)^2);
    end
end
end

<코드 내용>

 

k =
 
   3.0456e+05
 
 
a1 =
 
    0.0025
 
 
kb =
 
   5.9218e+03
 
 
cb =
 
  376.9911
 
 
zetas =
 
    0.0200
    0.0584
    0.0920
    0.1182
    0.1348
 
 
zetasb =
 
    0.0958
    0.0564
    0.0787
    0.1034
    0.1233
    0.1361
 
 
tablea =
 
  5×2 table
 
     Var1      Var2 
    ______    ______
 
       0.4      0.02
     0.137    0.0584
    0.0869     0.092
    0.0677    0.1182
    0.0593    0.1348
 
 
tableb =
 
  6×2 table
 
     Var1      Var2 
    ______    ______
 
    2.0298    0.0958
    0.2175    0.0564
    0.1136    0.0787
    0.0804    0.1034
    0.0657    0.1233
    0.0589    0.1361

<출력 내용1>

 

Phi = 
 
  6×5 string array
 
    "1 shape"    "2 shape"    "3 shape"    "4 shape"    "5 shape"
    "0.017"      "-0.0456"    "0.0597"     "0.0549"     "-0.0326"
    "0.0326"     "-0.0597"    "0.017"      "-0.0456"    "0.0549" 
    "0.0456"     "-0.0326"    "-0.0549"    "-0.017"     "-0.0597"
    "0.0549"     "0.017"      "-0.0326"    "0.0597"     "0.0456" 
    "0.0597"     "0.0549"     "0.0456"     "-0.0326"    "-0.017" 
 
 
Phib = 
 
  7×6 string array
 
    "1 shape"    "2 shape"    "3 shape"    "4 shape"    "5 shape"    "6 shape"
    "-0.0396"    "0.0559"     "-0.0504"    "0.0412"     "0.0291"     "-0.0151"
    "-0.0403"    "0.0417"     "-0.0007"    "-0.0405"    "-0.0577"    "0.0409" 
    "-0.0408"    "0.016"      "0.0496"     "-0.041"     "0.0287"     "-0.0558"
    "-0.0412"    "-0.014"     "0.0501"     "0.0406"     "0.029"      "0.0557" 
    "-0.0414"    "-0.0402"    "0.0002"     "0.0409"     "-0.0577"    "-0.0408"
    "-0.0416"    "-0.0554"    "-0.0499"    "-0.0408"    "0.0288"     "0.0149" 

snm =
 
   35.6271   30.0884   20.7694   10.6288    2.8863
   68.3680   39.4074    5.9116   -8.8307   -4.8562
   95.5701   21.5243  -19.0868   -3.2920    5.2843
  115.0296  -11.2165  -11.3442   11.5658   -4.0347
  125.1702  -36.2148   15.8578   -6.3173    1.5041
 
 
snmb =
 
   97.0608    2.2173    0.4907    0.1644    0.0549    0.0118
   98.6426    1.6529    0.0071   -0.1617   -0.1088   -0.0321
   99.9141    0.6355   -0.4836   -0.1639    0.0541    0.0438
  100.8712   -0.5560   -0.4884    0.1623    0.0546   -0.0437
  101.5110   -1.5952   -0.0024    0.1633   -0.1088    0.0320
  101.8314   -2.1972    0.4860   -0.1628    0.0544   -0.0117