Luận án Xấp xỉ tương đương cơ tính vật liệu tổ hợp có các cốt liệu phức hợp

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Ệ 
----------------------------- 
NGUYỄN THỊ HƯƠNG GIANG 
XẤP XỈ TƯƠNG ĐƯƠNG CƠ TÍNH VẬT LIỆU TỔ HỢP 
CÓ CÁC CỐT LIỆU PHỨC HỢP 
LUẬN ÁN TIẾN SĨ CƠ HỌC 
Hà Nội – Năm 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Ệ 
----------------------------- 
NGUYỄN THỊ HƯƠNG GIANG 
XẤP XỈ TƯƠNG ĐƯƠNG CƠ TÍNH VẬT LIỆU TỔ HỢP 
CÓ CÁC CỐT LIỆU PHỨC HỢP 
Chuyên ngành: Cơ học vật rắn 
Mã số: 9 44 01 07 
LUẬN ÁN TIẾN SĨ CƠ HỌC 
NGƯỜI HƯỚNG DẪN KHOA HỌC: 
 1. PGS. TSKH. Phạm Đức Chính 
 2. PGS.TS. Trần Bảo Việt 
Hà Nội – Năm 2020 
LỜI CAM ĐOAN 
 Tôi xin cam đoan đây là công trình nghiên cứu riêng của tôi, mọi số liệu 
và kết quả trong luận án là trung thực và cũng chưa có tác giả khác công bố ở 
bất cứ công trình nghiên cứu nào từ trước tới nay. 
 Tôi xin chịu hoàn toàn trách nhiệm về nội dung khoa học của công trình 
này. 
 Nghiên cứu sinh 
 Nguyễn Thị Hương Giang 
LỜI CẢM ƠN 
 Với lòng biết ơn sâu sắc, tôi xin chân thành cảm ơn PGS. TSKH. Phạm 
Đức Chính, TS. Trần Bảo Việt – những người thày đã tận tình hướng dẫn, 
động viên giúp đỡ và tạo mọi điều kiện cho tôi hoàn thành luận án này. 
 Tôi cũng xin bày tỏ lòng biết ơn đến quý Thày, Cô đã giảng dạy tôi trong 
thời gian học chuyên đề trong khuôn khổ chương trình đào tạo Tiến sĩ, các cán 
bộ của Học viện Khoa học & Công nghệ, nhóm nghiên cứu tại Viện Cơ học đã 
giúp đỡ hỗ trợ tôi tài liệu, kinh nghiệm để hoàn thiện luận án. 
 Các nghiên cứu trong luận án cũng được hỗ trợ bởi Quỹ phát triển Khoa 
học và Công nghệ quốc gia (NAFOSTED). 
 Xin gửi lời cảm ơn đến Trường Đại học Giao thông Vận tải, nơi tôi công 
tác, đã hỗ trợ học phí và tạo mọi điều kiện về thời gian cho tôi hoàn thành luận 
án. 
 Cuối cùng xin gửi lời cảm ơn đến gia đình nhỏ của tôi, những người luôn 
gần gũi và là động lực cho tôi trong suốt quá trình học tập và nghiên cứu để 
hoàn thành luận án này. 
MỤC LỤC 
Lời cam đoan .................................................................................................... 2 
Lời cảm ơn ...................................................................................................... 3 
Danh mục các ký hiệu, các chữ viết tắt ............................................................ 7 
Danh mục các bảng .......................................................................................... 8 
Danh mục các hình vẽ, đồ thị ........................................................................... 10 
MỞ ĐẦU .......................................................................................................... 13 
CHƯƠNG 1. TỔNG QUAN ........................................................................... 18 
1.1. Phân loại vật liệu Composite...................................................... 18 
1.2. Hệ số dẫn .................................................................................... 19 
1.3. Các mô đun đàn hồi .................................................................... 21 
1.4. Phần tử thể tích đặc trưng .......................................................... 23 
1.5. Các xấp xỉ và đánh giá xác định các giá trị hiệu dụng của 
vật liệu ........................................................................................ 24 
1.5.1. Phương pháp xấp xỉ trung bình ....................................... 24 
1.5.2. Đường bao của các giá trị hiệu dụng ............................. 29 
1.5.3. Phương pháp cốt tương đương ....................................... 30 
1.6. Phương pháp số .......................................................................... 34 
1.7. Kết luận ...................................................................................... 35 
CHƯƠNG 2. XẤP XỈ CỐT TƯƠNG ĐƯƠNG HỆ SỐ DẪN VĨ MÔ VẬT LIỆU 
 CÓ CỐT LIỆU PHỨC HỢP ........................................................................... 36 
 2.1. Vật liệu cốt sợi phức hợp đồng phương, lớp phủ quanh 
 cốt đẳng hướng ............................................................................ 37 
 2.1.1. Mô hình vật liệu ................................................................ 37 
 2.1.2. Các công thức đánh giá hệ số dẫn của vật liệu cốt tròn .. 38 
 2.1.3. Xấp xỉ tương đương với cốt tròn được phủ ...................... 43 
 2.1.4. So sánh với kết quả thực nghiệm ...................................... 52 
 2.2. Vật liệu cốt sợi phức hợp đồng phương, lớp phủ quanh 
 cốt dị hướng ................................................................................ 55 
 2.2.1. Mô hình vật liệu ................................................................ 55 
 2.2.2. Xấp xỉ tương đương với cốt tròn được phủ, lớp phủ 
 dị hướng ...................................................................................... 55 
 2.3. Kết luận ...................................................................................... 58 
