Luận án Chẩn đoán dầm cầu bằng phương pháp phân tích dao động trên mô hình số hoá kết cấu được cập nhật sử dụng thuật toán tối ưu hoá bầy đàn kết hợp mạng nơ ron nhân tạo

BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC GIAO THÔNG VẬN TẢI 
---------------------- 
HỒ KHẮC HẠNH 
CHẨN ĐOÁN DẦM CẦU BẰNG PHƯƠNG PHÁP 
PHÂN TÍCH DAO ĐỘNG TRÊN MÔ HÌNH SỐ HOÁ 
KẾT CẤU ĐƯỢC CẬP NHẬT SỬ DỤNG 
THUẬT TOÁN TỐI ƯU HOÁ BẦY ĐÀN KẾT HỢP 
MẠNG NƠ RON NHÂN TẠO 
LUẬN ÁN TIẾN SĨ KỸ THUẬT 
HÀ NỘI – 2021 
BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC GIAO THÔNG VẬN TẢI 
---------------------- 
HỒ KHẮC HẠNH 
CHẨN ĐOÁN DẦM CẦU BẰNG PHƯƠNG PHÁP 
PHÂN TÍCH DAO ĐỘNG TRÊN MÔ HÌNH SỐ HOÁ 
KẾT CẤU ĐƯỢC CẬP NHẬT SỬ DỤNG 
THUẬT TOÁN TỐI ƯU HOÁ BẦY ĐÀN KẾT HỢP 
MẠNG NƠ RON NHÂN TẠO 
Chuyên ngành : Kỹ thuật Xây dựng Cầu - Hầm 
Mã số : 9580205 
LUẬN ÁN TIẾN SĨ KỸ THUẬT 
NGƯỜI HƯỚNG DẪN KHOA HỌC: 
1. PGS.TS. Bùi Tiến Thành 
2. PGS.TS. Ngô Văn Minh 
HÀ NỘI – 2021 
LỜI CAM ĐOAN 
Tôi xin cam đoan đây là công trình nghiên cứu khoa học độc lập của riêng tôi. Các số 
liệu sử dụng phân tích trong luận án có nguồn gốc rõ ràng, đã công bố theo đúng quy định. 
Các kết quả nghiên cứu trong luận án do tôi tự tìm hiểu, phân tích một cách trung thực. Việc 
tham khảo các nguồn tài liệu đã được thực hiện trích dẫn và ghi nguồn tài liệu tham khảo 
đúng quy định. 
 Hà Nội, ngày 30 tháng 10 năm 2021 
