Luận án Đánh giá và mô phỏng các hệ số đàn hồi đa tinh thể hỗn độn

 BỘ GIÁO DỤC VÀ ĐÀO TẠO VIỆN HÀN LÂM KHOA HỌC 
VÀ CÔNG NGHỆ VIỆT NAM 
HỌC VIỆN KHOA HỌC VÀ CÔNG NGHỆ 
----------------------------- 
Vương Thị Mỹ Hạnh 
ĐÁNH GIÁ VÀ MÔ PHỎNG CÁC HỆ SỐ ĐÀN HỒI 
ĐA TINH THỂ HỖN ĐỘN 
LUẬN ÁN TIẾN SỸ NGÀNH CƠ HỌC 
Hà Nội-2020 
BỘ GIÁO DỤC VÀ ĐÀO TẠO VIỆN HÀN LÂM KHOA HỌC 
VÀ CÔNG NGHỆ VIỆT NAM 
HỌC VIỆN KHOA HỌC VÀ CÔNG NGHỆ 
----------------------------- 
Vương Thị Mỹ Hạnh 
ĐÁNH GIÁ VÀ MÔ PHỎNG CÁC HỆ SỐ ĐÀN HỒI 
ĐA TINH THỂ HỖN ĐỘN 
Chuyên ngành: Cơ học vật rắn 
Mã số: 9440107 
LUẬN ÁN TIẾN SỸ NGÀNH CƠ HỌC 
NGƯỜI HƯỚNG DẪN KHOA HỌC: 
 1. GS. TSKH. Phạm Đức Chính 
 2. TS. Lê Hoài Châu 
Hà Nội – 2020 
i 
LỜI CAM ĐOAN 
Tôi xin cam đoan đây là công trình nghiên cứu của bản thân tôi, mọi số liệu 
và kết quả trình bày trong luận án hoàn toàn trung thực và chưa có tác giả khác công 
bố ở bất kỳ công trình nào từ trước tới nay. 
Tôi xin chịu hoàn toàn trách nhiệm về các nội dung khoa học trong công 
trình nghiên cứu này. 
 Nghiên cứu sinh 
 VƢƠNG THỊ MỸ HẠNH
ii 
LỜI CẢM ƠN 
Để có thể hoàn thành luận án một cách hoàn chỉnh, bên cạnh sự nỗ lực cố 
gắng của bản thân còn có sự hướng dẫn nhiệt tình của quý Thầy Cô, cũng như sự 
động viên ủng hộ của gia đình và bạn bè trong suốt thời gian học tập nghiên cứu. 
Đầu tiên, tôi xin gửi lời cảm ơn chân thành nhất đến thầy GS. TSKH Phạm 
Đức Chính, người đã trực tiếp định hướng nghiên cứu, tận tình hướng dẫn và tạo 
những điều kiện tốt nhất để tôi có thể học hỏi và hoàn thành luận án này. 
Tôi cũng xin gửi lời cảm ơn đến TS. Lê Hoài Châu, mặc dù có khó khăn về 
khoảng cách địa lý nhưng thầy đã hướng dẫn và giúp đỡ tôi rất nhiều trong nội dung 
phương pháp số của luận án. 
Xin gửi lời cảm ơn đến các thầy, cô đã từng giảng dạy tôi trong lĩnh vực 
nghên cứu Cơ học, giúp tôi có thêm kiến thức để áp dụng trong luận án. 
Cảm ơn Học Viện Khoa học và Công Nghệ, cảm ơn Ban lãnh đạo Viện Cơ 
học đã tạo điều kiện để tôi tập trung nghiên cứu. 
Cảm ơn các anh chị trong nhóm Seminar và các anh chị em đồng nghiệp đã 
giúp đỡ tôi rất nhiều về tài liệu khoa học, thời gian và kinh nghiệm. 
Cuối cùng, xin gửi lời cảm ơn đến gia đình của tôi, những người đã luôn sát 
cánh cùng tôi trong mọi khó khăn, động viên, khích lệ để tôi có thể chuyên tâm 
nghiên cứu. 
iii 
MỤC LỤC 
LỜI CAM ĐOAN .......................................................................................................... i 
DANH MỤC CÁC KÝ HIỆU, CHỮ VIẾT TẮT ...................................................... v 
DANH MỤC CÁC BẢNG ........................................................................................ viii 
DANH MỤC CÁC HÌNH VẼ..................................................................................... ix 
MỞ ĐẦU ....................................................................................................................... 1 
1. Lý do lựa chọn đề tài .............................................................................................. 1 
2. Mục tiêu, phương pháp nghiên cứu của luận án ..................................................... 2 
3. Đối tượng, phạm vi nghiên cứu của luận án ........................................................... 3 
4. Những đóng góp mới của luận án ........................................................................... 3 
5. Bố cục luận án......................................................................................................... 3 
CHƢƠNG 1: TỔNG QUAN ........................................................................................ 6 
1.1. Tổng quan về vật liệu đa tinh thể......................................................................... 6 
1.2. Lịch sử nghiên cứu các hệ số đàn hồi vật liệu đa tinh thể ................................. 12 
1.3. Phương pháp nghiên cứu các hệ số đàn hồi vật liệu đa tinh thể ........................ 20 
1.4. Kết luận chương 1 .............................................................................................. 21 
CHƢƠNG 2: XÂY DỰNG CÁC ĐÁNH GIÁ CÁC MÔ ĐUN ĐÀN HỒI VẬT 
LIỆU ĐA TINH THỂ HỖN ĐỘN D CHIỀU .......................................................... 23 
2.1. Các công thức xuất phát..................................................................................... 23 
2.2. Mô đun đàn hồi khối vật liệu đa tinh thể hỗn độn d chiều ................................ 35 
2.3. Mô đun đàn hồi trượt vật liệu đa tinh thể hỗn độn d chiều ............................... 49 
2.4. Kết luận chương 2 .............................................................................................. 52 
CHƢƠNG 3: ĐÁNH GIÁ CÁC MÔ ĐUN ĐÀN HỒI VĨ MÔ CHO CÁC ĐA 
TINH THỂ HỖN ĐỘN TỪ CÁC LỚP ĐỐI XỨNG TINH THỂ CỤ THỂ .......... 53 
3.1. Các đa tinh thể 2 chiều ....................................................................................... 53 
iv 
3.2. Các đa tinh thể 3 chiều ....................................................................................... 69 
3.3. Kết luận chương 3 .............................................................................................. 83 
CHƢƠNG 4: ÁP DỤNG PHƢƠNG PHÁP PTHH VÀ SO SÁNH VỚI CÁC 
ĐÁNH GIÁ CHO MỘT SỐ MÔ HÌNH ĐA TINH THỂ CỤ THỂ ....................... 84 
4.1. Các công thức xuất phát..................................................................................... 84 
4.2. Quy trình tính toán PTHH ................................................................................. 85 
4.3. Áp dụng cho đối xứng tinh thể cụ thể ................................................................ 91 
4.4. Kết quả PTHH và so sánh .................................................................................. 93 
4.5. Kết luận chương 4 ............................................................................................ 102 
KẾT LUẬN VÀ HƢỚNG NGHIÊN CỨU TIẾP THEO ...................................... 103 
DANH MỤC CÁC CÔNG TRÌNH KHOA HỌC ĐÃ CÔNG BỐ ....................... 105 
TÀI LIỆU THAM KHẢO ....................................................................................... 106 
PHỤ LỤC .................................................................................................................. 115 
v 
DANH MỤC CÁC KÝ HIỆU, CHỮ VIẾT TẮT 
1. Ký hiệu thông thƣờng 
2D, 3D: không gian 2 chiều, 3 chiều. 
d: số chiều không gian 
E: mô đun Young 
G: mô đun trượt (mô đun cắt) 
n: véc tơ pháp tuyến 
U: trường chuyển vị 
e
U : trường chuyển vị trên biên 
 
