Luận án Phân tích động lực học kết cấu công trình biển hệ thanh cố định trên nền san hô chịu tác dụng của tải trọng sóng biển và gió theo mô hình bài toán không gian

BỘ GIÁO DỤC VÀ ĐÀO BỘ QUỐC PHÒNG 
HỌC VIỆN KỸ THUẬT QUÂN SỰ 
Lê Hoàng Anh 
PHÂN TÍCH ĐỘNG LỰC HỌC KẾT CẤU CÔNG TRÌNH BIỂN HỆ THANH 
CỐ ĐỊNH TRÊN NỀN SAN HÔ CHỊU TÁC DỤNG CỦA TẢI TRỌNG 
SÓNG BIỂN VÀ GIÓ THEO MÔ HÌNH BÀI TOÁN KHÔNG GIAN 
LUẬN ÁN TIẾN SĨ KỸ THUẬT 
Hà Nội - 2016 
BỘ GIÁO DỤC VÀ ĐÀO TẠO BỘ QUỐC PHÒNG 
HỌC VIỆN KỸ THUẬT QUÂN SỰ 
Lê Hoàng Anh 
PHÂN TÍCH ĐỘNG LỰC HỌC KẾT CẤU CÔNG TRÌNH BIỂN HỆ THANH 
CỐ ĐỊNH TRÊN NỀN SAN HÔ CHỊU TÁC DỤNG CỦA TẢI TRỌNG 
SÓNG BIỂN VÀ GIÓ THEO MÔ HÌNH BÀI TOÁN KHÔNG GIAN 
Chuyên ngành: Cơ kỹ thuật 
Mã số: 62.52.01.01 
LUẬN ÁN TIẾN SĨ KỸ THUẬT 
Người hướng dẫn khoa học: 
PGS.TS Nguyễn Thái Chung 
Hà Nội - 2016 
 i
LỜI CAM ĐOAN 
 Tôi là Lê Hoàng Anh, tôi xin cam đoan đây là công trình nghiên cứu 
của riêng tôi. Các số liệu, kết quả trong luận án là trung thực và chưa từng 
được ai công bố trong bất kỳ công trình nào. 
Tác giả 
Lê Hoàng Anh 
 ii
LêI C¶M ¥N 
 T¸c gi¶ luËn ¸n xin bµy tá lßng biÕt ¬n ch©n thµnh ®èi víi 
PGS.TS NguyÔn Th¸i Chung ®· tËn t×nh h−íng dÉn, gióp ®ì vµ cho 
nhiÒu chØ dÉn khoa häc cã gi¸ trÞ, gióp cho t¸c gi¶ hoµn thµnh luËn ¸n 
nµy. T¸c gi¶ tr©n träng sù ®éng viªn, khuyÕn khÝch vµ nh÷ng kiÕn thøc 
khoa häc còng nh− chuyªn m«n mµ ThÇy h−íng dÉn ®· chia sÎ cho t¸c 
gi¶ trong nhiÒu n¨m qua, gióp cho t¸c gi¶ n©ng cao n¨ng lùc, ph−¬ng 
ph¸p nghiªn cøu khoa häc. 
 T¸c gi¶ tr©n träng c¶m ¬n tËp thÓ Bé m«n C¬ häc vËt r¾n, Khoa 
C¬ khÝ, Phßng Sau ®¹i häc - Häc viÖn Kü thuËt qu©n sù vµ Tr−êng 
§¹i häc C«ng nghÖ giao th«ng vËn t¶i ®· t¹o mäi ®iÒu kiÖn thuËn lîi 
cho t¸c gi¶ trong qu¸ tr×nh nghiªn cøu. T¸c gi¶ xin tr©n träng c¶m ¬n 
GS.TS.NGND Hoµng Xu©n L−îng - Häc viÖn Kü thuËt qu©n sù, 
GS.TSKH NguyÔn TiÕn Khiªm - ViÖn C¬ häc ®· cung cÊp cho t¸c gi¶ 
nhiÒu tµi liÖu quý hiÕm, c¸c kiÕn thøc khoa häc hiÖn ®¹i vµ nhiÒu lêi 
khuyªn bæ Ých, cã gi¸ trÞ. 
 Cuèi cïng, t¸c gi¶ xin bµy tá lßng c¶m ¬n ®èi víi nh÷ng ng−êi 
th©n trong gia ®×nh ®· th«ng c¶m, ®éng viªn vµ chia sÎ nh÷ng khã 
kh¨n víi t¸c gi¶ trong suèt thêi gian lµm luËn ¸n. 
 T¸c gi¶ 
 iii