CHƯƠNG 3. XẤP XỈ CỐT TƯƠNG ĐƯƠNG MÔ ĐUN ĐÀN HỒI VĨ MÔ VẬT 
 LIỆU CÓ CỐT LIỆU PHỨC HỢP ................................................................. 59 
 3.1. Mô đun đàn hồi của vật liệu Composite cốt hạt hình cầu với 
 cốt phức hợp .......................................................................................... 59 
 3.1.1. Mô hình vật liệu ................................................................ 59 
 3.1.2. Mô đun đàn hồi thể tích.................................................... 59 
 3.1.3. Mô đun đàn hồi trượt ....................................................... 63 
 3.1.4. Công thức tổng quát ......................................................... 70 
 3.1.5. Kiểm tra và so sánh .......................................................... 71 
 3.2. Mô đun đàn hồi của vật liệu Composite cốt sợi phức hợp đồng 
 phương .................................................................................................. 79 
 3.2.1. Mô hình vật liệu ................................................................ 79 
 3.2.2. Mô đun đàn hồi diện tích hiệu dụng ................................. 80 
 3.2.3. Mô đun đàn hồi trượt doc hiệu dụng ................................ 81 
 3.2.4. Mô đun đàn hồi Young dọc trục và hệ số Poisson ........... 82 
 3.2.5. Mô đun đàn hồi trượt ngang ............................................ 84 
 3.3. Kết luận .......................................................................................... 85 
CHƯƠNG 4. MÔ PHỎNG SỐ PHẦN TỬ HỮU HẠN VẬT LIỆU PHỨC 
HỢP...86 
 4.1. Vật liệu tuần hoàn .......................................................................... 86 
 4.2. Các công thức xuất phát ................................................................. 86 
 4.2.1. Mô đun đàn hồi ................................................................ 87 
 4.2.2. Hệ số dẫn .......................................................................... 88 
 4.3. Phần mềm Cast3M ......................................................................... 89 
 4.4. Tính toán cho mô hình vật liệu cụ thể và so sánh kết quả ............. 91 
 4.4.1. Vật liệu Composite cốt sợi phức hợp đồng phương 
 4.4.1.1. Hệ số dẫn ngang khi các pha đồng nhất, đẳng 
 hướng .................................................................................. 91 
 4.4.1.2. Hệ số dẫn ngang khi lớp vỏ bọc dị hướng ............ 96 
 4.4.1.3. Mô đun đàn hồi ..................................................... 99 
 4.4.2. Vật liệu Composite với cốt hình cầu ................................ 106 
 4.5. Kết luận .......................................................................................... 118 
KẾT LUẬN VÀ KIẾN NGHỊ.......................................................................... 119 
DANH MỤC CÔNG TRÌNH CỦA TÁC GIẢ ................................................ 121 
TÀI LIỆU THAM KHẢO ................................................................................ 122 
PHỤ LỤC .................................................................................................... 130 
7 
DANH MỤC CÁC KÝ HIỆU, CÁC CHỮ VIẾT TẮT 
A Ten xơ mật độ biến dạng 
σ Trường ứng suất 
ij Toán tử Krӧnecker 
c, ceff Hệ số dẫn nhiệt, hệ số dẫn nhiệt hiệu dụng 
Ceff Ten xơ độ cứng hiệu dụng 
ε, E Trường biến dạng vi mô, vĩ mô 
Eeff Mô đun đàn hồi Young hiệu dụng 
q Véc tơ dòng nhiệt 
keff, Keff Mô đun đàn hồi thể tích, diện tích hiệu dụng 
T Nhiệt độ 
u Trường chuyển vị 
μeff Mô đun đàn hồi trượt hiệu dụng 
νeff Hệ số nở ngang hiệu dụng 
  Tỉ lệ thể tích pha α 
 Trung bình các đại lượng 
 Toán tử gradient 
