MATLAB Examples

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