Bối cảnh và vấn đề nghiên cứu

Báo cáo Đồ án Project 3 được thực hiện bởi nhóm sinh viên Khoa Kỹ thuật Giao thông, Trường Đại học Bách Khoa – Đại học Quốc gia TP. Hồ Chí Minh dưới sự hướng dẫn của TS. Nguyễn Song Thanh Thảo (hoàn thành ngày 10 tháng 12 năm 2021). Đề tài tập trung vào việc ứng dụng Phương pháp Phần tử Hữu hạn (Finite Element Method - FEM) để giải quyết các bài toán kỹ thuật nền tảng trong cơ học kết cấu, truyền nhiệt và động lực học.

Trong kỹ thuật giao thông và cơ khí hiện đại, việc phân tích trường nhiệt độ trên các kết cấu tản nhiệt, đánh giá ứng suất – biến dạng của vật liệu composite nhiều lớp và khảo sát dao động tự nhiên của các chi tiết máy đòi hỏi sự kết hợp chặt chẽ giữa lý thuyết giải tích số và công cụ mô phỏng chuyên dụng. Đề tài được xây dựng nhằm mục tiêu kiểm chứng độ chính xác giữa phương pháp tính toán giải tích số bằng tay (hand calculations) và mô phỏng số thông qua phần mềm ANSYS.

Nhiệm vụ nghiên cứu cụ thể của đề tài bao gồm ba bài toán kỹ thuật độc lập:

  1. Bài toán truyền nhiệt (Heat Transfer Problems - Problem 10.4): Xác định quy luật phân bố nhiệt độ và lượng tổn thất nhiệt trên cánh tản nhiệt dạng trụ tròn (pin fin) bằng phương pháp giải tay FEM với các mô hình 2 phần tử (1 bậc tự do), 3 phần tử (1 bậc tự do), 2 phần tử (2 bậc tự do) và đối chứng với mô hình mô phỏng dòng nhiệt trên phần mềm ANSYS Fluent.
  2. Bài toán kết cấu và vật liệu composite (Composite Materials and Structures): Khảo sát trường ứng suất (stress) và biến dạng (strain) trên tấm vật liệu composite nhiều lớp laminate cấu hình xếp lớp $[0/90]$ chịu tải trọng kéo dọc trục $P = 8\text{ kN}$ bằng phương pháp giải tích phần tử hữu hạn và mô phỏng ANSYS Static Structural.
  3. Bài toán phân tích động lực học kết cấu (Dynamic Analysis - Problem 11): Xây dựng ma trận độ cứng và ma trận khối lượng toàn cục cho bài toán dao động dọc trục của thanh đồng (copper bar); xác định tần số dao động tự nhiên thấp nhất (lowest natural frequency) và dạng dao động tương ứng (mode shape) bằng thuật toán lặp nghịch đảo (inverse iteration algorithm).

Đối tượng nghiên cứu của đồ án là các phần tử hình học mẫu: cánh tản nhiệt trụ tròn (đường kính $d = 0.5\text{ in}$, chiều dài $L = 6\text{ in}$), tấm laminate composite $[0/90]$, và thanh đồng chịu dao động dọc trục. Phạm vi nghiên cứu giới hạn ở trạng thái truyền nhiệt ổn định (steady-state), kết cấu đàn hồi tuyến tính tĩnh và dao động tự do không cản.

Cơ sở lý thuyết và phương pháp

Đồ án ứng dụng khung lý thuyết giải tích số phần tử hữu hạn kết hợp công cụ mô phỏng kỹ thuật số trên nền tảng ANSYS Workbench.