MỤC LỤC 
Lời cam đoan ............................................................................................................... i 
Mục lục.... ........................................................................................................... ii 
Danh mục các ký hiệu, các chữ viết tắt ..................................................................... vi 
Danh mục các bảng ................................................................................................... ix 
Danh mục các hình vẽ, đồ thị ................................................................................... xi 
MỞ ĐẦU.. ...................................................................................................... 1 
CHƯƠNG 1: TỔNG QUAN VẤN ĐỀ NGHIÊN CỨU............................................. 6 
1.1. Sơ lược về san hô và nền san hô .......................................................................... 6 
1.2. Công trình biển và tải trọng phổ biến tác dụng lên công trình biển ..................... 7 
1.2.1. Tổng quan về công trình biển ........................................................................... 7 
1.2.2. Tổng quan về tải trọng tác dụng lên công trình biển ...................................... 11 
1.2.2.1. Tải trọng sóng biển....................................................................................... 12 
1.2.2.2. Tải trọng gió ................................................................................................. 13 
1.3. Tổng quan về tính toán công trình biển ............................................................. 15 
1.3.1. Tình hình nghiên cứu trong nước .................................................................... 15 
1.3.2. Tình hình nghiên cứu ỏ nước ngoài ................................................................ 17 
1.4. Các kết quả nghiên cứu đạt được từ các công trình đã công bố ........................ 21 
1.5. Các vấn đề cần tiếp tục nghiên cứu .................................................................... 23 
1.6. Kết luận rút ra từ tổng quan ............................................................................... 23 
CHƯƠNG 2: PHÂN TÍCH ĐỘNG LỰC HỌC CÔNG TRÌNH BIỂN CỐ ĐỊNH 
CHỊU TÁC DỤNG CỦA TẢI TRỌNG SÓNG VÀ GIÓ  ................................ 25 
2.1. Đặt vấn đề .......................................................................................................... 25 
2.2. Giới thiệu bài toán và các giả thiết ..................................................................... 26 
2.3. Thiết lập các phương trình cơ bản của bài toán ................................................. 28 
2.3.1. Các quan hệ đối với phần tử thanh thuộc công trình ...................................... 28 
2.3.1.1. Trường chuyển vị ........................................................................................ 28 
2.3.1.2. Trường biến dạng ........................................................................................ 29 
 iv
2.3.1.3. Trường ứng suất .......................................................................................... 31 
2.3.1.4. Phương trình mô tả dao động của phần tử trong hệ tọa độ cục bộ .............. 31 
2.3.1.5. Phương trình mô tả dao động của phần tử trong hệ tọa độ tổng thể ............ 33 
2.3.2. Các quan hệ đối với phần tử thuộc các lớp nền san hô ................................... 34 
2.3.2.1. Các phương trình cơ bản của phần tử ......................................................... 34 
2.3.2.2. Phương trình mô tả dao động của phần tử .................................................. 37 
2.3.3. Quan hệ đối với phần tử tiếp xúc giữa thanh và nền san hô ........................... 38 
2.3.4. Tải trọng sóng và gió tác dụng lên công trình................................................. 42 
2.3.4.1. Tải trọng sóng tác dụng lên phần tử thanh ................................................... 42 
2.3.4.2. Tải trọng gió tác dụng lên công trình ........................................................... 44 
2.4. Xây dựng phương trình mô tả dao động của hệ ................................................. 46 
2.4.1. Tập hợp ma trận và véc tơ toàn hệ .................................................................. 46 
2.4.1.1. Tập hợp ma trận độ cứng tổng thể [K] – hàm assem() ................................ 46 
2.4.1.2. Tập hợp véc tơ tải trọng tổng thể {P} – hàm inset () ................................... 47 
2.4.2. Phương trình mô tả dao động của hệ............................................................... 48 
2.4.3. Điều kiện biên ................................................................................................. 50 
2.5. Thuật toán PTHH phân tích động lực học của hệ kết cấu công trình biển 
và nền san hô ............................................................................................................. 50 
2.6. Chương trình tính và kiểm tra độ tin cậy của chương trình tính ........................ 56 
2.6.1. Chương trình tính ............................................................................................ 56 
2.6.2. Kiểm tra độ tin cậy của chương trình .............................................................. 56 
2.6. Kết luận chương 2 .............................................................................................. 59 
CHƯƠNG 3: KHẢO SÁT ẢNH HƯỞNG CỦA MỘT SỐ YẾU TỐ ĐẾN 
PHẢN ỨNG ĐỘNG CỦA CÔNG TRÌNH BIỂN CỐ ĐỊNH CHỊU TÁC DỤNG 
CỦA TẢI TRỌNG SÓNG VÀ GIÓ  ................................................................ 60 
3.1. Đặt vấn đề .......................................................................................................... 60 
3.2. Bài toán xuất phát ............................................................................................... 60 
3.3. Khảo sát ảnh hưởng của một số yếu tố đến phản ứng động của hệ ................... 67 
3.3.1. Ảnh hưởng của mô hình tính........................................................................... 67 
 v
