Tổng quan nghiên cứu

Trong kỹ thuật xây dựng công trình dân dụng và công nghiệp hiện đại, sự xuất hiện của các dòng vật liệu composite và vật liệu đa pha đã tạo ra bước đột phá về hiệu quả chịu lực, giúp tăng tỷ lệ tối ưu trọng lượng kết cấu từ 30% đến 50% so với vật liệu đồng nhất truyền thống. Tuy nhiên, việc mô hình hóa chính xác ứng xử cơ học của các cấu trúc dị tính vi mô đối mặt với thách thức lớn về chi phí tính toán. Nếu sử dụng phương pháp phần tử hữu hạn tiêu chuẩn để chia lưới toàn bộ chi tiết vi mô cho một kết cấu công trình lớn, hệ phương trình số học có thể chạm ngưỡng hàng triệu bậc tự do, gây quá tải năng lực xử lý của máy tính tính toán thông thường.

Nhằm giải quyết bài toán hóc búa này, đề tài tập trung nghiên cứu phương pháp tính toán đồng nhất hóa vật liệu đa tỉ lệ (Multi-scale computational homogenization) kết hợp cùng kỹ thuật phần tử hữu hạn trơn dựa trên cạnh (Edge-based Smoothed Finite Element Method - ES-FEM). Đề tài được hoàn thành tại Trường Đại học Bách Khoa – Đại học Quốc gia TP. Hồ Chí Minh trong khoảng thời gian từ tháng 01/2013 đến tháng 06/2013.

Mục tiêu cốt lõi của nghiên cứu là xây dựng quy trình tính toán hai cấp độ vi mô - vĩ mô hiệu quả cho vật liệu đàn hồi tuyến tính biến dạng nhỏ. Điểm đột phá của công trình là đề xuất việc áp dụng linh hoạt ES-FEM tại cấp độ vĩ mô, giúp khắc phục triệt để hiện tượng quá cứng (overly stiff) của phần tử hữu hạn tam giác 3 nút cổ điển. Kết quả thực nghiệm số chỉ ra rằng phương pháp kết hợp tối ưu giúp cải thiện độ chính xác chuyển vị lên tới 28% so với phương pháp FE2 truyền thống, đồng thời tiết kiệm thời gian xử lý nhanh gấp 9 lần so với việc áp dụng kỹ thuật làm trơn đồng thời ở cả hai cấp độ.

Cơ sở lý thuyết và phương pháp nghiên cứu

Khung lý thuyết áp dụng

Nền tảng của nghiên cứu được xây dựng trên sự giao thoa giữa hai lý thuyết tính toán cơ học tiên tiến: Lý thuyết đồng nhất hóa tính toán đa tỉ lệ bậc nhất (First-order Multi-scale Homogenization) do Kouznetsova phát triển và Phương pháp phần tử hữu hạn trơn trên cạnh (ES-FEM) do G. Liu đề xuất năm 2008.

Mô hình nghiên cứu vận hành dựa trên 5 khái niệm then chốt:

  1. Phần tử thể tích đại diện (Representative Volume Element - RVE): Cấu trúc hình học ở cấp vi mô mang đầy đủ đặc trưng đại diện cho tính chất không đồng nhất của vật liệu.
  2. Tính tuần hoàn cục bộ (Local Periodicity): Giả định sát thực tế cho phép các vi cấu trúc biến đổi linh hoạt giữa các điểm vĩ mô khác nhau nhưng vẫn lặp lại trong một lân cận vi mô xác định.
  3. Điều kiện biên tuần hoàn (Periodic Boundary Conditions - PBC): Ràng buộc biến dạng trên các cạnh đối xứng của RVE nhằm đảm bảo tính liên tục của trường ứng suất và biến dạng.
  4. Điều kiện đồng nhất năng lượng Hill - Mandel: Nguyên lý cân bằng công cơ học, yêu cầu biến thiên mật độ công nội lực trung bình trên thể tích RVE vi mô phải tương đương với công biến dạng tại điểm tích phân vĩ mô.
  5. Miền làm trơn biến dạng dựa trên cạnh: Miền cục bộ được tạo thành bằng cách kết nối hai nút của một cạnh tam giác với hai trọng tâm của hai phần tử kề bên, cho phép tính toán biến dạng trung bình thông qua tích phân đường biên.