Tác giả 
Hồ Khắc Hạnh 
MỤC LỤC 
MỞ ĐẦU ......................................................................................................................... 1 
1. Đặt vấn đề nghiên cứu .............................................................................................. 1 
2. Mục tiêu nghiên cứu ................................................................................................. 4 
3. Phương pháp nghiên cứu .......................................................................................... 5 
4. Đối tượng và phạm vi nghiên cứu ............................................................................ 5 
5. Ý nghĩa khoa học và thực tiễn .................................................................................. 5 
6. Nội dung của luận án ................................................................................................ 6 
CHƯƠNG 1: TỔNG QUAN VỀ GIÁM SÁT SỨC KHỎE KẾT CẤU CÔNG 
TRÌNH DỰA TRÊN PHƯƠNG PHÁP ĐO NHẬN DẠNG DAO ĐỘNG................ 8 
1.1. Tổng quan về giám sát sức khỏe kết cấu công trình dựa vào kết quả đo các đặc 
trưng dao động ................................................................................................................. 8 
1.2. Tình hình nghiên cứu trên thế giới về giám sát sức khỏe kết cấu công trình dựa 
trên phương pháp đo nhận dạng dao động .................................................................... 11 
1.3. Tình hình nghiên cứu ở Việt Nam về giám sát sức khỏe kết cấu công trình dựa trên 
phương pháp đo nhận dạng dao động ............................................................................ 16 
Kết luận Chương 1 ......................................................................................................... 18 
CHƯƠNG 2: NGHIÊN CỨU VÀ ĐỀ XUẤT PHƯƠNG PHÁP CHẨN ĐOÁN HƯ 
HỎNG KẾT CẤU BẰNG PHƯƠNG PHÁP ĐO DAO ĐỘNG NGẪU NHIÊN ... 19 
2.1. Cơ sở lý thuyết về dao động kết cấu ....................................................................... 19 
2.2. Các phương pháp chẩn đoán hư hỏng kết cấu dựa trên dao động .......................... 23 
2.2.1. Phương pháp dựa trên sự thay đổi của tần số dao động riêng .................. 23 
2.2.2. Phương pháp dựa vào sự thay đổi của hình dạng dao động ...................... 25 
2.2.3. Phương pháp dựa trên hàm đáp ứng tần số ............................................... 34 
2.3. Chẩn đoán hư hỏng kết cấu dựa trên kết quả đo nhận dạng dao động sử dụng thuật 
toán bầy đàn PSO .......................................................................................................... 37 
2.4. Chẩn đoán hư hỏng kết cấu dựa trên kết quả đo nhận dạng dao động sử dụng mạng 
nơ ron nhân tạo. ............................................................................................................. 40 
2.5. Chẩn đoán hư hỏng kết cấu dựa trên kết quả đo nhận dạng dao động sử dụng thuật 
toán kết hợp PSO-ANN. ................................................................................................ 44 
Kết luận Chương 2 ......................................................................................................... 52 
CHƯƠNG 3: ÁP DỤNG PSO KẾT HỢP VỚI ANN ĐỂ CHẨN ĐOÁN CÁC HƯ 
HỎNG CHO MÔ HÌNH SỐ ....................................................................................... 53 
3.1. Cầu dầm giản đơn ................................................................................................... 53 
3.1.1. Trường hợp hư hỏng tại một phần tử ........................................................ 56 
3.1.2. Trường hợp hư hỏng tại nhiều phần tử ..................................................... 59 
3.2. Tấm composite ........................................................................................................ 61 
3.2.1. Trường hợp hư hỏng tại một phần tử ........................................................ 64 
3.2.2. Trường hợp hư hỏng tại nhiều phần tử ..................................................... 66 
3.3. Cầu dàn thép liên tục ............................................................................................... 68 
3.3.1. Trường hợp hư hỏng tại một phần tử ........................................................ 73 
3.3.2. Trường hợp hư hỏng tại nhiều phần tử ..................................................... 75 
Kết luận Chương 3 ......................................................................................................... 77 
CHƯƠNG 4: ÁP DỤNG PSO KẾT HỢP VỚI ANN ĐỂ CHẨN ĐOÁN CÁC HƯ 
HỎNG TRÊN MÔ HÌNH THỰC NGHIỆM ............................................................ 79 
4.1. Dầm giản đơn trong phòng thí nghiệm ................................................................... 79 
4.1.1. Mô hình bài toán ....................................................................................... 79 
4.1.2. Mô hình thực nghiệm ................................................................................ 80 
4.1.3. Phân tích, so sánh các kết quả thu được từ lý thuyết và thực nghiệm ...... 83 
4.2. Cầu Bến Quan ......................................................................................................... 88 
4.2.1. Giới thiệu về cầu Bến Quan ...................................................................... 88 
4.2.2. Mô hình thí nghiệm ................................................................................... 89 
4.2.3. Mô hình phần tử hữu hạn .......................................................................... 91 
4.2.4. Cập nhật mô hình và xác định hư hỏng trong kết cấu ............................... 92 
Kết luận Chương 4 ......................................................................................................... 96 
KẾT LUẬN VÀ KIẾN NGHỊ ..................................................................................... 97 
DANH MỤC HÌNH VẼ 
Hình 1.1. Cầu Tsing Ma được lắp đặt 600 đầu đo phục vụ cho việc giám sát sức khỏe 
kết cấu công trình 10 
Hình 2.1. Sơ đồ thuật toán tối ưu hoá bầy đàn – PSO ............................................................. 39 
Hình 2.2. Kiến trúc mạng ANN ................................................................................................ 40 
Hình 2.3. Sơ đồ thuật toán ANN .............................................................................................. 43 
Hình 2.4. Quá trình xác định giải pháp tối ưu dựa trên phương pháp GD của ANN (a) mạng 
với một điểm tối ưu toàn cục, (b) mạng với tối ưu toàn cục và cục bộ .................................... 45 
Hình 2.5. Mạng với nhiều tối ưu cục bộ ................................................................................... 46 
Hình 2.6. Áp dụng PSO để khắc phục vấn đề tối ưu cục bộ của ANN .................................... 47 
Hình 2.7. Sơ đồ thuật toán PSO kết hợp song song với ANN .................................................. 51 
Hình 3.1. Bố trí chung cầu ............................................................................................. 53 
Hình 3.2. Mặt cắt ngang cầu .......................................................................................... 54 
Hình 3.3. Kiến trúc mạng nơ ron cho bài toán xác định hư hỏng trong dầm giản đơn . 55 
Hình 3.4. Giá trị hồi quy (a) ANN và (b) ANN kết hợp với PSO ................................. 57 
Hình 3.5. Sai số phân bố (a) ANN và (b) ANN kết hợp với PSO ................................. 57 
Hình 3.6. Sai số thực (a) ANN và (b) ANN kết hợp với PSO ...................................... 58 
Hình 3.7. Kết quả phát hiện hư hỏng của PSO, ANN và ANNPSO của phần tử số 3: (a) 
7% hư hỏng; (b) 70% hư hỏng ...................................................................................... 59 
Hình 3.8. Giá trị hồi quy (a) ANN và (b) ANN kết hợp với PSO ................................ 60 
Hình 3.9. Sai số phân bố (a) ANN và (b) ANN kết hợp với PSO ................................ 60 
Hình 3.10. Sai số thực (a) ANN và (b) ANN kết hợp với PSO .................................... 60 
Hình 3.11. Kết quả phát hiện hư hỏng của PSO, ANN và ANNPSO trong dầm: 20% 
hư hỏng của phần tử 1 – và 60% hư hỏng của phần tử 3. ............................................. 61 
Hình 3.12. (a) Tấm composite 3 lớp ; (b) Mô hình phần tử hữu hạn ............................ 62 
Hình 3.13. Hình dạng dao động của 15 mode đầu tiên ................................................. 63 
Hình 3.14. Kiến trúc mạng dùng để phát hiện hư hỏng của các phần tử trong tấm 
composite ....................................................................................................................... 64 
Hình 3.15. Giá trị hồi quy (a) ANN và (b) ANN kết hợp với PSO ............................... 64 
Hình 3.16. Sai số thực (a) ANN và (b) ANN kết hợp với PSO .................................... 65 
Hình 3.17. Kết quả phát hiện hư hỏng của PSO, ANN và ANNPSO của phần tử số 9:
 ....................................................................................................................................... 66 