Δ Toán tử Laplace (Δ =  ) 
HSU, HSL Kết quả đường bao trên và dưới của Hashin – Strickman 
PTHH Phần tử hữu hạn 
RVE Phần tử thể tích đặc trưng 
TN Kết quả thí nghiệm 
CTĐ Cốt tương đương 
CTĐĐG Cốt tương đương đơn giản 
VP-CTĐ Xấp xỉ vi phân sử dụng cốt tương đương 
TĐLN Xấp xỉ theo mô hình đĩa tròn lồng nhau 
TT3Đ-CTĐ Xấp xỉ tương tác ba điểm sử dụng cốt tương đương. 
MA Xấp xỉ Maxwell 
8 
DANH MỤC CÁC BẢNG 
Bảng 1.1. Mối quan hệ giữa các mô đun đàn hồi 
Bảng2.1. Thông tin hình học bậc ba 2 của vật liệu có cốt dạng đĩa tròn phân bố ngẫu 
nhiên (không chồng lấn) 
Bảng 2.2. Hệ số dẫn hiệu dụng khi 32321 2051  ,c,c,c 
Bảng 2.3. Hệ số dẫn hiệu dụng khi 32321 1520  ,c,c,c 
Bảng 2.4. Kết quả thí nghiệm hệ số dẫn của composite sợi abaca 
Bảng 2.5. Hệ số dẫn của composite sợi abaca 
Bảng 3.1. Mô đun đàn hồi thể tích hiệu dụng ,5,1 1 GPaGPa IM  
21212 23025 IIMIII ,.,GPa  
Bảng 3.2. Mô đun đàn hồi trượt hiệu dụng ,GPa,GPa IM 51 1  
21212 23025 IIMIII ,.,GPa  
Bảng 3.3. Mô đun đàn hồi thể tích hiệu dụng ,GPa,GPa IM 525 1  
21212 2301 IIMIII ,.,GPa  
Bảng 3.4. Mô đun đàn hồi trượt hiệu dụng ,GPa,GPa IM 525 1  
21212 2301 IIMIII ,.,GPa  
Bảng 4.1. Hệ số dẫn hiệu dụng khi 32321 2051  ,c,c,c 
Bảng 4.2. Hệ số dẫn hiệu dụng khi 32321 1520  ,c,c,c 
Bảng 4.3. Hệ số dẫn hiệu dụng khi 3231 1050301001  .,c,c,c,c Tn Bảng 
4.4. Hệ số dẫn hiệu dụng khi 3231 1070501100  .,c,c,c,c Tn 
Bảng 4.5. Các giá trị của mô đun đàn hồi hiệu dụng khi ,E,.,E IMM 5201 1 
12221 401030 IIIII ,.,E,.  
Bảng 4.6. Các giá trị của mô đun đàn hồi hiệu dụng khi ,E,.,E IMM 54010 1  
12221 20130 IIIII ,.,E,.  
Bảng 4.7. Các giá trị của mô đun đàn hồi hiệu dụng khi ,.,E MM 201 
55805157 .,.E EIEI  
Bảng 4.8. Các giá trị của mô đun đàn hồi hiệu dụng khi ,.,E MM 4010 
27303 .,E EIEI  
9 
 Bảng 4.9. Các giá trị của mô đun đàn hồi hiệu dụng của nhân tử lập phương đơn giản 
khi 202024401 2211 IIIIMM ,K,,K,.,K  
 Bảng 4.10. Các giá trị của mô đun đàn hồi hiệu dụng của nhân tử lập phương tâm 
khối khi 202024401 2211 IIIIMM ,K,,K,.,K  
Bảng 4.11. Các giá trị của mô đun đàn hồi hiệu dụng của nhân tử lập phương tâm mặt 
khi 202024401 2211 IIIIMM ,K,,K,.,K  
Bảng 4.12. Các giá trị của mô đun đàn hồi hiệu dụng của nhân tử lập phương đơn 
giản khi 244011220 2211 IIIIMM ,K,.,K,,K  
Bảng 4.13. Các giá trị của mô đun đàn hồi hiệu dụng của nhân tử lập phương tâm 
khối khi 244011220 2211 IIIIMM ,K,.,K,,K  
Bảng 4.14. Các giá trị của mô đun đàn hồi hiệu dụng của nhân tử lập phương tâm mặt 
khi 244011220 2211 IIIIMM ,K,.,K,,K  
10 
DANH MỤC CÁC HÌNH VẼ 
Hình 0.1. Ứng dụng của vật liệu Composite trong máy bay Boeing 787 ...... 13 
Hình 0.2. Hình ảnh vật liệu Geopolymers có sợi alumina tráng Monazite ..... 14 
Hình 1.1. (a) - Bê tông cốt sợi, (b) – sợi polypropylene ................................. 18 
Hình 1.2. Bê tông nhựa ............................... .................................................... 19 
Hình 1.3. Phần tử thể tích đặc trưng RVE . .................................................... 23 
Hình 1.4. Mô hình cốt tương đương ................................................................ 30 
Hình 1.5. (a) – phần tử tam giác, (b) – phần tử tứ giác .................................... 34 
Hình 2.1. Vật liệu cốt sợi đồng phương ........................................................... 37 
Hình 2.2. Hệ số dẫn hiệu dụng khi 110 21 c,c ........................................... 42 
Hình 2.3. Hệ số dẫn hiệu dụng khi 101 21 c,c ........................................... 42 
Hình 2.4. Mô hình cốt tròn hai lớp đặt trong pha nền vô tận ........................... 43 
Hình 2.5. Hệ số dẫn hiệu dụng khi 2051 321 c,c,c ............................... 51 
Hình 2.6. Hệ số dẫn hiệu dụng khi 1520 321 c,c,c ............................... 51 
Hình 2.7. Mô hình vật liệu composite sợi abaca .............................................. 52 
Hình 2.8. Hệ số dẫn của composite sợi abaca .................................................. 53 
Hình 2.9. Vật liệu cốt sợi phức hợp đồng phương, lớp phủ quanh cốt dị 
hướng ............................................................................................................... 55 
Hình 3.1. (a) - Mô hình vật liệu composite với cốt hình cầu được phủ, 
(b) - Mô hình cốt tương đương ........................................................................ 59 
Hình 3.2. Mô đun đàn hồi thể tích hiệu dụng khi ,GPa,GPa IM 51 1  
21212 23025 IIMIII ,.,GPa  ............................................ 76 
Hình3.3. Mô đun đàn hồi trượt hiệu dụng khi ,GPa,GPa IM 51 1 
21212 23025 IIMIII ,.,GPa  ............................................ 76 
Hình 3.4. Mô đun đàn hồi thể tích hiệu dụng khi ,GPa,GPa IM 525 1  
21212 2301 IIMIII ,.,GPa  ............................................... 77 
Hình 3.5. Mô đun đàn hồi trượt hiệu dụng khi ,GPa,GPa IM 525 1  
21212 2301 IIMIII ,.,GPa  . .............................................. 77 
11 
Hình 3.6. Mô hình vật liệu cốt sợi dọc trục 3 pha ........................................... 79 
Hình 4.1. Ảnh vi mô có độ tương phản cao của mặt cắt ngang của vật liệu 
B/Al composite và RVE với phân bố tuần hoàn dạng lục giác ........................ 86 
Hình 4.2. Một số nhân tử tuần hoàn ................................................................. 87 
Hình 4.3. Sơ đồ khối chương trình tính toán bằng phương pháp PTHH ......... 90 
Hình 4.4. Một nửa nhân tử tuần hoàn mặt cắt ngang vật liệu cốt sợi đồng phương, 
(a) - mảng lục giác với cốt tròn, (b) - mảng hình vuông với cốt tròn .............. 91 
Hình 4.5. Trường nhiệt ..................................................................................... 92 
Hình 4.6. Hệ số dẫn hiệu dụng của mảng hình vuông, (a) - khi ,c,c 51 21 
;,c 323 20  (b) – khi ,c,c 520 21 323 1  ,c ............................... 95 
Hình 4.7. Hệ số dẫn hiệu dụng của mảng lục giác, (a) - khi ,c,c 51 21 
;,c 323 20  (b) – khi ,c,c 520 21 323 1  ,c .............................. 96 
Hình 4.8. (a) - nhân tử tuần hoàn với lớp vỏ bọc dị hướng, (b) – phần tử lớp 
 vỏ bọc ............................................................................................................. 97 