1
d
i
i
x
 x : vị trí trong không gian d chiều 
ε : trường biến dạng. 
σ : trường ứng suất. 
 : hệ số Poison. 
λ : hệ số Lame 
eτ : lực biên cho trước 
ij : toán tử Kronecker. 
 : hàm thế điều hòa 
  : hàm song điều hòa 
uV : biên cho chuyển vị 
V : biên cho lực 
2. Ký hiệu riêng trong lĩnh vực luận án 
A , B

 : Các tham số thống kê bậc ba về hình học pha vật liệu. 
eff
C : ten xơ hệ số đàn hồi (độ cứng) vĩ mô/ hiệu quả. 
C : ten xơ hệ số đàn hồi (độ cứng) thành phần pha α. 
I : hàm chỉ số hình học pha α 
effK : mô đun đàn hồi diện tích vĩ mô vật liệu 2D. 
effk : mô đun đàn hồi thể tích (khối) vĩ mô vật liệu 3D (d chiều). 
vi 
Vk : đánh giá mô đun đàn hồi thể tích (khối) của Voigt cho vật liệu 3D (d chiều) 
Rk : đánh giá mô đun đàn hồi thể tích (khối) của Ruess cho vật liệu 3D (d chiều) 
VK : đánh giá mô đun đàn hồi diện tích của Voigt cho vật liệu 2D 
RK : đánh giá mô đun đàn hồi diện tích của Ruess cho vật liệu 2D 
U
HSk : đánh giá trên mô đun đàn hồi thể tích (khối) của HS cho vật liệu 3D (d chiều) 
L
HSk : đánh giá dưới mô đun đàn hồi thể tích (khối) của HS cho vật liệu 3D (d chiều) 
U
HSK : đánh giá trên mô đun đàn hồi diện tích của HS cho vật liệu 2D 
L
HSK : đánh giá dưới mô đun đàn hồi diện tích của HS cho vật liệu 2D 
Uk : đánh giá trên mô đun đàn hồi thể tích (khối) của luận án cho vật liệu 3D (d 
chiều) 
Lk : đánh giá dưới mô đun đàn hồi thể tích/ khối của luận án cho vật liệu 3D (d 
chiều) 
UK : đánh giá trên mô đun đàn hồi diện tích của luận án cho vật liệu 2D 
LK : đánh giá dưới mô đun đàn hồi diện tích của luận án cho vật liệu 2D 
eff : mô đun đàn hồi trượt vĩ mô. 
V : đánh giá mô đun đàn hồi trượt của Voigt 
R : đánh giá mô đun đàn hồi trượt của Ruess 
U
HS : đánh giá trên mô đun đàn hồi trượt của HS 
L
HS : đánh giá dưới mô đun đàn hồi trượt của HS 
U : đánh giá trên mô đun đàn hồi trượt của luận án 
L : đánh giá dưới mô đun đàn hồi trượt của luận án 
1f , 1g , 2f , 2g : các tham số hình học pha vật liệu 
Ub , 1
Uf , 1
Ug : các giá trị các biến 1 1, ,b f g đạt được ứng với đánh giá trên của luận án 
Lb , 1
Lf , 1
Lg : các giá trị các biến 1 1, ,b f g đạt được ứng với đánh giá dưới của luận án 
: trung bình thể tích trên miền V 
: trung bình thể tích trên miền V 
v : tỷ lệ thể tích pha α 
vii 
α, β, γ: các pha của vật liệu (chỉ hướng của các tinh thể) 
S: ten xơ hệ số đàn hồi mềm 
S : ten xơ hệ số đàn hồi mềm thành phần pha α. 
kS , S : tương ứng là các tham số phân tán của hệ số đàn hồi khối/ thể tích/ diện tích 
và trượt vĩ mô của vật liệu đa tinh thể ngẫu nhiên 
LA
kS , 
LAS : các tham số phân tán cho mô đun đàn hồi diện tích/thể tích và trượt của 
Luận án 
irc
kS , 
ircS : các tham số phân tán cho mô đun đàn hồi diện tích/ thể tích và trượt của 
tinh thể dạng hạt tròn. 
VR
kS , 
VRS : các tham số phân tán cho mô đun đàn hồi diện tích/ thể tích và trượt của 
Voigt- Reuss 
HS
kS ,
HSS : các tham số phân tán cho mô đun đàn hồi diện tích/ thể tích và trượt của 
Hashin-Strickman 
3. Chữ viết tắt 
HS: Hashin-Strickman 
NCS: nghiên cứu sinh 
PĐC: Phạm Đức Chính 
PTHH: phần tử hữu hạn 
RVE: phần tử đặc trưng (representative volum element) 
SC: tự tương hợp (self-consistent) 
V-R: Voigt- Reuss 
viii 
DANH MỤC CÁC BẢNG 
Bảng 1.1: Các mạng Bravais ............. 9 
Bảng 1.2: Quan hệ giữa các cặp hệ số mô đun đàn hồi .................. 24 
Bảng 3.1: Các hệ số đàn hồi một số tinh thể orthorhombic 2D .............. 59 
Bảng 3.2: Kết quả mô đun đàn hồi diện tích orthorhombic 2D ............................. 61 
Bảng 3.3: Kết quả mô đun đàn hồi diện tích square ............................................... 66 
Bảng 3.4: Kết quả mô đun đàn hồi trượt square ..................................................... 66 
Bảng 3.5: Các hệ số đàn hồi một số tinh thể tetragonal 2D ................................... 69 
Bảng 3.6: Kết quả mô đun đàn hồi diện tích tetragonal 2D .................................... 69 
Bảng 3.7: Kết quả mô đun đàn hồi trượt tetragonal 2D .......................................... 70 
Bảng 3.8: Các hệ số đàn hồi một số tinh thể tetragonal 3D .................................... 80 
Bảng 3.9: Kết quả mô đun đàn hồi thể tích tetragonal 3D....................................... 81 
Bảng 3.10: Kết quả mô đun đàn hồi trượt tetragonal 3D ........................................ 82 
ix 
DANH MỤC CÁC HÌNH VẼ 
Hình 1.1: Ứng dụng của vật liệu đa tinh thể ............................................................. 6 
Hình 1.2: Mô hình vật liệu đa tinh thể hỗn độn ........................................................ 6 
Hình 3.1: Đối xứng tinh thể orthorhombic .. 55 
Hình 3.2: Quy trình tính các đánh giá bằng Matlab ................................................ 60 
Hình 3.3: Đối xứng tinh thể tetragonal  .. 67 
Hình 4.1 (a-e): Kích thước RVE ........................................................................ 87 
Hình 4.2: Quy trình tính toán bằng phương pháp PTHH ....................................... 88 
Hình 4.3: Kết quả PTHH cho mô đun trượt square Cu ........................................... 95 
Hình 4.4: Kết quả PTHH cho mô đun trượt square Pb ........................................... 96 
Hình 4.5: Kết quả PTHH cho mô đun trượt square Li ............................................ 97 
Hình 4.6: Kết quả PTHH cho mô đun diện tích orthorhombic 2D S(1) ................. 97 
Hình 4.7: Kết quả PTHH cho mô đun diện tích orthorhombic 2D TiO2 ................ 98 
Hình 4.8: Kết quả PTHH cho mô đun diện tích orthorhombic 2D U(1) ................. 98 
Hình 4.9: Kết quả PTHH cho mô đun trượt orthorhombic 2D S(1) ....................... 99 
Hình 4.10: Kết quả PTHH cho mô đun trượt orthorhombic 2D S(3) ..................... 99 
Hình 4.11: Kết quả PTHH cho mô đun diện tích tetragonal 2D BaTiO3 .............. 100 
Hình 4.12: Kết quả PTHH cho mô đun diện tích tetragonal 2D Hg2Cl2 ............... 100 
Hình 4.13: Kết quả PTHH cho mô đun diện tích tetragonal 2D In ....................... 101 
Hình 4.14: Kết quả PTHH cho mô đun trượt tetragonal 2D Hg2Cl2 ..................... 101 
Hình 4.15: Kết quả PTHH cho mô đun trượt tetragonal 2D In ............................. 102
1 
MỞ ĐẦU 
1. Lý do lựa chọn đề tài 
a. Nguyên nhân khách quan 
 Phần lớn trong số 90 nguyên tố hóa học và hơn 3000 khoáng chất trong tự 
