MATLAB Examples

# EXAMPLE 1: SINGLE WIND HISTORY

## Buffeting analysis: hypothesis and assumptions:

• The quasi-steady theory is applied.
• No modal coupling induced by the wind load is introduced.
• No structural modal coupling is introduced.
• Both the along and vertical wind component are taken into account for the buffeting load.
• The wind direction is assumed normal to the bridge deck
• The incidence angle of the wind is equal to 0 deg.
• The frequency domain approach is used.

## Mechanical properties of the bridge

clearvars;close all;clc; % phi and wn are matrix that I have built and store the mode shapes and % eigen frequencies of the main span. load('bridgeModalProperties.mat','wn','phi') % Lateral mode shapes: phiH = squeeze(phi(1,:,:)); % Vertical mode shapes: phiV = squeeze(phi(2,:,:)); % Torsional mode shapes: phiT = squeeze(phi(3,:,:)); % Lateral eigen frequencies: wnH = squeeze(wn(1,:)); % Vertical mode shapes: wnV = squeeze(wn(2,:)); % Torsional mode shapes: wnT = squeeze(wn(3,:)); [Nmodes,Nyy]=size(phiH); % Nyy: number of "nodes" % Structural modal damping ratio for the Bridge deck: % here, for simplicity, it is taken as 0.5 % for every modes Bridge.zetaStruct = 5e-3.*ones(1,Nmodes); % The mechanical properties of the single span suspension bridge are % computed. Bridge.B = 12.3 ; % width of the bridge deck (m) Bridge.D = 2.76 ; % height of the bridge deck (m) Bridge.L = 446 ; % Length of the main span (m) Bridge.x = linspace(0,1,Nyy) ; % normalized bridge main span Bridge.m = 5350 ; % lineic mass of girder (kg/m) Bridge.mc = 408 ; % lineic mass of cable (kg/m) Bridge.m_theta = 82430; % torsional stifness (kg.m^2/m) Bridge.k = 1/4; % bridge constant . cf aerodynamic of streamlined bridge % aerodynamic coefficient (quasi steady terms) Bridge.Cd = 1;% drag coefficient Bridge.dCd = 0;% first derivative of drag coefficient Bridge.Cl = 0.1;% lift coefficient Bridge.dCl = 3;% first derivative of lift coefficient Bridge.Cm = 0.02;% pitching moment coefficient Bridge.dCm = 1.12;% first derivative of pitching moment coefficient 

## Time and frequency definition

tmax = 600; % based on 600 seconds of records = 10 min fs = 10; % sampling frequency of 30 Hz is chosen fmin = 1/tmax; % minimal frequency recorded Nfreq = 512; f = logspace(log10(1/tmax),log10(fs/2),Nfreq); 

## Semi empirical wind spectrum (Von karman)

clear Wind Wind.U = 10; % mean wind speed normal to the bridge deck Wind.Cuy = 7; % decay coefficient for co-coherence for u-component along bridge span Wind.Cwy = 6; % decay coefficient for co-coherence for w-component along bridge span Wind.f = f; Lu = 100; Lw = Lu/10; stdU = 0.15*Wind.U; stdW = 0.55*stdU; [Su] = VonKarmanSpectrum(f,Wind.U,stdU,Lu,'u'); [Sw] = VonKarmanSpectrum(f,Wind.U,stdW,Lw,'w'); % Storage of the data into the structure variableWind . Wind.Su = Su; Wind.Sw = Sw; % Let's check the spectra: figure semilogx(f,f.*Su./stdU^2,f,f.*Sw./stdW^2) xlabel('f (Hz)') ylabel('Normalized PSD') legend('u','w') % Selection of position where the response is computed % We choose to compute the response at 1/3 of the bridge. % We recall that x is the normalized span length ranging from 0 to 1 [~,position] = min(abs(Bridge.x-1/3)); 

## Overview of the fundamental inputs

fprintf(' The varable Wind is : \n\n') disp(Wind) fprintf(' The varable Bridge is : \n\n') disp(Bridge) 
 The varable Wind is : U: 10 Cuy: 7 Cwy: 6 f: [1x512 double] Su: [1x512 double] Sw: [1x512 double] The varable Bridge is : zetaStruct: [0.0050 0.0050 0.0050 0.0050] B: 12.3000 D: 2.7600 L: 446 x: [1x30 double] m: 5350 mc: 408 m_theta: 82430 k: 0.2500 Cd: 1 dCd: 0 Cl: 0.1000 dCl: 3 Cm: 0.0200 dCm: 1.1200 

## bridge response

Bridge.DOF = 'lateral'; Bridge.wn = wnH; % eigen frequencies (rad/s) Bridge.phi = phiH; % mode shapes [Srx] = dynaRespFD3(Bridge,Wind,position); % compute PSD of response 

## Vertical bridge response

Bridge.DOF = 'vertical'; Bridge.wn = wnV; Bridge.phi = phiV; [Srz] = dynaRespFD3(Bridge,Wind,position); 

