Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set

Chia sẻ: Dạ Thiên Lăng | Ngày: | Loại File: PDF | Số trang:10

Thêm vào BST

Báo xấu

3
lượt xem 0
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Trong nghiên cứu "Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set", một phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất (LS-PFEM) mới được sử dụng để giải quyết bài toán level set (LS) đối lưu-khuếch tán không ổn định. Để ổn định các nghiệm số mà không có các dao động phi vật lý (gradient dốc), cả quá trình tiến hóa và tái khởi tạo của hàm LS đều được giải bằng cách tối thiểu hóa điều kiện ổn định bình phương nhỏ nhất. Mời các bạn cùng tham khảo!

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set

401 Tuyển tập công trình Hội nghị Cơ học toàn quốc lần thứ XI, Hà Nội, 02-03/12/2022 Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set Nguyễn Trần Bá Đình1, Nguyễn Hoàng Sơn2,3 và Phan Đức Huynh4,* 1 Công ty TNHH công nghệ và kỹ thuật Nhật Bản, Việt Nam 2 Viện Khoa học Tính toán, Trường Đại học Tôn Đức Thắng, Việt Nam 3 Khoa xây dựng, Trường Đại học Tôn Đức Thắng, Việt Nam 4 Khoa xây dựng, Trường Đại học Sư phạm Kỹ thuật TP.HCM, Việt Nam *Email: huynhpd@hcmute.edu.vn Tóm tắt. Trong nghiên cứu này, một phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất (LS-PFEM) mới được sử dụng để giải quyết bài toán level set (LS) đối lưu-khuếch tán không ổn định. Để ổn định các nghiệm số mà không có các dao động phi vật lý (gradient dốc), cả quá trình tiến hóa và tái khởi tạo của hàm LS đều được giải bằng cách tối thiểu hóa điều kiện ổn định bình phương nhỏ nhất. Phương pháp này cung cấp các tính chất toán học tốt như sự khuếch tán số tự nhiên và tính đối xứng xác định dương của các hệ phương trình đại số thu được. So với phần tử tam giác (T3) và tứ giác (Q4) thông thường, phần tử đa giác có khả năng cung cấp tính linh hoạt cao hơn trong việc tạo lưới cho các mô hình bài toán phức tạp cũng như mang lại kết quả tính toán chính xác hơn. Trong bài báo này, phương pháp đề xuất được áp dụng để khảo sát một số bài toán chuẩn như: vòng quay của đĩa Zalesak, dòng chảy của xoáy đơn đảo ngược theo thời gian. Kết quả chỉ ra rằng cách tiếp cận được đề xuất có hiệu quả cao với tỷ lệ hội tụ tuyệt vời so với các phần tử T3 và Q4. Từ khóa: Phần tử đa giác, Phương pháp level set, Phương pháp bình phương nhỏ nhất, Phương trình đối lưu-khuếch tán, Tái khởi tạo. 1. Mở đầu Trong những thập kỷ qua, phương pháp level set (LSM) đã được phát triển và áp dụng cho nhiều lĩnh vực kỹ thuật, chẳng hạn như hình học tính toán, tối ưu hóa, xử lý hình ảnh và động lực học Hamilton-Jacobi cho một hàm phụ 𝜙𝜙 có giá trị hàm LS bằng 0 xác định cho hình dạng của giao diện chất lỏng tính toán. Ý tưởng cơ bản của phương pháp này là giải một phương trình đạo hàm riêng LS (interface) tự do [1, 2]. Trong quá trình giải một bài toán LS điển hình bằng LSM, hàm LS ban đầu được khởi tạo như một hàm khoảng cách có dấu (SDF) sao cho nó thỏa mãn phương trình Eikonal [2]. Trong suốt quá trình phát triển của giao diện, hàm LS không đảm bảo đặc tính SDF và trở nên quá phẳng hoặc quá dốc do các sai số tích lũy. Vấn đề này có thể được giải quyết bằng cách sử dụng các kỹ thuật tái khởi tạo [3, 4, 5, 6]. Sussman và cộng sự (1994) [3] đã trình bày một phương pháp tái khởi tạo trong đó một phương trình đạo hàm riêng hypebol được giải đến trạng thái ổn định. Kết quả của họ cho thấy rằng hàm LS cần phải được tái khởi tạo sau mỗi bước thời gian để giữ cho nghiệm chính xác. Tuy nhiên, phương pháp này dẫn đến sự dịch chuyển đáng kể của giao diện và thời gian đạt đến trạng thái ổn định là quá lâu. Elias và cộng sự (2007) [6] đã đề xuất một phương pháp mới để tính toán các hàm khoảng cách trong lưới phi cấu trúc bởi áp đặt sự thỏa mãn của phương trình Eikonal ở cấp phần tử. Phương pháp này dễ thực hiện và có thể được sử dụng dễ dàng trong các bộ giải phần tử hữu hạn vì tất cả thông tin cần thiết đều có sẵn hoặc được xây dựng sẵn dưới dạng các cấu trúc dữ liệu. Lưu ý rằng hầu hết các nghiên cứu trên đều tập trung vào các phần tử T3 và Q4 để giải quyết các bài toán LS. Trong khi đó, cả hai phần tử này đều tồn tại những ưu, nhược điểm nhất định. Các phần tử T3 thích hợp để tạo lưới có dạng hình học phức tạp, nhưng độ chính xác và độ hội tụ của
402 Nguyễn Trần Bá Đình, Nguyễn Hoàng Sơn và Phan Đức Huynh nghiệm thấp. Ngược lại, độ chính xác và độ hội tụ nghiệm của phần tử Q4 cao nhưng khả năng tạo lưới cho các mô hình bài toán phức tạp của nó vẫn còn nhiều hạn chế. Trong những năm gần đây, phương pháp phần tử hữu hạn đa giác với rất nhiều tính năng tốt đã được ứng dụng rộng rãi trong các bài toán cơ học. Tuy nhiên, việc sử dụng nó vào LSM vẫn còn nhiều hạn chế. So với các phần tử T3 và Q4, phần tử đa giác có khả năng cung cấp tính linh hoạt cao hơn trong việc tạo lưới cho các mô hình bài toán phức tạp cũng như mang lại kết quả tính toán chính xác hơn [7, 8]. Trong khuôn khổ của phương pháp phần tử hữu hạn (FEM) cho LSM, sử dụng phương pháp Galerkin tiêu chuẩn để rời rạc miền không gian là phương pháp đơn giản và phổ biến nhất. Tuy nhiên, phương pháp này dễ gây ra hiện tượng dao dộng của nghiệm [4]. Để giải quyết vấn đề này, phương pháp bình phương nhỏ nhất tiêu chuẩn [1, 4] sẽ được áp dụng trong nghiên cứu này. Điều khác biệt ở đây là chúng ta không chỉ giới hạn trong việc sử dụng phần tử T3 và Q4 mà còn mở rộng ra cho cả phần tử đa giác. Hiệu quả của phương pháp này được khảo sát bởi một số bài toán chuẩn như: vòng quay của đĩa Zalesak, dòng chảy của xoáy đơn đảo ngược theo thời gian. Phần còn lại của bài viết này được tổ chức như sau: phần tiếp theo trình bày lý thuyết của phương pháp LS thông thường. Phần 3 tập trung vào việc xây dựng phương pháp phần tử hữu hạn đa giác (Poly-FEM) cho các bài toán LS. Các kết quả số được trình bày và thảo luận trong phần 4. Sau cùng, một số kết luận sẽ được rút ra trong phần 5. 2. Phương pháp level set Trong LSM [4], ta giả sử một miền tính toán 𝒟𝒟 ⊂ ℝ2 chứa một miền con 𝛺𝛺 sao cho 𝛺𝛺 ⊂ 𝒟𝒟. Với 2.1. Quá trình tiến hóa của hàm level set một điểm bất kỳ 𝒙𝒙 ∈ 𝒟𝒟, ta xác định một hàm LS 𝜙𝜙 = 𝜙𝜙( 𝒙𝒙) để biểu diễn cho một bề mặt, như được mô 𝜙𝜙( 𝒙𝒙) < 0 ⇔ ∀𝒙𝒙 ∈ 𝛺𝛺\𝜕𝜕𝜕𝜕 tả như trong Hình 1 sao cho � 𝜙𝜙( 𝒙𝒙) = 0 ⇔ ∀𝒙𝒙 ∈ 𝜕𝜕𝜕𝜕 𝜙𝜙( 𝒙𝒙) > 0 ⇔ ∀𝒙𝒙 ∈ 𝒟𝒟\( 𝛺𝛺 ∪ 𝜕𝜕𝜕𝜕) (1) Hình 1. Hình minh họa hàm level trong thiết kế 2D: (a) bề mặt level set, và (b) các miền và giao diện Phần giao nhau giữa 𝛺𝛺 và 𝒟𝒟 được gọi là giao diện, 𝜕𝜕𝜕𝜕. Nó được biểu thị một cách ngầm định tính toán bởi giá trị bằng 0 của hàm LS trong Eq. (1), tức là {∀𝒙𝒙 ∈ 𝒟𝒟 ∶ 𝜙𝜙( 𝒙𝒙) = 0}. Sự phát triển của giao diện
403 Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set trường vận tốc 𝓿𝓿 = 𝓿𝓿( 𝒙𝒙, 𝑡𝑡) sao cho được xác định bằng cách giải một phương trình đạo hàm riêng LS Hamilton-Jacobi [1] dưới một ∀𝒙𝒙 ∈ 𝒟𝒟, + ( 𝓿𝓿 ∙ ∇) 𝜙𝜙 = 0 với 𝜙𝜙( 𝒙𝒙, 0) = 𝜙𝜙0 ( 𝒙𝒙) 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 trong đó 𝜙𝜙0 ( 𝒙𝒙) là một hàm LS khởi tạo, 𝑡𝑡 là biến thời gian. Việc giải Eq. (2) thường cho kết quả có sự (2) một phương pháp độ nhớt nhân tạo được thực hiện bằng cách thêm một điều kiện khuếch tán 𝜀𝜀∇2 𝜙𝜙 dao động của nghiệm (nghiệm yếu). Để giải quyết vấn đề này, Osher và Fedkim (2003) [9] đã đề xuất vào vế bên phải của Eq. (2), với 𝜀𝜀 là hệ số nhớt rất nhỏ. Do đó, phương trình LS (2) có thể được viết + ( 𝓿𝓿 ∙ ∇) 𝜙𝜙 − 𝜀𝜀∇2 𝜙𝜙 = 0 trong 𝒟𝒟 × 𝒯𝒯 lại dưới dạng một phương trình đạo hàm riêng đối lưu-khuếch tán LS như sau: 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜙𝜙( 𝒙𝒙, 0) = 𝜙𝜙0 ( 𝒙𝒙) trong 𝒟𝒟 × {0} (3) trong đó 𝒯𝒯 là khoảng thời gian cho sự phát triển của hàm LS. Để giải Eq. (3), phương pháp Crank- (4) �� = �� 𝜙𝜙 (𝑛𝑛+1) + �𝓿𝓿(𝑛𝑛+1) ∙ ∇𝜙𝜙 (𝑛𝑛+1) − 𝜀𝜀∇2 𝜙𝜙 (𝑛𝑛+1) � 𝜙𝜙 (𝑛𝑛) − �𝓿𝓿(𝑛𝑛) ∙ ∇𝜙𝜙 (𝑛𝑛) − 𝜀𝜀∇2 𝜙𝜙 (𝑛𝑛) � Nicolson [10] được áp dụng để rời rạc biến thời gian của Eq. (3): ∆𝑡𝑡 ∆𝑡𝑡 2 2 (5) 𝔏𝔏�𝜙𝜙(𝑛𝑛+1) � 𝔣𝔣� (𝑛𝑛) � trong đó ℚ = { 𝜙𝜙 ∶ 𝜙𝜙( 𝒙𝒙, 𝑡𝑡) ∈ ℍ1 ( 𝒟𝒟), 𝑡𝑡 ∈ 𝒯𝒯 } là không gian vô hạn chiều của nghiệm phụ thuộc thời gian; ∆𝑡𝑡 là bước thời gian; 𝜙𝜙 (𝑛𝑛) , 𝜙𝜙 (𝑛𝑛+1) , 𝓿𝓿(𝑛𝑛) , và 𝓿𝓿(𝑛𝑛+1) là các giá trị của hàm LS và trường vận tốc đối lưu tại thời điểm 𝑡𝑡 (𝑛𝑛) = �𝑡𝑡 (0) + 𝑛𝑛∆𝑡𝑡� và 𝑡𝑡 (𝑛𝑛+1) = �𝑡𝑡 (𝑛𝑛) + ∆𝑡𝑡�, tương ứng. 𝜙𝜙 𝑡𝑡 2.2. Quá trình tái khởi tạo hàm level set Trong suốt quá trình phát triển của giao diện, hàm LS không đảm bảo đặc tính SDF như hàm LS khởi tạo do các sai số tích lũy. Để giải quyết vấn đề này, một quá trình tái khởi tạo nên được thực hiện bằng cách giải một phương trình Eikonal [3, 4, 9]. Trong nghiên cứu này, phương pháp độ nhớt nhân + ( 𝓬𝓬 ∙ ∇) 𝜓𝜓 − 𝜀𝜀∇2 𝜙𝜙 = 𝓈𝓈( 𝜙𝜙) trong 𝒟𝒟 × 𝒯𝒯 tạo cũng được áp dụng cho quá trình tái khởi tạo. Do đó, ta có: 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜓𝜓( 𝒙𝒙, 0) = 𝜙𝜙( 𝒙𝒙) trong 𝒟𝒟 × {0} (6) 𝜏𝜏 trong đó 𝜏𝜏 là biến thời gian ảo, 𝒯𝒯 là khoảng thời gian ảo của quá trình tái khởi tạo, 𝜓𝜓 = 𝜓𝜓( 𝒙𝒙, 𝜏𝜏) là một (7) hàm LS đã hiệu chỉnh, 𝓬𝓬( 𝜙𝜙; 𝒙𝒙, 𝜏𝜏) ∶= 𝓈𝓈 ( 𝜙𝜙) 𝒏𝒏( 𝒙𝒙, 𝜏𝜏) là trường vận tốc đặc trưng với 𝒏𝒏( 𝒙𝒙, 𝑡𝑡) = là vectơ pháp tuyến đơn vị hướng ra, 𝓈𝓈( 𝜙𝜙) ∶= 2 2 là hàm dấu làm mịn, và 𝜚𝜚 = ∇𝜓𝜓(𝒙𝒙,𝜏𝜏) 𝜏𝜏 max(10−8 ,|∇𝜓𝜓(𝒙𝒙,𝜏𝜏)|) �𝜙𝜙 +𝜚𝜚 2ℎ là kích thước của vùng làm mịn (ℎ là kích thước phần tử). Eq. (6) là một phương trình phi tuyến vì 𝜙𝜙 trường vận tốc đặc trưng 𝓬𝓬( 𝜙𝜙; 𝒙𝒙, 𝜏𝜏) phụ thuộc vào thời gian ảo 𝜏𝜏. Để tuyến tính hóa điều kiện đối lưu 𝓬𝓬( 𝜙𝜙; 𝒙𝒙, 𝜏𝜏), phương pháp bước phân số Euler lùi [10] được áp dụng để tách Eq. (6) thành hai phương trình rời rạc theo thời gian như sau: + �𝓬𝓬(𝑘𝑘) ∙ ∇�𝜓𝜓 (∗) − 𝜀𝜀∇2 𝜓𝜓 (∗) = 0 trong 𝒟𝒟 × �𝜏𝜏 (𝑘𝑘) , 𝜏𝜏 (𝑘𝑘+1) � (∗) −𝜓𝜓(𝑘𝑘) ∆𝜏𝜏 (8) = 𝓈𝓈 ( 𝜙𝜙) trong 𝒟𝒟 × �𝜏𝜏 (𝑘𝑘) , 𝜏𝜏 (𝑘𝑘+1) � 𝜓𝜓 (𝑘𝑘+1) −𝜓𝜓(∗) ∆𝜏𝜏 trong đó ∆𝜏𝜏 là bước thời gian ảo tăng dần; 𝜓𝜓 (𝑘𝑘) và 𝜓𝜓 (𝑘𝑘+1) lần lượt là các giá trị LS tại thời điểm ảo (9) 𝜓𝜓 𝜏𝜏 (𝑘𝑘) và 𝜏𝜏 (𝑘𝑘+1) ; 𝜓𝜓 (∗) là một giá trị LS trung gian tại 𝜏𝜏 (∗) ∈ �𝜏𝜏 (𝑘𝑘) , 𝜏𝜏 (𝑘𝑘+1) �.
404 Nguyễn Trần Bá Đình, Nguyễn Hoàng Sơn và Phan Đức Huynh 3. Phương pháp phần tử hữu hạn đa giác cho các bài toán level set 3.1. Hàm dạng trên một phần tử đa giác bất kỳ Trong nghiên cứu này, nhóm tác giả chỉ tập trung vào việc sử dụng hệ tọa độ Wachspress [11, cấu trúc chia lưới. Trong khuôn khổ của Poly-FEM, ta giả định rằng 𝛺𝛺 ⊂ ℝ2 đã được phân vùng thành 12] với ánh xạ đẳng tham số trên một phần tử tham chiếu [13] để xây dựng các xấp xỉ phù hợp trên các 𝑛𝑛 𝑒𝑒 phần tử đa giác không chồng chéo 𝛺𝛺 𝑒𝑒 sao cho ≈ 𝛺𝛺ℎ ≔ ∑ 𝑒𝑒=1 𝛺𝛺 𝑒𝑒 . Với mỗi phần tử đa giác 𝛺𝛺 𝑒𝑒 bất 𝑛𝑛 𝑒𝑒 kỳ, ta ký hiệu 𝒙𝒙 𝑎𝑎 ( 𝑎𝑎 = 1, … , 𝑛𝑛 𝑒𝑒 ) là các đỉnh của nó và được sắp xếp ngược chiều kim đồng hồ, với 𝑛𝑛 𝑒𝑒 𝑛𝑛 𝑛𝑛 là số nút của phần tử 𝛺𝛺 𝑒𝑒 , 𝒆𝒆 𝑎𝑎 là cạnh nối giữa 𝒙𝒙 𝑎𝑎 với 𝒙𝒙 𝑎𝑎+1 , và 𝒏𝒏 𝑎𝑎 là vectơ pháp tuyến đơn vị hướng ra cạnh 𝒆𝒆 𝑎𝑎 . Với bất kỳ điểm 𝒗𝒗 ∈ 𝛺𝛺 𝑒𝑒 , ta gọi 𝑑𝑑 𝑎𝑎 ( 𝒙𝒙) là khoảng cách vuông góc của 𝒗𝒗 tới cạnh 𝒆𝒆 𝑎𝑎 như trong 𝑁𝑁 𝑎𝑎 ( 𝒙𝒙) = ∑ 𝑛𝑛 với 𝓌𝓌𝑎𝑎 ( 𝒙𝒙) = det( 𝓹𝓹 𝑎𝑎−1 , 𝓹𝓹 𝑎𝑎 ) Hình 2. Các hàm dạng Wachspress có thể được xác định bởi 𝓌𝓌 𝑎𝑎 (𝒙𝒙) 𝑏𝑏=1 𝓌𝓌 𝑏𝑏 (𝒙𝒙) (10) 𝒙𝒙 𝛺𝛺 ℯ 𝒏𝒏 trong đó 𝓹𝓹 𝑎𝑎 = 𝒏𝒏 𝑎𝑎 /𝑑𝑑 𝑎𝑎 ( 𝒙𝒙). Hàm dạng Wachspress thỏa mãn các thuộc tính của một hàm dạng cơ bản Hình 2. Hình minh họa tọa độ Wachspress được xác định bởi các khoảng cách vuông góc như Kronecker-delta, phân vùng thống nhất, không âm và độ chính xác tuyến tính [7]. Gradient của nó ∇𝑁𝑁 𝑎𝑎 ( 𝒙𝒙) = 𝑁𝑁 𝑎𝑎 �𝓻𝓻 𝑎𝑎 − ∑ 𝑏𝑏=1 𝑁𝑁𝑏𝑏 𝓻𝓻 𝑏𝑏 � với 𝓻𝓻 𝑎𝑎 = 𝓹𝓹 𝑎𝑎−1 + 𝓹𝓹 𝑎𝑎 được tính bởi công thức: 𝑛𝑛 𝑛𝑛 𝑒𝑒 (11) 3.2. Poly-FEM cho phương pháp level set Chúng ta giả sử 𝕎𝕎 ∈ ℍ1 ( 𝒟𝒟) là một không gian của hàm kiểm tra không phụ thuộc thời gian. Ta 3.2.1. Poly-FEM cho quá trình tiến hóa của hàm level set gọi ℜ( 𝑤𝑤) = 𝔏𝔏( 𝑤𝑤) − 𝔣𝔣� 𝜙𝜙 (𝑛𝑛) � là một hàm dư sao cho ℜ( 𝑤𝑤) ∶ 𝕎𝕎 → ℝ2 . Với một hàm kiểm tra tùy ý 𝑤𝑤 ∈ 𝕎𝕎, bằng cách áp dụng phương pháp bình phương nhỏ nhất tiêu chuẩn, chúng ta thấy rằng ∀𝑤𝑤 ∈ 𝕎𝕎, ∫𝒟𝒟 𝔏𝔏( 𝑤𝑤) 𝔏𝔏�𝜙𝜙 (𝑛𝑛+1) �d𝒟𝒟 = ∫𝒟𝒟 𝔏𝔏( 𝑤𝑤) 𝔣𝔣� 𝜙𝜙 (𝑛𝑛) �d𝒟𝒟 (12) trong đó 𝔏𝔏( 𝑤𝑤) = � � 𝑤𝑤 �( 𝓿𝓿 ∙ ∇) 𝑤𝑤 − ∇ ∙ ( 𝜀𝜀∇𝑤𝑤)� + ��. Sử dụng tích phân từng phần và ∆𝑡𝑡 2 Điều kiện Galerkin Điều kiện ổn định định lý phân kỳ, chúng ta có thể viết lại dạng yếu (12) thành dạng thu gọn như sau: với 𝜙𝜙 (𝑛𝑛) ∈ ℚ cho trước, tìm 𝜙𝜙 (𝑛𝑛+1) ∈ ℚ sao cho ∀𝑤𝑤 ∈ 𝕎𝕎, 𝑡𝑡 𝑡𝑡
405 Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set �� + ℬ�𝑤𝑤, 𝜙𝜙 � � = 𝒜𝒜�𝑤𝑤, 𝜙𝜙 (𝑛𝑛+1) � � � �� 𝒜𝒜�𝑤𝑤, 𝜙𝜙 (𝑛𝑛) � �� + ℬ�𝑤𝑤, 𝜙𝜙 (𝑛𝑛) � �� (𝑛𝑛+1) Điều kiện Galerkin Điều kiện ổn định Điều kiện Galerkin Điều kiện ổn định (13) 𝒜𝒜�𝑤𝑤, 𝜙𝜙 (𝑡𝑡) � = ∫𝒟𝒟 𝑤𝑤𝜙𝜙 (𝑡𝑡) d𝒟𝒟 + 𝛼𝛼 ∫𝒟𝒟 𝑤𝑤 ( 𝓿𝓿 ∙ ∇) 𝜙𝜙 (𝑡𝑡) d𝒟𝒟 + 𝛼𝛼 ∫𝒟𝒟 𝜀𝜀∇𝑤𝑤 ∙ ∇𝜙𝜙 (𝑡𝑡) d𝒟𝒟 ∆𝑡𝑡 ∆𝑡𝑡 với 2 2 ℬ�𝑤𝑤, 𝜙𝜙 (𝑡𝑡) � = ∫𝒟𝒟 𝜙𝜙 (𝑡𝑡) ( 𝓿𝓿 ∙ ∇) 𝑤𝑤d𝒟𝒟 + ∫𝒟𝒟 ∇𝜙𝜙 (𝑡𝑡) ∙ ∇𝑤𝑤d𝒟𝒟 + (14) ∆𝑡𝑡 𝜀𝜀∆𝑡𝑡 2 2 𝛼𝛼 � � ∫𝒟𝒟 ( 𝓿𝓿 ∙ ∇) 𝑤𝑤 ( 𝓿𝓿 ∙ ∇) 𝜙𝜙 (𝑡𝑡) d𝒟𝒟 ∆𝑡𝑡 2 2 1 if 𝑡𝑡 = 𝑛𝑛 + 1 (15) trong đó 𝛼𝛼 = � −1 if 𝑡𝑡 = 𝑛𝑛 . Bằng cách sử dụng khuôn khổ của Poly-FEM, hệ thống đại số tạo ra từ sự rời rạc phần tử hữu hạn cho dạng yếu trong Eq. (13) có thể được viết theo công thức sau: ∆𝑡𝑡 2 �� 2 𝐂𝐂 +� 𝐃𝐃 𝐌𝐌 + �� + 𝐂𝐂 𝑇𝑇 + 𝐃𝐃 𝑇𝑇 + � � 𝐊𝐊 (𝑛𝑛+1) = ∆𝑡𝑡 ∆𝑡𝑡 ∆𝑡𝑡 ∆𝑡𝑡 � �2 �� 𝚽𝚽 2 2 2 Điều kiện Galerkin Điều kiện ổn định ∆𝑡𝑡 2 �� 2 𝐂𝐂 +� 𝐃𝐃 𝐌𝐌 − �� + 𝐂𝐂 𝑇𝑇 + 𝐃𝐃 𝑇𝑇 − � � 𝐊𝐊 (𝑛𝑛) ∆𝑡𝑡 ∆𝑡𝑡 ∆𝑡𝑡 ∆𝑡𝑡 � �2 �� 𝚽𝚽 2 2 2 Điều kiện Galerkin Điều kiện ổn định (16) trong đó 𝚽𝚽 = ⋁ 𝑒𝑒=1 𝝓𝝓 𝑒𝑒 (𝑛𝑛+1) 𝑛𝑛 𝑒𝑒 (𝑛𝑛+1) 𝑡𝑡 (𝑛𝑛+1) ; 𝚽𝚽 = 𝝓𝝓 là lắp ghép của vectơ giá trị nút đã cho của hàm LS tại thời điểm 𝑡𝑡 (𝑛𝑛) ; ⋁ 𝑒𝑒=1 𝑒𝑒 (𝑛𝑛) (𝑛𝑛) 𝑛𝑛 𝑒𝑒 là lắp ghép của vectơ giá trị nút chưa biết của hàm LS tại thời điểm 𝐌𝐌 = ⋁ 𝑒𝑒=1 𝑴𝑴 , 𝐂𝐂 = 𝐂𝐂( 𝓿𝓿) = ⋁ 𝑒𝑒=1 𝑪𝑪 𝑒𝑒 ( 𝓿𝓿), 𝐃𝐃 = 𝐃𝐃( 𝜀𝜀 ) = ⋁ 𝑒𝑒=1 𝑫𝑫 𝑒𝑒 ( 𝜀𝜀 ) và 𝐊𝐊 = 𝐊𝐊( 𝓿𝓿) = ⋁ 𝑒𝑒=1 𝑲𝑲 𝑒𝑒 ( 𝓿𝓿) là 𝑛𝑛 𝑒𝑒 𝑒𝑒 𝑛𝑛 𝑒𝑒 𝑛𝑛 𝑒𝑒 𝑛𝑛 𝑒𝑒 lắp ghép của các ma trận khối lượng, đối lưu, khuếch tán và đối lưu cộng thêm, tương ứng; trong đó, các tenxơ bậc hai cơ sở của chúng được định nghĩa bởi: 𝑴𝑴 𝑒𝑒 = ∫𝐾𝐾 𝑒𝑒 𝑵𝑵 𝑒𝑒 ⨂𝑵𝑵 𝑒𝑒 d𝛺𝛺 hoặc 𝑀𝑀 𝑎𝑎𝑎𝑎 = ∫𝐾𝐾 𝑒𝑒 𝑁𝑁 𝑎𝑎 𝑁𝑁𝑏𝑏 d𝛺𝛺 𝑒𝑒 𝑒𝑒 𝑒𝑒 𝑪𝑪 𝑒𝑒 ( 𝓿𝓿) = ∫𝐾𝐾 𝑒𝑒 ( 𝑵𝑵 𝑒𝑒 ⨂𝓿𝓿) ∙ ( 𝑵𝑵 𝑒𝑒 ⨂∇)d𝛺𝛺 hoặc 𝐶𝐶 𝑎𝑎𝑎𝑎 ( 𝓿𝓿) = ∫𝐾𝐾 𝑒𝑒 𝑁𝑁 𝑎𝑎 𝓋𝓋𝑖𝑖 d𝛺𝛺 (17) 𝑒𝑒 𝑒𝑒 𝜕𝜕𝑁𝑁 𝑏𝑏 𝑒𝑒 𝜕𝜕𝑥𝑥 𝑖𝑖 (18) 𝑫𝑫 𝑒𝑒 ( 𝜀𝜀 ) = 𝜀𝜀 ∫𝐾𝐾 𝑒𝑒 (∇⨂𝑵𝑵 𝑒𝑒 ) ∙ ( 𝑵𝑵 𝑒𝑒 ⨂∇)d𝛺𝛺 hoặc 𝐷𝐷 𝑎𝑎𝑎𝑎 ( 𝜀𝜀 ) = 𝜀𝜀 ∫𝐾𝐾 𝑒𝑒 d𝛺𝛺 𝑒𝑒 𝜕𝜕𝑁𝑁 𝑎𝑎 𝜕𝜕𝑁𝑁 𝑏𝑏 𝑒𝑒 𝑒𝑒 𝜕𝜕𝑥𝑥 𝑖𝑖 𝜕𝜕𝑥𝑥 𝑖𝑖 (19) 𝑲𝑲 𝑒𝑒 ( 𝓿𝓿) = ∫𝐾𝐾 𝑒𝑒(∇⨂𝑵𝑵 𝑒𝑒 ) ∙ ( 𝓿𝓿⨂𝓿𝓿) ∙ ( 𝑵𝑵 𝑒𝑒 ⨂∇)d𝛺𝛺 hoặc 𝐾𝐾 𝑎𝑎𝑎𝑎 ( 𝓿𝓿) = ∫𝐾𝐾 𝑒𝑒 𝓋𝓋𝑖𝑖 𝓋𝓋𝑗𝑗 d𝛺𝛺 𝑒𝑒 𝜕𝜕𝑁𝑁 𝑎𝑎 𝑒𝑒 𝜕𝜕𝑁𝑁 𝑏𝑏 𝑒𝑒 𝜕𝜕𝑥𝑥 𝑖𝑖 𝜕𝜕𝑥𝑥 𝑗𝑗 (20) 3.2.2. Poly-FEM cho quá trình tái khởi tạo hàm level set Tương tự như tiểu mục phía trên, ta cũng áp dụng phương pháp bình phương nhỏ nhất tiêu chuẩn cho sự rời rạc không gian phần tử hữu hạn để thu được nghiệm ổn định trong Eq. (8). Trong khi đó, phương pháp Galerkin tiêu chuẩn được áp dụng cho bài toán đơn giản trong Eq. (9). Do đó, hệ �� +� 𝐌𝐌 + ∆𝜏𝜏𝐂𝐂 �∆𝜏𝜏𝐃𝐃 � �� + �� ∆𝜏𝜏𝐂𝐂 𝑇𝑇 + ∆𝜏𝜏𝐃𝐃 𝑇𝑇 + ∆𝜏𝜏 2 𝐊𝐊� 𝚿𝚿 (∗) = thống đại số tạo ra từ sự rời rạc phần tử hữu hạn cho các Eqs. (8) và (9) có thể được tạo thành như sau: Điều kiện Galerkin Điều kiện ổn định � ⏟ 𝐌𝐌 + �� 𝑇𝑇 + ∆𝜏𝜏𝐃𝐃 𝑇𝑇 ∆𝜏𝜏𝐂𝐂 �� 𝚿𝚿 (𝑘𝑘) � � Điều kiện Galerkin Điều kiện ổn định (21) 𝐌𝐌𝚿𝚿 = 𝐌𝐌𝚿𝚿 + ∆𝜏𝜏𝐅𝐅 (𝑘𝑘+1) (∗) trong đó 𝚿𝚿 = ⋁ 𝑒𝑒=1 𝝍𝝍 𝑒𝑒 là lắp ghép của vectơ giá trị nút chưa biết của hàm LS tại 𝜏𝜏 (𝑘𝑘+1) ; (𝑘𝑘+1) (22) (𝑘𝑘+1) 𝑒𝑒𝑛𝑛 𝚿𝚿 = ⋁ 𝑒𝑒=1 𝝍𝝍 𝑒𝑒 là lắp ghép của vectơ giá trị nút đã cho của hàm LS tại 𝜏𝜏 (𝑘𝑘) ; 𝚿𝚿 = ⋁ 𝑒𝑒=1 𝝍𝝍 𝑒𝑒 (𝑘𝑘) 𝑛𝑛 𝑒𝑒 (𝑘𝑘) (∗) 𝑒𝑒 𝑛𝑛 (∗)
406 Nguyễn Trần Bá Đình, Nguyễn Hoàng Sơn và Phan Đức Huynh là lắp ghép của véc tơ giá trị nút đã cho của hàm LS tại một thời điểm trung gian 𝜏𝜏 (∗) ∈ �𝜏𝜏 (𝑘𝑘) , 𝜏𝜏 (𝑘𝑘+1) �; 𝐌𝐌, 𝐂𝐂 = 𝐂𝐂�𝓬𝓬(𝑘𝑘) �, 𝐃𝐃 = 𝐃𝐃( 𝜀𝜀 ), và 𝐊𝐊 = 𝐊𝐊�𝓬𝓬(𝑘𝑘) � là lắp ghép của các ma trận khối lượng, đối lưu, khuếch tán và đối lưu cộng thêm như được định nghĩa trong Eq. (16), và 𝐅𝐅 = 𝐅𝐅�𝓈𝓈( 𝜙𝜙)� = ⋁ 𝑒𝑒=1 𝑭𝑭 𝑒𝑒 �𝓈𝓈 ( 𝜙𝜙)� là 𝑛𝑛 𝑒𝑒 𝑭𝑭 𝑒𝑒 �𝓈𝓈( 𝜙𝜙)� = ∫𝐾𝐾 𝑒𝑒 𝑵𝑵 𝑒𝑒 𝓈𝓈( 𝜙𝜙)d𝛺𝛺 𝐹𝐹𝑎𝑎𝑒𝑒 �𝓈𝓈 ( 𝜙𝜙)� = ∫𝐾𝐾 𝑒𝑒 𝑁𝑁 𝑎𝑎 𝓈𝓈( 𝜙𝜙)d𝛺𝛺 lắp ghép của vectơ lực với vectơ cơ sở của nó được xác định bởi: hoặc 𝑒𝑒 (23) 4. Kết quả số Trong phần này, hai bài toán chuẩn được khảo sát là: vòng quay của đĩa Zalesak và dòng chảy các phần tử T3, Q4 và đa giác. Trong các ví dụ số này, hệ số nhớt 𝜀𝜀 được chọn là 10−3 cho cả quá của xoáy đơn đảo ngược theo thời gian. Các ví dụ số này sẽ được tính toán dựa trên lưới cấu trúc của trình tiến hóa và tái khởi tạo của hàm LS. Các hàm LS thường được tái khởi tạo sau mỗi năm lần lặp lại của quá trình phát triển giao diện; trong khi đó, năm bước thời gian ảo sẽ được áp dụng cho quá trình tái khởi tạo của hàm LS. 4.1. Vòng quay của đĩa Zalesak một đĩa có rãnh có tâm đặt tại điểm có tọa độ (50, 75) nằm trong một miền (0, 100) × (0, 100) với Bài toán về vòng quay của đĩa Zalesak được xem xét như là ví dụ số đầu tiên. Trong ví dụ này, đồng hồ với tâm quay tại điểm có tọa độ (50, 50) dưới tác dụng của một trường vận tốc đối lưu 𝓿𝓿 = bán kính là 15, chiều dài là 25 và chiều rộng là 5. Đĩa Zalesak quay theo quỹ đạo tròn ngược chiều kim ( 𝑢𝑢, 𝑣𝑣 ), trong đó 𝑢𝑢 = ( 𝜋𝜋/314)(50 − 𝑦𝑦) � 𝑣𝑣 = ( 𝜋𝜋/314)( 𝑥𝑥 − 50) (24) Đĩa hoàn thành một vòng quay và trở lại vị trí ban đầu sau mỗi 628s. Chúng ta có thể xác định Error = � các sai số phân tán của giao diện bằng công thức sau: ∑|𝜙𝜙|
407 Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set -0.6 -0.7 Đa giác -0.8 -0.9 -1 (Error) -1.1 1 0 log -1.2 -1.3 -1.4 -0.4 -0.3 -0.2 -0.1 0 0.1 log (h ) 1 0 Hình 3. Đồ thị biểu diễn sai số phân tán theo kích thước phần tử (ℎ) cho bài toán vòng quay của đĩa Trong ví dụ này, ta cũng khảo sát hình dạng cũa đĩa Zalesak tại 𝑡𝑡 = 0, 157, 314, 471, và 628s Zalesak bằng cách sử dụng LS-PFEM với kích thước phần tử là 0.378 như trong Hình 4. Vì sai số luôn được và hình dạng ban đầu của nó cũng sẽ tăng lên. Tuy nhiên, hình dạng của đĩa Zalesak tại 𝑡𝑡 = 628s vẫn tích lũy trong quá trình tính toán, khi thời gian tăng lên thì sự khác biệt về hình dạng giữa đĩa Zalesak gần như tương tự với hình dạng ban đầu. Điều này đã phần nào cho thầy được tính hiệu quả của phương pháp được đề xuất trong nghiên cứu hiện tại. 100 t = 0s t = 628s 80 60 t = 157s t = 471s 40 y 20 t = 314s 0 0 20 40 60 80 100 x Hình 4. Hình dạng của đĩa Zalesak tại 𝑡𝑡 = 0, 157, 314, 471, và 628s bằng cách sử dụng LS-PFEM với kích thước phần tử là 0.378
408 Nguyễn Trần Bá Đình, Nguyễn Hoàng Sơn và Phan Đức Huynh 4.2. Mô phỏng dòng chảy của xoáy đơn đảo ngược theo thời gian này, một phần tử chất lỏng hình tròn đặt tại điểm có tọa độ (50, 75) nằm trong một miền (0, 100) × Ví dụ số thứ hai là bài toán về dòng chảy của xoáy đơn đảo ngược theo thời gian. Trong ví dụ (0, 100) với bán kính bằng 15 và được tác động bởi trường vận tốc đối lưu 𝓿𝓿 = ( 𝑢𝑢, 𝑣𝑣), trong đó 𝑢𝑢 = −sin2 ( 𝜋𝜋𝜋𝜋/100) ∙ sin( 𝜋𝜋𝜋𝜋/50) ∙ cos( 𝜋𝜋𝜋𝜋/𝑇𝑇) � 𝑣𝑣 = sin( 𝜋𝜋𝜋𝜋/50) ∙ sin2 ( 𝜋𝜋𝜋𝜋/100) ∙ cos( 𝜋𝜋𝜋𝜋/𝑇𝑇) (26) trong đó 𝑇𝑇 = 800 là thời gian mà xoáy đơn quay trở lại hình dạng tròn ban đầu của nó. Chúng ta có thể xác định sai số phân tán của giao diện bằng công thức sau Error = � ∑|𝜙𝜙|
409 Phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất mới cho bài toán level set Hình 6 cho thấy giao diện của xoáy đơn ở 𝑡𝑡 = 200, 400, 600, và 800s bằng cách sử dụng LS- PFEM với kích thước phần tử là 0.283. Chúng ta thấy rằng ở thời điểm 𝑡𝑡 là 200s và 600s thì hình dạng của xoáy đơn gần giống nhau, và khi 𝑡𝑡 là 800s, hình dạng của xoáy đơn tương tự như hình tròn ban đầu. 100 100 80 80 60 60 y y 40 40 20 20 0 0 0 20 40 60 80 100 0 20 40 60 80 100 x x (a) (b) 100 100 80 80 60 60 y y 40 40 20 20 0 0 0 20 40 60 80 100 0 20 40 60 80 100 x x (c) (d) PFEM với kích thước phần tử là 0.283: (a) 𝑡𝑡 = 200𝑠𝑠, (b) 𝑡𝑡 = 400𝑠𝑠, (c) 𝑡𝑡 = 600𝑠𝑠, và (d) 𝑡𝑡 = 800𝑠𝑠. Hình 6. Sự phát triển về hình dạng của phần tử chất lỏng tròn theo thời gian bằng cách sử dụng LS- 5. Kết luận Trong nghiên cứu này, một phương pháp phần tử hữu hạn đa giác bình phương nhỏ nhất (LS- PFEM) mới được sử dụng để giải quyết bài toán level set (LS) đối lưu-khuếch tán không ổn định. Để ổn định các nghiệm số mà không có các dao động phi vật lý (gradient dốc), cả quá trình tiến hóa và tái khởi tạo của hàm LS đều được giải bằng cách tối thiểu hóa điều kiện bình phương nhỏ nhất ổn định. Bên cạnh đó, hệ tọa độ Wachspress được sử dụng với ánh xạ đẳng tham số trên một phần tử tham chiếu để xây dựng các xấp xỉ phù hợp trên các cấu trúc chia lưới nhằm tính toán cho các tích phân số của các công thức dạng yếu. Phương pháp đề xuất được áp dụng để khảo sát một vài bài toán chuẩn chẳng hạn như: vòng quay của đĩa Zalesak, dòng xoáy đơn đảo ngược theo thời gian. Kết quả cho thấy rằng: so với phần tử T3 và Q4 thông thường, phần tử đa giác có khả năng cung cấp tính linh hoạt cao hơn trong việc tạo lưới cho các mô hình bài toán phức tạp cũng như mang lại kết quả tính toán chính xác hơn.
410 Nguyễn Trần Bá Đình, Nguyễn Hoàng Sơn và Phan Đức Huynh Trong tương lai, phương pháp được đề xuất sẽ phát triển thêm bằng cách sử dụng lưới thích ứng xung quanh giao diện để cải thiện độ chính xác cũng như mở rộng sang các vấn đề 3D và các vấn đề khác nhau của lĩnh vực kỹ thuật. Tài liệu tham khảo [1] S.H. Nguyen, H.-G. Kim. Stress-constrained shape and topology optimization with the level set method using trimmed hexahedral meshes. Methods Appl. Mech. Engrg., 366, (2020), pp. 61–113. [2] C. Basting, D. Kuzmin. A minimization-based finite element formulation for interface-preserving level set reinitialization. Computing, 95, (1), (2012), pp. 13–25. [3] M. Sussman, P. Smereka, S. Osher. A level set approach for computing solutions to incompressible two- phase flow. J. Comput. Phys., 144, (1), (1994), pp. 146–159. [4] H.G. Choi. A least-square weighted residual method for level set simulation. Int. J. Numer. Meth. Fluids, 68, (2012), pp. 887–904. [5] J.A. Sethian. A fast marching level set method for monotonically advancing fronts. Proc. Natl. Acad. Sci. USA, 93, (4), (1996), pp. 1591–1595. [6] R.N. Elias, M.A.D. Martins, A.L.G.A. Coutinho. Simple finite element-based computation of distance functions in unstructured grids. Int. J. Numer. Methods Engrg, 72, (2007), pp. 1095–1110. [7] H. Nguyen-Xuan, S. Nguyen-Hoang, T. Rabczuk, K. Hackl. A polytree-based adaptive approach to limit analysis of cracked structures. Computer Methods in Applied Mechanics and Engineering, 313, (2016), pp. 1006–1039. [8] H. Nguyen-Xuan. A polygonal finite element method for plate analysis. Computers & Structures, 188, (2017), pp. 45–62. [9] S. Osher, R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces, Springer, New York, (2003). [10] J. Donea, A. Huerta. Finite element methods for flow problems, Wiley, Chichester, (2003). [11] E.L. Wachspress. A rational finite element basis, Academic Press, New York, (1975). [12] M. S. Floater, K. Hormann, G. Kós. A general construction of barycentric coordinates over convex polygons. Advances in Computational Mathematics, 24, (1-4), (2006), pp. 311–331. [13] N. Sukumar, A. Tabarraei. Conforming polygonal finite elements. Int. J. Numer. Meth. Engng, 61, (2004), pp. 2045-2066.