Hình 3.18. Giá trị hồi quy (a) ANN và (b) ANN kết hợp với PSO ............................... 67 
Hình 3.19. Sai số thực (a) ANN và (b) ANN kết hợp với PSO .................................... 67 
Hình 3.20. Kết quả phát hiện hư hỏng của PSO, ANN và ANNPSO trong dầm: 50% 
hư hỏng của phần tử 3 và 70% hư hỏng của phần tử 5. ................................................ 68 
Hình 3.21. Cầu Bến Thủy; (a) Kết cấu nhịp; (b) Mặt cắt ngang của cầu ...................... 69 
Hình 3.22. Mô hình phần tử hữu hạn cầu Bến Thủy ..................................................... 70 
Hình 3.23. Hình dạng dao động của một số mode tính toán ......................................... 72 
Hình 3.24. Kiến trúc  ... 08;109;308;309 
 120;320 
 132;332 
 143;144;343;344 
 151]+0.02]; 
DOF=removedof(DOF,seldof); 
% Assembly of stiffness matrix K 
[K,M]=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); 
% bvib manual 
nMode=20; 
[phi,omega]=eigfem(K,M,nMode); 
% Display eigenfrequenties 
disp('Lowest eigenfrequencies [Hz]'); 
a=omega/2/pi; 
idx=(1:15); 
% Display eigenfrequenties 
% disp('Lowest eigenfrequencies [Hz]'); 
Show_animate = 0; 
for show_01 = 1:10 
 if Show_animate == 0 
 % Disp Max 
 figure; 
 plotdisp(Nodes,Elements,Types,DOF,phi(:,show_01),'DispMax','off') 
 title(['Eigenmode', num2str(show_01), '**** Frequency: ' , num2str(a(show_01))]) 
 elseif Show_animate == 1 
 % Animate eigenmodes 
 figure; 
 animdisp(Nodes,Elements,Types,DOF,phi(:,show_01)) 
 title(['Eigenmode', num2str(show_01)]) 
 elseif Show_animate == 2 
 close all 
 continue 
18 
 end 
end 
figure; 
subplot(2,2,1),plotdisp(Nodes,Elements,Types,DOF,phi(:,1),'DispMax','off'); 
subplot(2,2,2),plotdisp(Nodes,Elements,Types,DOF,phi(:,2),'DispMax','off'); 
subplot(2,2,3),plotdisp(Nodes,Elements,Types,DOF,phi(:,3),'DispMax','off'); 
subplot(2,2,4),plotdisp(Nodes,Elements,Types,DOF,phi(:,4),'DispMax','off'); 
 % Animate eigenmodes 