3.3.2. Ảnh hưởng của dạng kết cấu ........................................................................... 71 
3.3.3. Ảnh hưởng của vật liệu kết cấu ....................................................................... 74 
3.3.3.1. Ảnh hưởng của mô đun đàn hồi cọc chính .................................................. 74 
3.3.3.2. Ảnh hưởng của mô đun đàn hồi cọc phụ ..................................................... 78 
3.3.4. Ảnh hưởng của đường kính ngoài cọc chính .................................................. 81 
3.3.5. Ảnh hưởng của đường kính ngoài cọc phụ ..................................................... 85 
3.3.6. Ảnh hưởng của lớp nền san hô ........................................................................ 88 
3.3.7. Ảnh hưởng của tải trọng .................................................................................. 91 
3.3.7.1. Ảnh hưởng của giản đồ vận tốc gió ............................................................. 91 
3.3.7.2. Ảnh hưởng của phương tải trọng gió ........................................................... 95 
3.3.7.3. Ảnh hưởng của chiều cao sóng biển ............................................................ 99 
3.4. Kết luận chương 3 ............................................................................................ 102 
CHƯƠNG 4: NGHIÊN CỨU PHẢN ỨNG ĐỘNG CỦA KẾT CẤU HỆ THANH 
MÔ PHỎNG CÔNG TRÌNH BIỂN BẰNG THỰC NGHIỆM ............................. 104 
4.1. Mục đích thí nghiệm ........................................................................................ 104 
4.2. Mô hình và các thiết bị thí nghiệm .................................................................. 105 
4.2.1. Mô hình thí nghiệm ....................................................................................... 105 
4.2.2. Thiết bị thí nghiệm ........................................................................................ 106 
4.2.2.1. Các thiết bị gây tải...................................................................................... 106 
4.2.2.2. Thiết bị cảm biến gia tốc biến dạng ........................................................... 107 
4.2.2.3. Máy đo dao động ........................................................................................ 108 
4.3. Phương pháp xác định gia tốc, biến dạng của kết cấu ..................................... 109 
4.4. Cơ sở phân tích và xử lý số liệu thí nghiệm ..................................................... 111 
4.5. Thí nghiệm và kết quả thí nghiệm ................................................................... 112 
4.5.1. Thí nghiệm xác định đáp ứng động của hệ liên hợp giàn thép - bể chứa trên 
nền san hô bãi cạn ven đảo Song Tử Tây chịu tác dụng của xung lực va chạm .... 112 
4.5.1.1. Mô tả thí nghiệm ........................................................................................ 112 
4.5.1.2. Thí nghiệm và kết quả ................................................................................ 113 
 vi
4.5.2. Thí nghiệm xác định đáp ứng động của hệ giàn thép không gian mô phỏng 
một dạng công trình DKI tại bể tạo sóng ................................................................ 120 
4.5.2.1. Mô tả thí nghiệm ........................................................................................ 120 
4.5.2.2. Thí nghiệm và kết quả ................................................................................ 121 
4.6. Kết luận chương 4 ............................................................................................ 125 
KẾT LUẬN VÀ KIẾN NGHỊ................................................................................. 126 
DANH MỤC CÁC CÔNG TRÌNH CÔNG BỐ CỦA TÁC GIẢ ........................... 129 
TÀI LIỆU THAM KHẢO ....................................................................................... 131 
PHỤ LỤC. ........................................................................................................... 141 
 vii
DANH MỤC CÁC KÝ HIỆU, CÁC CHỮ VIẾT TẮT 
1. Danh mục các ký hiệu 
1.1. Các ký hiệu bằng chữ La tinh 
ax, ay, az Gia tốc hạt nước theo các phương x, y, z, 
[B], [ ]iB Ma trận đạo hàm, ma trận nội suy của phần tử, 
[B]e, [B]se Ma trận quan hệ biến dạng – chuyển vị của phần tử, PTTX 
CD, C1, Cp Hệ số lực cản, hệ số lực quán tính, hệ số áp lực gió, 
Cw, cx, cy Tốc độ truyền sóng và các cosin chỉ phương 
[C], [C]e, Ma trận cản, ma trận cản phần tử, 
t tC +∆   Ma trận cản phụ thuộc thời gian, 
Dch, Dph Đường kính ngoài cọc chính, đường kính ngoài cọc phụ, 
[D], [Dse] Ma trận quan hệ ứng suất – biến dạng của phần tử, PTTX, 
E, Ech, Eph, Ef Môđun đàn hồi của vật liệu, cọc chính, cọc phụ, nền san hô 
{F}e Véc tơ tải trọng quy nút của phần tử, 
fx, fy, fz Tải trọng tác dụng lên phẩn tử thanh theo phương x, y và z, 
{f}, {f0}, {fs}eVéc tơ tải trọng tổng thể, ngoại lực, tải trọng nút phần tử 
{ }
w
e
f , { }win
e
f Véc tơ tải trọng nút do sóng, gió tác dụng lên phần tử thanh, 
G, {g} Mô đun đàn hồi biến dạng trượt, véc tơ lực thể tích, 
Jz Mômen quán tính trục của mặt cắt ngang phần tử thanh, 
knz, kres Độ cứng pháp tuyến theo phương z, chống trượt tới hạn 
ksx, ksy Độ cứng tiếp tuyến theo phương x và phương y, 
kη, kξ Độ cứng pháp tuyến, tiếp tuyến của PTTX, 
[K], [K]e Ma trận độ cứng tổng thể, phần tử 
[Kb]e, [Ks]e Ma trận độ cứng của phần tử thanh, PTTX 
[K*] Ma trận độ cứng hiệu quả, 
kw, Tw, Hw, Lw, Số sóng, chu kỳ sóng, chiều cao sóng, chiều dài bước sóng, 
lx, mx, nx, ly, my, ny, lz, mz, nz Cosin chỉ phương của trục x, y, z 
 viii 
