Tổng quan nghiên cứu

Trong kỹ thuật công trình hiện đại, kết cấu dạng tấm chiếm hơn 65% tỷ trọng các cấu kiện chịu lực chính trong các công trình công nghiệp, cầu thép và kết cấu hàng không. Sự xuất hiện của các vết nứt tế vi hoặc khuyết tật hình học trong quá trình gia công và khai thác có thể làm suy giảm nghiêm trọng độ cứng uốn, khiến kết cấu bị mất ổn định đàn hồi dưới tải trọng nén thấp hơn từ 15% đến 45% so với tải trọng thiết kế ban đầu. Bài toán ổn định tấm uốn đàn hồi đặt ra thách thức lớn khi các phương pháp giải tích cổ điển không thể giải quyết trọn vẹn cho các miền có hình học phức tạp hoặc vết nứt bất liên tục, trong khi phương pháp phần tử hữu hạn truyền thống (FEM) đòi hỏi phải tái lập lưới phần tử liên tục tại lân cận đỉnh nứt, làm tiêu tốn từ 70% đến 80% thời gian xử lý dữ liệu.

Nghiên cứu tập trung phát triển giải pháp số mô phỏng bài toán mất ổn định đàn hồi của tấm Reissner - Mindlin có vết nứt bằng phương pháp phần tử hữu hạn mở rộng (XFEM) kết hợp phần tử đẳng tham số tứ giác 8 nút (Q8). Mục tiêu cụ thể là xác định chính xác hệ số tải trọng tới hạn $k$, khảo sát sự biến thiên của các dạng mất ổn định (buckling modes) dưới tác động của nhiều điều kiện biên, vị trí, kích thước và góc nghiêng của vết nứt. Nghiên cứu được thực hiện tại Trường Đại học Bách khoa – ĐHQG TP. Hồ Chí Minh trong giai đoạn từ tháng 01/2014 đến tháng 01/2015, tập trung vào lĩnh vực Xây dựng công trình Dân dụng và Công nghiệp. Kết quả đạt được mang ý nghĩa khoa học và thực tiễn sâu sắc khi giảm thiểu thời gian tạo lưới về 0 đối với các bài toán nứt phức tạp, đồng thời kiểm soát sai số tương đối dưới 1.2% so với nghiệm giải tích và phần mềm thương mại ABAQUS.

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

Khung lý thuyết áp dụng

Nghiên cứu được xây dựng dựa trên 2 nền tảng lý thuyết chủ đạo: Lý thuyết tấm Reissner - Mindlin (lý thuyết biến dạng trượt bậc nhất) và Lý thuyết phân chia đơn vị trong phương pháp phần tử hữu hạn mở rộng (XFEM). Không giống như lý thuyết tấm mỏng cổ điển Kirchhoff bỏ qua biến dạng trượt, mô hình Reissner - Mindlin xét đến góc xoay của đường pháp tuyến trung hòa quanh trục $x$ và $y$ kết hợp thành phần trượt trung bình, áp dụng hệ số hiệu chỉnh lực cắt $\alpha = 5/6$ (hoặc $0.833$) để biểu diễn trạng thái ứng suất tiếp phân bố bậc hai theo chiều dày tấm $h$.

Khung lý thuyết vận hành dựa trên 4 khái niệm cốt lõi:

  1. Hàm bước nhảy Heaviside $H(x)$: Nhận giá trị $+1$ hoặc $-1$ để mô tả sự gián đoạn dịch chuyển của vật liệu qua bề mặt thân vết nứt.
  2. Bộ 4 hàm làm giàu tiệm cận đỉnh nứt $B_\alpha(r, \theta)$: Dẫn xuất từ nghiệm giải tích trường ứng suất kỳ dị bậc $1/\sqrt{r}$ quanh đầu nứt trong tọa độ cực.
  3. Kỹ thuật Level Set và phương pháp hình học: Xác định chính xác vị trí vết nứt thông qua 2 hàm khoảng cách $\psi(x)$ và $\Phi(x)$ để tự động kích hoạt các nút làm giàu mà không cần chỉnh sửa lưới gốc.
  4. Hiện tượng khóa cắt (Shear Locking): Trạng thái ma trận độ cứng phần tử bị cứng hóa giả tạo khi tỷ số chiều dày trên nhịp $h/a$ tiến dần về giới hạn tấm mỏng $0.001$.

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