nhiên đều có dạng tinh thể. Trong chế tạo, các vật liệu có cấu trúc tinh thể được sử 
dụng nhiều và cho hiệu quả cao. Ngoài ra, vật liệu đa tinh thể còn được ứng dụng 
rất nhiều trong sản xuất các sản phẩm phục vụ mọi mặt của đời sống con người. Với 
những ưu điểm là có sẵn trong tự nhiên, dễ dàng trong chế tạo và có tính ứng dụng 
cao, trong tương lai vật liệu đa tinh thể sẽ ngày càng được sử dụng rộng rãi. 
 Trong quá trình chế tạo và sử dụng, ta cần quan tâm lớn đến các tính chất đàn 
hồi của vật liệu này. Mô đun đàn hồi vĩ mô của vật liệu đa tinh thể phụ thuộc vào các 
tính chất của từng đơn tinh thể thành phần cấu thành, hình học và tương tác giữa 
chúng. Do đó việc nghiên cứu các tính chất này là cần thiết và có ý nghĩa khoa học 
lẫn thực tiễn, giúp giải thích được mối quan hệ giữa tính chất vĩ mô của vật liệu với 
các tính chất và hình học vi mô của các đơn tinh thể, giúp thiết kế vật liệu với các 
tính chất theo yêu cầu. 
 Hướng nghiên cứu “đồng nhất hóa vật liệu” hình thành từ thế kỷ 19, trải qua 