[M], [M]e Ma trận khối lượng, ma trận khối lượng phần tử, 
[N], [N]e Ma trận hàm dạng, ma trận các hàm dạng của phần tử, 
[Nu(ξ)],[Nv(ξ)],Nw(ξ)] Véc tơ hàng của hàm dạng chuyển vị ngang theo 
phương x, y, z 
[Nθx(ξ)], [Nθy(ξ)], [Nθz(ξ)] Véc tơ hàng của hàm dạng chuyển vị xoay 
quanh trục x, y, z 
pwin(t) Áp lực gió tác động lên diện tích của công trình theo thời gian, 
{P}, {P*} Véc tơ tải trọng tổng thể, véctơ tải trọng hiệu quả, 
( )winq t Lực gió phân bố theo chiều dài của thanh, 
{q}T, {q}e Véc tơ chuyển vị nút phần tử 
{ }t tq +∆ Luỹ tích của véc tơ chuyển vị nút 
tch, tph Chiều dày thành ống cọc chính, cọc phụ 
[T]e Ma trận chuyển hệ trục tọa độ, 
Uwin(t) Hàm vận tốc gió theo thời gian, 
Ux, xU
 , xU

 Đáp ứng chuyển vị ngang, vận tốc, gia tốc tại đỉnh giàn, 
Vx, Vy, Vz Chuyển vị của hạt nước theo các phương x, y, z, 
1.2. Các ký hiệu bằng chữ Hy Lạp 
α, δ Các tham số trong tích phân Newmark, 
αr, βr Các hằng số cản Rayleigh, 
β: Góc nghiên của cọc chính, 
∆σ, ∆εz Số gia ứng suất và số gia biến dạng theo phương pháp tuyến z, 
∆τzx, ∆τzy Số gia ứng suất và số gia biến dạng trong mặt phẳng xoz, yoz 
∆γzx, ∆γzy Số gia biến dạng trong mặt phẳng xoz, yoz
{ }seU∆ Véc tơ số gia chuyển vị nút của PTTX 
{∆εse}(i) Số gia biến dạng của phần tử tiếp xúc, 
 ... end 
%Stress and internal forces at Slip2D nodes 
 DoIt=True; 
 if DoIt==True 
 CalcSlip2DStress_disk(ElementList,nTime,... 
 T2L_SL.TempDirectory,... 
 T2L_SL.fileNodeDisp,... 
 T2L_SL.fileSlip2DStress); 
 ElementNumber=160; 
 NodeList=NodeToPrint; 
 switch etype 
 case 1 
 Nodes=FS1array(ElementNumber).Node; 
 case 2 
 Nodes=FS2array(ElementNumber).Node; 
 end 
 for ii=2:5 
 switch ii 
 case 1 
 flag=[14 15 16 17]; %[Xicmax Xicmay Xicmaz Tauxy] 
 case 2 
 flag=14; %Xicmax 
 case 3 
 flag=15; %Xicmay 
 case 4 
 flag=16; %Xicmaz 
 case 5 
 flag=17; %Tauxy 
 end 
WindowTitle=[estring'Ung suat phan tu'... 
num2str(ElementNumber)',nut'num2str(NodeToPrint)]; 
 DrawSlip2DStress_disk(ElementNumber,... 
 Nodes,... 
 NodeList,... 
 nTime,dt,... 
 T2L_SL.TempDirectory,... 
 T2L_SL.fileSlip2DStress,... 
 WindowTitle,... 
 flag); 
 end 
 end 
 158 
%Strain and displate at Plate2D nodes 
 DoIt=True; 
 if DoIt==True 
 CalcPlate2DStrain_disk(ElementList,nTime,... 
 T2L_SL.TempDirectory,... 
 T2L_SL.fileNodeDisp,... 
 T2L_SL.filePlate2DStrain); 
 ElementNumber=280; 
 NodeList=NodeToPrint; 
 switch etype 
 case 1 
 Nodes=FS1array(ElementNumber).Node; 
 case 2 
 Nodes=FS2array(ElementNumber).Node; 
 end 
 for ii=2:5 
 switch ii 
 case 1 
 flag=[14 15 16]; %[Strainx Strainy Tauxy] 
 case 2 
 flag=14; %Strainx 
 case 3 
 flag=15; %Strainy 
 case 4 
 flag=16; %Tauxy 
 end 
 WindowTitle=[estring'Bien dang phan tu'... 
num2str(ElementNumber)',nut'num2str(NodeToPrint)]; 
 DrawPlate2DStrain_disk(ElementNumber,... 
 Nodes,... 
 NodeList,... 
 nTime,dt,... 
 T2L_SL.TempDirectory,... 
 T2L_SL.filePlate2DStrain,... 
 WindowTitle,... 
 flag); 
 end 
 end 
%-----------------------CalculateElementParameters.m--------------------- 
%1.For Slip2D 
if T2L_SP.NumberOfSlip2D > 0 
 CalculateElementParametersSlip2D; 
end 
%2.For Plate84 
 159 
if T2L_SP.NumberOfPlate84 > 0 
 CalculateElementParametersPlate84; 
end 
%-------------------CalculateElementParametersSlip2D.m------------------- 
%Calculate element parameters for Slip2D 
global SLV SLR SLTG ax; 
PublicDeclarations; 
MohinhPTHH; 
TsPhantu=T2L_SL.NuberOfSlip2D; 
% Xay dung mang bac tu do nut NDF 
 NDF = zeros(TsNut,6); 
 Bactudo=0; 
 for i=1:TsNut 
 if JF(i,1)== 0 
 Bactudo = Bactudo+1; 
 NDF(i,1)= Bactudo; 
 end 
 if JF(i,2)== 0 
 Bactudo = Bactudo+1; 
 NDF(i,2)= Bactudo; 
 end 
 if JF(i,3)== 0 
 Bactudo = Bactudo+1; 
 NDF(i,3)= Bactudo; 
 end 
 end % KET THUC for i=1:TsNut 
 Neq = Bactudo; % Neq - tong so phuong trinh 