Quy trình tính toán chuyển đổi hai chiều diễn ra khép kín: ten-sơ biến dạng tại mỗi điểm vĩ mô đóng vai trò điều kiện biên động học áp đặt lên phần tử RVE vi mô; sau khi giải bài toán vi mô với kỹ thuật tích phân chỉ cần 1 điểm Gauss trên miền trơn, ten-sơ ứng suất vĩ mô và ma trận độ cứng tương đương 6x6 được đồng nhất hóa để cập nhật ngược lại cho toàn bộ kết cấu công trình.

Phương pháp nghiên cứu

Nghiên cứu sử dụng phương pháp mô phỏng số học (numerical simulation) kết hợp phương pháp đối chuẩn (benchmarking) với các nghiệm giải tích chính xác. Nguồn dữ liệu kiểm chứng được xây dựng từ các mô hình RVE 2D chuẩn và các bài toán kết cấu thực tế trong cơ học công trình.

Mẫu nghiên cứu vi mô bao gồm:

  • Mẫu RVE 2 pha vật liệu: Pha nền có mô đun đàn hồi E1 = 70 GPa, hệ số Poisson 0.2; pha gia cường có E2 = 700 GPa, hệ số Poisson 0.3, khảo sát dưới 3 dạng bố trí hình học (phân bố theo phương x, phương y và phân bố xen kẽ).
  • Mẫu RVE composite thực tế: Ma trận nhựa Epoxy (E = 3.13 GPa, hệ số Poisson 0.34) kết hợp sợi thủy tinh (E = 73 GPa, hệ số Poisson 0.2) với tỷ lệ thể tích sợi thay đổi từ 5% đến 35%.
  • Mẫu RVE chứa lỗ rỗng: Khảo sát 1 lỗ rỗng tại tâm, 2 lỗ rỗng và 4 lỗ rỗng đối xứng với tỷ lệ thể tích rỗng chiếm từ 10%, 30% đến 50% thể tích RVE.

Mẫu kết cấu vĩ mô bao gồm:

  • Dầm công-xôn ngàm một đầu chịu lực tập trung P = 10 kN tại đầu tự do, kích thước chiều dài L = 4m, chiều cao D = 1m, bề rộng B = 0.1m.
  • Tấm phẳng vô hạn có lỗ khoét tròn bán kính a = 1m chịu lực kéo vô cực q = 1 N/m, khai thác tính đối xứng để mô phỏng một phần tư tấm có kích thước 5m x 5m với chiều dày t = 1m.

Phương pháp chia lưới sử dụng phần tử tam giác 3 nút (T3) linh hoạt, sau đó tự động tái tạo hệ thống miền trơn không chồng lấn theo số cạnh của lưới. Toàn bộ thuật toán đồng nhất hóa, khử bậc tự do phụ thuộc và giải hệ phương trình ma trận độ cứng được lập trình độc lập trên môi trường phần mềm tính toán MATLAB trong suốt 6 tháng thực hiện đề tài. Lý do lựa chọn phân tích ES-FEM kết hợp FEM là nhằm tận dụng khả năng tính toán nhanh của FEM ở vi mô và tính mềm dẻo khử khóa chuyển vị của ES-FEM ở cấp độ kết cấu vĩ mô.

Kết quả nghiên cứu và thảo luận

Những phát hiện chính

Quá trình mô phỏng số và phân tích đối chuẩn đã đem lại các phát hiện khoa học quan trọng:

Thứ nhất, quy luật phân bố pha vi mô quyết định mạnh mẽ đến độ cứng tương đương của vật liệu. Với cùng tỷ lệ thể tích pha gia cường có độ cứng gấp 10 lần pha nền, kiểu phân bố xen kẽ nâng mô đun đàn hồi tương đương lên mức 1.24x10^8 kPa (tăng 77% so với pha nền ban đầu) và mô đun cắt đạt 5.12x10^7 kPa (tăng 75%). Trong khi đó, việc bố trí theo phương x hoặc phương y chỉ giúp tăng mô đun đàn hồi khoảng 64% (đạt 1.15x10^8 kPa), chứng minh hình thái vi mô xen kẽ tạo ra liên kết truyền lực tối ưu nhất.

