MATLAB Examples

Contents

Introduction

This example is only a verification of the function eigenBridge2.m based on the works of [1]. For this reason, I use a slightly different function whose only purpose is to verifiy the numerical implementation. The function I use here is : eigenBridge_Verification.m

First we will start with the computations of the eigen-frequencies in a table in a similar fashion as table 1 and 2 from [1].

Then we will plot the evolution of the eigen-frequencies as a function of mu and lambda (fig. 2 from [1])

Finally we will plot the mode shapes for different values of lambda and mu (fig. 3 and 44 from [1])

[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.

clear all;close all;clc;

Part 1: eigen frequencies calculation

Approximatively 2 minutes of computation is required for this part.
Mode-shapes and eigen-frequencies are calculated for every combination
of _mu_ and _lambda_ . Separation between anti-symmetric and symmetric mode
shapes and eigen frequencies is done.
% Following Table 1 in [1], we have:
lambda = sqrt([0,100,225,1e8]); % 1e8 instead of  + infinity
mu = sqrt([0,1e-5,1e-4,5e-4,1e-3,2.5e-3,5e-3,1e-2,2.5e-2,5e-2,0.1]);
Nmodes = 4; % I only want 4 modes here
Nyy=200; % Number of integration points for deck span
% Non-dimensional deck span x/L, where L = span length
x = linspace(0,1,Nyy);


% Preallocation
w_s = zeros(numel(lambda),numel(mu),Nmodes);
w_a = zeros(numel(lambda),numel(mu),Nmodes);
phi_s = zeros(numel(lambda),numel(mu),Nmodes,Nyy);
phi_a = zeros(numel(lambda),numel(mu),Nmodes,Nyy);

% Calculation of eigen frequencies and mode shapes
tic
for ii=1:numel(lambda),
    for jj=1:numel(mu),
        [w_a(ii,jj,:),w_s(ii,jj,:),phi_a(ii,jj,:,:),phi_s(ii,jj,:,:)] = ...
            eigBridge_Verification(mu(jj),lambda(ii),Nmodes,x);
    end
end
toc

% We display the eigen frequencies in the workspace. A comparison with
% Table 1 in [1] shows a good agreement, which suggests that the
% implementation of the code is fine.

testVerif_s = w_s./pi;
testVerif_a = w_a./pi;
fprintf('1st mode \n')
disp([mu(:).^2,testVerif_s(1,:,1)',testVerif_s(2,:,1)',...
    testVerif_s(3,:,1)',testVerif_s(4,:,1)',testVerif_a(4,:,1)'])
fprintf('2nd mode \n')
disp([mu(:).^2,testVerif_s(1,:,2)',testVerif_s(2,:,2)',...
    testVerif_s(3,:,2)',testVerif_s(4,:,2)',testVerif_a(4,:,2)'])
fprintf('3rd mode \n')
disp([mu(:).^2,testVerif_s(1,:,3)',testVerif_s(2,:,3)',...
    testVerif_s(3,:,3)',testVerif_s(4,:,3)',testVerif_a(4,:,3)'])
fprintf('4th mode \n')
disp([mu(:).^2,testVerif_s(1,:,4)',testVerif_s(2,:,4)',...
    testVerif_s(3,:,4)',testVerif_s(4,:,4)',testVerif_a(4,:,4)'])
Elapsed time is 119.530244 seconds.
1st mode 
         0    1.0000    2.5971    2.7878    2.8607    2.0000
    0.0000    1.0000    2.5980    2.7887    2.8619    2.0004
    0.0001    1.0005    2.6041    2.7989    2.8734    2.0039
    0.0005    1.0025    2.6308    2.8435    2.9230    2.0196
    0.0010    1.0049    2.6608    2.8973    2.9842    2.0391
    0.0025    1.0123    2.7352    3.0488    3.1605    2.0964
    0.0050    1.0244    2.8205    3.2713    3.4336    2.1885
    0.0100    1.0482    2.9125    3.6198    3.9232    2.3620
    0.0250    1.1166    3.0138    4.1180    5.1178    2.8192
    0.0500    1.2221    3.0841    4.3303    6.6479    3.4490
    0.1000    1.4096    3.1783    4.4617    8.9553    4.4487

2nd mode 
         0    3.0000    3.4798    4.4621    4.9182    4.0000
    0.0000    3.0013    3.4804    4.4646    4.9243    4.0032
    0.0001    3.0133    3.4880    4.4853    4.9784    4.0315
    0.0005    3.0659    3.5202    4.5588    5.2123    4.1549
    0.0010    3.1304    3.5600    4.6215    5.4905    4.3043
    0.0025    3.3164    3.6800    4.7285    6.2516    4.7240
    0.0050    3.6052    3.8907    4.8335    7.3466    5.3510
    0.0100    4.1224    4.3230    5.0248    9.1514    6.4239
    0.0250    5.3839    5.5004    5.8044   13.1564    8.8975
    0.0500    6.9980    7.0760    7.2240   17.9441   11.9303
    0.1000    9.4310    9.4840    9.5668   24.8957   16.3909

3rd mode 
         0    5.0000    5.0525    5.3024    6.9417    6.0000
    0.0000    5.0062    5.0586    5.3069    6.9586    6.0106
    0.0001    5.0613    5.1127    5.3473    7.1075    6.1057
    0.0005    5.2995    5.3460    5.5306    7.7352    6.5112
    0.0010    5.5829    5.6252    5.7678    8.4546    6.9851
    0.0025    6.3578    6.3913    6.4754   10.3155    8.2448
    0.0050    7.4728    7.4991    7.5513   12.8308    9.9978
    0.0100    9.3105    9.3303    9.3624   16.7651   12.8027
    0.0250   13.3870   13.3999   13.4184   25.1073   18.8620
    0.0500   18.2599   18.2691   18.2815   34.8218   25.9914
    0.1000   25.3348   25.3413   25.3499   48.7536   36.2643

4th mode 
         0    7.0000    7.0149    7.0487    8.9547    8.0000
    0.0000    7.0169    7.0318    7.0652    8.9907    8.0252
    0.0001    7.1673    7.1817    7.2132    9.3058    8.2488
    0.0005    7.8005    7.8132    7.8384   10.5940    9.1768
    0.0010    8.5262    8.5377    8.5578   12.0114   10.2189
    0.0025   10.4040   10.4129   10.4266   15.5049   12.8478
    0.0050   12.9416   12.9485   12.9584   20.0153   16.3135
    0.0100   16.9106   16.9156   16.9226   26.8520   21.6393
    0.0250   25.3263   25.3298   25.3340   41.0155   32.7818
    0.0500   35.1261   35.1287   35.1315   57.3095   45.6650
    0.1000   49.1802   49.1821   49.1840   80.5515   64.0827

Part 2: Evolution of eigen frequencies with lambda and mu

Here, I am reproducing Figure 2 from [1], where the x-axis is logarithmic. The figures I obtain are very close to what [1] get. Hurra !

% 3 figures are plotted
IndTitle = {'First','Second','Third'};
indLim = [6,12,20];
clf;close all;
for ii=1:3,
    figure
    hold on; box on;
    plot(mu.^2,w_s(:,:,ii)./pi);
    plot(mu.^2,w_a(:,:,ii)./pi,'k--');
    set(gca,'xscale','log');
    legend('\lambda^2 =0','\lambda^2 =100','\lambda^2 =225','\lambda^2 =1e6','antisymm.','location','NorthWest')
    xlabel('\mu^2')
    ylabel('\omega/\pi')
    title([IndTitle{ii},' mode shape'])
    ylim([0,indLim(ii)]);
    set(gcf,'color','w');
end

Part 3: Evolution of mode shapes with lambda and mu

Here, I am reproducing Figure 3 from [1], where the mode-shapes for different values og mu and lambda are calculated. It seems that the mode shapes I obtain agree well with [1] Hurra !

clf;close all;
figure('units','normalized','outerposition',[0,0,0.7,1]);
jj=1;
IndTitle = {'First','Second','First','Second',};
legLambda = [100,100,225,225];
indLambda = [2,2,3,3];
indMode = [1,2,1,2];
for ii=1:4,
    subplot(2,2,ii)
    hold on; box on;
    plot(x,squeeze(phi_s(indLambda(ii),2:2:end,indMode(ii),:)));
    if ii==1,
        legend('\mu^2 =1e-5','\mu^2 =5e-4','\mu^2 =2.5e-3','\mu^2 =1e-2','\mu^2 =5e-2','location','South')
    end
    xlabel('x/L')
    ylabel('\phi')
    title([IndTitle{ii},' mode shape with \lambda^2 =',num2str(legLambda(ii))])
    set(gcf,'color','w');
end

Conclusions:

The numerical implementation seems fine. However, numerical issues may occur if the options "tolFun" and "tolX" in the fsolve function are higher than 1e-8. The user is welcome to try 1e-6 or higher and see the differences.