Hình 4.9. Hệ số dẫn  ...  stresses, Acta. Mech. , 2007, 188, 39-
54. 
[59] Vũ Lâm Đông, Đánh giá và mô phỏng mô đun đàn hồi vật liệu nhiều thành phần, 
Luận án tiến sĩ cơ kỹ thuật, 2016, Hà Nội. 
[60] Nguyễn Tiến Luật, Đánh giá các tính chất dẫn của vật liệu nhiều thành phần và 
đa tinh thể hỗn độn, Luận án tiến sĩ cơ kỹ thuật, 2017, Hà Nội. 
[61] Vương 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ĩ cơ học vật rắn, 2019, Hà Nội. 
[62] Đỗ Quốc Hoàng, Một cách tiếp cận xấp xỉ và mô hình hóa phần tử hữu hạn hệ 
số dẫn và mô đun đàn hồi của vật liệu nhiều thành phần, Luận án tiến sĩ cơ kỹ thuật, 
2019, Hà Nội. 
[63] K. J. Bathe, Finite element procedures, 1996, Prentice-Hall. 
[64] E. B. Becker, G. F. Carey and J. T. Oden, Finite elements: an introduction, 
1980, Prentice-Hall. 
[65] P. Wriggers, Nichtlineare Finite-Element-Methoden, 2001, Springer-Verlag. 
[66] T. J. R. Hughes, The finite element method, 1989, Prentice Hall. 
[67] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and 
nonlinear properties of composites, CR. Acad. Sc. Paris II, 1994, 318, 1417-1423. 
[68] Nguyen Van Luat, Nguyen Trung Kien, FFT-simulation and multi-coated 
inclusion model for macroscopic conductivity of 2D suspensions of compound 
inclusions, Vietnam Journal of Mechanics, 2015, 37, 169-176. 
126 
[69] Nguyễn Văn Luật, Nguyễn Trung Kiên, Phạm Đức Chính, Các đánh giá bậc ba 
và mô phỏng số FFT cho hệ số dẫn một số vật liệu nhiều thành phần, Hội nghị Khoa 
học toàn quốc Cơ học vật rắn biến dạng lần thứ XI Thành phố Hồ Chí Minh, 2013. 
[70] D.C. Pham and S.Torquato, Strong-contrast expansions and approximations for 
the effective conductivity of isotropic multiphase composites, J. Appl. Phys., 2003, 
94, 6591-6602. 
[71] S. Torquato, Random Heterogeneous Materials, Springer-Verlag, New York, 
2002. 
[72] H. Le-Quang, D.C. Pham, G. Bonnet, Q.C. He, Estimation of the effective 
conductivity of anisotropic multiphase composites with imperfect interfaces, Int. J. 
Heat. Mass Transf., 2013, 58, 175-187. 
[73] D.S. Liu, D.Y. Chiou, Modeling of inclusions with interphases in heterogeneous 
material using the infinite element method, Comput. Mater. Sci., 2004, 31, 405-420. 
[74] R. M. Christensen, Mechanics of Composite Materials, 1979, New York: Wiley. 
[75] D. C. Pham, A. B. Tran, Q. H. Do, On the effective medium approximations for 
the properties of isotropic multicomponent matrix-based composites, Int. J. Eng. Sci, 
2013, 68, 75-85. 
[76] A. S. Sarvestani, On the overall elastic moduli of composites with spherical 
coated fillers, Int. J. Solids Struct., 2003, 7553-7566. 
[77] S. Nemat-Nasser, M. Hori, Micromechanics: Overall Properties of 
Heterogeneous Materials, 1993, North-Holland, Amsterdam. 
[78] Z. Hashin, Viscoelastic fiber reinforced materials, 1966, AIAA Journal, 4(8), 
1411-1417. 
[79] George J. Dvorak, Micromechanics of Composite Materials, 2013, New York. 
[80] M. L. Accorsi, S. Nemat-Nasser, Bounds on the overall elastic and instantaneous 
elastoplastic moduli of periodic composites, Mechanics of Materials, 1986, 5, 209-
220. 
[81] J. L. Teply, J. N. Reddy, A unified formulation of micromechanics models of 
fiber-reinforced composites, In G. J. Dvorak (Ed.), Inelastic deformation of 
composite materials, 1990, 341-372. New York: Springer. 
127 
[82] J. Aboudi, Mechanics of composite materials – A unified micromechanical 
approach, 1991, Amsterdam: Elsevier. 
[83] S. Nemat – Nasser, M. Hori, Micromechanics: Overall properties of hetero-
geneous materials (2nd ed.), 1999, Amsterdam: Elsevier. 
[84] J. Fish, K. Shek, M. S. Shephard, M. Pandheeradi, Computational plasticity for 
composite structures based on mathematical homogenization: Theory and practice, 
Computer Methods in Applied Mechanics and Engineering, 1997, 157, 69-94. 
[85] L. Dormieux, D. Kondo, F. J. Ulm: Microporomechanics, 2006, Wiley, London. 
[86]  
PHỤ LỤC 
1. Tính toán mô đun đàn hồi trượt của cốt tương đương theo công thức (3.32) bằng phần mềm 
MATPLE. 
> restart; 
> ur:=A*R-6*(v/(1-2*v))*B*R^3+3*C/R^4+(5-4*v)/(1-2*v)*D/R^2: 
128 
up:=A*R-((7-4*v)/(1-2*v))*B*R^3-2*C/R^4+2*D/R^2: 
>sr:=2*(A+3*(v/(1-2*v))*B*R^2-12*C/R^5-2*(5-v)/(1-2*v)*D/R^3)*m: 
sp:=(A-((7+2*v)/(1-2*v))*B*R^2+8*C/R^5+2*(1+v)/(1-2*v)*D/R^3)*m: 
> ur11:=subs({R=R1,A=A1,B=B1,C=C1,D=D1,v=v1},ur): 
up11:=subs({R=R1,A=A1,B=B1,C=C1,D=D1,v=v1},up): 
ur21:=subs({R=R1,A=A2,B=B2,C=C2,D=D2,v=v2},ur): 
up21:=subs({R=R1,A=A2,B=B2,C=C2,D=D2,v=v2},up): 
> ur22:=subs({R=R2,A=A2,B=B2,C=C2,D=D2,v=v2},ur): 
up22:=subs({R=R2,A=A2,B=B2,C=C2,D=D2,v=v2},up): 
ur32:=subs({R=R2,A=A3,B=B3,C=C3,D=D3,v=v3},ur): 
up32:=subs({R=R2,A=A3,B=B3,C=C3,D=D3,v=v3},up): 
> sr11:=subs({R=R1,A=A1,B=B1,C=C1,D=D1,v=v1,m=m1},sr): 
sp11:=subs({R=R1,A=A1,B=B1,C=C1,D=D1,v=v1,m=m1},sp): 
sr21:=subs({R=R1,A=A2,B=B2,C=C2,D=D2,v=v2,m=m2},sr): 
sp21:=subs({R=R1,A=A2,B=B2,C=C2,D=D2,v=v2,m=m2},sp): 
> sr22:=subs({R=R2,A=A2,B=B2,C=C2,D=D2,v=v2,m=m2},sr): 
sp22:=subs({R=R2,A=A2,B=B2,C=C2,D=D2,v=v2,m=m2},sp): 
sr32:=subs({R=R2,A=A3,B=B3,C=C3,D=D3,v=v3,m=m3},sr): 
sp32:=subs({R=R2,A=A3,B=B3,C=C3,D=D3,v=v3,m=m3},sp): 
> ex1:=ur11-ur21: 
ex2:=up11-up21: 
ex3:=ur22-ur32: 
ex4:=up22-up32: 
ex5:=sr11-sr21: 
ex6:=sp11-sp21: 
ex7:=sr22-sr32: 
ex8:=sp22-sp32: 
> C1:=0;D1:=0;A3:=E0;B3:=0; 
>solus:=solve({ex1=0,ex2=0,ex3=0,ex4=0,ex5=0,ex6=0,ex7=0,ex8=0},{A1,B1,A2,B2,C2,D2,C3,
D3}): 
>A1:=eval(A1,solus);B1:=eval(B1,solus);A2:=eval(A2,solus);B2:=eval(B2,solus); 
C2:=eval(C2,solus);D2:=eval(D2,solus);C3:=eval(C3,solus);D3:=eval(D3,solus); 
Equivalent 
> ur2e:=subs({R=R2,A=Ae2,B=Be2,C=Ce2,D=De2,v=v2e},ur): 
up2e:=subs({R=R2,A=Ae2,B=Be2,C=Ce2,D=De2,v=v2e},up): 
ur3e:=subs({R=R2,A=Ae3,B=Be3,C=Ce3,D=De3,v=v3},ur): 
up3e:=subs({R=R2,A=Ae3,B=Be3,C=Ce3,D=De3,v=v3},up): 
129 
> sr2e:=subs({R=R2,A=Ae2,B=Be2,C=Ce2,D=De2,v=v2e,m=m2e},sr): 
sp2e:=subs({R=R2,A=Ae2,B=Be2,C=Ce2,D=De2,v=v2e,m=m2e},sp): 
sr3e:=subs({R=R2,A=Ae3,B=Be3,C=Ce3,D=De3,v=v3,m=m3},sr): 
sp3e:=subs({R=R2,A=Ae3,B=Be3,C=Ce3,D=De3,v=v3,m=m3},sp): 
> exe1:=ur2e-ur3e: 
exe2:=up2e-up3e: 
exe3:=sr2e-sr3e: 
exe4:=sp2e-sp3e: 
Ce2:=0;De2:=0;Ae3:=E0;Be3:=0: 
solue:=solve({exe1=0,exe2=0,exe3=0,exe4=0},{Ae2,Be2,Ce3,De3}): 
Equivallent 2 - 3 
> Meff2:=m3*(1+fe*Aed*(me-m3)/m3): 
> Meff1:=m3*(1+f1*A1d*(m1-m3)/m3+f2*A2d*(m2-m3)/m3): 
> Mtess:=m3*(1+fe*(me-m3)/(m3+beta*(me-m3))): beta:=6*(k3+2*m3)/(5*(3*k3+4*m3)): 
> simplify(solve(Mtess-Meff1=0,me)): 
> Aed:=(Ae2-21/5/(1-2*ve)*Be2*R2^2+4*(4-5*ve)/5/(1-2*ve)*De2/R2^3)/E0: 
> A2d:=simplify((A12d-(R1/R2)^3*A1d)/(1-(R1/R2)^3)): 
> A1d:=((A1-21/5/(1-2*v1)*B1*R1^2+4*(4-5*v1)/5/(1-2*v1)*D1/R1^3)/E0): 
> A12d:=((A2-21/5/(1-2*v2)*B2*R2^2+4*(4-5*v2)/5/(1-2*v2)*D2/R2^3)/E0): 
>Ae2:=eval(Ae2,solue);Be2:=eval(Be2,solue):Ce3:=eval(Ce3,solue):De3:=eval(De3,solue): 
> A1:=eval(A1,solus):B1:=eval(B1,solus):A2:=eval(A2,solus):B2:=eval(B2,solus): 
C2:=eval(C2,solus):D2:=eval(D2,solus):C3:=eval(C3,solus):D3:=eval(D3,solus): 
> Meq:=simplify(solve(Meff1-Meff2=0,me)); 
2. Mô hình tính hệ số dẫn ngang vật liệu cốt sợi dọc trục với cốt phức hợp (khi các pha đồng 
nhất đẳng hướng) – lập trình trên Cast3m 
*Periodic-Homogeneization of Conductivity Properties 
*2D Problem (2 spherical inclusion) 
option dime 2 elem tri3 echo 1; 
conm=1; N=0.05; conia=5; conib=20; 
m=1.; l=(3.)**(1./2.); ll=(-1)*l; 
listC=prog conm; listP=prog 0; nn=20; 
j=0; 
repere b1 nn; 
j=j+1; 
phii=0.05*j; phiib=(10./11.)*phii; phiia=phii-phiib; 
130 
a=((m*2.*l*phii/3.14)**(1./2.)); aa=(-1)*a; 
b=((m*2.*l*phiib/3.14)**(1./2.)); bb=(-1)*b; vt=2.*l*m; vib=(3.14*b*b); 
via=(3.14*a*a)-vib; vm=vt-via-vib; phiiaa=via/vt; phiibb=vib/vt; phim=vm/vt; 
p1=ll 0; p2=l 0; p3=l m; p4=ll m; p5=aa 0; p6=0 a; p7=a 0; p8=bb 0; p9=0 b; 
p10=b 0; p2a=l (m-a); p2b=l (m-b); p1a=ll (m-a); p1b=ll (m-b); p3a=(l-a) m; 
p3b=(l-b) m; p4a=(ll+a) m; p4b=(ll+b) m; 
 l15=droi p1 p5 dini (N) dfin (N); l72=droi p7 p2 dini (N) dfin (N); 