Thứ hai, mức độ suy giảm cơ tính phi tuyến tính theo thể tích và hình thái lỗ rỗng. Khi tỷ lệ lỗ rỗng đạt mức 50%, mẫu RVE chứa 1 lỗ rỗng tại tâm ghi nhận sự sụt giảm nghiêm trọng nhất: mô đun đàn hồi giảm tới 83% (chỉ còn 1.18x10^7 kPa) và mô đun cắt giảm 83% (còn 4.88x10^6 kPa). Ngược lại, nếu chia tách độ rỗng 50% thành 2 lỗ rỗng thì mức giảm là 69%, và nếu phân tán thành 4 lỗ rỗng nhỏ thì độ cứng chỉ giảm khoảng 28.8%. Điều này khẳng định sự tập trung khuyết tật tại trung tâm gây phá hủy độ cứng nguy hiểm hơn nhiều so với khuyết tật phân tán.

Thứ ba, tính chuẩn xác của điều kiện biên tuần hoàn trong dự báo composite. Đối với vật liệu Epoxy gia cường sợi thủy tinh, khi tỷ lệ thể tích sợi dưới ngưỡng 20%, kết quả mô đun cắt tương đương tính toán theo mô hình hoàn toàn trùng khớp với đường cong nghiệm giải tích của Nemat-Nasser với sai số kiểm định dưới 2%.

Thứ tư, xác lập cấu hình tính toán đa tỉ lệ tối ưu vượt trội. Trong bài toán dầm công-xôn chịu uốn, phương pháp ESFEM-M (áp dụng ES-FEM ở vĩ mô và FEM ở vi mô) cho độ võng đầu dầm lớn hơn 28% so với phương pháp FE2 truyền thống, tiệm cận chính xác với nghiệm giải tích Timoshenko. Về mặt tài nguyên máy tính, thời gian chạy của ESFEM-M chỉ tăng 1.6 lần so với FE2, trong khi phương pháp áp dụng ES-FEM ở cả hai cấp độ (ESFEM-M-m) làm tăng thời gian tính toán lên tới 14.4 lần nhưng độ lệch kết quả chuyển vị chỉ chênh lệch khoảng 4%.

Thảo luận kết quả

Nguyên nhân cốt lõi giúp ESFEM-M đạt được độ chính xác cao nằm ở cơ chế làm mềm ma trận độ cứng của kỹ thuật làm trơn cạnh. Phương pháp FEM tiêu chuẩn với phần tử tam giác bậc thấp thường áp đặt trường biến dạng tương thích quá ngặt nghèo, khiến kết cấu bị cứng giả tạo. Khi chuyển sang ES-FEM ở cấp vĩ mô, việc tính toán biến dạng dựa trên tích phân đường biên của các miền trơn cục bộ quanh cạnh đã giải phóng năng lượng biến dạng, giúp mô hình phản ánh đúng độ võng thực tế của kết cấu dầm và sự tập trung ứng suất quanh lỗ tròn tấm phẳng.

Dữ liệu nghiên cứu được biểu diễn trực quan thông qua các đồ thị phân bố chuyển vị dọc trục trung hòa của dầm và đường cong chuyển vị các nút biên x(y=0) của tấm phẳng vô hạn. Đồ thị so sánh cho thấy đường cong chuyển vị của ESFEM-M bám sát hoàn hảo với nghiệm giải tích của Timoshenko và Liu, trong khi đường kết quả của FEM cổ điển bị lệch đáng kể về phía dưới do hiện tượng khóa cứng.

So sánh với các nghiên cứu của Lê Văn Cảnh và cộng sự năm 2011 về phương pháp FE2 thuần túy, mô hình cải tiến trong luận văn này đã khắc phục hoàn toàn nhược điểm tính toán chậm mà vẫn giữ được độ chính xác tương đương các phương pháp phần tử cấp cao. Bằng chứng số học từ bảng so sánh thời gian khẳng định việc duy trì FEM ở cấp độ vi mô (nơi chỉ cần giải các bài toán RVE đơn giản với ma trận nhỏ) kết hợp ES-FEM ở cấp độ vĩ mô là sự đánh đổi hoàn hảo giữa độ chính xác cơ học và hiệu năng tính toán.

Đề xuất và khuyến nghị