Khung lý thuyết áp dụng

  • Phương pháp Phần tử Hữu hạn (FEM): Rời rạc hóa miền liên tục thành các phần tử 1D/2D; xây dựng các hàm dạng (shape functions), thiết lập ma trận độ cứng/ma trận dẫn nhiệt phần tử và vector tải phần tử; ghép nối (assembly) thành hệ phương trình đại số toàn cục $[K]{U} = {F}$ hoặc $[K_T]{T} = {Q}$.

  • Lý thuyết truyền nhiệt 1D trên cánh tản nhiệt: Phương trình vi phân truyền nhiệt kết hợp dẫn nhiệt dọc trục và đối lưu nhiệt qua bề mặt xung quanh:

    $$[K_c + K_h]{T} = {F_Q + F_h}$$

    Trong đó $[K_c]$ là ma trận dẫn nhiệt, $[K_h]$ là ma trận tỏa nhiệt đối lưu qua bề mặt, ${F_h}$ là vector nguồn nhiệt do đối lưu với môi trường xung quanh.

  • Lý thuyết cơ học vật liệu composite nhiều lớp: Áp dụng ma trận chuyển đổi tọa độ $[T]$ và ma trận nghịch đảo chuyển vị $[T]^{-T}$ để biến đổi độ cứng của từng lớp đơn hướng (lamina) ở các góc định hướng sợi $0^\circ$ và $90^\circ$ về hệ tọa độ chung của tấm laminate.

  • Lý thuyết dao động kết cấu và thuật toán lặp: Phương trình dao động tự do không cản $[K]{X} = \omega^2 [M]{X}$. Áp dụng thuật toán lặp nghịch đảo (Inverse Iteration Algorithm) để tìm trị riêng nhỏ nhất $\lambda = \omega_1^2$ và vector riêng tương ứng ${X_1}$.

Phương pháp nghiên cứu và nguồn dữ liệu

  • Phương pháp giải tích số giải tay (Hand calculations): Thiết lập thủ công các ma trận phần tử, thực hiện khử bậc tự do theo điều kiện biên (Boundary Conditions) và giải hệ đại số tuyến tính.
  • Phương pháp mô phỏng số (Numerical Simulation): Sử dụng các module chuyên biệt trong hệ sinh thái ANSYS:
    • ANSYS Fluent / Mesh: Mô phỏng trường phân bố nhiệt độ và dòng nhiệt đối lưu.
    • ANSYS Composite PrepPost (ACP) & Static Structural: Mô phỏng kết cấu tấm composite, định nghĩa đặc tính lớp ép (fabric, stackup, rosette, oriented selection set).
  • Nguồn dữ liệu thông số đầu vào:
    • Bài toán truyền nhiệt: Đường kính fin $d = 0.5\text{ in}$, chiều dài $L = 6\text{ in}$, nhiệt độ gốc $T_{\text{root}} = 147^\circ\text{F}$, nhiệt độ môi trường $T_\infty = 80^\circ\text{F}$, đầu đỉnh cánh cách nhiệt ($q_{\text{tip}} = 0$).
    • Bài toán composite: Lực kéo dọc trục $P = 8\text{ kN}$, cấu hình $[0/90]$, kích thước mô hình theo đề bài.
    • Bài toán dao động: Đặc tính vật liệu thanh đồng, ngàm cố định một đầu (chuyển vị tại nút gốc bằng 0).

Thiết kế và triển khai mô phỏng kỹ thuật

Quá trình mô phỏng trên phần mềm ANSYS được nhóm tác giả chuẩn hóa qua các bước thiết lập hình học, chia lưới và gán điều kiện biên cụ thể:

Bài toán kỹ thuật Module ANSYS sử dụng Kích thước phần tử lưới (Element Size) Thiết lập điều kiện biên và tải trọng
Truyền nhiệt cánh tản nhiệt ANSYS Mesh & Fluid Flow (Fluent) $1.0\text{ mm}$ (chia lưới tự động) - Gốc (base): Nhiệt độ cố định $T = 147^\circ\text{F}$
- Bề mặt xung quanh (Surface) & đỉnh (End-face): Đối lưu nhiệt với không khí (air) ở $T_\infty = 80^\circ\text{F}$
Kết cấu tấm composite ANSYS ACP & Static Structural $0.4\text{ mm}$ - Gốc dưới tấm: Cố định hoàn toàn (Fixed Support)
- Mặt biên trên: Tác dụng lực kéo $P = 8\text{ kN}$ dọc trục Y
Động lực học thanh đồng ANSYS Modal / Dynamic Analysis Chia lưới tự động theo cấu trúc thanh - Gốc nút 1: Cố định chuyển vị dọc trục ($u_1 = 0$)
- Cố định biến dạng theo phương Y và phương Z

Nội dung chính theo từng phần

Phần 1: Phân tích bài toán truyền nhiệt (Heat Transfer Problems)

