Add files via upload

eigBridge_Verification.m

@@ -0,0 +1,173 @@
+function [w_a,w_s,phi_a,phi_s] = eigBridge_Verification(mu,lambda,Nmodes,x)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% GOAL: Compute the non-dimensional eigen frequencies and mode shapes
+% following [1], as a verification procedure
+
+% [1] Luco, J. E., & Turmo, J. (2010).
+% Linear vertical vibrations of suspension bridges:
+% A review of continuum models and some new results.
+% Soil Dynamics and Earthquake Engineering,
+% 30(9), 769-781.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%  INPUTS:
+
+% mu: positive scalar: relative bending stifness of the girder
+% lambda: positive scalar: Irvine-Caughy vable parameter
+% Nmodes: positive natural: Number of modes to be computed (< 6 )
+% x: vector [1 x Nyy] Non-dimensional bridge span discretized in Nyy points
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% OUTPUTS
+
+
+% w_s: non dimensional symmetric eigen frequencies
+% w_a: non dimensional asymmetric eigen frequencies
+% phi_s: normalized symmetric mode shapes
+% phi_a: normalized asymmetric mode shapes
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Author: Etienne Cheynet  -- last modified: 26/12/2015
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% INITIALISATION
+% Number of integration points
+Nyy = numel(x);
+% non-dimensionla pulsation
+CST = (2*pi*[1:1:Nmodes]);
+
+%% Antisymmetric eigen frequencies and mode shapes
+% egien-frequencies
+w_a = CST.*sqrt(1+(CST.*mu).^2);
+% mode shapes
+phi_a = sin(2*pi.*[1:1:Nmodes]'*x);
+
+%% Symmetric eigen frequencies and mode shapes
+
+% Here there are 2 steps:
+% 1) Solve a characteristic equation to find w
+% 2) Report omega into the analytical expression of the mode shape
+
+% 1) Solve the characteristic equation
+% Nmber of solutions Nsol: be careful not to choose too many modes !
+Nsol = max(200,Nmodes^3);
+dummy = zeros(1,Nsol);
+% check matlab version
+checkMatlabVersion = version;
+
+% First, we need too check if lambda is not too close to zero,
+% i.e. below 0.1. Otherwise it is assumed to be an inelastic cable.
+if lambda >0.1,
+    if str2double(checkMatlabVersion(end-5:end-2))<=2012, % R2012b or earlier
+        for ii=1:Nsol,
+            dummy(ii) =fsolve(@(w) charFun(w,mu,lambda),ii,...
+                optimset('Display','off','TolFun', 1e-8, 'TolX', 1e-8));
+        end
+    elseif str2double(checkMatlabVersion(end-5:end-2))>=2013, % R2013a or later
+        for ii=1:Nsol,
+            dummy(ii) =fsolve(@(w) charFun(w,mu,lambda),ii,...
+                optimoptions('fsolve','Display','off',...
+                'TolFun', 1e-8, 'TolX', 1e-8));
+        end
+    end
+
+    % Analytically, many solutions are identical to each other,
+    % but numerically, it is not the case.
+    % Therefore I need to limit the prceision of the solutions to 1e-3.
+    dummy = unique(round(dummy*1e3).*1e-3);
+    % remove the first solution that is found to be meaningless
+    dummy(dummy<1e-1)=[];
+    % Check if any  issues:
+    if isempty(dummy), % no solutions found
+        warning('no mode shapes found');
+        w_a =NaN(1,Nmodes);
+        w_s =NaN(1,Nmodes);
+        phi_a=NaN(Nmodes,Nyy);
+        phi_s=NaN(Nmodes,Nyy);
+        return
+    elseif numel(dummy)<Nmodes, % less solutions than Nmodes found
+        warning('not enough mode shapes found');
+        w_a =NaN(1,Nmodes);
+        w_s =NaN(1,Nmodes);
+        phi_a=NaN(Nmodes,Nyy);
+        phi_s=NaN(Nmodes,Nyy);
+        return
+    end
+    % Limit the number of solutions to Nmodes
+    w_s=dummy(1:Nmodes);
+    % Get the mode shapes
+    phi_s = zeros(Nmodes,Nyy);
+    for ii=1:Nmodes
+        [phi_s(ii,:)] = modeShapeFun(x,mu,lambda,w_s(ii));
+    end
+else % case of an inelastic cable, i.e. lambda = 0
+    N = 2.*[1:1:Nmodes]-1;
+    dummy = N.*pi.*sqrt(1+(N.*pi.*mu).^2);
+    w_s=dummy(1:Nmodes); % eigen-frequencies
+    phi_s = sin(N'*x.*pi); % mode shapes
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Nested functions
+    function H = charFun(w,mu,lambda)
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % GOAL: choose the corresponding characteristic equation
+        % INPUT
+        %   w: Non dimensional eigen frequency (unknown)
+        %   mu: positive scalar: relative bending stifness of the girder
+        %   lambda: positive scalar: Irvine-Caughy vable parameter
+        % OUTPUT: H: Characteristic equation to be solved for w
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        if mu<=1e-6, % case of a perfectly flexible deck (limit case)
+            H = 1-(w/lambda).^2-tan(w/2)/(w/2);
+        elseif mu>=1e6, % very stiff deck ( limit case)
+            A = sqrt((sqrt(1+(2*mu*w).^2)-1)/(2*mu^2));
+            H = 1-(A*mu/lambda).^2-(1/2)*tan(A/2)/(A/2)-...
+                (1/2)*tanh(A/2)/(A/2);
+        else % case of deck usually met
+            A = sqrt((sqrt(1+(2*mu*w).^2)-1)/(2*mu^2));
+            B = sqrt((sqrt(1+(2*mu*w).^2)+1)/(2*mu^2));
+            H = 1-(w/lambda).^2-(B^2/(A^2+B^2))*tan(A/2)/(A/2)-...
+                (A^2/(A^2+B^2))*tanh(B/2)/(B/2);
+        end
+    end
+
+    function [phi] = modeShapeFun(x,mu,lambda,omega)
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % GOAL: calculate the mode shapes
+        % INPUT
+        %   x:  [1 x Nyy] vector of non-dimensional deck span
+        %   mu: positive scalar: relative bending stifness of the girder
+        %   lambda: positive scalar: Irvine-Caughy vable parameter
+        %   omega: Non dimensional eigen frequency (it is now known)
+        % OUTPUT: phi: [1 x Nyy] vector of mode shape
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        if mu<=1e-6, % case of a perfectly flexible deck (limit case)
+            dummyPhi = (lambda.*omega).^2.*(1-cos(omega*(x-1/2))/cos(omega/2));
+            phi = dummyPhi./max(abs(dummyPhi));% normalization
+        elseif mu>=1e6, % very stiff deck ( limit case)
+            A = sqrt((sqrt(1+(2*mu*omega).^2)-1)/(2*mu^2));
+            dummyPhi = (1-0.5*cos(A*(x-1/2))/cos(A/2)-...
+                0.5*cosh(A*(x-1/2))/cosh(A/2));
+            phi = dummyPhi./max(abs(dummyPhi));% normalization
+        else % case of deck usually met
+            A = sqrt((sqrt(1+(2*mu*omega).^2)-1)/(2*mu^2));
+            B = sqrt((sqrt(1+(2*mu*omega).^2)+1)/(2*mu^2));
+            dummyPhi = (lambda.*omega).^2.*(1-(B.^2/(A.^2+B.^2))*cos(A*(x-1/2))/cos(A/2)-...
+                (A.^2/(A.^2+B.^2))*cosh(B*(x-1/2))/cosh(B/2));
+            phi = dummyPhi./max(abs(dummyPhi)); % normalization
+        end
+        if abs(min(phi))>max(phi),
+            phi=-phi; % positiv maximum
+        end
+    end
+
+
+end
+

eigenBridge.m

@@ -0,0 +1,223 @@
+function [wn,phi,phi_cables] = eigenBridge(Bridge,Nmodes)
+% function [wn,phi,phi_cables] = eigenBridge(Bridge) 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
+
+if isfield(Bridge,'ec') && isfield(Bridge,'sag')
+    Bridge.sag = Bridge.ec;
+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.sag).*trapz(x.*Bridge.L,sin(pp.*pi.*x).*sin(nn.*pi.*x)./...
+                (1+Bridge.hm/Bridge.sag-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.sag/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.sag/(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.sag*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
+