Add files via upload

eigenBridge.m

@@ -0,0 +1,219 @@
+function [wn,phi,phi_cables] = eigenBridge(Bridge,Nmodes)
+% function [wn,phi,phi_cables] = eigenBridge(Bridge,Nyy,Nmodes) computes
+% the mode shape and eigen-frequency of a single-span suspension bridge.
+% It is designed to be fast and straitghforward, although it might slightly
+% less accurate than a Finite element model.
+%
+% INPUT:
+% Bridge: type: structure (see file studyCase.m)
+% Nyy: type: float [1 x 1] : number of "nodes" consituting the deck
+% Nmodes: type: float [1 x 1] : number of "modes" to be computed
+%
+% OUTPUT
+% wn: type: 2D matrix [3 x Nmodes] : eigen-modes of the suspension bridge
+% wn(1,:) --> all the eigen-frequencies for the lateral displacement (x axis)
+% wn(2,:) --> all the eigen-frequencies for the vertical displacement (z axis)
+% wn(3,:) --> all the eigen-frequencies for the torsional angle (around y axis)
+%
+% phi: type: 3D matrix [3 x Nmodes x Nyy] : modes shapes of the bridge girder ( =deck)
+% phi(1,i,:) --> mode shape for the i-th eigen frequency for the lateral
+% bridge displacement ( x axis)
+% phi(2,i,:) --> mode shape for the i-th eigen frequency for the vertical
+% bridge displacement ( z axis)
+% phi(3,i,:) --> mode shape for the i-th eigen frequency for the torsional
+% bridge angle (around y axis)
+%
+% phi_cables: type: 2D matrix [Nmodes x Nyy]:
+% modes shapes of the bridge main cables along x axis. No mode
+% shapes for vertical and torsional motion of the cable are calculated.
+%
+% Author info:
+% E. Cheynet , University of Stavanger. last modified: 31/12/2016
+%
+% see also LysefjordBridge hardangerBridge
+%
+
+%%
+% preallocation
+if isfield(Bridge,'Nyy'),
+ Nyy = Bridge.Nyy;
+else
+ error(' ''Nyy'' is not a field of the structure ''Bridge'' ')
+end
+
+wn = zeros(3,Nmodes); % eigen frequency matrix
+phi = zeros(3,Nmodes,Nyy); % mode shapes of bridge girder
+phi_cables = zeros(Nmodes,Nyy); % mode shapes of cables
+
+% discretisation of bridge span
+x = linspace(0,Bridge.L,Nyy); % vector bridge span
+x = x./Bridge.L ;% reduced length
+
+
+%% LATERAL MOTION
+% preallocation
+alpha = zeros(2*Nmodes,2*Nmodes); % reduced variable
+beta = zeros(2*Nmodes,2*Nmodes); % reduced variable
+gamma = zeros(2*Nmodes,2*Nmodes); % reduced variable
+
+m_tilt= Bridge.m;
+mc_tilt = 2.*Bridge.mc;
+
+% calculation of alpha, beta and gamma
+% alpha and beta
+% cf E.N. Strømmen "STRUCTURAL DYNAMICS" for explanations
+for nn=1:2*Nmodes,
+ alpha(nn,nn) = (Bridge.E.*Bridge.Iz).*(nn*pi./Bridge.L).^4;
+ beta(nn,nn) = 2*Bridge.H_cable.*(nn*pi./Bridge.L).^2;
+end
+
+% gamma
+% cf E.N. Strømmen "STRUCTURAL DYNAMICS" for explanations
+for pp=1:2*Nmodes,
+ for nn=1:2*Nmodes,
+ if and(rem(pp,2)==1,rem(nn,2)==0)||and(rem(pp,2)==0,rem(nn,2)==1)
+ gamma(pp,nn)=0;
+ else
+ gamma(pp,nn)= 2*Bridge.m*Bridge.g/(Bridge.L*Bridge.ec).*trapz(x.*Bridge.L,sin(pp.*pi.*x).*sin(nn.*pi.*x)./...
+ (1+Bridge.hm/Bridge.ec-4.*(x).*(1-x)));
+ end
+ end
+end
+
+% matrix mass
+M = diag(repmat([m_tilt;mc_tilt],[Nmodes,1]));
+
+% Matrix stiwness K
+for p=1:2:2*Nmodes,
+ for nn=1:2:2*Nmodes,
+ if p==nn,
+ clear Omega
+ Omega(1,1) = alpha((p+1)/2,(nn+1)/2)+gamma((p+1)/2,(nn+1)/2);
+ Omega(1,2) = -gamma((p+1)/2,(nn+1)/2);
+ Omega(2,1) = -gamma((p+1)/2,(nn+1)/2);
+ Omega(2,2) = beta((p+1)/2,(nn+1)/2)+gamma((p+1)/2,(nn+1)/2);
+ K(p:p+1,nn:nn+1) = Omega;
+ else
+ clear V
+ V = gamma((p+1)/2,(nn+1)/2).*[1,-1;-1,1];
+ K(p:p+1,nn:nn+1) = V;
+ end
+ end
+end
+
+% eigen-value problem solved for non-trivial solutions
+[vector,lambda]=eig(K,M,'chol');
+
+wn(1,:) = sqrt(diag(lambda(1:Nmodes,1:Nmodes))); % filling the matrix wn
+
+% Normalization
+for ii=1:Nmodes,
+ % deck mode shape construction using series expansion
+ phi(1,ii,:) = vector(1:2:end,ii)'*sin([1:1:Nmodes]'.*pi*x);
+ % cables mode shape construction using series expansion
+ phi_cables(ii,:) = vector(2:2:end,ii)'*sin([1:1:Nmodes]'.*pi*x); 
+
+ modeMax = max([max(abs(phi_cables(ii,:))),max(abs(phi(1,ii,:)))]);
+ phi(1,ii,:) = phi(1,ii,:)./modeMax; % normalisation
+ phi_cables(ii,:) = phi_cables(ii,:)./modeMax; % normalisation
+end
+
+
+%% VERTICAL MOTION
+
+clear K M
+% INITIALISATION
+kappa = zeros(Nmodes,Nmodes); %reduced variable
+lambda = zeros(Nmodes,Nmodes); %reduced variable
+mu = zeros(Nmodes,Nmodes); %reduced variable
+
+le = Bridge.L*(1+8*(Bridge.ec/Bridge.L)^2); % effective length
+% cf page 124 "STRUCTURAL DYNAMICS" of E.N. Strømmen
+
+for nn=1:Nmodes,
+ kappa(nn,nn) = Bridge.E*Bridge.Iy.*(nn*pi./Bridge.L).^4;
+ lambda(nn,nn) = 2*Bridge.H_cable.*(nn*pi./Bridge.L).^2;
+end
+
+for p=1:Nmodes,
+ for nn=1:Nmodes,
+ if and(rem(p,2)==1,rem(nn,2)==1) % sont impaires
+ mu(p,nn)=(32*Bridge.ec/(pi*Bridge.L)).^2*(Bridge.Ec*Bridge.Ac)/(Bridge.L*le)/(p*nn);
+ else
+ mu(p,nn)= 0;
+ end
+ end
+end
+
+M = (2*Bridge.mc+Bridge.m).*eye(Nmodes); % mass matrix
+K= kappa+lambda +mu; % stiwness matrix
+
+
+% eigen-value problem solved for non-trivial solutions
+[vector,lambda]=eig(K,M,'chol');
+
+wn(2,:) = sqrt(diag(lambda(1:Nmodes,1:Nmodes))); % filling of matrix wn
+
+for ii=1:Nmodes,
+ phi(2,ii,:) = vector(1:Nmodes,ii)'*sin([1:Nmodes]'.*pi*x); % mode shape construction using series expansion
+ phi(2,ii,:) = phi(2,ii,:)./max(abs(phi(2,ii,:))); % normalisation
+end
+
+%% TORSIONAL MOTION
+
+clear K M
+omega = zeros(Nmodes,Nmodes);
+v = zeros(Nmodes,Nmodes);
+V = zeros(Nmodes,Nmodes);
+xi = zeros(Nmodes,Nmodes);
+m_tilt_theta_tot = zeros(Nmodes,Nmodes);
+m_tilt = 2*Bridge.mc+Bridge.m;
+
+
+m_theta_tot = Bridge.m_theta + Bridge.mc*(Bridge.bc^2/2);
+
+for nn=1:Nmodes,
+ omega(nn,nn) = (nn*pi./Bridge.L).^2.*(Bridge.GIt+(nn*pi/Bridge.L)^2*(Bridge.E*Bridge.Iw));
+ V(nn,nn) = Bridge.H_cable * (Bridge.bc^2/2)*(nn*pi./Bridge.L).^2;
+ v(nn,nn) = (m_tilt)*Bridge.g*Bridge.hr;
+ m_tilt_theta_tot(nn,nn) = m_theta_tot;
+end
+
+for p=1:Nmodes,
+ for nn=1:Nmodes,
+ if and(rem(p,2)==1,rem(nn,2)==1) % sont impaires
+ xi(p,nn)=(16*Bridge.ec*Bridge.bc/(pi*Bridge.L)).^2*Bridge.Ec*Bridge.Ac/(Bridge.L*le)*1/(p*nn);
+ else
+ xi(p,nn)= 0;
+ end
+ end
+end
+
+
+K= omega+V+v+xi; % stiwness matrix
+M = diag(diag(m_tilt_theta_tot)); % mass matrix
+
+% eigen-value problem solved for non-trivial solutions
+[vector,lambda]=eig(K,M,'chol');
+wn(3,:) = sqrt(diag(lambda(1:Nmodes,1:Nmodes))); % filling the wn matrix
+
+for ii=1:Nmodes,
+ phi(3,ii,:) = vector(1:Nmodes,ii)'*sin([1:Nmodes]'.*pi*x); % mode shape construction using series expansion
+ phi(3,ii,:) = phi(3,ii,:)./max(abs(phi(3,ii,:))); % normalisation
+end
+
+
+% positiv max value
+for ii=1:3,
+ for jj=1:Nmodes,
+ if abs(min(squeeze(phi(ii,jj,:))))> abs(max(squeeze(phi(ii,jj,:)))),
+ phi(ii,jj,:)=-phi(ii,jj,:);
+ end
+ end
+end
+
+
+
+
+end
+