Nhóm tác giả thực hiện giải bài toán truyền nhiệt cánh tản nhiệt Problem 10.4 qua 4 phương pháp tính toán khác nhau nhằm khảo sát ảnh hưởng của số lượng phần tử và bậc tự do:

                  [Pin Fin: d = 0.5 in, L = 6 in]
     Root (x=0)                                        Tip (x=6 in)
       |==================================================|
    T = 147°F             Đối lưu xung quanh            Cách nhiệt
                        Môi trường T_inf = 80°F
  • Phương pháp giải tay 2 phần tử (1 bậc tự do/nút): Hệ được chia thành 2 phần tử chiều dài $L_e = 3\text{ in}$. Thiết lập ma trận dẫn nhiệt $[K_c]$, ma trận đối lưu $[K_h]$ và vector nhiệt ${F}$. Sau khi khử điều kiện biên tại gốc ($T_1 = 147^\circ\text{F}$), nhiệt độ tính được tại vị trí nút giữa ($x = 3\text{ in}$) đạt $94.3^\circ\text{F}$. Tổng năng lượng nhiệt tổn thất ra môi trường được tổng hợp từ hai phần tử.
  • Phương pháp giải tay 3 phần tử (1 bậc tự do/nút): Chiều dài mỗi phần tử giảm xuống $L_e = 2\text{ in}$. Ma trận toàn cục mở rộng thành cấp $4 \times 4$. Sau khi áp dụng điều kiện biên tại nút 1, kết quả nhiệt độ tại vị trí $x = 2\text{ in}$ đạt $105.5^\circ\text{F}$.
  • Phương pháp giải tay 2 phần tử (2 bậc tự do/nút): Áp dụng hàm dạng bậc cao để tăng độ chính xác phân bố nhiệt dọc theo chiều dài phần tử.
  • Mô phỏng bằng ANSYS Fluent: Dựng hình học 3D, đặt tên các mặt biên (base, Surface, end-face), chia lưới kích thước $1\text{ mm}$. Kết quả mô phỏng trích xuất trường nhiệt độ liên tục dọc theo 20 điểm trên chiều dài cánh và đối chứng tại các vị trí nút rời rạc.

Phần 2: Vật liệu và kết cấu composite (Composite Materials and Structures)

Nghiên cứu tập trung vào phân tích trạng thái ứng suất và biến dạng của tấm laminate $[0/90]$ chịu lực kéo $P = 8\text{ kN}$:

  • Phương pháp giải tay FEM:
    • Rời rạc hóa kết cấu thành 2 phần tử tam giác/tứ giác liên kết qua 4 nút tọa độ.
    • Xây dựng ma trận chuyển đổi tọa độ $[T]$ cho góc xếp lớp $0^\circ$ và $90^\circ$.
    • Thiết lập ma trận độ cứng phần tử thứ nhất $K^{(1)}$ và phần tử thứ hai $K^{(2)}$, sau đó ghép nối thành ma trận độ cứng toàn cục $[K]$.
    • Thiết lập vector lực toàn cục ${F}$, giải hệ phương trình tìm trường chuyển vị tại các nút tự do.
    • Tính toán ứng suất ($\sigma_x, \sigma_y, \tau_{xy}$) và biến dạng ($\varepsilon_x, \varepsilon_y, \gamma_{xy}$) riêng biệt cho từng lớp đặt $0^\circ$ và $90^\circ$.
  • Mô phỏng bằng ANSYS:
    • Thiết lập dữ liệu vật liệu trong Engineering Data.
    • Tạo hình dạng tấm và chia lưới phần tử với kích thước $0.4\text{ mm}$.
    • Định nghĩa cấu trúc composite: tạo lớp vải (fabric), thiết lập góc đan (stackup), gán hệ tọa độ cục bộ (rosette), xác định hướng sợi (oriented selection set) cho từng lớp ép $0^\circ$ và $90^\circ$.
    • Áp dụng điều kiện biên trong Static Structural: ngàm cố định đáy tấm, đặt lực kéo $8\text{ kN}$ theo phương Y tại mép trên. Trích xuất biểu đồ phân bố biến dạng (Strain) và ứng suất (Stress).
                      Lực kéo P = 8 kN (Phương Y)
                         +-----------------+
                         |  Lớp ép [0°]    |
                         |  Lớp ép [90°]   | (Lưới: 0.4 mm)
                         +-----------------+
                          Ngàm cố định gốc

Phần 3: Phân tích động lực học kết cấu (Dynamic Analysis)