% Xay dung mang bac tu do phan tu ND 
 NED = 12 ; 
 ND = zeros(TsPhantu,NED); 
 for i=1:TsPhantu 
 Nut1 = LNC(i,1); ND(i,1) = NDF(Nut1,1); 
ND(i,2)= NDF(Nut1,2); 
Nut2 = LNC(i,2); ND(i,3) = NDF(Nut2,3); 
ND(i,4)= NDF(Nut2,4); 
 Nut3 = LNC(i,3); ND(i,5) = NDF(Nut3,5); 
ND(i,6)= NDF(Nut3,6); 
 Nut4 = LNC(i,4); ND(i,7)= NDF(Nut4,7); 
ND(i,8)= NDF(Nut4,8); 
 Nut5 = LNC(i,5); ND(i,9)= NDF(Nut5,9); 
ND(i,10)= NDF(Nut5,10); 
 Nut6 = LNC(i,6); ND(i,11)= NDF(Nut6,11); 
ND(i,12)= NDF(Nut6,12); 
 end 
 160 
% Xay dung cau truc ma tran do cung 
% Xay dung mang chieu cao cot CHT 
 knz=AM.E*(1-AM.Nuy)/(1+AM.Nuy)/(1-2*AM.Nuy); 
ksx=ksy=AM.E/(1+AM.Nuy)/2; 
CHT = zeros(Neq,1); 
 for i=1:TsPhantu 
 MinCs=10^6 ; 
 for j=1:NED 
 Chiso = ND(i,j); 
 if Chiso>0 
 if Chiso < MinCs 
 MinCs = Chiso; 
 end 
 end 
 end % ket thuc vong lap for j=1:NED 
 for j=1:NED 
 Chiso = ND(i,j); 
 if Chiso > 0 
 Caocot = Chiso - MinCs; 
 if CHT(Chiso,1)< Caocot 
 CHT(Chiso,1)= Caocot; 
 end 
 end 
 end % ket thuc vong lap for j=1:NED 
 end % ket thuc vong lap for i=1:TsPhantu 
 Nuagiai = max(CHT); 
% Tinh mang luu dia chi cac phan tu tren duong cheo NDS 
 NDS = zeros(Neq+1,1); 
 NDS(1,1)=1; 
 for i=2:(Neq+1) 
 NDS(i,1)=NDS(i-1)+CHT(i-1)+1; 
 end 
 Nsky = NDS(Neq+1,1)-1; % Tong so phan tu trong SK 
% Xay dung tinh chat phan tu 
 SK = zeros(Nsky,1); % Ma tran do cung 
 SM = zeros(Nsky,1); % Ma tran khoi luong 
 SC = zeros(Nsky,1); % Ma tran can nhot 
 PT = zeros(Neq,1); % Vec to tai 
% Vong lap chinh xay dung tinh chat phan tu 
 for Phantu = 1:TsPhantu 
 X8 = zeros(8,1); Y8 = zeros(8,1); Z8 = zeros(8,1); 
 for i=1:8 
 X8(i,1)= XYZ(LNC(Phantu,i),1); 
 Y8(i,1)= XYZ(LNC(Phantu,i),2); 
 161 
 Z8(i,1)= XYZ(LNC(Phantu,i),3); 
 end 
 ChisoVl = SLV.Pt.ChisoVl(Phantu,1); 
 VATLIEU = SLV.Vl.VATLIEU; 
 E = VATLIEU(ChisoVl,1); 
 Mu = VATLIEU(ChisoVl,2); 
 Ro=VATLIEU(ChisoVl,3); 
 [EK,EQ,EM] = ld8(X8,Y8,Z8,E,Mu,Ro); 
% EQ tai trong nut tinh theo trong luong ban than 
% Gui EK vao SK; EQ vao PT; EM vao SM 
 for i=1:NED 
 m = ND(Phantu,i); 
 if m>0 PT(m,1)= PT(m,1)+ EQ(i,1); end 
 for j=1:NED 
 n = ND(Phantu,j); 
 if (m>0)&(n>=m) 
 Chiso = NDS(n,1)+ n-m ; 
 SK(Chiso,1)= SK(Chiso,1)+ EK(i,j); 
 SM(Chiso,1)= SM(Chiso,1)+ EM(i,j); 
 end 
 end; % ket thuc for j=1:NED 
 end; % ket thuc for i=1:NED 
 end ; % ket thuc for Phantu = 1:TsPhantu 
 COMBO = PT; 
 PT = zeros(Neq,1); 
% Nhan vec to tai tu cac tai trong nut 
 TT6 = SLV.Nut.TT6; 
 for i=1:TsNut 
 if NDF(i,1)>0 PT(NDF(i,1),1)= TT6(i,1); end 
 if NDF(i,2)>0 PT(NDF(i,2),1)= TT6(i,2); end 
 if NDF(i,3)>0 PT(NDF(i,3),1)= TT6(i,3); end 
 end; 
 COMBO = [COMBO PT]; 
 SLV.COMBO = COMBO; 
 SLTG.Neq = Neq; 
 SLTG.NDF = NDF; 
 SLTG.ND = ND; 
 SLTG.SK = SK; 
 SLTG.SM = SM; 
 SLTG.PT = PT; 
 SLTG.NDS = NDS; 
 SLTG.maDECOM = 0; 