figure; 
animdisp(Nodes,Elements,Types,DOF,phi(:,1)) 
title('Eigenmode 1') 
for i = 1:nMode 
 fre_FEM=a(i); 
 disp(['Mode ',num2str(i),' - Frequency: ', num2str(a(i))]) 
end 
✓ ANN 
clc; 
clear; 
close all; 
rng default 
tic 
%% 
% filename = 'ANN_Benthuy_04'; 
filename = 'Case2_Benthuy.xlsx'; 
sheetname1 = 'Sheet1'; 
sheetname2 = 'Sheet2'; 
input = xlsread(filename,sheetname1,'A1:Z10000'); %call datas from sheetname1 
target = xlsread(filename,sheetname2,'A1:Z10000'); %call datas from sheetname2 
inputs=input'; 
targets=target'; 
% x=inputs; 
x = awgn(inputs,90,'measured'); 
t = targets; 
19 
trainFcn = 'trainlm'; % Levenberg-Marquardt backpropagation. 
% Create a Fitting Network 
hiddenLayerSize = 10; 
net = fitnet(hiddenLayerSize,trainFcn); 
net=feedforwardnet(hiddenLayerSize,trainFcn); 
net.divideParam.trainRatio =70/100; 
net.divideParam.valRatio = 15/100; 
net.divideParam.testRatio = 15/100; 
net.trainParam.epochs = 1000; 
[net] = train(net,x,t) 
y = net(x); 
performance = perform(net,t,y) 
%% 
nntraintool('plot','plotregression'); 
set(findall(gcf,'-property','FontSize'),'FontSize',10); 
for i1 = 1:4 
 subplot(2,2,i1) 
 grid on 
