Mô phỏng dòng chảy nhớt không nén được trong một miền vuông chứa vật cản trụ tròn ở tâm miền tính toán

Chia sẻ: _ _ | Ngày: | Loại File: PDF | Số trang:8

Thêm vào BST

Báo xấu

6
lượt xem 1
download

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

Bài viết này được tổ chức như sau, phần 2 trình bày hệ phương trình chuyển động của dòng chảy nhớt khôn nén qua một vật cản không đàn hồi. Trong phần 3, phương pháp số để giải quyết bài toán được dẫn ra. Kết quả mô phỏng được minh họa ở phần 4. Sau cùng, một số kết luận được trình bày ở phần 5.

Chủ đề:

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

Lưu

Nội dung Text: Mô phỏng dòng chảy nhớt không nén được trong một miền vuông chứa vật cản trụ tròn ở tâm miền tính toán

MÔ PHỎNG DÒNG CHẢY NHỚT KHÔNG NÉN ĐƯỢC TRONG MỘT MIỀN VUÔNG CHỨA VẬT CẢN TRỤ TRÒN Ở TÂM MIỀN TÍNH TOÁN Lê Quốc Cường1, Nguyễn Bá Duy2,* 1. Viện Kỹ thuật-Công nghệ, Trường Đại học Thủ Dầu Một 2. Khoa Kiến trúc, Trường Đại học Thủ Dầu Một * Email: duynb@tdmu.edu.vn TÓM TẮT Phương pháp biên nhúng được áp dụng để mô phỏng dòng chảy nhớt không nén qua một miền vuông có trụ tròn cố định ở tâm của miền tính toán. Với phương pháp biên nhúng, sự ảnh hưởng của vật cản lên dòng lưu chất được xử lý bằng cách đưa một thành phần lực cưỡng bức vào hệ phương trình Navier-Stokes của dòng chảy. Trong phương trình chuyển động của dòng chảy, sự kết hợp của vận tốc và áp xuất được xử lý bằng phương pháp chiếu. Kết quả mô phỏng số cho bài toán dòng chảy nhớt không nén trong một miền vuông chứa vật cản trụ tròn ở tâm miền tính toán được thực hiện ở hệ số Reynolds và được so sánh với các kết quả tính toán đã được công bố. Từ khóa: dòng chảy nhớt không nén qua vật cản cố định, động lực học chất lỏng tính toán, hệ phương trình Navier-Stokes, phương pháp biên nhúng, phương pháp sai phân hữu hạn. 1. ĐẶT VẤN ĐỀ Trong những năm gần đây, phương pháp biên nhúng ngày càng trở nên phổ biến trong việc mô phỏng các bài toán tương tác rắn – lỏng do tính đơn giản, hiệu quả và linh hoạt, cũng như khả năng xử lý dòng chảy phức tạp và kết cấu biến dạng lớn. Phương pháp biên nhúng (IBM-Immersed boundary method) ban đầu được phát triển bởi (C. S. Peskin, 1977) để nghiên cứu lưu lượng máu qua tim, và từ đó được nghiên cứu rộng rãi và áp dụng cho hàng loạt các vấn đề tương tương tác lưu chất – kết cấu (K. Goncharuk and nnk., 2023; W.-X. Huang and F.-B. Tian, 2019; W. Kim and H. Choi, 2019; F. Sotiropoulos and X. Yang, 2014; S. Tschisgale and J. Fröhlich, 2020; J. Yang, 2016). Phương pháp biên nhúng giải quyết các phương trình chất lỏng với một thành phần bổ sung đó là lực tương tác lưu chất – kết cấu, đại diện cho tác động của biên nhúng lên sự chuyển động của dòng lưu chất. Lực tương tác lưu chất – kết cấu được tính toán từ cấu trúc biên nhúng, sau đó được sử dụng để tính toán cho vận tốc và áp suất chất lỏng. Về cơ bản, các phương trình chất lỏng được giải quyết trong toàn bộ miền lưu chất với một lưới Euler cố định, biên kết cấu di chuyển được tính toán trên một lưới Lagrangian. Với phương pháp này, việc cập nhật lưới được loại bỏ hoàn toàn. Phân tích chi tiết và các ứng dụng khác nhau của phương pháp biên nhúng được trình bày bởi (R. Mittal and G. Iaccarino, 2005). Bài báo này được tổ chức như sau, phần 2 trình bày hệ phương trình chuyển động của dòng chảy nhớt khôn nén qua một vật cản không đàn hồi. Trong phần 3, phương pháp số để giải quyết bài toán được dẫn ra. Kết quả mô phỏng được minh họa ở phần 4. Sau cùng, một số kết luận được trình bày ở phần 5. 917
2. HỆ PHƯƠNG TRÌNH CHUYỂN ĐỘNG Xét bài toán dòng chảy nhớt không nén trong miền chữ nhật hai chiều  f =  0, lx   0, l y    chứa một biên nhúng không khối lượng ở dạng một đường cong khép kín  như trình bày ở hình 1. Hình 1. Hệ lưu chất-kết cấu đơn giản và lưới rời rạc Euler (đánh dấu sáng) và lưới Lagrange (đánh dấu tối) Cấu trúc của biên nhúng được cho ở dạng tham số: X(s,t), 0  s  Lb , X ( 0, t ) = X ( Lb , t ) , ở đây Lb là chiều dài biên  , s là chiều dài cung và t là thời gian. Ảnh hưởng của biên nhúng lên lưu chất được trình bày bởi thành phần lực cưỡng bức f tác dụng lên lưu chất. Vì vậy, sự chuyển động của dòng lưu chất được mô tả bằng hệ phương trình Navier-Sokes như sau u  +  ( u  ) u + p = u + f (1) t  u = 0 (2) Với x = ( x, y ) là tọa độ trên lưới Euler và X = ( X , Y ) là điểm biên trên lưới Lagrange, u ( x, t ) = ( u ( x, t ) , v ( x, t ) ) là vận tốc của lưu chất và p ( x, t ) là áp suất lưu chất. Các hệ số  và  lần lượt là khối lượng riêng và độ nhớt của lưu chất. Thành phần lực khối tác dụng lên lưu chất là f ( x, t ) = ( f x ( x, t ) , f y ( x, t ) ) có dạng công thức toán học là f ( x, t ) =  F ( s, t )  ( x − X ( s, t ) ) ds (3)  Trong đó:  ( x ) =  ( x )  ( y ) là hàm rời rạc Delta và F ( s, t ) = ( Fx ( s, t ) , Fy ( s, t ) ) là lực khối tại các điểm biên được xác định theo đề xuất của Lai & Peskin [9] như sau F ( s, t ) =  ( X e ( s ) − X ( s, t ) ) (4) 918
Trong đó:  là một hằng số dương và  1, X e ( s ) là vị trí cân bằng cố định, và X ( s, t ) là các điểm biên nhúng ở thời điểm t . Sự chuyển động của biên nhúng được tính theo công thức sau X ( s, t ) = U ( s, t ) = u ( X ( s, t ) , t ) =  u ( x, t )  ( x − X ( s, t ) ) dx (5) t  Phương trình (3) và (5) thể hiện sự tương tác giữa biên nhúng và lưu chất. Trong phương trình (3) là thành phần lực khối tác dụng đến lưu chất gây ra bởi biên nhúng, trong khi đó ở phương trình (5) biên nhúng được di chuyển cùng với lưu chất. 3. HỆ PHƯƠNG TRÌNH CHUYỂN ĐỘNG 3.1. Phương pháp chiếu Khó khăn chủ yếu trong việc giải hệ phương trình Navier-Stokes đó là sự kết hợp của vận tốc – áp suất, để giải quyết vấn đề này, một phương pháp chiều được đề xuất bởi (A. J. Chorin, 1968) đã được sử dụng. Xét phương trình Navier-Stokes trong không gian hai chiều ở bước thời gian thứ n + 1như sau un+1 − un  = −pn+1 −  ( un  ) un + un + f n (6) t  un+1 = 0 (7) với điều kiện biên u n+1  = ub +1 n (8) Hệ phương trình (6) – (8) được giải theo trình tự như sau: Bước 1: tính trực tiếp vận tốc trung gian u * từ phương trình động lượng (4.8) bỏ qua thành phần gradient áp suất u − u n  = − ( u n  ) u n + u n + f n 1 (9) t   Trong đó: u n là vận tốc ở bước thời gian thứ n . Ở bước thời gian tiếp theo, ta có u n +1 − u 1 = − p n +1 (10) t  Bước 2: Hiệu chỉnh áp suất  p n +1 =   u (11) t Đây là phương trình Poisson cho áp suất, giải phương trình này chúng ta sẽ tìm được áp suất ở bước thời gian n + 1. Bước 3: Cập nhật vận tốc ở bước thời gian n + 1 Với áp suất vừa tìm được từ phương trình (11), thay vào phương trình (10) ta có vận tốc ở bước thời gian kế tiếp n + 1 được tính như sau 919
t un+1 = u − p n+1 (12)  3.2. Xác định lực cưỡng bức f Để giải hệ phương trình Navier – Stokes, thành phần lực cững bức f n phải được xác định. Thành phần lực khối sẽ sẽ được tính toán tại các điểm trên biên nhúng và sau đó sẽ được phân bố đến các điểm lưới Đề các xung quanh trên toàn miền lưu chất bằng một biểu diễn rời rạc của hàm Dirac delta. Hình 2 minh họa quá trình nội suy vận tốc tại các điểm trên biên nhúng và quá trình phân bố lực khối trên biên nhúng đến các điểm lưới xung quanh trên toàn miền lưu chất. Hình 2. Phân bố lực cưỡng bức từ một điểm trên biên nhúng đến các điểm lưới lân cận và nội suy vận tốc ở một điểm khác trên biên nhúng Lực cưỡng bức tại các điểm trên biên nhúng được tính như sau Fkn =  ( Xe − X k ) k n (13) Trong đó:  là một hằng số dương và  e 1, Xk là vị trí cân bằng cố định của điểm biên nhúng thứ k , và Xn là các điểm biên nhúng thứ k ở thời điểm t = nt . k Sau khi xác định được lực cưỡng bức tại các điểm biên nhúng ở bước thời gian thứ n , lực này sẽ được phân bố đến các điểm lân cận theo công thức sau Nb fin j =  Fkn h ( xin, j − Xk ) s , n (14) k =1 Trong đó: x i , j và fi , j là tọa độ điểm lưới ( i, j ) và lực khối tương ứng tại điểm đó.  h ( x ) là một hàm rời rạc delta hai chiều được tính theo công thức sau 1  x  y h ( x) =      (15) h2  h   h  920
Trong đó: h là kích thước lưới, x và y là hai thành phần tọa độ Đề Các của x và  ( r ) là hàm rời rạc delta một chiều và  ( r ) được chọn như trong (C. S. Peskin, 2002) như sau 1 (  8 3 − 2 r + 1 + 4 r − 4r , 2 ) 0  r 1  1 ( )  ( r ) =  5 − 2 r − −7 + 12 r − 4r 2 , 1  r  2 8 (16) 0 , 2 r   Vì lực khối được tính ở các điểm trên biên nhúng và mở rộng đến toàn bộ các điểm lưới, hệ phương trình Navier – Stokes sẽ được giải để tìm vận tốc và áp suất ở các điểm lưới ở bước thời gian thứ n + 1. Trường vận tốc sau đó được nội suy để tìm vận tốc tại các điểm biên nhúng như sau Un+1 ( Xk ) =  u n+1 ( xi , j )  h ( xin,+j1 − Xn+1 ) h2 k k (17) i, j Từ vận tốc tại các trên biên nhúng, vị trí của các điểm biên được tính theo công thức sau Xn+1 = Xn + tUn+1 ( Xk ) k k k (18) 3.3. Giải thuật tổng quát Giả sử các biến của bài toán ở bước thời gian thứ n đã biết. Để giải quyết bài toán ở bước thời gian thứ n + 1, chúng ta tiến hành các bước sau: Bước 1: Tính lực cưỡng bức tại các điểm trên biên nhúng F n theo công thức (13). Bước 2: Phân bố lực cưỡng bức từ các điểm trên biên nhúng đến các điểm lưới xung quanh theo công thức (14). Bước 3: Cập nhật vận tốc trung gian u * từ phương trình (9). Bước 4: Giải phương trình Poisson áp suất (11) để tìm p n +1 . Bước 5: Sử dụng công thức (12) để tìm u n +1 . Bước 6: Nội suy vận tốc tại các điểm trên biên nhúng U k từ phương trình (17). Bước 7: Từ vận tốc vừa tìm được ở bước 6, cập nhật vị trí của biên nhúng Xn theo phương trình (18). 4. KẾT QUẢ MÔ PHỎNG SỐ Trong phần này, phương pháp IB sẽ được áp dụng để mô phỏng bài toán dòng chảy trong một miền vuông với vật cản là một trụ tròn cố định ở giữa miền tính toán. Miền tính toán và điều kiện biên của bài toán được cho như hình 4. Dòng chảy được dẫn bởi sự di chuyển của biên trên ở một hằng số vận tốc nhất định u0 . Hệ số Reynolds của bài toán được định nghĩa như sau  u0 D Re = (21)  Trong đó: D là đường kính trụ tròn, u0 là vận tốc dòng chảy theo phương ngang ở biên trên của miền tính toán,  và  lần lượt là khối lượng riêng và độ nhớt của lưu chất. 921
Hình 4. Miền tính toán và điều kiện biên của bài toán Chúng ta sẽ tiến hành mô phỏng bài toán với các thông số như sau: L = 1 , u0 = 1 , Re = 1000 và D = 0.4L . Lưới đều 301 301 điểm lưới được sử dụng để rời rạc miền lưu chất, bước thời gian t = 10−4 s được sử dụng để mô phỏng bài toán. Hình 5 trình bày đường dòng của bài toán ở hệ số Re = 1000 . Từ kết quả ở hình 5, có ba xoáy nước xuất hiện. Một xoáy ở vị trí phía trên bên phải của trụ tròn và hai xoáy ở phía dưới gần góc phải và góc trái. Có thể nhận thấy rằng xoáy phía trên được tạo ra do sự xuất hiện của trụ tròn cố định ở tâm của miền tính toán. Thành phần vận tốc theo phương ngang u ở vị trí x = 0.5 và thành phần vận tốc theo phương đứng v ở vị trí y = 0.5 được trình bày ở hình 6. Kết quả cho thấy sự đồng thuận khá tốt của phương pháp đề xuất với kết quả của (S.-G. Cai and nnk., 2017. Để đánh giá ảnh hưởng của bước lưới đến độ chính xác của lời giải, chúng ta sẽ tiến hành khảo sát bài toán với các bước lưới khác nhau h = (0.02500, 0.01250, 0.00625, 0.003125). Vì bài toán không có lời giải chính xác nên chúng ta sử dụng lời giải ở bước lưới h = 0.0015625 như lời giải tham khảo để tính toán sai số. Bài toán được khảo sát đến thời điểm t = 0.2s với bước thời gian Δt = 5x10-5s để đảm bảo sự ổn định của bài toán ở các bước lưới nhỏ. Hình 7 trình bày các giá trị sai số của thành phần vận tốc theo phương ngang eu 2 . Từ kết quả ở hình 7 cho thấy sai số có bậc hội tụ khoảng 1.32. Hình 5. Đường dòng của bài toán Lid-driven cavity 922
với trụ tròn ở tâm miền tính toán Hình 6. Thành phần vận tốc theo phương ngang u ở vị trí x =0.5 (trái) và vận tốc theo phương đứng ν ở vị trí y = 0.5 (phải) của bài toán Hình 7. Sai số của thành phần vận tốc theo phương ngang ở các bước lưới khác nhau. 5. KẾT LUẬN Bài báo đã trình bày kết quả mô phỏng dòng chảy nhớt không nén qua miền vuông có trụ tròn cố định ở tâm của miền tính toán. Các kết quả mô phỏng số về đường dòng cũng như các thành phần vận tốc của dòng chảy được thực hiện ở hệ số Reynolds Re = 1000 . Các kết quả mô phỏng cho thấy tính hiệu quả và độ chính xác của phương pháp đề xuất với sai số có bậc hội tụ khoảng 1.32. 923
TÀI LIỆU THAM KHẢO 1. C. S. Peskin, "Numerical analysis of blood flow in the heart," Journal of Computational Physics, vol. 25, no. 3, pp. 220-252, 1977/11/01/ 1977. 2. K. Goncharuk, O. Oshri, and Y. Feldman, "The immersed boundary method: A SIMPLE approach," Journal of Computational Physics, vol. 487, p. 112148, 2023/08/15/ 2023. 3. W.-X. Huang and F.-B. Tian, "Recent trends and progress in the immersed boundary method," Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 233, no. 23-24, pp. 7617-7636, 2019/12/01 2019. 4. W. Kim and H. Choi, "Immersed boundary methods for fluid-structure interaction: A review," International Journal of Heat and Fluid Flow, vol. 75, pp. 301-309, 2019/02/01/ 2019. 5. F. Sotiropoulos and X. Yang, "Immersed boundary methods for simulating fluid–structure interaction," Progress in Aerospace Sciences, vol. 65, pp. 1-21, 2014/02/01/ 2014. 6. S. Tschisgale and J. Fröhlich, "An immersed boundary method for the fluid-structure interaction of slender flexible structures in viscous fluid," Journal of Computational Physics, vol. 423, p. 109801, 2020/12/15/ 2020. 7. J. Yang, "Sharp interface direct forcing immersed boundary methods: A summary of some algorithms and applications," Journal of Hydrodynamics, Ser. B, vol. 28, no. 5, pp. 713-730, 2016/10/01/ 2016. 8. R. Mittal and G. Iaccarino, "Immersed Boundary Methods," Annual Review of Fluid Mechanics, vol. 37, no. 1, pp. 239-261, 2005. 9. M.-C. Lai and C. S. Peskin, "An Immersed Boundary Method with Formal Second-Order Accuracy and Reduced Numerical Viscosity," Journal of Computational Physics, vol. 160, no. 2, pp. 705- 719, 2000/05/20/ 2000. 10. A. J. Chorin, "Numerical Solution of the Navier-Stokes Equations," Mathematics of Computation, vol. 22, no. 104, pp. 745-762, 1968. 11. C. S. Peskin, "The immersed boundary method," in Acta Numerica 2002, vol. 11, A. Iserles, Ed. (Acta Numerica, Cambridge: Cambridge University Press, 2002, pp. 479-518. 12. S.-G. Cai, A. Ouahsine, J. Favier, and Y. Hoarau, "Moving immersed boundary method," International Journal for Numerical Methods in Fluids, vol. 85, no. 5, pp. 288-323, 2017/10/20 2017. 924