Phần này giải quyết bài toán dao động dọc trục của thanh đồng (Problem 11):

  • Thiết lập hệ ma trận toàn cục: Xây dựng ma trận độ cứng $[K]$ và ma trận khối lượng $[M]$ của hệ thanh chịu dao động dọc trục.
  • Khử điều kiện biên ngàm: Do chuyển vị tại nút ngàm số 1 bằng 0 ($u_1 = 0$), nhóm tác giả tiến hành khử hàng 1 và cột 1 trên cả ma trận độ cứng và ma trận khối lượng, đưa về hệ phương trình rút gọn cho các bậc tự do chuyển vị còn lại ($u_2, u_3$).
  • Áp dụng thuật toán lặp nghịch đảo (Inverse Iteration Algorithm):
    • Vòng lặp 1: Khởi tạo vector thử nghiệm (guess vector) ${X_0}$, giải hệ phương trình tìm vector chuyển vị mới, tính toán xấp xỉ tần số riêng $\omega_1^{(1)}$ và chuẩn hóa vector dạng dao động.
    • Vòng lặp 2: Sử dụng vector đã chuẩn hóa từ vòng 1 để tiếp tục quá trình lặp, tính toán độ hội tụ của giá trị riêng $\omega_1^{(2)}$ và dạng dao động chuẩn (mode shape).
  • Kiểm chứng bằng ANSYS: Thiết lập mô hình phân tích dao động, cố định biến dạng phương Y và phương Z, trích xuất kết quả chuyển vị và tần số dao động tự nhiên tại các nút.

Kết quả và đóng góp

Kết quả tính toán và mô phỏng

Nhóm tác giả đã tổng hợp và so sánh chi tiết kết quả phân bố nhiệt độ tại các vị trí nút giữa các phương pháp:

Vị trí khảo sát trên cánh tản nhiệt Giải tay FEM 2 phần tử (1 bậc tự do) Giải tay FEM 3 phần tử (1 bậc tự do) Giải tay FEM 2 phần tử (2 bậc tự do) Mô phỏng ANSYS Fluent
Gốc cánh ($x = 0\text{ in}$) $147.0^\circ\text{F}$ $147.0^\circ\text{F}$ $147.0^\circ\text{F}$ $147.0^\circ\text{F}$
Vị trí $x = 2\text{ in}$ — $105.5^\circ\text{F}$ — Tương thích đồ thị
Vị trí $x = 3\text{ in}$ $94.3^\circ\text{F}$ — — Tương thích đồ thị
Đỉnh cánh ($x = 6\text{ in}$) — — — $85.35^\circ\text{F}$

Trong bài toán composite và động lực học, kết quả giải tay các thành phần ứng suất, biến dạng và tần số dao động tự nhiên bằng thuật toán lặp nghịch đảo cho thấy sự trùng khớp cao với các trường dữ liệu đồ họa trích xuất từ phần mềm ANSYS.

Đóng góp kỹ thuật

  • Xây dựng hoàn chỉnh quy trình giải tích số từng bước cho ba dạng bài toán cơ bản trong kỹ thuật: truyền nhiệt, composite nhiều lớp và động lực học kết cấu.
  • Chứng minh tính khả thi và độ chính xác của mô hình phần tử hữu hạn giải tay bậc thấp khi đối chiếu với công cụ mô phỏng 3D chuyên sâu (độ sai lệch giữa các phương pháp là rất nhỏ).
  • Xác định rõ nguyên nhân sai số phát sinh chủ yếu do quá trình rút gọn, làm tròn số học trong các phép toán ma trận giải tay và giả thiết đơn giản hóa miền hình học.

Hạn chế và hướng nghiên cứu tiếp

  • Hạn chế:
    • Các bài toán giải tay sử dụng số lượng phần tử còn ít (2 đến 3 phần tử), giả định trường biến dạng và truyền nhiệt 1D/2D đơn giản.
    • Bài toán truyền nhiệt giả định điều kiện đỉnh cánh cách nhiệt lý tưởng, chưa xét đến bức xạ nhiệt và biến thiên phi tuyến của hệ số đối lưu.
    • Quá trình lặp nghịch đảo trong bài toán động lực học mới dừng lại ở các vòng lặp đầu tiên để tìm mode dao động cơ bản, chưa mở rộng cho các mode bậc cao.
  • Hướng nghiên cứu tiếp:
    • Mở rộng mô hình giải tích số cho kết cấu composite có góc xếp lớp phức tạp hơn (ví dụ: $[0/45/-45/90]_s$) và khảo sát phá hủy lớp (delamination).
    • Tăng mật độ chia lưới trong mô phỏng ANSYS để đánh giá độ hội tụ của lưới (mesh convergence study).
    • Khảo sát dao động cưỡng bức có cản và phản ứng động lực học phi tuyến của kết cấu.