Dựa trên các kết quả đạt được, luận văn đưa ra 4 khuyến nghị và giải pháp kỹ thuật ứng dụng thực tiễn:

  1. Tối ưu hóa cấu trúc vật liệu xốp trong xây dựng: Các doanh nghiệp sản xuất bê tông bọt khí và vật liệu cách nhiệt nhẹ nên điều chỉnh quy trình công nghệ để phân tán các bọt khí thành nhiều lỗ rỗng siêu nhỏ đường kính dưới 2mm thay vì để hình thành các lỗ rỗng lớn tập trung. Giải pháp này giúp duy trì mô đun đàn hồi của cấu kiện cao hơn từ 40% đến 54% dù cùng đạt tỷ lệ thể tích rỗng 50%.
  2. Chuẩn hóa quy trình mô phỏng hybrid ESFEM-M trong phần mềm kỹ thuật: Các kỹ sư kết cấu nên áp dụng mô hình lai ESFEM-M khi phân tích các kết cấu composite nhiều pha. Việc tích hợp này giúp tiết kiệm khoảng 85% thời gian chạy mô phỏng so với mô hình ES-FEM toàn phần mà vẫn đảm bảo sai số độ võng dưới ngưỡng 4% cho phép của tiêu chuẩn thiết kế.
  3. Tự động hóa thuật toán khử bậc tự do điều kiện biên tuần hoàn: Đề xuất các nhóm nghiên cứu phát triển các module lập trình mở rộng trên ngôn ngữ Python hoặc C++ nhằm tự động nhận dạng cặp nút biên đối xứng trên lưới RVE bất kỳ, giúp rút ngắn thời gian thiết lập ma trận độ cứng tương đương xuống dưới 5 giây cho mỗi phần tử đại diện.
  4. Mở rộng khung tính toán sang miền phi tuyến và nứt gãy: Các cơ quan nghiên cứu và trường đại học cần triển khai lộ trình trong 12 đến 24 tháng tới để mở rộng thuật toán ESFEM-M cho các bài toán phi tuyến hình học, biến dạng dẻo lớn và mô phỏng lan truyền vết nứt vi mô trong bê tông cốt sợi.

Đối tượng nên tham khảo luận văn

Luận văn là tài liệu tham khảo học thuật và kỹ thuật giá trị cao cho 4 nhóm đối tượng chuyên môn:

  1. Học viên cao học và nghiên cứu sinh chuyên ngành Cơ kỹ thuật, Kỹ thuật Xây dựng: Cung cấp tài liệu hoàn chỉnh về lý thuyết đồng nhất hóa đa tỉ lệ bậc nhất, nguyên lý bảo toàn Hill - Mandel và các bước triển khai chi tiết phương pháp phần tử hữu hạn trơn ES-FEM.
  2. Kỹ sư R&D trong ngành công nghiệp vật liệu mới và composite: Cung cấp công cụ mô phỏng số đáng tin cậy để dự báo trước các chỉ tiêu cơ lý (mô đun đàn hồi, mô đun cắt) của vật liệu nhiều pha mà không cần thực hiện hàng loạt thí nghiệm phá hủy tốn kém chi phí.
  3. Chuyên gia phát triển phần mềm tính toán kết cấu (CAE/FEA): Nắm bắt kiến trúc thuật toán rút gọn ma trận độ cứng 6x6, kỹ thuật tích phân biên miền trơn 1 điểm Gauss và phương pháp gán điều kiện biên tuần hoàn phục vụ phát triển các solver thương mại hoặc mã nguồn mở.
  4. Giảng viên các trường đại học khối kỹ thuật: Sử dụng làm bài giảng chuyên đề, tài liệu tham khảo cho các môn học Cơ học tính toán nâng cao, Phương pháp phần tử hữu hạn mở rộng và Cơ học vật liệu dị tính.