return 
%-------------------CalculateElementParametersPlate84.m------------------ 
 162 
%Calculate element parameters for Plate84 
function CalculateElementParametersPlate84 
PublicDeclarations; 
MohinhPTHH; 
if symPlate84.Built==False 
 symPlate84.Built=True; 
 [symPlate84.DisplacementShape,... 
 symPlate84.BendingStrainStress,... 
 symPlate84.ShearStrainStress,... 
 symPlate84.BendingMaterialMatrix,... 
 symPlate84.ShearMaterialMatrix]=... 
BuildSymbolicPlate84(symFS1.HermitianFunctions,symFS1.dHdxy); 
end 
for ENumber=1:T2L_SP.NumberOfPlate84 
 Nd1=Plate84array(ENumber).Node(1); 
 Nd2=Plate84array(ENumber).Node(2); 
 Nd3=Plate84array(ENumber).Node(3); 
 Nd4=Plate84array(ENumber).Node(4); 
n1=[aNodalData(Nd1).X,aNodalData(Nd1).Y,aNodalData(Nd1).Z]; 
n2=[aNodalData(Nd2).X,aNodalData(Nd2).Y,aNodalData(Nd2).Z]; 
n3=[aNodalData(Nd3).X,aNodalData(Nd3).Y,aNodalData(Nd3).Z]; 
n4=[aNodalData(Nd4).X,aNodalData(Nd4).Y,aNodalData(Nd4).Z]; 
 imt=Plate84array(ENumber).MaterialType; 
 E=aMaterialType(imt).E; 
 Nuy=aMaterialType(imt).Nuy; 
 Ro=aMaterialType(imt).Ro; 
 th=aMaterialType(imt).thickness; 
 Plate84array(ENumber).BendingMaterialMatrix=... 
 eval(symPlate84.BendingMaterialMatrix); 
 Plate84array(ENumber).ShearMaterialMatrix=... 
 eval(symPlate84.ShearMaterialMatrix); 
 [Plate84array(ENumber).TM12,... 
Plate84array(ENumber).TM3]=CoordinateTransformFS1(n1,n2,n3,n4); 
 ne=(n1+n2+n3+n4)/4; 
 XYZ=[n1-ne n2-ne n3-ne n4-ne]; 
 XYZ=Platearray(ENumber).TM12*XYZ; 
 en1=XYZ(1:3); 
 en2=XYZ(4:6); 
 en3=XYZ(7:9); 
 en4=XYZ(10:12); 
 flag=0; 
 [Plate84array(ENumber).Stiffness,M]=... 
 StiffMassPlate84(en1,en2,en3,en4,... 
 Plate84array(ENumber).thickness,... 
 163 
 aMaterialType(imt).Ro,... 
 symFS1.JacobianMatrix,... 
 symPlate84.DisplacementShape,... 
 symPlate84.BendingStrainStress,... 
 symPlate84.ShearStrainStress,... 
 Plate84array(ENumber).BendingMaterialMatrix,... 
 Plate84array(ENumber).ShearMaterialMatrix,... 
 GaussIntegrationConstant.Gp1,... 
 GaussIntegrationConstant.Gw1,... 
 GaussIntegrationConstant.Gip,... 
 GaussIntegrationConstant.Gwf,... 
 flag); 
 if SolveOpt.DynamicCalculus==True 
 Plate84(ENumber).Mass=M; 
 Plate84(ENumber).Damping=C; 
 end 
if~isempty(Plate84array(ENumber).StaticUniformLoadSurfaceGCS) 
Plate84array(ENumber).StaticUniformLoadOnTheSurface=... 
Plate84array(ENumber).StaticUniformLoadOnTheSurface+... 
 Plate84array(ENumber).TM12*... 
Plate84array(ENumber).StaticUniformLoadSurfaceGCS; 
 end 
if~isempty(Plate84array(ENumber).DynamicUniformLoadSurfaceGCS) 
Plate84array(ENumber).DynamicUniformLoadOnTheSurface=... 
Plate84array(ENumber).DynamicUniformLoadOnTheSurface+... 
 Plate84array(ENumber).TM12*... 
Plate84array(ENumber).DynamicUniformLoadSurfaceGCS; 
 end 
if~isempty(Plate84array(ENumber).StaticUniformLoadOnTheSurface) 
 Rs=SurfaceLoadPlate84(... 
 en1,en2,en3,en4,... 
Plate84array(ENumber).StaticUniformLoadOnTheSurface,... 
 symFS1.H3x12,... 
 symFS1.JacobianMatrix,... 
 GaussIntegrationConstant.Gip,... 
 GaussIntegrationConstant.Gwf); 
 Plate84array(ENumber).CalcElemLoadStatic=... 
 Plate84array(ENumber).CalcElemLoadStatic+Rs; 
 end 
 if SolveOpt.GravityLoad==True 
 Rg=GravityLoadPlate84(... 
 en1,en2,en3,en4,... 
 Plate84array(ENumber).Thickness,... 
 aMaterialType(imt).Ro,... 
 164 
 T2L_SP.Gravity,... 
 Plate84array(ENumber).TM3,... 
 symFS1.H3x12,... 
 symFS1.JacobianMatrix,... 
 GaussIntegrationConstant.Gip,... 
 GaussIntegrationConstant.Gwf); 
 Plate84array(ENumber).CalcElemLoadStatic=... 
 Plate84array(ENumber).CalcElemLoadStatic+Rg; 
 end 
 if SolveOpt.DynamicCalculus==True 
