Luận án Nghiên cứu chế tạo và tính toán mô phỏng màng mỏng với nền cacbon giống kim cương có tính chống mòn cao và ma sát thấp

BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC SƯ PHẠM KỸ THUẬT 
THÀNH PHỐ HỒ CHÍ MINH 
LÊ VĂN AN 
NGHIÊN CỨU CHẾ TẠO VÀ TÍNH TOÁN 
MÔ PHỎNG MÀNG MỎNG VỚI NỀN CACBON 
GIỐNG KIM CƯƠNG CÓ TÍNH 
CHỐNG MÒN CAO VÀ MA SÁT THẤP 
LUẬN ÁN TIẾN SĨ 
NGÀNH: KỸ THUẬT CƠ KHÍ 
TP. Hồ Chí Minh, tháng 7/2020 
 BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC SƯ PHẠM KỸ THUẬT 
THÀNH PHỐ HỒ CHÍ MINH 
LÊ VĂN AN 
NGHIÊN CỨU CHẾ TẠO VÀ TÍNH TOÁN 
MÔ PHỎNG MÀNG MỎNG VỚI NỀN CACBON 
GIỐNG KIM CƯƠNG CÓ TÍNH 
CHỐNG MÒN CAO VÀ MA SÁT THẤP 
NGÀNH: KỸ THUẬT CƠ KHÍ – 62520103 
Hướng dẫn khoa học: 
Phản biện 1: ................................................................................................ 
Phản biện 2: ................................................................................................ 
Phản biện 3: ................................................................................................ 
ii 
LÝ LỊCH CÁ NHÂN 
I LÝ LỊCH SƠ LƯỢC 
Họ và tên : LÊ VĂN AN Giới tính: Nam 
Này, tháng, năm sinh : 07/3/1978 Nơi sinh: Nghệ An 
Quê quán : Nghệ An Dân tộc: Kinh 
Chỗ ở riêng hoặc địa chỉ liên hệ: P2-02.12, Vinhomes Central Park, 720A Điện Biên 
Phủ, P.22, Q. Bình Thạnh, TP. HCM 
Điện thoại cơ quan: 
Điện thoại nhà riêng: 0989.235.449 
E-mail: anlevan78@gmail.com 
II QUÁ TRÌNH ĐÀO TẠO 
1. Trung học chuyên nghiệp: 
- Hệ đào tạo: Thời gian đào tạo từ ../ đến ../ 
- Nơi học (trường, thành phố): 
- Ngành học: 
2. Đại học 
- Hệ đào tạo: chính quy Thời gian đào tạo từ 10/1999 đến 03/2003 
- Nơi học (trường, thành phố): Trường Đại học Sư phạm Kỹ thuật TP. HCM 
- Ngành học: Chế tạo máy 
- Tên đồ án, luận văn hoặc môn thi tốt nghiệp: Khai thác và sử dụng phần mềm 
lập trình MTS 
3. Cao học 
- Hệ đào tạo: chính quy Thời gian đào tạo từ 10/2004 đến 10/2006 
- Nơi học (trường, thành phố): Trường Đại học Sư phạm Kỹ thuật TP. HCM 
- Ngành học: Công nghệ Chế tạo máy 
- Tên đồ án, luận văn hoặc môn thi tốt nghiệp: Nghiên cứu các loại màng chống 
mòn có cấu trúc nano được phủ bằng phương pháp plasma 
iii 
III QUÁ TRÌNH CÔNG TÁC CHUYÊN MÔN KỂ TỪ KHI TỐT NGHIỆP ĐẠI 
HỌC 
Thời gian Nơi công tác Công việc đảm nhận 
2003 – 2006 Trường ĐH GTVT TP. HCM Giảng viên 
2006 – 2019 Trường ĐH GTVT TP. HCM Trưởng Bộ môn 
2019 – nay Trường ĐH GTVT TP. HCM Giảng viên 
IV CÁC ĐỀ TÀI, DỰ ÁN, NHIỆM VỤ KHÁC ĐÃ CHỦ TRÌ HOẶC THAM GIA 
Tên đề tài, dự án, nhiệm vụ 
Thời gian 
(bắt đầu – 
kết thúc 
Thuộc 
chương trình 
(nếu có) 
Tình trạng đề tài 
(đã nghiệm thu, 
chưa nghiệm thu) 
Xây dựng bài giảng điện tử môn 
học Công nghệ gia công KTS – 
chủ trì 
01/2008-
05/2009 
 Đã nghiệm thu 
Thiết kế và mô phỏng bộ truyền 
bánh răng vi sai phẳng ba cấp – 
chủ trì 
10/2008-
09/2009 
 Đã nghiệm thu 
 TP. HCM, ngày 19 tháng 7 năm 2020 
Nghiên cứu sinh 
Lê Văn An 
iv 
LỜI CAM ĐOAN 
 Tôi cam đoan đây là công trình nghiên cứu của tôi. 
 Các số liệu, kết quả nêu 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 khác. Các hình và thông tin bắt nguồn từ nghiên cứu 
khác, tôi cam kết đã trích dẫn đầy đủ trong luận án. 
 TP. HCM, ngày 19 tháng 7 năm 2020 
Nghiên cứu sinh 
Lê Văn An 
v 
LỜI CẢM ƠN 
 Tôi xin gửi lời cảm ơn chân thành đến cán bộ hướng dẫn, PGS. TS. Bùi Xuân 
Lâm và PGS. TS. Lê Hiếu Giang, các thầy đã tận tình giúp đỡ tôi trong suốt thời gian 
chuẩn bị, tiến hành và hoàn thành luận án này. 
 Tôi xin chân thành cảm ơn GS. TS. Sam Zhang, nguyên giáo sư Trường Kỹ 
thuật Cơ khí và Hàng không, ĐH Công nghệ Nanyang, Singapore vì sự giúp đỡ tận 
tình khi tôi tiến hành các thí nghiệm tại đây. 
 Tôi xin chân thành cảm ơn các thầy, cô đã tận tâm giảng dạy trong suốt quá 
trình tôi theo học NCS tại Khoa Cơ khí Chế tạo máy, Trường Đại học Sư phạm Kỹ 
thuật Tp. HCM. 
 Tôi xin chân thành cảm ơn Khoa Cơ khí Chế tạo máy, Phòng Đào tạo – Bộ phận 
Sau Đại học và Ban Giám hiệu Trường Đại học Sư phạm Kỹ thuật TP. HCM đã quan 
tâm, giúp đỡ và tạo điều kiện tốt nhất để tôi hoàn thành luận án này. 
 Tôi xin chân thành cảm ơn bố, mẹ, các anh, em đã động viên tôi hoàn thành tốt 
việc học tập. 
 Cuối cùng, tôi xin cảm ơn vợ và các con, những người luôn theo sát, động viên 
tôi trong quá trình thực hiện luận án này. 
vi 
BM13.1-ĐT-BVCS 
BỘ GIÁO DỤC & ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC SƯ PHẠM KỸ THUẬT 
THÀNH PHỐ HỒ CHÍ MINH 
CỘNG HÒA XÃ HỘI CHỦ NGHĨA VIỆT NAM 
Độc lập – Tự do – Hạnh phúc 
TÓM TẮT NHỮNG ĐÓNG GÓP MỚI CỦA LUẬN ÁN 
Họ & tên NCS : Lê Văn An MS NCS: 13252010302 
Thuộc chuyên ngành : Kỹ thuật cơ khí Khoá: 2013-2016 
Tên luận án : Nghiên cứu chế tạo và tính toán mô phỏng màng mỏng 
với nền cacbon giống kim cương có tính chống mòn cao và ma sát thấp 
Người hướng dẫn chính : PGS. TS. Bùi Xuân Lâm 
Người hướng dẫn phụ : PGS. TS. Lê Hiếu Giang 
 Tóm tắt những đóng góp mới về lý luận và học thuật của luận án: 
 - Việc nghiên cứu, chế tạo màng mỏng có nền cacbon giống kim cương (DLC) 
bằng phương pháp kết hợp phún xạ magnetron và hóa học đã được thực hiện. Phương 
pháp kết hợp này cho tốc độ phủ màng cao hơn nhiều (gấp 2 lần) so với phún xạ 
magnetron thuần túy. 
 - Thông qua việc thiết kế các thí nghiệm, tác giả đã có được các thông số phù 
hợp để phủ màng có cấu trúc nanocomposite với các hạt tinh thể nano TiC (kích thước 
từ 2-4 nm) phân tán trong nền vô định hình DLC. Cấu trúc tế vi của màng dễ dàng 
thay đổi thông qua việc thay đổi mật độ năng lượng của bia Ti. 
 - Thành phần, cấu trúc tế vi, cơ tính và tính chất ma sát của màng được khảo sát 
chi tiết bằng các thiết bị chẩn đoán chuyên biệt dùng cho vật liệu có cấu trúc nano. 
Các màng nanocomposite có độ cứng cao (20-30 GPa), ứng suất nội thấp (1-2 GPa), 
bám dính lên bề mặt chi tiết phủ rất tốt và hệ số ma sát trong điều kiện không bôi trơn 
khi trượt với thép rất nhỏ (0,17-0,19). Các thông số này cho thấy màng 
nanocomposite có nền DLC được chế tạp bằng phương pháp mới vượt trội so với các 
màng chống mòn thông thường đang được sử dụng như carbide hoặc nitride. 
vii 
 - Trên cơ sở lý thuyết về đo độ cứng nano, lý thuyết biến dạng và phương pháp 
PTHH, mô hình toán và lưu đồ giải thuật của màng mỏng chịu tải bởi đầu đâm nano 
đã được thiết lập. Chương trình tính toán được lập trình bằng Matlab. Các màng 
nanocomposite có nền DLC được chế tạo trong phần thực nghiệm đã được sử dụng 
làm đối tượng nghiên cứu để khảo sát độ cứng, phân bố ứng suất, biến dạng và hư 
hỏng khi chịu tải bởi đầu đâm nano. Độ chính xác của tính toán và mô phỏng được 
đánh giá dựa vào sai số theo chuẩn năng lượng và tốc độ hội tụ. Kết quả có được từ 
luận án là cơ sở để phát triển các nghiên cứu chuyên sâu về khả năng chịu tải và hư 
hỏng của màng mỏng ở các chế độ làm việc phức tạp. 
 Tp. Hồ Chí Minh, ngày 07 tháng 5 năm 2020 
Nghiên cứu sinh 
(Ký và ghi rõ họ tên) 
Lê Văn An 
Người hướng dẫn chính 
(Ký và ghi rõ họ tên) 
Người hướng dẫn phụ 
(Ký và ghi rõ họ tên) 
viii 
BM13.2-ĐT-BVCS 
MINISTRY OF EDUCATION AND TRAINING 
HO CHI MINH CITY UNIVERSITY OF 
TECHNOLOGY AND EDUCATION 
SOCIALIST REPUBLIC OF VIETNAM 
Independence – Freedom - Happiness 
SUMMARY OF CONTRIBUTIONS OF THE THESIS 
PhD candidate : Le Van An Candidate code: 13252010302 
Major : Mechanical Engineering Intake: 2013-2016 
Thesis title : Deposition, calculation and simulation of DLC-based thin films 
with good wear resistance and low friction 
Supervisor : Assoc. Prof. Dr. Bui Xuan Lam 
Co-supervisor : Assoc. Prof. Dr. Le Hieu Giang 
 Summary of theoretical and academic contribution of the thesis: 
 - The research and producing of diamond-like carbon (DLC) based thin films 
deposited via a combination method of magnetron sputtering and chemical vapor 
deposition were carried out. This deposition method gives much higher depostion 
rates compared to that of magnetron sputtering 
 - Through experiment design, the deposition parematers, which result in 
nanocomposite structure with TiC nanocrysrtals (2-4 nm) embedded in the 
amorphous matrix of DLC, were determined. The microstructure of the thin films can 
be modified easily via changing the power density of Ti target. 
 - Composition, microstructure, mechanical and tribological properties of the 
deposited thin films were investigated using advanced characterization techniques. 
Nanocomposite thin films have high hardness (20-30 GPa), low residual stress (1-2 
GPa), very good adhesion to the substrate, and low coefficient of friction (0.17-0.19) 
when sliding against steel counterpart. These excellent properties prove that the DLC 
nanocomposite thin films deposited via this new method are much better than carbide 
or nitride thin films, which are being used in tribological applications. 
 - Basing on the theory of nanoindentation, theory of deformation, and finite 
ix 
element analysis, the mathematical model and algorithm of a thin film loaded by the 
nanoindenter were built. The programming was carried out with Matlab code. 
Hardness, stresses distribution, deformation and failure of DLC based nanocomposite 
thin films deposited in the experiment part were investigated. The accuracy of the 
calculation and simulation was estimated using the energy norm and the rate of 
convergence. The results of this research can be used to develop studies on load-
bearing capacity and failure of thin films under complicated working conditions. 
 Ho Chi Minh City, 07 May 2020 
Candidate 
Le Van An 
Supervisor Co-supervisor 
x 
MỤC LỤC 
LÝ LỊCH CÁ NHÂN ................................................................................................. ii 
LỜI CAM ĐOAN ..................................................................................................... iv 
LỜI CẢM ƠN ............................................................................................................. v 
MỤC LỤC ................................................................................................................... x 
DANH SÁCH CÁC CHỮ VIẾT TẮT .................................................................... xiii 
DANH MỤC CÁC ĐỒ THỊ, HÌNH VẼ ................................................................. xvi 
DANH MỤC CÁC BẢNG ...................................................................................... xix 
MỞ ĐẦU ..................................................................................................................... 1 
1. Lý do chọn đề tài ................................................................................................. 1 
2. Mục đích nghiên cứu ........................................................................................... 2 
3. Nhiệm vụ nghiên cứu .......................................................................................... 2 
4. Phạm vi nghiên cứu ............................................................................................. 3 
5. Phương pháp nghiên cứu ..................................................................................... 3 
6. Ý nghĩa khoa học và thực tiễn ............................................................................. 3 
7. Cấu trúc của luận án ............................................................................................ 3 
 .................................................................................................................... 5 
TỔNG QUAN VỀ MÀNG MỎNG VÀ CÁC PHƯƠNG PHÁP CHẾ TẠO ............. 5 
1.1. Tổng quan về màng mỏng ................................................................................ 5 
1.2. Các phương pháp chế tạo màng mỏng ............................................................. 7 
1.2.1. Phủ hoá học (CVD) ................................................................................... 7 
1.2.2. Phủ vật lý (PVD) ..................................................................................... 10 
1.3. Cơ chế hình thành và phát triển của màng phủ bằng phương pháp phún xạ . 13 
1.4. Ứng suất dư của màng mỏng .......................................................................... 15 
1.5. Các kết quả đạt được từ các công trình đã công bố ........................................ 16 
1.6. Các vấn đề cần được tiếp tục nghiên cứu ....................................................... 17 
1.7. Kết luận rút ra từ tổng quan ............................................................................ 17 
 .................................................................................................................. 18 
xi 
CHẾ TẠO MÀNG MỎNG CÓ NỀN CACBON GIỐNG KIM CƯƠNG BẰNG 
PHƯƠNG PHÁP PHÚN XẠ MAGNETRON KẾT HỢP HÓA HỌC .................... 18 
2.1. Mô tả hệ thống phún xạ magnetron ................................................................ 18 
2.2. Màng mỏng cacbon giống kim cương ............................................................ 20 
2.3. Chế tạo màng mỏng có nền DLC ................................................................... 24 
2.3.1. Thành phần của màng .............................................................................. 24 
2.3.2. Cấu trúc của màng ................................................................................... 27 
2.3.3. Hình ảnh bề mặt của màng ...................................................................... 29 
2.4. Kết luận chương 2 .......................................................................................... 30 
 .................................................................................................................. 32 
CÁC TÍNH CHẤT CƠ HỌC VÀ MA SÁT CỦA MÀNG MỎNG CÓ NỀN 
CACBON GIỐNG KIM CƯƠNG PHỦ BẰNG PHƯƠNG PHÁP PHÚN XẠ 
MAGNETRON KẾT HỢP HÓA HỌC .................................................................... 32 
3.1. Ứng suất dư và tính bám dính ........................................................................ 32 
3.2. Độ cứng và mô đun đàn hồi của màng ........................................................... 35 
3.3. Tính chất ma sát và mòn của màng ................................................................ 36 
3.4. Kết luận chương 3 .......................................................................................... 40 
 .................................................................................................................. 41 
TÍNH TOÁN MÔ PHỎNG ỨNG XỬ TĨNH CỦA MÀNG NANOCOMPOSITE CÓ 
NỀN CACBON GIỐNG KIM CƯƠNG ................................................................... 41 
4.1. Mở đầu ............................................................................................................ 41 
4.2. Phương pháp đo độ cứng tiếp xúc .................................................................. 42 
4.3. Cơ sở phương pháp phần tử hữu hạn (PTHH) xác định tính chất màng chịu tải 
trọng tĩnh................................................................................................................ 46 
4.3.1. Phương trình PTHH .......................................................... ...  sd; 
 t(6) = ( 729.0 / 4.0 ) * ra * rc * rd * sa * sc * sd; 
 t(7) = - ( 729.0 / 4.0 ) * ra * rb * rd * sa * sc * sd; 
41 
 t(8) = ( 243.0 / 4.0 ) * ra * rb * rc * sa * sc * sd; 
 t(9) = ( 243.0 / 4.0 ) * rb * rc * rd * sa * sb * sd; 
 t(10) = - ( 729.0 / 4.0 ) * ra * rc * rd * sa * sb * sd; 
 t(11) = ( 729.0 / 4.0 ) * ra * rb * rd * sa * sb * sd; 
 t(12) = - ( 243.0 / 4.0 ) * ra * rb * rc * sa * sb * sd; 
 t(13) = - ( 81.0 / 4.0 ) * rb * rc * rd * sa * sb * sc; 
 t(14) = ( 243.0 / 4.0 ) * ra * rc * rd * sa * sb * sc; 
 t(15) = - ( 243.0 / 4.0 ) * ra * rb * rd * sa * sb * sc; 
 t(16) = ( 81.0 / 4.0 ) * ra * rb * rc * sa * sb * sc; 
 dbcd = 3.0 * r * r - 4.0 * r + 11.0 / 9.0; 
 dacd = 3.0 * r * r - 10.0 * r / 3.0 + 2.0 / 3.0; 
 dabd = 3.0 * r * r - 8.0 * r / 3.0 + 1.0 / 3.0; 
 dabc = 3.0 * r * r - 2.0 * r + 2.0 / 9.0; 
 dtdr(1) = ( 81.0 / 4.0 ) * dbcd * sb * sc * sd; 
 dtdr(2) = - ( 243.0 / 4.0 ) * dacd * sb * sc * sd; 
 dtdr(3) = ( 243.0 / 4.0 ) * dabd * sb * sc * sd; 
 dtdr(4) = - ( 81.0 / 4.0 ) * dabc * sb * sc * sd; 
 dtdr(5) = - ( 243.0 / 4.0 ) * dbcd * sa * sc * sd; 
 dtdr(6) = ( 729.0 / 4.0 ) * dacd * sa * sc * sd; 
 dtdr(7) = - ( 729.0 / 4.0 ) * dabd * sa * sc * sd; 
 dtdr(8) = ( 243.0 / 4.0 ) * dabc * sa * sc * sd; 
 dtdr(9) = ( 243.0 / 4.0 ) * dbcd * sa * sb * sd; 
 dtdr(10) = - ( 729.0 / 4.0 ) * dacd * sa * sb * sd; 
 dtdr(11) = ( 729.0 / 4.0 ) * dabd * sa * sb * sd; 
 dtdr(12) = - ( 243.0 / 4.0 ) * dabc * sa * sb * sd; 
 dtdr(13) = - ( 81.0 / 4.0 ) * dbcd * sa * sb * sc; 
 dtdr(14) = ( 243.0 / 4.0 ) * dacd * sa * sb * sc; 
 dtdr(15) = - ( 243.0 / 4.0 ) * dabd * sa * sb * sc; 
 dtdr(16) = ( 81.0 / 4.0 ) * dabc * sa * sb * sc; 
42 
 dbcd = 3.0 * s * s - 4.0 * s + 11.0 / 9.0; 
 dacd = 3.0 * s * s - 10.0 * s / 3.0 + 2.0 / 3.0; 
 dabd = 3.0 * s * s - 8.0 * s / 3.0 + 1.0 / 3.0; 
 dabc = 3.0 * s * s - 2.0 * s + 2.0 / 9.0; 
 dtds(1) = ( 81.0 / 4.0 ) * rb * rc * rd * dbcd; 
 dtds(2) = - ( 243.0 / 4.0 ) * ra * rc * rd * dbcd; 
 dtds(3) = ( 243.0 / 4.0 ) * ra * rb * rd * dbcd; 
 dtds(4) = - ( 81.0 / 4.0 ) * ra * rb * rc * dbcd; 
 dtds(5) = - ( 243.0 / 4.0 ) * rb * rc * rd * dacd; 
 dtds(6) = ( 729.0 / 4.0 ) * ra * rc * rd * dacd; 
 dtds(7) = - ( 729.0 / 4.0 ) * ra * rb * rd * dacd; 
 dtds(8) = ( 243.0 / 4.0 ) * ra * rb * rc * dacd; 
 dtds(9) = ( 243.0 / 4.0 ) * rb * rc * rd * dabd; 
 dtds(10) = - ( 729.0 / 4.0 ) * ra * rc * rd * dabd; 
 dtds(11) = ( 729.0 / 4.0 ) * ra * rb * rd * dabd; 
 dtds(12) = - ( 243.0 / 4.0 ) * ra * rb * rc * dabd; 
 dtds(13) = - ( 81.0 / 4.0 ) * rb * rc * rd * dabc; 
 dtds(14) = ( 243.0 / 4.0 ) * ra * rc * rd * dabc; 
 dtds(15) = - ( 243.0 / 4.0 ) * ra * rb * rd * dabc; 
 dtds(16) = ( 81.0 / 4.0 ) * ra * rb * rc * dabc; 
 return 
end 
function node=square_node_array(pt1,pt2,pt3,pt4 ,nnx,nny) 
if ( nargin < 6 ) 
 disp('Not enough parameters specified for quare_node_array function') 
end 
% get node spacing along u direction 
 xi_pts=linspace(-1,1,nnx); 
% get node spacing along v direction 
43 
 eta_pts=linspace(-1,1,nny); 
x_pts=[pt1(1),pt2(1),pt3(1),pt4(1)]; 
y_pts=[pt1(2),pt2(2),pt3(2),pt4(2)]; 
for r=1:nny 
 eta=eta_pts(r); 
 for c=1:nnx 
 xi=xi_pts(c); 
 % get interpolation basis at xi, eta 
 N=shape_func('Q4',[xi,eta]); 
 N=N(:,1); 
 node((r-1)*nnx+c,:)=[x_pts*N,y_pts*N]; 
 end 
end 
end 
function [node,element]=structured_q8_mesh(pt1,pt2,pt3,pt4,numx,numy) 
%Forms the element and node matrices for Q8 element, 
%with the element and node matrices arranged in a nodal 
%counterclockwise order (primary nodes then secondary nodes) 
nnx=numx+1; 
nny=numy+1; 
node=square_node_array(pt1,pt2,pt3,pt4,nnx,nny); 
inc_u=1; 
inc_v=nnx; 
node_pattern=[ 1 2 nnx+2 nnx+1 ]; 
element=make_elem(node_pattern,numx,numy,inc_u,inc_v); 
[element,node]=q4totq8(element,node,numx,numy); 
end 
function [node,element]=structured_q9_mesh(pt1,pt2,pt3,pt4,numx,numy) 
%Forms the element and node matrices for Q9 element, with the element 
44 
%and node matrices arranged in a nodal counterclockwise order 
%(primary nodes-> secondary nodes then internal nodes) 
nnx=numx+1; 
nny=numy+1; 
node=square_node_array(pt1,pt2,pt3,pt4,nnx,nny); 
inc_u=1; 
inc_v=nnx; 
node_pattern=[ 1 2 nnx+2 nnx+1 ]; 
element=make_elem(node_pattern,numx,numy,inc_u,inc_v); 
[element,node]=q4totq9(element,node,numx,numy); 
end 
function element_node = grid_q12_element ( nelemx, nelemy ) 
 element = 0; 
 for j = 1 : nelemy 
 for i = 1 : nelemx 
 base = ( j - 1 ) * ( 5 * nelemx + 3 ) + 1; 
 element = element + 1; 
 element_node( 1,element) = base + ( i - 1 ) * 3; 
 element_node( 2,element) = base + ( i - 1 ) * 3 + 1; 
 element_node( 3,element) = base + ( i - 1 ) * 3 + 2; 
 element_node( 4,element) = base + ( i - 1 ) * 3 + 3; 
 element_node( 5,element) = base + 3 * nelemx + i; 
 element_node( 6,element) = base + 3 * nelemx + i + 1; 
 element_node( 7,element) = base + 4 * nelemx + i + 1; 
 element_node( 8,element) = base + 4 * nelemx + i + 2; 
 element_node( 9,element) = base + 5 * nelemx + 3 * i; 
 element_node(10,element) = base + 5 * nelemx + 3 * i + 1; 
 element_node(11,element) = base + 5 * nelemx + 3 * i + 2; 
45 
 element_node(12,element) = base + 5 * nelemx + 3 * i + 3; 
 end 
 end 
 return 
end 
 element = 0; 
 for j = 1 : nelemy 
 for i = 1 : nelemx 
 base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2; 
 element = element + 1; 
 element_node( 1,element) = base; 
 element_node( 2,element) = base + 1; 
 element_node( 3,element) = base + 2; 
 element_node( 4,element) = base + 3; 
 element_node( 5,element) = base + ( 3 * nelemx + 1 ); 
 element_node( 6,element) = base + ( 3 * nelemx + 1 ) + 1; 
 element_node( 7,element) = base + ( 3 * nelemx + 1 ) + 2; 
 element_node( 8,element) = base + ( 3 * nelemx + 1 ) + 3; 
 element_node( 9,element) = base + 2 * ( 3 * nelemx + 1 ); 
 element_node(10,element) = base + 2 * ( 3 * nelemx + 1 ) + 1; 
 element_node(11,element) = base + 2 * ( 3 * nelemx + 1 ) + 2; 
 element_node(12,element) = base + 2 * ( 3 * nelemx + 1 ) + 3; 
 element_node(13,element) = base + 3 * ( 3 * nelemx + 1 ); 
 element_node(14,element) = base + 3 * ( 3 * nelemx + 1 ) + 1; 
 element_node(15,element) = base + 3 * ( 3 * nelemx + 1 ) + 2; 
 element_node(16,element) = base + 3 * ( 3 * nelemx + 1 ) + 3; 
 end 
 end 
 return 
46 
end 
INTERGRATION 
function Bfem =Bmatrix(pt,elemType,iel) 
% Gives the strain displacement matrix (B matrix of size 3x8)of each element 
global node element 
sctr = element(iel,:); 
nn = length(sctr); 
[N,dNdxi] = shape_func(elemType,pt); % element shape functions 
J0 = node(sctr,:)'*dNdxi; % element Jacobian matrix 
invJ0 = inv(J0); 
dNdx = dNdxi*invJ0; % derivatives of N w.r.t XY 
Gpt = N'*node(sctr,:); % GP in global coord, used 
 Bfem = zeros(3,2*nn); 
 Bfem(1,1:2:2*nn) = dNdx(:,1)' ; 
 Bfem(2,2:2:2*nn) = dNdx(:,2)' ; 
 Bfem(3,1:2:2*nn) = dNdx(:,2)' ; 
 Bfem(3,2:2:2*nn) = dNdx(:,1)' ; 
end % end of function 
function Bfem4 =Bmatrix4(pt,elemType,iel) 
% Gives the strain displacement matrix (B matrix of size 4x8)of each element 
global node element 
sctr = element(iel,:); 
nn = length(sctr); 
[N,dNdxi] = shape_func(elemType,pt); % element shape functions 
J0 = node(sctr,:)'*dNdxi; % element Jacobian matrix 
invJ0 = inv(J0); 
dNdx = dNdxi*invJ0; % derivatives of N w.r.t XY 
Gpt = N'*node(sctr,1); % GP in global coord, used 
 Bfem4 = zeros(4,2*nn); 
47 
 Bfem4(1,1:2:2*nn) = dNdx(:,1)' ; 
 Bfem4(2,2:2:2*nn) = dNdx(:,2)' ; 
 Bfem4(3,1:2:2*nn) = dNdx(:,2)' ; 
 Bfem4(3,2:2:2*nn) = dNdx(:,1)' ; 
end % end of function 
function [W,Q] = gauss_pt_wt( quadorder, qt, sdim ) 
% Returns the weights and coordinates of the gauss integration points 
 if ( nargin < 3 ) % set default arguments 
 if ( strcmp(qt,'GAUSS') == 1 ) 
 dim = 1; 
 else 
 dim = 2; 
 end 
 end 
 if ( nargin < 2 ) 
 type = 'GAUSS'; 
 end 
 if ( strcmp(qt,'GAUSS') == 1 ) 
 quadpoint=zeros(quadorder^sdim ,sdim); 
 quadweight=zeros(quadorder^sdim,1); 
 r1pt=zeros(quadorder,1); r1wt=zeros(quadorder,1); 
 switch ( quadorder ) 
 case 1 
 r1pt(1) = 0.000000000000000; 
 r1wt(1) = 2.000000000000000; 
 case 2 
 r1pt(1) = 0.577350269189626; 
 r1pt(2) =-0.577350269189626; 
 r1wt(1) = 1.000000000000000; 
48 
 r1wt(2) = 1.000000000000000; 
 case 3 
 r1pt(1) = 0.774596669241483; 
 r1pt(2) =-0.774596669241483; 
 r1pt(3) = 0.000000000000000; 
 r1wt(1) = 0.555555555555556; 
 r1wt(2) = 0.555555555555556; 
 r1wt(3) = 0.888888888888889; 
 case 4 
 r1pt(1) = 0.861134311594053; 
 r1pt(2) =-0.861134311594053; 
 r1pt(3) = 0.339981043584856; 
 r1pt(4) =-0.339981043584856; 
 r1wt(1) = 0.347854845137454; 
 r1wt(2) = 0.347854845137454; 
 r1wt(3) = 0.652145154862546; 
 r1wt(4) = 0.652145154862546; 
 case 5 
 r1pt(1) = 0.906179845938664; 
 r1pt(2) =-0.906179845938664; 
 r1pt(3) = 0.538469310105683; 
 r1pt(4) =-0.538469310105683; 
 r1pt(5) = 0.000000000000000; 
 r1wt(1) = 0.236926885056189; 
 r1wt(2) = 0.236926885056189; 
 r1wt(3) = 0.478628670499366; 
 r1wt(4) = 0.478628670499366; 
 r1wt(5) = 0.568888888888889; 
49 
 case 6 
 r1pt(1) = 0.932469514203152; 
 r1pt(2) =-0.932469514203152; 
 r1pt(3) = 0.661209386466265; 
 r1pt(4) =-0.661209386466265; 
 r1pt(5) = 0.238619186003152; 
 r1pt(6) =-0.238619186003152; 
 r1wt(1) = 0.171324492379170; 
 r1wt(2) = 0.171324492379170; 
 r1wt(3) = 0.360761573048139; 
 r1wt(4) = 0.360761573048139; 
 r1wt(5) = 0.467913934572691; 
 r1wt(6) = 0.467913934572691; 
 case 7 
 r1pt(1) = 0.949107912342759; 
 r1pt(2) = -0.949107912342759; 
 r1pt(3) = 0.741531185599394; 
 r1pt(4) = -0.741531185599394; 
 r1pt(5) = 0.405845151377397; 
 r1pt(6) = -0.405845151377397; 
 r1pt(7) = 0.000000000000000; 
 r1wt(1) = 0.129484966168870; 
 r1wt(2) = 0.129484966168870; 
 r1wt(3) = 0.279705391489277; 
 r1wt(4) = 0.279705391489277; 
 r1wt(5) = 0.381830050505119; 
 r1wt(6) = 0.381830050505119; 
 r1wt(7) = 0.417959183673469; 
50 
 otherwise 
 disp('unsupported integration order') 
 end % end of quadorder switch 
 n=1; 
 if ( sdim == 1 ) 
 for i = 1:quadorder 
 quadpoint(n,:) = [ r1pt(i) ]; 
 quadweight(n) = r1wt(i); 
 n = n+1; 
 end 
 elseif ( sdim == 2 ) 
 for i = 1:quadorder 
 for j = 1:quadorder 
 quadpoint(n,:) = [ r1pt(i), r1pt(j)]; 
 quadweight(n) = r1wt(i)*r1wt(j); 
 n = n+1; 
 end 
 end 
 end 
 Q=quadpoint; 
 W=quadweight; 
 % END OF GAUSSIAN QUADRATURE DEFINITION FOR RECTANGULAR 
ELEMENTS 
 elseif ( strcmp(qt,'TRIANGULAR') == 1) 
 if ( quadorder > 7 ) % check for valid quadrature order 
 disp('Quadrature order too high for triangular quadrature'); 
 quadorder = 1; 
 end 
 if ( quadorder == 1 ) % set quad points and quadweights 
51 
 quadpoint = [ 0.3333333333333, 0.3333333333333 ]; 
 quadweight = 1; 
 elseif ( quadorder == 2 ) 
 quadpoint = zeros( 3, 2 ); 
 quadweight = zeros( 3, 1 ); 
 quadpoint(1,:) = [ 0.1666666666667, 0.1666666666667 ]; 
 quadpoint(2,:) = [ 0.6666666666667, 0.1666666666667 ]; 
 quadpoint(3,:) = [ 0.1666666666667, 0.6666666666667 ]; 
 quadweight(1) = 0.3333333333333; 
 quadweight(2) = 0.3333333333333; 
 quadweight(3) = 0.3333333333333; 
 elseif ( quadorder <= 5 ) 
 quadpoint = zeros( 7, 2 ); 
 quadweight = zeros( 7, 1 ); 
 quadpoint(1,:) = [ 0.1012865073235, 0.1012865073235 ]; 
 quadpoint(2,:) = [ 0.7974269853531, 0.1012865073235 ]; 
 quadpoint(3,:) = [ 0.1012865073235, 0.7974269853531 ]; 
 quadpoint(4,:) = [ 0.4701420641051, 0.0597158717898 ]; 
 quadpoint(5,:) = [ 0.4701420641051, 0.4701420641051 ]; 
 quadpoint(6,:) = [ 0.0597158717898, 0.4701420641051 ]; 
 quadpoint(7,:) = [ 0.3333333333333, 0.3333333333333 ]; 
 quadweight(1) = 0.1259391805448; 
 quadweight(2) = 0.1259391805448; 
 quadweight(3) = 0.1259391805448; 
 quadweight(4) = 0.1323941527885; 
 quadweight(5) = 0.1323941527885; 
 quadweight(6) = 0.1323941527885; 
 quadweight(7) = 0.2250000000000; 
52 
 else 
 quadpoint = zeros( 13, 2 ); 
 quadweight = zeros( 13, 1 ); 
 quadpoint(1 ,:) = [ 0.0651301029022, 0.0651301029022 ]; 
 quadpoint(2 ,:) = [ 0.8697397941956, 0.0651301029022 ]; 
 quadpoint(3 ,:) = [ 0.0651301029022, 0.8697397941956 ]; 
 quadpoint(4 ,:) = [ 0.3128654960049, 0.0486903154253 ]; 
 quadpoint(5 ,:) = [ 0.6384441885698, 0.3128654960049 ]; 
 quadpoint(6 ,:) = [ 0.0486903154253, 0.6384441885698 ]; 
 quadpoint(7 ,:) = [ 0.6384441885698, 0.0486903154253 ]; 
 quadpoint(8 ,:) = [ 0.3128654960049, 0.6384441885698 ]; 
 quadpoint(9 ,:) = [ 0.0486903154253, 0.3128654960049 ]; 
 quadpoint(10,:) = [ 0.2603459660790, 0.2603459660790 ]; 
 quadpoint(11,:) = [ 0.4793080678419, 0.2603459660790 ]; 
 quadpoint(12,:) = [ 0.2603459660790, 0.4793080678419 ]; 
 quadpoint(13,:) = [ 0.3333333333333, 0.3333333333333 ]; 
 quadweight(1 ) = 0.0533472356088; 
 quadweight(2 ) = 0.0533472356088; 
 quadweight(3 ) = 0.0533472356088; 
 quadweight(4 ) = 0.0771137608903; 
 quadweight(5 ) = 0.0771137608903; 
 quadweight(6 ) = 0.0771137608903; 
 quadweight(7 ) = 0.0771137608903; 
 quadweight(8 ) = 0.0771137608903; 
 quadweight(9 ) = 0.0771137608903; 
 quadweight(10) = 0.1756152576332; 
 quadweight(11) = 0.1756152576332; 
 quadweight(12) = 0.1756152576332; 
 quadweight(13) =-0.1495700444677; 
53 
 end 
 Q=quadpoint; 
 W=quadweight/2; 
 end 
 end % end of function 
 function [W,Q] = gauss_rule(iel,elemType,normal_order) 
% Provides the weight vector and gauss point coordinate matrix of 
% the element based on the integration order selected. 
global node element 
sctr = element(iel,:); % element connectivity 
 if ((elemType == 'Q4') & (normal_order <8)) 
 [W,Q] = gauss_pt_wt(normal_order,'GAUSS',2); 
 elseif ((elemType == 'Q8') & (normal_order < 8)) 
 [W,Q] = gauss_pt_wt(normal_order,'GAUSS',2); 
 elseif ((elemType == 'Q9') & (normal_order < 8)) 
 [W,Q] = gauss_pt_wt(normal_order,'GAUSS',2); 
 elseif elemType == 'T3' 
 [W,Q] = gauss_pt_wt(normal_order,'TRIANGULAR',2); 
 elseif elemType == 'T6' 
 [W,Q] = gauss_pt_wt(normal_order,'TRIANGULAR',2); 
 end 
end % end of function 
 function [N,dNdxi]=shape_func(type,coord,di) 
% Gives the shape function and its derivatives with respect to x and y 
 if ( nargin == 2 ) 
 dim=1; 
 end 
 switch type 
 case 'L2' 
54 
 if size(coord,2) < 1 
 disp('Error coordinate needed for the L2 element') 
 else 
 xi=coord(1); 
 N=([1-xi,1+xi]/2)'; 
 dNdxi=[-1;1]/2; 
 end 
 if size(coord,2) < 2 
 disp('Error two coordinates needed for the T3 element') 
 else 
 xi=coord(1); eta=coord(2); 
 N=[1-xi-eta;xi;eta]; 
 dNdxi=[-1,-1;1,0;0,1]; 
 end 
 case 'Q4' 
 if size(coord,2) < 2 
 disp('Error two coordinates needed for the Q4 element') 
 else 
 xi=coord(1); eta=coord(2); 
 N=1/4*[ (1-xi)*(1-eta); 
 (1+xi)*(1-eta); 
 (1+xi)*(1+eta); 
 (1-xi)*(1+eta)]; 
 dNdxi=1/4*[-(1-eta), -(1-xi); 
 1-eta, -(1+xi); 
 1+eta, 1+xi; 
 -(1+eta), 1-xi]; 
 end 
end % end of function
1