Nghiên cứu sử dụng nguồn dữ liệu mô phỏng số thông qua chương trình tính toán tự phát triển trên ngôn ngữ lập trình MATLAB, được đối chuẩn trực tiếp với nghiệm giải tích kinh điển của Euler/Mindlin và kết quả mô phỏng từ phần mềm phần tử hữu hạn ABAQUS. Không gian bài toán khảo sát tổng cộng 32 bộ tham số thử nghiệm khác nhau, bao gồm các cấu hình tấm vuông ($a/b = 1.0$), tấm chữ nhật ($a/b = 2.0$), chiều dày tương đối biến thiên từ $h/a = 0.001$ đến $0.2$, cùng 3 nhóm điều kiện biên điển hình (4 cạnh tựa đơn SSSS, 4 cạnh ngàm CCCC, 2 cạnh ngàm 2 cạnh tựa SCSC).

Thuật toán phân tích lựa chọn phần tử đẳng tham số tứ giác 8 nút (Q8) với 24 bậc tự do cơ sở cho mỗi phần tử chuẩn (3 bậc tự do $w, \theta_x, \theta_y$ trên mỗi nút). Khi vết nứt cắt qua, phần tử thân nứt được làm giàu thêm 1 bậc tự do tại mỗi nút (tổng cộng 6 bậc tự do/nút), trong khi phần tử chứa đỉnh nứt được bổ sung 4 hàm nhánh tiệm cận (tổng cộng 15 bậc tự do/nút). Để triệt tiêu hoàn toàn hiện tượng khóa cắt, nghiên cứu áp dụng kỹ thuật tích phân suy giảm có chọn lọc theo sơ đồ Gauss: tích phân ma trận độ cứng uốn dùng lưới $3 \times 3$ điểm Gauss (9 điểm), còn tích phân thành phần biến dạng cắt dùng lưới $2 \times 2$ điểm Gauss (4 điểm). Ma trận độ cứng hình học được tích hợp qua sơ đồ $3 \times 3$ điểm Gauss để giải bài toán giá trị riêng tìm hệ số lực tới hạn $k_{cr}$. Toàn bộ quy trình tính toán được thực hiện theo lộ trình 10 tháng với các bài toán kiểm tra độ hội tụ lưới nghiêm ngặt từ $10 \times 10$ đến $40 \times 40$ phần tử.

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

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

Nghiên cứu đã đưa ra các kết quả định lượng chi tiết về hành vi mất ổn định của tấm nứt thông qua 4 phát hiện quan trọng:

  1. Ảnh hưởng của vị trí vết nứt: Vết nứt xuất hiện ở biên làm suy giảm tải trọng tới hạn nghiêm trọng hơn so với vết nứt nằm tại tâm tấm. Cụ thể, khi tỷ lệ chiều dài vết nứt $c/a$ tăng từ $0.1$ lên $0.5$ ở tấm tựa đơn chịu nén 1 phương ($h/a = 0.01$), hệ số lực tới hạn $k$ đối với vết nứt ở cạnh giảm từ $3.82$ xuống $2.35$ (mức sụt giảm đạt $38.48%$), trong khi vết nứt ở tâm chỉ giảm từ $3.91$ xuống $2.94$ (sụt giảm $24.81%$).
  2. Ảnh hưởng của phương vết nứt so với hướng tác dụng lực: Vết nứt nằm ngang (vuông góc với phương lực nén) gây nguy hiểm cục bộ cao hơn vết nứt đứng (song song với phương lực nén). Khi $c/a = 0.4$, hệ số lực tới hạn của tấm có nứt ngang thấp hơn tấm có nứt đứng từ $18.5%$ đến $27.3%$ tùy thuộc vào điều kiện biên.
  3. Tác động của chiều dày tương đối $h/a$: Khi tỷ số $h/a$ tăng từ $0.001$ (tấm rất mỏng) lên $0.2$ (tấm dày) trên tấm nguyên vẹn tựa đơn, hệ số lực tới hạn $k$ giảm từ giá trị giải tích $4.000$ xuống còn $3.002$ (giảm $24.95%$) do biến dạng trượt ngang chiếm ưu thế, khẳng định tính tất yếu của việc áp dụng mô hình Reissner - Mindlin.
  4. Mức độ kiểm soát của điều kiện biên: Tấm có 4 biên ngàm (CCCC) cho khả năng duy trì ổn định cao nhất, với hệ số $k$ đạt khoảng $9.45$ ở trạng thái nứt biên $c/a = 0.2$, cao gấp $2.68$ lần so với tấm 4 biên tựa đơn ($k = 3.52$) trong cùng điều kiện hình học và vết nứt.

Thảo luận kết quả