Câu hỏi thường gặp

  1. Tại sao điều kiện biên tuần hoàn (PBC) lại vượt trội hơn điều kiện biên chuyển vị và lực kéo? Điều kiện biên tuần hoàn đảm bảo tính tương thích biến dạng và cân bằng ứng suất trên các mặt đối xứng của phần tử RVE mà không gây ra hiện tượng bó cứng nhân tạo như điều kiện biên chuyển vị, hay gây mất ổn định số học như điều kiện biên lực kéo. Kết quả kiểm chứng trên mẫu composite Epoxy - sợi thủy tinh cho thấy sai số của PBC so với nghiệm giải tích chỉ dưới 2%.

  2. Cơ chế nào giúp phương pháp ES-FEM khắc phục hiện tượng quá cứng của phần tử tam giác cổ điển? Thay vì tính đạo hàm hàm dạng trực tiếp trên từng phần tử tam giác gây cứng lưới, ES-FEM chia miền bài toán thành các miền trơn cục bộ dựa trên cạnh và áp dụng kỹ thuật làm trơn biến dạng. Quá trình này biến đổi tích phân diện tích thành tích phân đường biên, làm mềm ma trận độ cứng K và giúp nghiệm chuyển vị hội tụ tiệm cận chính xác nghiệm giải tích.

  3. Vì sao cấu hình kết hợp ESFEM-M lại được khuyến nghị thay vì chạy ES-FEM ở cả hai cấp độ? Việc chạy ES-FEM ở cả hai cấp độ vi mô và vĩ mô đòi hỏi tái tạo miền trơn lặp lại hàng nghìn lần tại mọi điểm Gauss, khiến thời gian tính toán tăng vọt gấp 14.4 lần so với phương pháp FE2 thông thường. Cấu hình ESFEM-M tận dụng FEM nhanh gọn ở vi mô và ES-FEM mềm dẻo ở vĩ mô, cho kết quả chỉ lệch 4% so với mô hình toàn phần nhưng tốc độ xử lý nhanh hơn gấp 9 lần.

  4. Sự suy giảm độ cứng khi vật liệu xuất hiện lỗ rỗng phụ thuộc vào những yếu tố nào? Sự suy giảm cơ tính phụ thuộc đồng thời vào phần trăm thể tích rỗng và cấu trúc hình học của lỗ rỗng. Khi thể tích rỗng chiếm 50%, một lỗ rỗng đơn lẻ tại trung tâm làm giảm tới 83% mô đun đàn hồi, trong khi việc phân tán thành 4 lỗ rỗng đối xứng chỉ làm giảm 28.8% độ cứng nhờ giảm thiểu mức độ tập trung ứng suất cục bộ.

  5. Phương pháp đồng nhất hóa đa tỉ lệ trong nghiên cứu có thể mở rộng cho bài toán 3D không? Hoàn toàn có thể mở rộng. Nguyên lý chuyển đổi vi mô - vĩ mô và điều kiện biên tuần hoàn vẫn giữ nguyên bản chất toán học. Khi chuyển sang không gian 3 chiều, các miền trơn cạnh sẽ được nâng cấp thành các miền trơn mặt (Face-based Smoothed Finite Element Method - FS-FEM) trên lưới tứ diện, giải quyết hiệu quả các kết cấu không gian phức tạp.

Kết luận

  • Đã tích hợp thành công phương pháp phần tử hữu hạn trơn dựa trên cạnh (ES-FEM) vào bài toán đồng nhất hóa vật liệu đa tỉ lệ bậc nhất, mở ra hướng tiếp cận mới trong phân tích vật liệu dị tính.
  • Xác lập mô hình lai ESFEM-M là cấu hình tối ưu nhất, giúp tăng 28% độ chuẩn xác chuyển vị so với FE2 tiêu chuẩn và rút ngắn thời gian tính toán gấp 9 lần so với mô hình ES-FEM toàn phần.
  • Định lượng chi tiết quy luật suy giảm độ bền của cấu trúc có lỗ rỗng: khuyết tật phân tán 4 điểm giúp bảo toàn độ cứng cao hơn 54% so với khuyết tật tập trung tại tâm ở cùng độ rỗng 50%.
  • Kiểm chứng thành công tính chính xác của thuật toán lập trình trên MATLAB qua các đối chuẩn nghiệm giải tích kinh điển của Timoshenko, Nemat-Nasser và Liu.
  • Đặt nền móng vững chắc cho lộ trình 12 tháng tiếp theo nhằm mở rộng công cụ mô phỏng cho bài toán nứt gãy và vật liệu phi tuyến 3D trong thực tế xây dựng.

Các đơn vị tư vấn thiết kế và kỹ sư kết cấu có thể ứng dụng ngay giải pháp mô phỏng ESFEM-M để tối ưu hóa thiết kế vật liệu composite, cắt giảm chi phí thử nghiệm và nâng cao tuổi thọ công trình.