l22a=droi p2 p2a dini (N) dfin (N); l2a2b=droi p2a p2b dini (N) dfin (N); 
l2b3=droi p2b p3 dini (N) dfin (N); 
l23=l22a et l2a2b et l2b3; 
l11a=droi p1 p1a dini (N) dfin (N); l1a1b=droi p1a p1b dini (N) dfin (N); 
l1b4=droi p1b p4 dini (N) dfin (N); l14=l11a et l1a1b et l1b4; 
l33b=droi p3 p3b dini (N) dfin (N); l3b3a=droi p3b p3a dini (N) dfin (N); 
l3a4a=droi p3a p4a dini (N) dfin (N); l4b4a=droi p4b p4a dini (N) dfin (N); 
l44b=droi p4 p4b dini (N) dfin (N); l58=droi p5 p8 dini (N) dfin (N); 
l810=droi p8 p10 dini (N) dfin (N); l107=droi p10 p7 dini (N) dfin (N); 
l567=cer3 p5 p6 p7 dini (N) dfin (N); l8910=cer3 p8 p9 p10 dini (N) dfin (N); 
l3b=cerc p2b p3 p3b dini (N) dfin (N); l3a=cerc p2a p3 p3a dini (N) dfin (N); 
l4b=cerc p1b p4 p4b dini (N) dfin (N); l4a=cerc p1a p4 p4a dini (N) dfin (N); 
dm=l15 et l72 et l567 et l22a et l3a et l3a4a et l4a et l11a; 
sm= surf dm; dia12=l567 et l107 et l8910 et l58; 
sia12=surf dia12; dia3=l2a2b et l3b et l3b3a et l3a; sia3=surf dia3; 
dia4=l1a1b et l4b et l4b4a et l4a; sia4=surf dia4; sia=sia12 et sia3 et sia4; 
dib12=l8910 et l810; sib12=surf dib12; dib3=l2b3 et l33b et l3b; sib3=surf dib3; 
dib4=l1b4 et l44b et l4b; sib4=surf dib4; 
sib=sib12 et sib3 et sib4; s=sia et sm et sib; 
*trac s; 
modia=mode sia thermique; matia=mate modia k conia; modib=mode sib thermique; 
matib=mate modib k conib; modm=mode sm thermique; matm=mate modm k conm; 
cia=cond modia matia; cib=cond modib matib; cm=cond modm matm;cim=cia et cib et cm; 
*blocage des periodicites 
n1 = nbno l22a; 
i = 1; 
k22a = rela T (point 1 l22a) - T (point 1 l11a); 
dk22a=depi k22a 2.*((3.)**(1./2.)); 
repeter bouc (n1-1); 
i = i + 1; 
tata=(rela T (point i l22a) - T (point i l11a)); toto=depi tata 2.*((3.)**(1./2.)); 
131 
k22a=k22a et tata; dk22a=dk22a et toto; 
fin bouc; 
n1 = nbno l2a2b; 
i = 1; 
k2a2b = rela T (point 1 l2a2b) - T (point 1 l1a1b); 
dk2a2b=depi k2a2b 2.*((3.)**(1./2.)); 
repeter bouc (n1-1); 
i = i + 1; 
tata=(rela T (point i l2a2b) - T (point i l1a1b)); toto=depi tata 2.*((3.)**(1./2.)); 
k2a2b=k2a2b et tata; dk2a2b=dk2a2b et toto; 
fin bouc; 
 n1 = nbno l2b3; 