## Torsional bridge response

Bridge.DOF = 'torsional'; Bridge.wn = wnT; Bridge.phi = phiT; [Srt] = dynaRespFD3(Bridge,Wind,position); 

## Data plot

figure subplot(221) loglog(f,Srx); xlabel(' frequency (Hz)'); ylabel('S_{rx} (m^2/s)'); axis tight grid on subplot(222) loglog(f,Srz); xlabel(' frequency (Hz)'); ylabel('S_{rz} (m^2/s)'); axis tight grid on subplot(223) loglog(f,Srt); xlabel(' frequency (Hz)'); ylabel('S_{rt} (rad^2/Hz)'); axis tight grid on 

## Multiple coherence models

Bridge.DOF = 'vertical'; Bridge.wn = wnV; Bridge.phi = phiV; Wind.Cuy = 7 ; % large decay coeff Wind.Cwy = 7; % large decay coeff [Srz1] = dynaRespFD3(Bridge,Wind,position); [Srz2] = dynaRespFD3(Bridge,Wind,position,'cohType','ESDU'); figure plot(f,Srz1,'k',f,Srz2,'r--'); legend('empirical coherence model','ESDU coherence model','location','SouthWest') set(gca,'Xscale','log'); set(gca,'Yscale','log'); xlabel(' frequency (Hz)'); ylabel('S_{rz} (m^2/s)'); axis tight 

## Influence of mean wind velocity

myU = [5,10,20,30,40]; Srz_test = zeros(numel(myU),Nfreq); for ii=1:numel(myU) Wind.U = myU(ii); [Srz_test(ii,:)] = dynaRespFD3(Bridge,Wind,position); end Wind.U = 10; figure plot(f,Srz_test); legend('U = 5 m/s','U = 10 m/s','U = 20 m/s','U = 30 m/s','U = 40 m/s',... 'location','SouthWest') set(gca,'Xscale','log'); set(gca,'Yscale','log'); xlabel(' frequency (Hz)'); ylabel('S_{rz} (m^2/s)'); axis tight 

## Influence of modal damping ratios

Bridge.DOF = 'vertical'; Bridge.wn = wn(2,:); Bridge.phi = squeeze(phi(2,:,:)); Wind.Cuy = 7 ; Wind.Cwy = 7; myZeta = [1,2,5,7,10].*1e-3; Srz_test = zeros(numel(myZeta),Nfreq); for ii=1:numel(myZeta) Bridge.zetaStruct = myZeta(ii).*ones(1,Nmodes); [Srz_test(ii,:)] = dynaRespFD3(Bridge,Wind,position); end figure plot(f,Srz_test); legend('\zeta = 1e-3','\zeta = 2e-3','\zeta = 5e-3','\zeta = 7e-3',... '\zeta = 1e-2','location','SouthWest') set(gca,'Xscale','log'); set(gca,'Yscale','log'); xlabel(' frequency (Hz)'); ylabel('S_{rz} (m^2/s)'); xlim([0.1,0.7]) ylim([1e-6,1e-2]) 

## Influence of wind turbulence

Bridge.zetaStruct = 5e-3.*ones(1,Nmodes); myIu = [0.05,0.1,0.2,0.3]; Srz_test = zeros(numel(myIu),Nfreq); for ii=1:numel(myIu) stdU = myIu(ii)*Wind.U; stdW = 0.55*stdU; Su = VonKarmanSpectrum(f,Wind.U,stdU,Lu,'u'); Sw = VonKarmanSpectrum(f,Wind.U,stdW,Lw,'w'); Wind.Su = Su; Wind.Sw = Sw; [Srz_test(ii,:)] = dynaRespFD3(Bridge,Wind,position); end figure plot(f,Srz_test); legend('I_u = 0.05','I_u = 0.1','I_u = 0.2','I_u = 0.3',... 'location','SouthWest') set(gca,'Xscale','log'); set(gca,'Yscale','log'); xlabel(' frequency (Hz)'); ylabel('S_{rz} (m^2/s)'); axis tight Wind.Iu = 0.1; % TI for u-component Wind.Iw = Wind.Iu*0.6; % TI for w-component 

## Influence of turbulence length scales

Lu = [50,100,150,200]; Srz = zeros(numel(Lu),Nfreq); for ii=1:numel(Lu) Su = VonKarmanSpectrum(f,Wind.U,stdU,Lu(ii),'u'); Sw = VonKarmanSpectrum(f,Wind.U,stdW,0.1*Lu(ii),'w'); % Lw = 0.1*Lu Wind.Su = Su; Wind.Sw = Sw; [Srz(ii,:)] = dynaRespFD3(Bridge,Wind,position); end figure plot(f,Srz); legend('Lw = 5 m','Lwx = 10 m','Lwx = 15 m','Lwx = 20 m',... 'location','SouthWest') set(gca,'Xscale','log'); set(gca,'Yscale','log'); xlabel(' frequency (Hz)'); ylabel('S_{rz} (m^2/s)'); axis tight