nhiều năm phát triển đã đạt được nhiều thành tựu nổi bật với các kết quả ngày càng 
tốt hơn. Những kết quả giải tích cho các hệ số đàn hồi của vật liệu đa tinh thể tiêu 
biểu như: Voigt, Ruess, Hill, Paul, Hashin-Strikman, Phạm Đức Chính... Tuy nhiên 
những kết quả mô phỏng số cụ thể chưa nhiều và mới có công bố của Lê Hoài Châu -
áp dụng phương pháp phần tử hữu hạn (PTHH) trong nghiên cứu mô đun trượt của 
đa tinh thể cấu tạo từ vô số các đơn tinh thể có cấu trúc square liên kết hỗn độn với 
nhau (gọi tắt là đối xứng tinh thể square). 
 Câu hỏi đặt ra là: những đánh giá trên có phải là những kết quả tốt nhất, có tồn 
tại hay không các đánh giá tốt hơn, làm thế nào để xây dựng được các đánh giá có 
thể đáp ứng được những yêu cầu ngày càng cao của cuộc sống con người, những tính 
toán số cho mô đun đàn hồi như thế nào, liệu có khác biệt nhiều so với các kết quả 
giải tích, ? 
 Do đó, việc nghiên cứu các hệ số đàn hồi vật liệu đa tinh thể là một vấn đề có 
ý nghĩa lớn trong lĩnh vực khoa học vật liệu nói riêng và cơ học nói chung. 
b. Nguyên nhân chủ qu ... 1,4); C(3,5) = C(1,4); C(3,6) = C(1,4); 
% C(4,1) = C(1,4); C(4,2) = C(1,4); C(4,3) = C(1,4); C(4,5) = 
C(1,4); C(4,6) = C(1,4); 
% C(5,1) = C(1,4); C(5,2) = C(1,4); C(5,3) = C(1,4); C(5,4) = 
C(1,4); C(5,6) = C(1,4); 
% C(6,3) = C(1,4); C(6,4) = C(1,4); C(6,5) = C(1,4); 
% %----------------------------------------------------------------------
---- 
% 
% % Du lieu dau vao Tetragonal CaMoO4 
% C(1,1) = 144; C(2,2) = C(1,1); 
% C(3,3) = 127; 
% C(1,2) = 65; C(2,1) = C(1,2); 
% C(1,3) = 47; C(3,1) = C(1,3); C(2,3) = C(1,3); C(3,2) = 
C(1,3); 
% C(4,4) = 36.8; C(5,5) = C(4,4); 
% C(6,6) = 45.8; 
% C(1,6) = -13.5; C(6,1) = C(1,6); C(2,6) = -C(1,6); C(6,2) = 
-C(1,6); 
% C(1,4) = 0; C(1,5) = C(1,4); 
% C(2,4) = C(1,4); C(2,5) = C(1,4); 
% C(3,4) = C(1,4); C(3,5) = C(1,4); C(3,6) = C(1,4); 
% C(4,1) = C(1,4); C(4,2) = C(1,4); C(4,3) = C(1,4); C(4,5) = 
C(1,4); C(4,6) = C(1,4); 
% C(5,1) = C(1,4); C(5,2) = C(1,4); C(5,3) = C(1,4); C(5,4) = 
C(1,4); C(5,6) = C(1,4); 
% C(6,3) = C(1,4); C(6,4) = C(1,4); C(6,5) = C(1,4); 
122 
% %----------------------------------------------------------------------
---- 
% 
f1 = a(1); 
g1 = a(2); 
b = x; 
f3=2/3-f1; 
g3=4/5-g1; 
% Tinh S = C^-1 
S(1,1) = C(2,2)/(C(1,1)*C(2,2) - (C(1,2))^2); 
S(1,2) = -C(1,2)/(C(1,1)*C(2,2) - (C(1,2))^2); 
S(2,2) = C(1,1)/(C(1,1)*C(2,2) - (C(1,2))^2); 
S(2,1) = S(1,2); 
S(3,3) = 1/C(3,3); 
S(1,3) = 0; 
S(2,3) = S(1,3); S(3,2) = S(1,3); 
% S 
% Tinh KV, MV, KR, MR: 
KV = 1/4*(C(1,1) + 2*C(1,2) + C(2,2)); 
MV = 1/8*(C(1,1) + C(2,2)) - 1/4*C(1,2) + 1/2*C(3,3); 
KR = (S(1,1) + 2*S(1,2) + S(2,2))^-1; 
MR = 2* (S(1,1) + S(2,2) - 2*C(1,2) + S(3,3))^-1; 
% Nhap cac thong so Fi, Gi: 
F1 = -1/15 - 8*b/105 -13*b^2/315; G1 = 43*b^2/1890; 
F2 = 1/10 + 4*b/35 + 16*b^2/315; G2 = -4*b^2/189; G4 = G2; G6 = G2; 
F3 = -4*b/105 - 4*b^2/315; G3 = -b^2/945; 
F4 = 2*b/35 + 16*b^2/315; F6 = F4; F5 = 1/10 + 4*b/35 - 4*b^2/315; G5 = 
10*b^2/189; 
F7 = -4*b^2*KV/105 + 4*b^2*MV/63; G7 = b^2*KV/630 - 5*b^2*MV/189; 
F8 = b^2*KV/105 - 10*b^2*MV/63; G8 = 2*b^2*KV/63 + 5*b^2*MV/27; 
% Tinh D: 
D1 = (KV - MV)*(f2*F1 + g2*G1) + MV*(f2*F2 + g2*G2) + 2*KV*(f2*F3 + 
g2*G3) +... 
 (KV + 2*MV)*(f2*F4 + g2*G2) + 1/2*(f1 + f2)*F7 + 5/8*(g1 + g2)*G7 