i = 1; 
k2b3 = rela T (point 1 l2b3) - T (point 1 l1b4); dk2b3=depi k2b3 2.*((3.)**(1./2.)); 
repeter bouc (n1-1); 
i = i + 1; 
tata=(rela T (point i l2b3) - T (point i l1b4)); toto=depi tata 2.*((3.)**(1./2.)); 
k2b3=k2b3 et tata; dk2b3=dk2b3 et toto; 
fin bouc; 
 k11=k2a2b et k22a et k2b3; dk11=dk2a2b et dk22a et dk2b3; 
 c1=bloq p9 T; dc1=depi c1 0; ctot=cia et cib et cm et c1 et k11; 
qtot=dk11 et dc1; r=reso ctot qtot; 
*trac r s; 
gria=grad modia r; grib=grad modib r; grm=grad modm r; 
xiia=intg modia gria t,x; xiib=intg modib grib t,x; xim=intg modm grm t,x; 
gria1=change chpo gria modia; gria2=vect gria1 0.2 t,x t,y rouge; 
grib1=change chpo grib modib; grib2=vect grib1 0.2 t,x t,y vert; 
grm1=change chpo grm modm; grm2=vect grm1 0.2 t,x t,y jaune; 
*trac (grm2 et gria2 et grib2) s; 
resm=(phim*conm*((xim/vm))); resa=(phiia*conia*((xiia/via))); 
resb=(phiib*conib*((xiib/vib))); 
res=resm+resa+resb; 
mess"result"; 
mess (phiiaa+phiibb) res; 
listC = listC 'ET' ('PROG' res); 
listP=listP et (prog phii); 
fin b1; 
132 
evo1 = 'EVOL' 'MANUEL' listP listC; dess evo1; list evo1; 
fin; 
3. Mô hình tính hệ số dẫn ngang vật liệu cốt sợi dọc trục với cốt phức hợp (khi lớp vỏ bọc dị 
hướng) – lập trình trên Cast3m 
*Periodic-Homogeneization of Conductivity Properties 
*2D Problem (2 spherical inclusion with anisotropic coated) 
option dime 2 mode plan elem tri3 echo 1; 
conm=100; 
conib=1; 
*kt:he so dan tiep, kn: he so dan vuong goc cua inlusion a. 
kn = 10; kt = 60; N=0.1; 
l=1.;ll=(-1)*l; 
phii= 0.5; phiib=(10./11.)*phii; phiia=phii-phiib; 
a=((4.*phii/3.14)**(1./2.)); aa=(-1)*a; b=((4.*phiib/3.14)**(1./2.)); 
bb=(-1)*b; vt=2*l*l; vib=(3.14*b*b)/2; via=((3.14*a*a)/2)-vib; 
vm=vt-via-vib; phiiaa=via/vt; phiibb=vib/vt; phim=vm/vt; 
p0=0 0; p1=ll 0; p2=l 0; p3=l l; p4=ll l; p5=aa 0; p6=0 a; p7=a 0; p8=bb 0; 
p9=0 b; p10=b 0; 
pe = table; po = table; loe = table; co = table; ce = table; da = table; sa = table; 
beta = table; modia = table; matia = table; k11d = table; k21d = table; k22d = table; 
cia = table; 
po.0 = p10; pe.0 = p7; loe.0 = droi po.0 pe.0 dini (N) dfin (N); 
*kk: so diem chia, kk tang gia tri se chinh xac, de kiem tra code, de kk gia tri nho. 
kk = 45; 
i = 0; 
repere b2 (2*kk); 
i = i+1; 
alpha = 90.0/kk; beta.i = (i*alpha)-(alpha/2.); 
po.i=(b*(cos(i*alpha))) (b*(sin(i*alpha))); pe.i=(a*(cos(i*alpha))) (a*(sin(i*alpha))); 
loe.i = droi po.i pe.i dini (N) dfin (N); co.i = cerc po.(i-1) p0 po.i dini (N) dfin (N); 
ce.i = cerc pe.(i-1) p0 pe.i dini (N) dfin (N); da.i = loe.(i-1) et ce.i et loe.i et co.i; 
sa.i = surf da.i; k11d.i = (kn*(cos(beta.i)**2))+(kt*(sin(beta.i)**2)); 
k22d.i = (kn*(sin(beta.i)**2))+(kt*(cos(beta.i)**2)); 
k21d.i = ((sin(beta.i))*(cos(beta.i)))*(kn-kt); 
toto = (kn*(cos(beta.i)**2))+(kt*(sin(beta.i)**2)); 
tata = (kn*(sin(beta.i)**2))+(kt*(cos(beta.i)**2)); 
tete = ((sin(beta.i))*(cos(beta.i)))*(kn-kt); 
modia.i=mode sa.i thermique anisotrope; 
133 
matia.i=mate modia.i direction (1. 0.) para k11 toto k21 tete k22 tata; 
cia.i=cond modia.i matia.i; 
fin b2; 
cot = co.1; cet = ce.1; sia = sa.1; ciat = cia.1; 
i = 1; 
repere b3 (2*kk-1); 
i = i+1; 
cot = cot et co.i; cet = cet et ce.i; sia = sia et sa.i; ciat=ciat et cia.i; 
fin b3; 
l15 = droit pe.(2*kk) p1 dini (N) dfin (N); l72 = droit pe.0 p2 dini (N) dfin (N); 
l810 = droit po.0 po.(2*kk) dini (N) dfin (N); 
dm=l15 et l72 et cet et l23 et l34 et l41; 
sm = surf dm; dib=cot et l810; sib=surf dib; s=sia et sm et sib; 
*trac s; 
modib=mode sib thermique; matib=mate modib k conib;modm=mode sm thermique; 
matm=mate modm k conm;cib=cond modib matib; cm=cond modm matm; 
cim=ciat et cib et cm; 
*blocage des periodicites 
n1 = nbno l23; i = 1; 
k11 = rela T (point 1 l23) - T (point 1 l41); 
dk11=depi k11 2; 
repeter bouc (n1-1); 
i = i + 1; 
tata=(rela T (point i l23) - T (point i l41)); 
toto=depi tata 2; k11=k11 et tata; dk11=dk11 et toto; 
fin bouc; 
c1=bloq p9 T; dc1=depi c1 0; 
ctot=ciat et cib et cm et c1 et k11; 
qtot=dk11 et dc1; r=reso ctot qtot; 
*trac r s; 
grib=grad modib r; grm=grad modm r; xiib=intg modib grib t,x; 
xim=intg modm grm t,x; resm=(phim*conm*((xim/vm))); 
resb=(phiib*conib*((xiib/vib))); 
gria=table; xiiax=table; xiiay=table; resaa=table; 
resa=0; 
i=0; 
repere b4 (2*kk); 
i=i+1; 
134 
gria.i=grad modia.i r; xiiax.i=intg modia.i gria.i t,x; xiiay.i=intg modia.i gria.i t,y; 
resai=(phiia/(2*kk))*(((k11d.i*xiiax.i)+(k21d.i*xiiay.i))/(via/(2*kk))); 
resai=(phiia/(2*kk))*(((k11d.i*xiiax.i)+(k21d.i*xiiay.i))/(via/(2*kk))); 
resa=resa+resai; resaa.i=resai; 
fin b4; 
res=resm+resb+resa; 
mess"result"; 
mess (phiiaa+phiibb) res; 
-fin;