MATLAB Examples

EXAMPLE 1: SINGLE WIND HISTORY

Contents

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