if~isempty(Plate84array(ENumber).DynamicUniformLoadOnTheSurface) 
 Rs=SurfaceLoadPlate84(... 
 en1,en2,en3,en4,... 
Plate84array(ENumber).DynamicUniformLoadOnTheSurface,... 
 symFS1.H3x12,... 
 symFS1.JacobianMatrix,... 
 GaussIntegrationConstant.Gip,... 
 GaussIntegrationConstant.Gwf); 
 Plate84array(ENumber).CalcElemLoadDynamic=... 
 Plate84array(ENumber).CalcElemLoadDynamic+Rs; 
 end 
 end 
end 
%---------------------------Matranchiso.m-------------------------------- 
% Xay dung mang chi so phan tu LNC va so do Skyline 
MohinhPTHH; %Goi ket qua chia mo hinh PTHH 
%tu file MohinhPTHH.inp cua ANSYS 12.1 
CsPhantu = 1:TsPhantu; 
LNC = zeros(TsPhantu,4) ; 
Phantu = 0; 
 for j=1:nY 
 for i=1:nX 
 Phantu = Phantu + 1; 
 Nut1 = (j-1)*(nX+1)+i ; Nut2 = Nut1+1; 
 Nut4 = (j-1)*(nX+1)+i ; Nut3 = Nut4+1; 
 LNC(Phantu,1)= Nut2; LNC(Phantu,4)= Nut1; 
 LNC(Phantu,5)= Nut3; LNC(Phantu,8)= Nut4; 
 end % Ket thuc for i=1:nX 
 end % Ket thuc for j=1:nY 
 ChisoHh = ones(TsPhantu,1); 
 ChisoVl = ones(TsPhantu,1); 
ChisoTt = zeros(TsPhantu,1); 
ChonPhantu = zeros(TsPhantu,1); % chua chon phan tu 
 SLV.Pt.TsPhantu = TsPhantu; 
 165 
 SLV.Pt.LNC = LNC; 
 SLV.Chon.Phantu = ChonPhantu; 
 HT = [0 0 0 0 0 0 0]; % chua hien thi thong tin 
 SLV.Hienthi = HT; 
 COMBO = []; % to hop tai trong 
 SLV.COMBO = COMBO; 
%------------------------------Newmark.m--------------------------------- 
%============================FOR FRAME_W1_2012=========================== 
function [Cvi,Vtoc,Giatoc,Thgian]=TichPhanNewmark(nBuocTg,BuocTg) 
global alfa delta F0 Tanso SoBacTdo X Y NDF 
% Cac he so Ai 
A0=1/(alfa*BuocTg*BuocTg); A1=delta/(alfa*BuocTg); 
A2=1/(alfa*BuocTg); A3=0.5/alfa-1; 
A4=delta/alfa-1; A5=0.5*BuocTg*(delta/alfa-2); 
A6=BuocTg*(1-delta); A7=delta*BuocTg; 
[SK,SM]=MangPtuMaTrDoCung; SC=0.5*SM; 
K_hqua=SK+A0*SM+A1*SC; 
Fnut=zeros(SoBacTdo,1); 
Cvi=zeros(SoBacTdo,nBuocTg); 
Vtoc=zeros(SoBacTdo,nBuocTg); 
Giatoc=zeros(SoBacTdo,nBuocTg); 
Thgian=zeros(1,nBuocTg); 
t=0; 
for i=2:nBuocTg 
 t=t+BuocTg; 
 Thgian(i)=t; 
 Luc=F0*sin(2*pi*Tanso*t); 
 Fnut(475)=Luc; 
% Cach 1 
% FM1=A0*TichMaTrKhLg_VectoCvi(SM,Cvi(:,i-1)); 
% FM2=A2*TichMaTrKhLg_VectoCvi(SM,Vtoc(:,i-1)); 
% FM3=A3*TichMaTrKhLg_VectoCvi(SM,Giatoc(:,i-1)); 
% FC1=A1*TichMaTrKhLg_VectoCvi(SC,Cvi(:,i-1)); 
% FC2=A4*TichMaTrKhLg_VectoCvi(SC,Vtoc(:,i-1)); 
% FC3=A5*TichMaTrKhLg_VectoCvi(SC,Giatoc(:,i-1)); 
% FnutTdt=Fnut+FM1+FM2+FM3+FC1+FC2+FC3; 
% Cvi(:,i)=Decompostion(K_hqua,FnutTdt); 
% Vtoc(:,i)=A1*(Cvi(:,i)-Cvi(:,i-1))-A4*Vtoc(:,i-1)-A5*Giatoc(:,i-1); 
% Giatoc(:,i)=A0*(Cvi(:,i)-Cvi(:,i-1))-A2*Vtoc(:,i-1)-A3*Giatoc(:,i-1); 
% CvDauCoc(i)=DauCOC(Cvi(:,i)); 
% VtDauCoc(i)=DauCOC(Vtoc(:,i)); 
% GtDauCoc(i)=DauCOC(Giatoc(:,i)); 
% Cach 2 
 QT=TichMaTrKhLg_VectoCvi(SK,Cvi(:,i-1)); 
 166 
 FM2=A2*TichMaTrKhLg_VectoCvi(SM,Vtoc(:,i-1)); 
 FM3=A3*TichMaTrKhLg_VectoCvi(SM,Giatoc(:,i-1)); 
 FC2=A4*TichMaTrKhLg_VectoCvi(SC,Vtoc(:,i-1)); 
 FC3=A5*TichMaTrKhLg_VectoCvi(SC,Giatoc(:,i-1)); 
 FnutTdt=Fnut-QT+FM2+FM3+FC2+FC3; 
 deltaU=Decompostion(K_hqua,FnutTdt); 
 Cvi(:,i)=Cvi(:,i-1)+deltaU; 
 Vtoc(:,i)=A1*(Cvi(:,i)-Cvi(:,i-1))-A4*Vtoc(:,i-1)-A5*Giatoc(:,i-1); 
 Giatoc(:,i)=A0*(Cvi(:,i)-Cvi(:,i-1))-A2*Vtoc(:,i-1)-A3*Giatoc(:,i-1); 