Giá trị tham khảo

Báo cáo đồ án là tài liệu tham khảo trực quan cho:

  • Sinh viên ngành Kỹ thuật Giao thông, Kỹ thuật Cơ khí, Kỹ thuật Hàng không và Cơ kỹ thuật: Tham khảo phương pháp thiết lập ma trận phần tử hữu hạn, cách áp dụng điều kiện biên và quy trình khử bậc tự do ma trận.
  • Người học và nghiên cứu phần mềm ANSYS: Cung cấp hướng dẫn từng bước thiết lập bài toán truyền nhiệt trong ANSYS Fluent, cấu hình vật liệu nhiều lớp trong ANSYS Composite PrepPost (ACP) và phân tích động lực học kết cấu.
  • Học phần Phương pháp Phần tử Hữu hạn và Đồ án chuyên ngành: Làm tài liệu mẫu cho việc trình bày bài toán đối chứng giữa giải tích lý thuyết và mô phỏng số.

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

1. Đồ án khảo sát bài toán truyền nhiệt trên đối tượng nào và điều kiện biên ra sao?

Đồ án khảo sát cánh tản nhiệt dạng trụ tròn (pin fin) có đường kính $d = 0.5\text{ in}$ và chiều dài $L = 6\text{ in}$. Điều kiện biên gồm: nhiệt độ tại gốc $T_{\text{root}} = 147^\circ\text{F}$, nhiệt độ môi trường đối lưu $T_\infty = 80^\circ\text{F}$, và đỉnh cánh tản nhiệt được giả thiết cách nhiệt hoàn toàn ($q_{\text{tip}} = 0$).

2. Kích thước lưới phần tử được thiết lập như thế nào trong các mô phỏng ANSYS?

Trong mô phỏng bài toán truyền nhiệt cánh tản nhiệt (ANSYS Fluent), kích thước phần tử lưới (Element size) được đặt là $1.0\text{ mm}$. Trong bài toán phân tích kéo nén tấm composite (ANSYS Static Structural), kích thước phần tử lưới được chia mịn ở mức $0.4\text{ mm}$.

3. Trong bài toán composite, các lớp vật liệu được định hướng như thế nào và chịu tải trọng gì?

Tấm composite được cấu tạo theo dạng nhiều lớp laminate với góc xếp sợi $[0/90]$. Tấm chịu tác dụng của tải trọng kéo dọc trục $P = 8\text{ kN}$ tác dụng theo phương Y tại biên trên, trong khi phần gốc đáy được cố định hoàn toàn.

4. Thuật toán lặp nghịch đảo (Inverse Iteration Algorithm) được sử dụng nhằm mục đích gì trong đồ án?

Thuật toán lặp nghịch đảo được sử dụng trong Phần 3 để giải bài toán dao động dọc trục của thanh đồng. Sau khi khử bậc tự do tại nút ngàm cố định, thuật toán sử dụng vector thử nghiệm (guess vector) để lặp ma trận qua các vòng, từ đó xác định tần số dao động tự nhiên thấp nhất (lowest natural frequency) và dạng dao động (mode shape).

5. Tại sao có sự chênh lệch kết quả giữa tính toán giải tay FEM và mô phỏng phần mềm ANSYS?

Nhóm tác giả chỉ ra rằng sự chênh lệch (rất nhỏ) xuất phát từ các sai số công thức giải tích, việc làm tròn số thập phân trong các bước biến đổi ma trận giải tay, cũng như sự khác biệt giữa mô hình rời rạc hóa 1D giải tay và mô hình phần tử khối 3D được chia lưới mịn trong ANSYS.

Kết luận

Báo cáo Đồ án Project 3 của nhóm sinh viên Trường Đại học Bách Khoa – ĐHQG TP.HCM đã giải quyết trọn vẹn ba bài toán kỹ thuật trọng tâm: truyền nhiệt trên cánh pin fin, ứng suất – biến dạng trên tấm composite $[0/90]$ và dao động dọc trục của thanh đồng. Nghiên cứu đã kết hợp thành công phương pháp tính toán giải tích số phần tử hữu hạn giải tay với công cụ mô phỏng số chuyên dụng ANSYS (Fluent, ACP, Static Structural). Kết quả đối chứng giữa các phương pháp cho thấy độ tin cậy cao của các mô hình tính toán, cung cấp cơ sở tài liệu học thuật giá trị cho việc học tập và nghiên cứu phương pháp phần tử hữu hạn trong kỹ thuật.