end 
nntraintool('plot','plotperform'); 
set(findall(gcf,'-property','FontSize'),'FontSize',10); 
grid on; 
nntraintool('plot','plottrainstate'); 
set(findall(gcf,'-property','FontSize'),'FontSize',10); 
grid on; 
nntraintool('plot','ploterrhist'); 
set(findall(gcf,'-property','FontSize'),'FontSize',10); 
grid on; 
timeElapsed = toc 
20 
Giao diện của ANN truyền thống trên nền tảng chương trình MATLAB 
✓ ANNPSO 
%% Empty space 
clc 
clear 
rng default 
tic; 
%% Read data 
% Read the input file 
filename = 'Case2_Benthuy.xlsx'; 
sheet_input = 'Sheet1'; 
sheet_output = 'Sheet2'; 
input = xlsread(filename,sheet_input,'A1:Z10000'); 
target = xlsread(filename,sheet_output,'A1:Z10000'); 
inputs = awgn(input',100,'measured'); 
% inputs=input'; 
targets=target'; 
% Select number neuron of input, output and hidden layer 
21 
number_input=length(inputs(:,1)); 
number_output=length(targets(:,1)); 
number_hidden=100; 
 %% STEP 1: TRAIN DATA WITH ANN 
% Chose your train function in ANN 
trainFcn = 'trainlm'; 
net = feedforwardnet(number_hidden,trainFcn); 
net = configure(net,inputs,targets); 
net = fitnet(number_hidden,trainFcn); 
 % Close train window 
net.trainParam.showWindow = 0; 
 % Setup Division of Data for Training, Validation, Testing 
net.divideParam.trainRatio = 70/100; 
net.divideParam.valRatio = 15/100; 
net.divideParam.testRatio = 15/100; 
% Train 
[net] = train(net,inputs,targets) 
 % error NMSE PSO optimized NN 
outputs=net(inputs); 
error_1=outputs - targets; 
err_1st = mean(error_1.^2)/mean(var(targets',1));%.^ means square. 
%% 
x=inputs; 
t = targets; 
y = net(x); 
perf = perform(net,t,y) 
%% 
% mean is used to calculate mean values. For example. A= [1 2 3 4] 
% mean (A)=(1+2+3+4)/4; 
%var returns the variance of the elements of A along the first array 
%dimension; see equation of var on Mathwork. 
err_1st=norm(err_1st) 
disp (['the differen between output and target is: ' num2str(err_1st)]) 
22 
% get the normal NN weights and bias 
weight_bias = getwb(net); 
% setting up parameters 
parameters.MaxIt = 50; 
parameters.nPop = 100; 
parameters.Lower_Boundary = min(weight_bias); 
parameters.Upper_Boundary = max(weight_bias); 
% parameters.Lower_Boundary = -1.5; 
% parameters.Upper_Boundary = 1.5; 
parameters.c1 = 2; 
parameters.c2 = 2; 
parameters.w = 1; 
parameters.w_damp = 0.99; 
parameters.nn=err_1st; 
 % Max epoch 
max_epoch = 1000; 
for epoch = 1 : max_epoch 
 % numWeightElements(nWE) 
nWE=number_hidden*number_input+number_output*number_hidden+number_hidden+number_outp
ut; 
 %% STEP 2: TRAIN WITH PSO 
% disp('STEP 2: TRAIN WITH PSO') 
 [FinalCost,FinalPosition,BestCosts] = PSO (parameters,weight_bias,nWE, 
number_hidden,number_input,number_output,net,inputs,targets); 
 xo = FinalPosition(1,:); 
 xbest = xo; 
 ybest = FinalCost; 
 %% STEP 3: TRAIN NN SECOND TIME 
 disp('STEP 3: TRAIN NN SECOND TIME') 
 % Setup Division of Data for Training, Validation, Testing 
 net_2nd.divideParam.trainRatio = 70/100; 
 net_2nd.divideParam.valRatio = 15/100; 
 net_2nd.divideParam.testRatio = 15/100; 
 % Chose your train function in ANN 
23 
 trainFcn_2nd = 'trainlm'; 
 net_2nd=feedforwardnet(number_hidden,trainFcn_2nd); 
 net_2nd=configure(net_2nd,inputs,targets); 
 [a, b]=min(ybest); 
 xo=xbest(b,:); 
 k=0; 
 % Setup global best for weight and bias in neural network 
 % inputlayers*hiddenlayers 
 xi=zeros(number_hidden,number_input); 
 for i=1:number_hidden 
 for j=1:number_input 
 k=k+1; 
 xi(i,j)=xo(k); 
 end 
 end 
 % outputlayers*hiddenlayers 
 xl=zeros(number_output,number_hidden); 
 for i=1:number_hidden 
 for j=1:number_output 
 k=k+1; 
 xl(j,i)=xo(k); 
 end 
 end 
 % hiddenlayers 
 xb1=zeros(number_hidden,1); 
 for i=1:number_hidden 
 k=k+1; 
 xb1(i,1)=xo(k); 
 end 
 % outputlayers 
 xb2=zeros(number_output,1); 
 for i=1:number_output 
 k=k+1; 
24 
 xb2(i,1)=xo(k); 
 end 
 % Chose your train function in ANN 
 trainFcn = 'trainlm'; 
 net_2nd = feedforwardnet(number_hidden,trainFcn); 
 net_2nd = configure(net,inputs,targets); 
 net_2nd = fitnet(number_hidden,trainFcn); 
 % Close train window 
 net_2nd.trainParam.showWindow = 1; 
 % Train 
 [net_2nd] = train(net_2nd,inputs,targets); 
 outputs_2=net_2nd(inputs); 
 error_2=outputs_2 - targets; 
 err_2nd = mean(error_2.^2)/mean(var(targets',1)) 
 x=inputs; 
t = targets; 
y = net_2nd(x); 
perf3 = perform(net_2nd,t,y) 
 % This line will caculate for 2 outputs 
 err_2nd=norm(err_2nd); 
 disp (['the differen between output and target is: ' num2str(err_2nd)]) 
 disp (['Epoch --------------------------- ' num2str(epoch)]) 
 if err_2nd < err_1st 
 disp ('complete model') 
 break %Terminate execution of for or while loop 
 else 
 weight_bias = getwb(net_2nd); 
 end 
end 
disp('Complete Final model net_f') 
toc; 
%% Test 
disp('Testing out put....') 
25 
% Test for all value in used data 
Test_ouputs1=net_2nd(inputs); 
Test_ouputs1 = [Test_ouputs1; targets]; 
 nntraintool('plot','plotregression') 
set(findall(gcf,'-property','FontSize'),'FontSize',10) 
for i1 = 1:4 
 subplot(2,2,i1) 
 grid on 
end 
nntraintool('plot','plotperform') 
set(findall(gcf,'-property','FontSize'),'FontSize',10) 
grid on; 
nntraintool('plot','plottrainstate') 
set(findall(gcf,'-property','FontSize'),'FontSize',10) 
grid on; 
nntraintool('plot','ploterrhist') 
set(findall(gcf,'-property','FontSize'),'FontSize',10) 
grid on; 
timeElapsed = toc 
✓ Dầm trong phòng thí nghiệm 
clear all; 
% materials 
E = 1.96E11; poisson = 0.30; 
% E = 1.82095E11; poisson = 0.30; 
% matriz C 
C=E/(1-poisson^2)*[1 poisson 0;poisson 1 0;0 0 (1-poisson)/2]; 
 % % load 
% P = 1e6; 
rho=7810; 
thickness=15e-3; 
% thickness=6e-3; 
%Mesh generation 
Lx=800e-3; 
26 
Ly=50e-3; 
% Lx=0.6; % hauteur de la section droite de la poutre 
% Ly=38e-3; % largeur de la section droite de la poutre 
L=Lx; 
B=Ly; 
numberElementsX=101; 
numberElementsY=1; 
Nx=numberElementsX; 
Ny=numberElementsY; 
numberElements=numberElementsX*numberElementsY; 
[nodeCoordinates, elementNodes] = MeshRectangular(L,B,Nx,Ny); 
% rectangularMesh(Lx,Ly,numberElementsX,numberElementsY); 
xx=nodeCoordinates(:,1); 
yy=nodeCoordinates(:,2); 
% drawingMesh(nodeCoordinates,elementNodes,'Q4','k-'); 
numberNodes=size(xx,1); 
% GDof: global number of degrees of freedom 
GDof=2*numberNodes; 
[stiffness,mass]=formStiffness2D(GDof,numberElements,... 
 elementNodes,numberNodes,nodeCoordinates,C,rho,thickness); 
 [Val_prop_Range,Indice] = sort(diag(Val_prop)); % définition de l'ordre des indices 
Vect_prop_range = Vect_prop(:,Indice); % Matrice ordonnée (Vecteurs propres) 
K_range=stiffness; 
M_range=mass; 
Puls_prop = Val_prop_Range.^(1/2) ;% Pulsation propre (en rad/s) 
freqHz = Puls_prop/(2*pi) ; % UNITS : Hertz 
freqHzz=freqHz(4:10,:) 
➢ Cầu Bến Quan 
clc; 
close all; 
% bvib manual 
% Tutorial: static analysis 
% Units: m, kN 
27 
% Nodes=[NodID X Y Z] 
lspan = 12; 
%dy = 0.0001; %m, distance between nodes at a truss joint 
b=0; 
% reference node 
% Check the node coordinates as follows: 
figure 
plotnodes(Nodes); 
Element types -> {EltTypID EltName} 
% Sections=[SecID A ky kz Ixx Iyy Izz yt yb zt zb] 
Sections=[1 A inf inf 0 0 Izz yt yb b/2 b/2]; 
% 2 pi*r^2 NaN NaN NaN NaN NaN NaN NaN NaN NaN]; 
% Materials=[MatID E nu]; 
E=2E11; 
Materials= [1 E 0.3 7850] 
Code=1; 
% Elements=[EltID(ten phan tu) TypID SecID MatID n1 n2 n3)] 
 Elements=Elements_Sources_beam_Benquan(Code); 
figure; 
plotnodes(Nodes,'numbering','off'); 
hold('on') 
plotelem(Nodes,Elements,Types,'numbering','on'); 
DOF=getdof(Elements,Types); 
seldof = [0.01;0.04 
 [Edge4]+0.02 
 [Edge4]+0.03 
 [Edge4]+0.05 
 [Edge1;Edge3]+0.05 
 Edge2+0.03 
 Edge2+0.05]; 
DOF=removedof(DOF,seldof); 
 [K,M]=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); 
% bvib manual 
% Tutorial: dynamic analysis: eigenvalue problem 
28 
% Units: m, N 
% Assembly of M and K 
% Eigenvalue problem 
nMode=20; 
[phi,omega]=eigfem(K,M,nMode); 
% Display eigenfrequenties 
% disp('Lowest eigenfrequencies [Hz]'); 
a=omega/2/pi; 
idx=[1 2 3 4]; 
figure; 
plotdisp(Nodes,Elements,Types,DOF,phi(:,1),'DispMax','off') 
figure; 
plotdisp(Nodes,Elements,Types,DOF,phi(:,2),'DispMax','off') 
figure; 
plotdisp(Nodes,Elements,Types,DOF,phi(:,3),'DispMax','off') 
figure; 
plotdisp(Nodes,Elements,Types,DOF,phi(:,4),'DispMax','off') 
fre_FEM=a(idx)