end 
Return 
%------------------------Newmark_NewtonRaphson.m------------------------- 
%============================FOR FRAME_W2_2012=========================== 
%To solve nonlinear equations M*U2+C*U1+K*U=P*F(t) 
%Using NewtonRaphson iteration method & 
%Newmark direct integration 
function [NodeDisplacement,NodeVelocity,NodeAcceleration]=... 
 Newmark_NewtonRaphson(TimeFunc,ETol,... 
 Write2Disk,TempDirectory,FName2Write) 
PublicDeclarations; 
Ut=T2L_SL.NodalDisplacement; 
Vt=T2L_SL.NodalVelocity; 
At=T2L_SL.NodalAcceleration; 
abc=T2L_SL.Index; 
nTime=length(TimeFunc); 
dt=TimeFunc(1,2)-TimeFunc(1,1); 
dU=Ut*0; 
Ft=T2L_SL.Stiffness*Ut; 
flag=1; 
%Integration constant 
c1=4/dt; 
c2=4/dt^2; 
c3=2/dt; 
switch Write2Disk 
 case 0 
 NodeDisplacement=Ut; 
 NodeVelocity=Vt; 
 NodeAcceleration=At; 
 case 1 
 NodeDisplacement=[]; 
 NodeVelocity=[]; 
 NodeAcceleration=[]; 
end 
for tt=1:nTime 
 167 
 if Write2Disk==1 
 outfile=[TempDirectory num2str(tt-1) FName2Write]; 
 NodalDisplacement=Ut; 
 save(outfile,'-mat','NodalDisplacement'); 
 end 
 Utdt=Ut; 
T2L_SL.Damping=T2L_SL.Mass*AlphaR+T2L_SL.Stiffness*BetaR; 
EffectStiff=T2L_SL.Stiffness+c3*T2L_SL.Damping+c2*T2L_SL.Mass; 
Rtdt=T2L_SL.NodalLoadStatic+T2L_SL.NodalLoadDynamic*TimeFunc(2,tt); 
 ii=1 
 while 1 
 Ftdt=T2L_SL.Stiffness*Utdt; 
 Temp1=Rtdt-Ftdt-... 
 T2L_SL.Damping*(c3*(Utdt-Ut)-Vt)-... 
 T2L_SL.Mass*(c2*(Utdt-Ut)-c1*Vt-At); 
 dU(abc)=EffectStiff(abc,abc)\Temp1(abc); 
 Atdt=c2*(Utdt-Ut+dU)-c1*Vt-At; 
 Vtdt=Vt+(dt/2)*(At+Atdt); 
 Utdt=Ut+(dt/2)*(Vt+Vtdt); 
 Temp1=dU(abc).'*... 
 (Rtdt(abc)-Ftdt(abc)-T2L_SL.Mass(abc,abc)*Atdt(abc)); 
 Temp2=(Utdt(abc)-Ut(abc)).'*... 
 (Rtdt(abc)-Ft(abc)-T2L_SL.Mass(abc,abc)*At(abc)); 
 if abs(Temp1/Temp2)<=ETol 
 break; 
 else 
 ii=ii+1 
 end 
 end 
 Ut=Utdt; 
 Vt=Vtdt; 
 At=Atdt; 
 Ft=Ftdt; 
 if Write2Disk==0 
 NodeDisplacement=[NodeDisplacement Ut]; 
 NodeVelocity=[NodeVelocity Vt]; 
 NodeAcceleration=[NodeAcceleration At]; 
 end 
%Lap tren phan tu tiep xuc-Kiem tra tach, truot cuc bo 
 for ii=1:T2L_SL.NumberOfSlip2D 
 ProgramBd2; 
 if xicmaz<=0 
 knz=knz*1E-4;knx=knx*1E-4;kny=knx*1E-4; 
 CalculateElementParametersSlip2D; 
 168 
 T2L_SL.Stiffness=T2L_SL.Stiffness; 
 end 
 if xicmaz>0 
 if tauxy<=taugh 
 T2L_SL.Stiffness=T2L_SL.Stiffness; 
 if tauxy>taugh 
 knz=knz;knx=knx*1E-4;kny=knx*1E-4; 
 CalculateElementParametersSlip2D; 
 T2L_SL.Stiffness=T2L_SL.Stiffness; 
 end 
 end 
 end 
 end 
return