Căn nguyên vật lý của sự suy giảm tải trọng tới hạn bắt nguồn từ việc vết nứt làm đứt gãy tính liên tục của trường chuyển vị, dẫn đến sự tập trung ứng suất rất lớn tại vùng lân cận đỉnh nứt và làm giảm mômen kháng uốn hữu hiệu của tiết diện mặt cắt ngang. Sự biến thiên của các hệ số lực tới hạn $k$ và các đường cong phân nhánh mất ổn định được biểu diễn trực quan thông qua các đồ thị quan hệ giữa $k$ và $c/a$, kết hợp cùng biểu đồ phân bố trường chuyển vị của các mode mất ổn định thứ nhất đến thứ tư.

So sánh với các nghiên cứu công bố trước đây sử dụng phần tử Q4R, MITC4 hoặc phương pháp không lưới CS-DSG3, phương pháp XFEM kết hợp phần tử Q8 tích phân suy giảm có chọn lọc trong luận văn này cho thấy tốc độ hội tụ vượt trội. Sai số số học so với phần mềm ABAQUS dao động ổn định trong khoảng $0.45%$ đến $1.15%$, nhưng loại bỏ hoàn toàn yêu cầu phải chia mịn hàng nghìn phần tử tam giác xung quanh đỉnh nứt. Bảng dữ liệu kiểm tra độ hội tụ trên các lưới $10 \times 10$, $20 \times 20$ và $30 \times 30$ chứng minh rằng lưới $20 \times 20$ phần tử Q8 đã đạt độ chính xác tương đương với lưới $60 \times 60$ phần tử bậc một truyền thống, giúp tiết kiệm hơn $60%$ dung lượng bộ nhớ giải ma trận thưa.

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

Nhằm chuyển hóa các kết quả nghiên cứu thành ứng dụng thực tiễn trong ngành kỹ thuật xây dựng và cơ học kết cấu, các giải pháp cụ thể được kiến nghị như sau:

  1. Xây dựng và tích hợp mô-đun XFEM-Q8 vào các phần mềm tính toán kết cấu kỹ thuật: Các kỹ sư phát triển phần mềm CAE cần lập trình gói thư viện mở rộng cho MATLAB hoặc Python trong vòng 12 tháng, hướng tới mục tiêu giảm $75%$ thời gian thiết lập mô hình nứt kết cấu dạng tấm.
  2. Chuẩn hóa quy trình thẩm tra độ ổn định kết cấu tấm có khuyết tật: Viện Khoa học Công nghệ Xây dựng và các hội nghề nghiệp cần bổ sung các chỉ dẫn kỹ thuật đánh giá hệ số suy giảm tải trọng $k$ cho các tấm thép và bê tông cốt thép trong thời hạn 18 tháng, nâng cao độ tin cậy an toàn công trình thêm ít nhất $20%$.
  3. Áp dụng quy tắc tích phân suy giảm $3 \times 3$ và $2 \times 2$ Gauss trong phân tích tấm dày: Các kỹ sư tư vấn thiết kế kết cấu cần áp dụng ngay sơ đồ tích phân này khi mô phỏng các tấm có tỷ số $h/a \ge 0.1$, kiểm soát sai số khóa cắt dưới ngưỡng $1.5%$ cho mọi bài toán kiểm tra độ ổn định cục bộ.
  4. Mở rộng khung thuật toán sang bài toán vết nứt động và vật liệu FGM: Các nhóm nghiên cứu sau đại học tại các trường đại học khối kỹ thuật cần phát triển tiếp mô hình cho kết cấu vật liệu biến tính chức năng FGM chịu tải trọng va đập hoặc dao động điều hòa trong lộ trình 24 tháng tới, giúp tăng khả năng dự báo phá hủy mỏi sớm hơn $30%$ chu kỳ tải.

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

Nội dung và kết quả của luận văn mang lại giá trị tham khảo chuyên môn cao cho 4 nhóm đối tượng chính:

  1. Kỹ sư kết cấu công trình dân dụng, công nghiệp và cầu đường: Ứng dụng trực tiếp các biểu đồ tra cứu hệ số tải trọng tới hạn $k$ để đánh giá nhanh độ an toàn chịu nén của các vách thép, bản sàn đáy và tấm dầm khi phát hiện có vết nứt bề mặt trong thực tế kiểm định.
  2. Giảng viên, nghiên cứu viên và học viên cao học chuyên ngành Cơ học kỹ thuật: Sử dụng toàn bộ quy trình thiết lập ma trận độ cứng $K$ và ma trận độ cứng hình học $K_G$ cho phần tử Q8 kết hợp hàm Heaviside và hàm nhánh làm tài liệu giảng dạy và mở rộng nghiên cứu bài toán phi tuyến hình học.
  3. Các đơn vị tư vấn giám định và kiểm định chất lượng công trình: Sử dụng phương pháp định lượng suy giảm cường độ tới hạn (với mức giảm từ $15%$ đến $40%$) để đưa ra quyết định gia cường kết cấu chính xác, tránh các biện pháp sửa chữa tốn kém không cần thiết.
  4. Kỹ sư phát triển thuật toán phần mềm mô phỏng số (CAE Developers): Khai thác thuật toán phát hiện vết nứt tự động bằng phương pháp hình học và kỹ thuật Level Set để tối ưu hóa mô-đun chia lưới tự động trong các phần mềm tính toán kết cấu thương mại.

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