+ (b^2*KV)/16; 
D2 = 2*MV*(f2*F1 + g2*G1) + KV*(f2*F2 + g2*G2) + (KV + 2*MV)*(f2*F5 + 
g2*G5) + ... 
 2*KV*(f2*F6 + g2*G6) + 1/2*(f1 + f2)*F8 + 5/8*(g1 + g2)*G8; 
D(1,1) = D1 + D2; D(2,2) = D(1,1); 
D(3,3) = 1/2*D2; 
D(1,2) = D1; D(2,1) = D(1,2); 
D(1,3) = 0 ; D(3,1) = D(1,3); D(2,3) = D(1,3); D(3,2) = D(1,3); 
D; 
% Tinh CK: 
CK(1,1) = 1/9*(C(1,1) + C(1,2) + C(1,3))*(1+2*b/5) + b*KV/5; 
CK(2,2) = CK(1,1); 
CK(3,3) = 1/9*(C(3,3) + 2*C(1,3))*(1+2*b/5) + b*KV/5; 
123 
CK(1,2)= 0; CK(1,3)= CK(1,2); CK(2,3)= CK(1,2); 
% CK 
B1 = 1/4*(1+b/2)^2 + f1*F1 + g1*G1; 
B2 = f1*F2 + g1*G2; 
B3 = b/8 + b^2/16 + f1*F3 + g1*G3; 
B4 = f1*F4 + g1*G4; 
B5 = f1*F5 + g1*G5; 
B6 = f1*F6 + g1*G6; 
% Tinh CA: 
CA(1,1) = C(1,1)*(B1 + B2) + (C(1,1) + C(1,2) + C(1,3))*(B3 + B6) + 
(C(1,1) + C(4,4) + C(6,6))*(B4 + B5); 
CA(2,2) = CA(1,1); 
CA(3,3) = C(3,3)*(B1 + B2) + (C(3,3) + 2*C(1,3))*(B3 + B6) + (C(3,3) + 
2*C(4,4))*(B4 + B5); 
CA(1,2) = C(1,2)*B1 + C(6,6)*B2 + (C(1,1) + C(1,2) + C(1,3))*B3 + 
(C(1,1) + C(4,4) + C(6,6))*B4; 
CA(2,1) = CA(1,2); 
CA(1,3) = C(1,3)*B1 + C(4,4)*B2 + (C(1,1) + C(3,3) + C(1,2) + 
3*C(1,3))*1/2*B3 + (C(1,1) + C(3,3) + 3*C(4,4) + C(6,6))*1/2*B4; 
CA(3,1) = CA(1,3); 
CA(2,3) = CA(1,3); CA(3,2) = CA(1,3); 
CA(4,4) = C(4,4)*B1 + (C(1,3) + C(4,4))*1/2*B2 + (C(1,1) + C(3,3) + 
3*C(4,4) + C(6,6))*1/4*B5 + (C(1,1) + C(3,3) + C(1,2) + 3*C(1,3))*1/4*B6; 
CA(5,5) = CA(4,4); 
CA(6,6) = C(6,6)*B1 + (C(1,2) + C(6,6))*1/2*B2 + (C(1,1) + C(4,4) + 
C(6,6))*1/2*B5 + (C(1,1) + C(1,2) + C(1,3))*1/2*B6; 
CA(1,6) = C(1,6)*(B1 + B2); 
CA(6,1) = CA(1,6); 
CA(2,6) = -CA(1,6); 
CA(6,2) = CA(2,6); 
A = CA + D; 
SA = inv(A); 
% Tinh C_AC: 
C_AC(1,1) = (SA(1,1)+ SA(1,2))*CK(1,1) + SA(1,3)*CK(3,3); 
C_AC(2,2) = C_AC(1,1); 
C_AC(3,3) = SA(3,3)*CK(3,3) + 2*SA(1,3)*CK(1,1); 
C_AC(1,2) = 0; C_AC(2,1) = C_AC(1,2); C_AC(1,3) = C_AC(1,2); C_AC(3,1) = 
C_AC(1,2); C_AC(2,3) = C_AC(1,2); C_AC(3,2) = C_AC(1,2); 
% Tinh C_CAC: 
C_CAC = 2*CK(1,1)*C_AC(1,1) + CK(3,3)*C_AC(3,3); 
KKR = (2*SA(1,1) + SA(3,3) + 2*SA(1,2) + 4*SA(1,3))^-1; 
% Cong thuc cuoi: 
K_U = KV + (2*C_AC(1,1)+C_AC(3,3))^2*KKR - C_CAC; 
end 
124 
b. Chương trình tính cho kết quả trong chương 4 của luận án 
Chương trình tính cho mô đun đàn hồi trượt đối xứng tinh thể square: 
%Calculating HS bounds 
% %Copper 
% c11=169; 
% c12=122; 
% c33=75.3; 
% folder = 'matlab'; 
% [Eh]=homogenization2D('unitCell8.dat',c11,c12,c33) 
% Eh 
% return 
%Lead 
%c11=48.8; 
%c12=41.4; 
%c33=14.8; 
%folder = 'lead'; 
% Lithium 
c11=13.6; 
c12=11.4; 
c33=9.8; 
folder = 'Lithium'; 
% mu4 = []; 
% for i = 1:10 
% [Eh]=homogenization2D('unitCell4.dat',c11,c12,c33) 
% mu4(i)=Eh(3,3); 
% end 
% save('mu4') 
% 
% mu8 = []; 
% for i = 1:10 
% [Eh]=homogenization2D('unitCell8.dat',c11,c12,c33) 
% mu8(i)=Eh(3,3); 
% end 
% save('mu8') 
% 
% mu16 = []; 
% for i = 1:10 
% [Eh]=homogenization2D('unitCell16.dat',c11,c12,c33) 
% mu16(i)=Eh(3,3)/2; 
% end 
% save('mu16') 
% 
% mu16x16 = []; 
% for i = 1:10 
% [Eh]=homogenization2D('unitCell16x16.dat',c11,c12,c33) 
% mu16x16(i)=Eh(3,3)/2; 
% end 
% save('mu16x16') 
% 
% mu32 = []; 
% for i = 1:10 
125 
% [Eh]=homogenization2D('unitCell32.dat',c11,c12,c33) 
% mu32(i)=Eh(3,3)/2; 
% end 
% save('mu32') 
% mu64 = []; 
% for i = 1:10 
% [Eh]=homogenization2D('unitCell64.dat',c11,c12,c33) 
% mu64(i+2)=Eh(3,3)/2; 
% end 
% save('mu64'); 
% 
% mu64 = []; 
% for i = 1:10 
% [Eh]=homogenization2D('unitCell128.dat',c11,c12,c33) 
% mu64(i+2)=Eh(3,3)/2; 
% end 
% save('mu128'); 
%Plots 
% figure; 
% hold on; 
% 
disp('VRH') 
mu_v=(c11-c12+2*c33)/4 
mu_r=2*(c11-c12)*c33/(2*c33+c11-c12) 
% H1=plot([0,130],[mu_v,mu_v],':r', 'LineWidth',2); 
% plot([0,130],[mu_r,mu_r],':r', 'LineWidth',2); 
% 
disp('HS') 
muo_u = max((c11-c12)/2,c33) 
muo_l = min((c11-c12)/2,c33) 
Kv=(c11+c12)/2 
[P_u]=HS(Kv,muo_u,c11,c12,c33) 
[P_l]=HS(Kv,muo_l,c11,c12,c33) 
% H2=plot([0,130],[P_u,P_u],'--b', 'LineWidth',2); 
% plot([0,130],[P_l,P_l],'--b', 'LineWidth',2); 
disp('PDC') 
[P_u]=HS(Kv,mu_v,c11,c12,c33) 
[P_l]=HS(Kv,mu_r,c11,c12,c33) 
% H3=plot([0,130],[P_u,P_u],'-.m', 'LineWidth',2); 
% plot([0,130],[P_l,P_l],'-.m', 'LineWidth',2); 
disp('Self-consistent') 
[P_sc]=HS(Kv,2.9042,c11,c12,c33) %copper 40.5055, lead 7.1517; Lithium 
2.9042 
% H4=plot([0,130],[P_sc,P_sc],'k', 'LineWidth',2); 
% legend([H1 H2 H3 H4],'VRH bounds', 'HS bounds', 'New bounds', 'Self-
consistent value') 
% load(['..\',folder,'\mu4']); 
% for i=1:10 
% plot(4+(i-5)*0.02,mu4(i),'ok', 'LineWidth',1.1); 
% end 
126 
% load(['..\',folder,'\mu8']); 
% for i=1:10 
% plot(8+(i-5)*0.04,mu8(i),'sk', 'LineWidth',1.1); 
% end 
% load(['..\',folder,'\mu16']); 
% for i=1:10 
% plot(16+(i-5)*0.08,mu16(i)*2,'*k', 'LineWidth',1.1); 
% end 
% load([folder,'\..\mu16x16']); 
% for i=1:10 
% plot(16+1+(i-5)*0.08,mu16x16(i)*2,'*k', 'LineWidth',1.1); 
% end 
% 
% load(['..\',folder,'\mu32']); 
% for i=1:10 
% plot(32+1+(i-5)*0.08,mu32(i)*2,'.k', 'LineWidth',1.1); 
% end 
% 
% load(['..\',folder,'\mu64']); 
% for i=1:10 
% plot(64+1+(i-5)*0.08,mu64(i)*2,'.k', 'LineWidth',1.1); 
% end 
% load([folder,'\..\mu128']); 
% for i=1:10 
% plot(128+1+(i-5)*0.08,mu64(i)*2,'.b', 'LineWidth',1.1); 
% end 
% grid on 
% 
% xlabel('Number of hexagons per unit cell side') 
% ylabel('Effective shear moduli ( \mu^{eff} )') 
% box on 
% xlim([0,70]); 
% ylim([35,50]); 
% set(gca,'YTick',[32 34 36 38 40 42 44 46 48 50]) 
% %xlim([0,70]); 
%ylim([1.5,6.0]); 
%set(gca,'YTick',[32 34 36 38 40 42 44 46 48 50]) 
Các hàm con 
function [mpcNum, mpcDofs1, 
mpcDofs2]=createPeriodicMpcs(boundaryNodeSets) 
set1=cell2mat(boundaryNodeSets(1)); 
set2=cell2mat(boundaryNodeSets(2)); 
set3=cell2mat(boundaryNodeSets(3)); 
set4=cell2mat(boundaryNodeSets(4)); 
set5=cell2mat(boundaryNodeSets(5)); 
tmp = length(set2)*2 + length(set4)*2 + 6; 
mpcNum = [tmp,tmp,tmp]; 
mpcDofs1=zeros(mpcNum(1),3); 
mpcDofs2=zeros(mpcNum(1),3); 
d = 1; 
for i=1:length(set2) 
 node1 = set2(i); 
 node2 = set3(i); 
127 
 dof1 = node1*2-1; 
 dof2 = node2*2-1; 
 mpcDofs1(d,:) = [dof1,dof1,dof1]; 
 mpcDofs2(d,:) = [dof2,dof2,dof2]; 
 d=d+1; 
 dof1 = node1*2; 
 dof2 = node2*2; 
 mpcDofs1(d,:) = [dof1,dof1,dof1]; 
 mpcDofs2(d,:) = [dof2,dof2,dof2]; 
 d=d+1; 
end 
for i=1:length(cell2mat(boundaryNodeSets(4))) 
 node1 = set4(i); 
 node2 = set5(i); 
 dof1 = node1*2-1; 
 dof2 = node2*2-1; 
 mpcDofs1(d,:) = [dof1,dof1,dof1]; 
 mpcDofs2(d,:) = [dof2,dof2,dof2]; 
 d=d+1; 
 dof1 = node1*2; 
 dof2 = node2*2; 
 mpcDofs1(d,:) = [dof1,dof1,dof1]; 
 mpcDofs2(d,:) = [dof2,dof2,dof2]; 
 d=d+1; 
end 
for i=2:4 
 node1 = set1(1); 
 node2 = set1(i); 
 dof1 = node1*2-1; 
 dof2 = node2*2-1; 
 mpcDofs1(d,:) = [dof1,dof1,dof1]; 
 mpcDofs2(d,:) = [dof2,dof2,dof2]; 
 d=d+1; 
 dof1 = node1*2; 
 dof2 = node2*2; 
 mpcDofs1(d,:) = [dof1,dof1,dof1]; 
 mpcDofs2(d,:) = [dof2,dof2,dof2]; 
 d=d+1; 
end 
function [KE]=elem_K(coord, E) 
M1=[1 0 0 0; 
 0 0 0 1; 
 0 1 1 0]; 
% Gauss Quadrature 
n=2; 
g_coord=[-1/sqrt(3), 1/sqrt(3)]; 
KE=zeros(8,8); 
for i=1:n 
 for j=1:n 
 x=g_coord(i); y=g_coord(j); 
 temp=[-(1+y), (1+y), (1-y), -(1-y); 
128 
 (1-x), (1+x), -(1+x), -(1-x)]; 
 temp=temp/4; 
 J_mat=temp*coord; 
 J=det(J_mat); 
 % M2 relates derivatives w.r.t. physical coordinates to that 
w.r.t. 
 % natural coordinates 
 M2(1:2,1:2)=inv(J_mat); M2(3:4,3:4)=M2(1:2,1:2); 
 M3=[-(1+y), 0 , (1+y), 0 , (1-y), 0 , -(1-y) , 0; 
 (1-x), 0 , (1+x), 0 , -(1+x), 0 , -(1-x) , 0; 
 0 , -(1+y), 0 , (1+y), 0 , (1-y), 0 , -(1-y); 
 0 , (1-x), 0 , (1+x), 0 , -(1+x), 0 , -(1-x)]; 
 M3=M3/4; 
 B=M1*M2*M3; 
 KE=KE+B'*E*B*J; 
 end 
end 
%%%%%%%%%% FE-ANALYSIS --------------------------------------------------
------------------------ 
function 
[U]=fea_static(coord,elemCon,E,F,spcDofs,spcNum,mpcDofs1,mpcDofs2,mpcNum) 
%global variable required: pC - the power parameter for interpolation of 
stiffness 
%only isoparametric quadrilateral is implemented for now 
[elemNum,~]=size(elemCon); 
[nodeNum,~]=size(coord); 
dofNum=nodeNum*2; 
[~,caseNum]=size(F); 
% assemble K 
Kval=zeros(elemNum*64,1); i_ind=zeros(elemNum*64,1); 
j_ind=zeros(elemNum*64,1); 
for elem=1:elemNum 
 ind=(elem-1)*64+1; 
 nodes=elemCon(elem,:); 
 eCoord=coord(nodes,:); 
 edofs=[nodes(1)*2-1; nodes(1)*2; nodes(2)*2-1; nodes(2)*2; 
nodes(3)*2-1; nodes(3)*2; nodes(4)*2-1; nodes(4)*2]; 
 i_ind(ind:ind+63)=[edofs; edofs; edofs; edofs; edofs; edofs; edofs; 
edofs]; 
 j_ind(ind:ind+63)=reshape([edofs'; edofs'; edofs'; edofs'; edofs'; 
edofs'; edofs'; edofs'],64,1); 
 KE=elem_K(eCoord,E(:,:,elem)); 
 Kval(ind:ind+63)=reshape(KE,64,1); 
end 
K=sparse(i_ind,j_ind,Kval,dofNum,dofNum); 
alldofs = [1:1:dofNum]; 
for i=1:caseNum 
 removedDofs=union(mpcDofs2(1:mpcNum(i),i),spcDofs(1:spcNum(i),i)); 
 freedofs = setdiff(alldofs,removedDofs); 
 % transformation matrix for mpcs 
 T=speye(dofNum,dofNum); 
 for j=1:mpcNum(i) 
 T(mpcDofs2(j,i),mpcDofs1(j,i))=1; 
 T(mpcDofs2(j,i),mpcDofs2(j,i))=0; 
129 
 end 
 Kt=T'*K*T; 
 Ft=T'*F(:,i); 
 U(freedofs,i) = Kt(freedofs,freedofs) \ Ft(freedofs); 
 U(spcDofs(1:spcNum(i),i),i)= 0.0; 
 U(mpcDofs2(1:mpcNum(i),i),i)= U(mpcDofs1(1:mpcNum(i),i),i); 
end 
function [Eh]=homogenization2D(meshFile,c11,c12,c33) 
[coord,elemCon,elemProp,~,~,boundaryNodeSets] = readNastran(meshFile); 
%mpcs 
[mpcNum, mpcDofs1, mpcDofs2]=createPeriodicMpcs(boundaryNodeSets); 
%spcs 
nodeSet1=cell2mat(boundaryNodeSets(1)); 
node1 = nodeSet1(1); 
spcDofs = [node1*2-1, node1*2-1, node1*2-1; 
 node1*2, node1*2, node1*2]; 
spcNum = [2,2,2]; 
% 
[elemNum,~]=size(elemCon); 
[nodeNum,~]=size(coord); 
%Get Random Elasticity Tensor 
[E]=randomCrystalOrientation(c11,c12,c33,elemNum,elemProp); 
dofNum=nodeNum*2; 
% Three test strain cases 
epsilon=[1 0 0; 
 0 1 0; 
 0 0 1]; 
%Material forces: 3 load cases 
F=zeros(dofNum,3); 
% 
for elem=1:elemNum 
 nodes=elemCon(elem,:); 
 eDofs=[nodes(1)*2-1; nodes(1)*2; nodes(2)*2-1; nodes(2)*2; 
nodes(3)*2-1; nodes(3)*2; nodes(4)*2-1; nodes(4)*2]; 
 eCoord=coord(nodes,:); 
 Ba=int_B(eCoord); 
 F(eDofs,:)=F(eDofs,:)-Ba'*E(:,:,elem)*epsilon; 
end 
% FE Analysis to get the characteristic displacements 
%mpcNum = mpcNum-6; 
%mpcDofs1 = mpcDofs1(1:116,:); 
%mpcDofs2 = mpcDofs2(1:116,:); 
[chi]=fea_static(coord,elemCon,E,F,spcDofs,spcNum,mpcDofs1,mpcDofs2,mpcNu
m); 
% while 1 
% cla; 
% 
plot_results_Q4(coord+(reshape(chi(:,3),2,length(chi(:,3))/2))'*0.1,elemC
on); 
% ylim([-1100,100]); 
% xlim([-100,1100]); 
130 
% pause(1); 
% cla; 
% 
plot_results_Q4(coord+(reshape(chi(:,3),2,length(chi(:,3))/2))'*0.0,elemC
on); 
% ylim([-1100,100]); 
% xlim([-100,1100]); 
% pause(1); 
% end 
% Calculate homogenized stiffness tensor and sensitivity 
Eh=zeros(3,3); 
%[1 2 3; 
% 0 4 5; 
% 0 0 6]; 
Area=0; 
for elem=1:elemNum 
 nodes=elemCon(elem,:); 
 eDofs=[nodes(1)*2-1; nodes(1)*2; nodes(2)*2-1; nodes(2)*2; 
nodes(3)*2-1; nodes(3)*2; nodes(4)*2-1; nodes(4)*2]; 
 eCoord=coord(nodes,:); 
 eChi=chi(eDofs,:); 
 for i=1:3 
 for j=i:3 
 Eh(i,j)=Eh(i,j)+int_hom(eCoord,epsilon,eChi,E(:,:,elem),i,j); 
 end 
 end 
 Ae=int_A(eCoord); 
 Area=Area+Ae; 
end 
Eh(2,1)=Eh(1,2); 
Eh(3,1)=Eh(1,3); 
Eh(3,2)=Eh(2,3); 
Eh=Eh/Area; 
function [P]=HS(k,mu,c11,c12,c33) 
mu_star=k*mu/(k+2*mu); 
c11_ps=c11+mu+mu_star; 
c33_ps=c33+mu_star; 
c12_ps=c12+mu-mu_star; 
P=2*(c11_ps-c12_ps)*c33_ps/(2*c33_ps+c11_ps-c12_ps)-mu_star; 
function [A]=int_A(coord) 
%integration of B over isoparametric the element 
% Gauss Quadrature 
n=2; 
g_coord=[-1/sqrt(3), 1/sqrt(3)]; 
A=0; 
for i=1:n 
 for j=1:n 
 x=g_coord(i); y=g_coord(j); 
 temp=[-(1+y), (1+y), (1-y), -(1-y); 
 (1-x), (1+x), -(1+x), -(1-x)]/4; 
 J_mat=temp*coord; 
 J=det(J_mat); 
 A=A+J; 
131 
 end 
end 
function [B]=int_B(coord) 
%integration of B over isoparametric the element 
M1=[1 0 0 0; 
 0 0 0 1; 
 0 1 1 0]; 
% Gauss Quadrature 
n=2; 
g_coord=[-1/sqrt(3), 1/sqrt(3)]; 
B=zeros(3,8); 
for i=1:n 
 for j=1:n 
 x=g_coord(i); y=g_coord(j); 
 temp=[-(1+y), (1+y), (1-y), -(1-y); 
 (1-x), (1+x), -(1+x), -(1-x)]; 
 temp=temp/4; 
 J_mat=temp*coord; 
 J=det(J_mat); 
 % M2 relates derivatives w.r.t. physical coordinates to that 
w.r.t. 
 % natural coordinates 
 M2(1:2,1:2)=inv(J_mat); M2(3:4,3:4)=M2(1:2,1:2); 
 M3=[-(1+y), 0 , (1+y), 0 , (1-y), 0 , -(1-y) , 0; 
 (1-x), 0 , (1+x), 0 , -(1+x), 0 , -(1-x) , 0; 
 0 , -(1+y), 0 , (1+y), 0 , (1-y), 0 , -(1-y); 
 0 , (1-x), 0 , (1+x), 0 , -(1+x), 0 , -(1-x)]; 
 M3=M3/4; 
 B=B+M1*M2*M3*J; 
 end 
end