Ưu điểm nổi bật nhất của phương pháp XFEM so với FEM truyền thống là gì? XFEM cho phép đưa các hàm gián đoạn Heaviside và hàm tiệm cận đỉnh nứt vào trường chuyển vị thông qua phân chia đơn vị. Nhờ đó, vết nứt hoàn toàn độc lập với hệ lưới phần tử, giúp loại bỏ $100%$ công đoạn chia lại lưới phức tạp khi vết nứt xuất hiện hoặc phát triển trong tấm.

Tại sao nghiên cứu lại chọn phần tử tứ giác 8 nút Q8 kết hợp lý thuyết Reissner - Mindlin? Lý thuyết Reissner - Mindlin xét đến biến dạng trượt ngang với hệ số hiệu chỉnh $5/6$, phù hợp cho cả tấm mỏng và tấm dày. Phần tử Q8 với phép nội suy bậc hai giúp mô tả chính xác trường ứng suất uốn phức tạp mà không đòi hỏi mật độ lưới quá dày đặc như phần tử 4 nút Q4.

Làm thế nào để luận văn khắc phục hoàn toàn hiện tượng khóa cắt trong bài toán tấm mỏng? Nghiên cứu ứng dụng kỹ thuật tích phân suy giảm có chọn lọc theo phương pháp cầu phương Gauss. Cụ thể, ma trận biến dạng uốn được tích phân chính xác với sơ đồ $3 \times 3$ điểm Gauss, trong khi ma trận biến dạng cắt được giảm bậc tích phân xuống $2 \times 2$ điểm Gauss, loại trừ sai số cứng hóa giả tạo.

Vết nứt ở cạnh hay ở tâm tấm gây nguy hiểm lớn hơn cho độ ổn định của kết cấu? Vết nứt ở cạnh nguy hiểm hơn đáng kể. Kết quả khảo sát định lượng cho thấy với cùng tỷ lệ chiều dài nứt $c/a = 0.5$, hệ số lực tới hạn $k$ của tấm có nứt cạnh giảm đến $38.48%$, trong khi nứt tâm chỉ làm giảm $24.81%$ khả năng chịu tải do biên tự do làm tăng độ linh động uốn cục bộ.

Độ chính xác của chương trình tính toán viết bằng MATLAB được kiểm chứng như thế nào? Kết quả tính toán của chương trình được đối chứng trực tiếp với 2 nguồn tin cậy: nghiệm giải tích kinh điển của Euler/Mindlin (sai số đạt dưới $0.2%$ trên tấm nguyên vẹn) và mô hình phần tử hữu hạn trên phần mềm ABAQUS, với sai số tương đối tối đa không vượt quá $1.2%$.

Kết luận

  • Luận văn đã xây dựng thành công thuật toán và chương trình phần tử hữu hạn mở rộng (XFEM) kết hợp phần tử Q8 cho bài toán phân tích ổn định tấm Reissner - Mindlin có vết nứt.
  • Khắc phục triệt để hiện tượng khóa cắt nhờ kỹ thuật tích phân suy giảm có chọn lọc ($3 \times 3$ cho uốn và $2 \times 2$ cho cắt), đảm bảo độ chính xác cao cho cả tấm mỏng ($h/a = 0.001$) và tấm dày ($h/a = 0.2$).
  • Xác lập bộ dữ liệu số chi tiết về 32 kịch bản phân tích ổn định, chứng minh vết nứt ở cạnh và vết nứt vuông góc với lực nén làm suy giảm tải trọng tới hạn mạnh nhất, lên tới $38.48%$.
  • Khẳng định tính ưu việt của công cụ tính toán tự phát triển trên MATLAB với sai số dưới $1.2%$ so với ABAQUS mà không cần tái tạo lưới phần tử.
  • Đóng góp giải pháp tính toán hiệu quả, mở đường cho việc tích hợp mô-đun XFEM vào quy trình kiểm tra an toàn kết cấu công trình trong lộ trình 12 đến 24 tháng tới.

Hãy ứng dụng phương pháp phần tử hữu hạn mở rộng ngay hôm nay để nâng cao độ chính xác và tối ưu hóa thời gian tính toán phân tích ổn định kết cấu công trình của bạn!