Mô phỏng dòng chảy lưu chất nhớt không nén đƣợc qua trụ tròn dao động bằng phương pháp biên nhúng

Nội dung text: Mô phỏng dòng chảy lưu chất nhớt không nén đƣợc qua trụ tròn dao động bằng phương pháp biên nhúng

1. MÔ PHỎNG DÒNG CHẢY LƢU CHẤT NHỚT KHÔNG NÉN ĐƢỢC QUA TRỤ TRÒN DAO ĐỘNG BẰNG PHƢƠNG PHÁP BIÊN NHÚNG Võ Trung Nam1, Phan Đức Huynh2 1 Học viên cao học 2 Trường Đại học Sư Phạm Kỹ Thuật Thành Phố Hồ Chí Minh Tóm tắt Trong bài báo này, tác giả sử dụng phương pháp biên nhúng (IBM) để mô phỏng dòng chảy lưu chất nhớt không nén được qua trụ tròn dao động. Trong phương pháp này thì phương trình điều khiển Navier – Stokes cho dòng không nén được giải bằng phương pháp sai phân hữu hạn trên hệ lưới xen kẽ Đềcác và lưới Euler. Trụ tròn được rời rạc trên lưới Lagrange hình thành các điểm điều khiển. Thành phần lực khối sẽ sẽ được tính toán tại các điểm điều khiển và sau đó sẽ tính toán lan rộng đến các điểm lưới Đềcác xung quanh điểm điều khiển 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. Kết quả mô phỏng số cho hệ số Reynold tại Re=200 sẽ được so sánh với các phương pháp số khác. Từ khóa: Phương pháp biên nhúng; Phương pháp sai phân hữu hạn, Trụ tròn dao động; Phương trình Navier – Stokes; Phương pháp lưới Đềcác. I. Giới thiệu Phương pháp biên nhúng được đề xuất bởi C.S. Peskin [1,2] gần đây đã được đưa ra để giải quyết đối với các bài toán có những dạng hình học phức tạp trong khi yêu cầu tính toán ít hơn, tạo lưới dễ hơn so với các phương pháp khác mà vẫn đảm bảo được sự chính xác. Phương pháp này quy định một lực khối để thay thế sự hiện diện của một bề mặt mà không làm thay đổi lưới tính toán.Đặc biệt đối với trường hợp biên di chuyển, cụ thể hơn là trong sự tương tác giữa lưu chất và kết cấu làm biên dịch chuyển. Lưới bám theo vật thể phải được chia lại ở mỗi bước thời gian khi tính toán, việc chia lưới ảnh hưởng trực tiếp đến chi phí tính toán, độ chính xác và sự ổn định của lời giải. Mục đích của bài báo này là nhằm khảo sát những ưu điểm của phương pháp cũng như là đánh giá sự hợp lý của phương pháp thông qua bài toán mô phỏng dòng chảy lưu chất nhớt không nén được qua trụ tròn dao động.
2. Một hệ lưu chất – kết cấu đơn giản được trình bày trong hình 1.1 a) b) Hình 1.1 a) Hệ lưu chất-kết cấu đơn giản. b) Lưới rời rạc Euler (đánh dấu sáng) và lưới Lagrange (đánh dấu tối) II. Mô hình bài toán Ta xét mô hình bài toán cho dòng chảy lưu chất nhớt không nén được trong miền không gian hai chiều  f chứa biên nhúng dưới dạng là một đường cong khép kín đơn giản b (Hình 2.1a), hình dạng đường biên nhúng sẽ được cho dưới dạng tham số như sau: X s,t , 0 s Lb , X 0,t X Lb ,t , trong đó Lb là chiều dài của đường cong kín b , và X st, là hàm vector cho biết vị trí của các điểm trên biên nhúng tại thời gian t và chiều dài cung là s. Biên được mô hình bởi một lực đơn mà nó thì được kết hợp lại thành lực khối, f, trong phương trình Navier-Stokes. Phương trình Navier-Stokes sau đó được giải để xác định vận tốc của dòng lưu chất trên toàn bộ miền lưu chất  f . Kể từ khi biên nhúng tiếp xúc với dòng chất lưu chất xung quanh, vận tốc của nó phải phù hợp với điều kiện biên không trượt. Phương trình chuyển động của hệ như sau [1]: u u u p  u f (2.1) t u 0 (2.2) Với x x, y là điểm lưới 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 mà dạng công thức toán học là
3. f x,t F s,t  x X s,t ds (2.3)  Trong đó F s,t Fx s,t ,Fy s,t là lực khối tại điểm biên và x  x  y là hàm Dirac delta. Chuyển động của biên nhúng là X s,t U s,t u X s,t ,t u x,t  x X s,t dx (2.4) t  Phương trình (2.3) và (2.4) thể hiện sự tương tác giữa biên nhúng và lưu chất. Trong phương trình (2.3) thành phần lực khối tác dụng đến lưu chất bởi biên nhúng, trong khi đó phương trình (2.4) biên nhúng được di chuyển cùng với lưu chất. III. Phƣơng pháp số 1. Sự rời rạc hóa theo không gian và thời gian Phương pháp biên nhúng là phương pháp sai phân hữu hạn hỗn hợp Euler-Lagrange cho việc tính toán dòng lưu chất tương tác với một biên nhúng. Một cặp lưới tính toán bao gồm: một ô lưới Đềcác trung tâm cho biến Euler và một tập các điểm rời rạc cho biến Lagrange. Cho miền lưu chất  f 0,lx  0,ly  và có N x N y ô lưới Euler, với h hx hy lx / N x l y / N y là kích thước lưới Euler. Một cặp chỉ số trên biến sẽ chỉ rõ vị trí mà tại đó biến Euler đang được xem xét, đánh giá, chẳng hạn như uij chỉ rõ giá trị của biến u tại điểm lưới thứ ij. Chúng ta sẽ sử dụng một tập các điểm lưới Lagrange Nb (biên lưới được chia với s Lb / Nb ) và các điểm lưới Lagrange này được định nghĩa bởi một chỉ số duy nhất, với các biến tại nhiều điểm lưới thì được định nghĩa bằng các chỉ số tương ứng, do đó Fk biểu thị cho giá trị của biến F tại điểm lưới thứ k. Vị trí của điểm lưới Lagrange thứ k thì được theo dõi một cách rõ ràng trong Xk. Chúng ta sử dụng các chỉ số để biểu thị giá trị của biến tại một bước thời gian nhất định, do đó un x u x,n t và Xn s X s,n t . 2. Phƣơng pháp giải (giải vật thể) Thành phần lực khối sẽ sẽ được tính toán tại các điểm điều khiển và sau đó sẽ tính toán lan rộng đến các điểm lưới Đềcác xung quanh điểm điều khiển 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.
4. Nb n n n n fi, j Fk t  h xi, j Xi, j sk (2.5) k 1 Trong đó h x là một hàm delta Dirac hai chiều, 1 x x  h x   (2.6) h2 h h Ở đây là một hàm liên tục 1 3 2r 1 4 r 4 r2 , 0 r 1 8 1 2  r 5 2r 7 12 r 4 r , 1 r 2 8 0 2 r (2.7) Với x/h và y/h ký hiệu là r. Hàm liên tục 휙 sẽ thu được từ [1]. Phương trình Navier-Stokes với các thành phần cưỡng bức sau đó được giải cho trường n 1 n 1 áp suất pi, j và trường vận tốc ui, j tại các điểm lưới Đềcác bằng cách sử dụng phương pháp sai phân hữu hạn trên lưới Đềcác xen kẽ [3]. Trường vận tốc sau đó được nội suy để tìm ra vận tốc tại các điểm điều khiển, n 1 n 1 n 1 n 1 2 dXk / dt Uk ui, j  h xi, j Xk h (2.8) i, j 3. Giải phƣơng trình Navier - Stokes Xét phương trình Navier-Stokes cho dòng chảy lưu chất nhớt không nén được trong không gian hai chiều, phương trình (2.1)-(2.2): n 1 n n n 1 n 1 n [ ui, j ui, j / t  u ui, j ] pi, j  u i, j fi, j (2.9) n 1 ui, j 0 (2.10) Ta sẽ tìm lời giải tại bước thời gian thứ n 1 st bằng cách tiếp cận ba bước như sau:
5. Xử lý các thành phần phi tuyến, độ nhớt và thành phần lực khối u un / t un  un  / un 1/ f n (2.11) un 1 u / t pn 1 / (2.12) Sự hiệu chỉnh áp suất : Chúng ta hiệu chỉnh trường vận tốc trung bình u bằng cách lấy gradient của áp suất p n 1 , tiến hành lấy divergence ( ) cho cả hai vế của phương trình (2.12) ta có hệ sau pn 1 / u / t (2.13) Phương trình này gọi là phương trìnhPoissoncho áp suất p n 1 tại thời điểm thứ n 1 st . Cập nhật trường vận tốc mới : Tính toán trường vận tốc mới un 1 sử dụng mối liên hệ un 1 u t pn 1 / (2.14) Với giá trị áp suất p n 1 được tính ở bước trên (Giải phương trìnhPoisson). 4. Phƣơng pháp biên nhúng cho biên cứng Mô phỏng dòng chảy quanh một biên cứng, ta nên cho biên di chuyển một chút hơn là cố định chúng. Miễn sao biên nhúng nằm gần bề mặt vật thể, ta cần phải giải phương trình Navier- Stokes hai chiều với một số ràng buộc vận tốc biên như sau, u u u p  u F s,t  x X s,t ds (2.15) t  u 0 (2.16) X s,t U s,t u X s,t ,t u x,t  x X s,t dx (2.17) t  Khi biên nhúng là biên cứng, vấn đề chính là định luật liên tục ( ví dụ như định luật Húc cho lò xo) nhìn chung thì không được phù hợp trong giới hạn cứng, nghĩa là biến dạng nhỏ của một biên rất cứng. Giải pháp cho vấn đề này là ta giả sử vật thể đàn hồi nhưng rất cứng. Thành phần
6. cưỡng bức trong phương trình (2.15) đã nói lên rằng các điểm biên sẽ ở sát bề mặt vật thể, theo định luật Húc cho lò xo ta có: F s,t  X s,t Xe s or F s,t  Xe s X s,t (2.18) Với  là hệ số đàn hồi, Xe s là vị trí cân bằng, và X s,t là các điểm biên nhúng mà chúng sẽ tương tác với lưu chất. Áp đặt điều kiện biên chính xác trên lưới biên nhúng đòi hỏi κ lớn và sẽ chọn nó đủ lớn để sự chuyển động là không đáng kể.Vì thế, nếu các điểm trên biên di chuyển ra khỏi vị trí mong muốn thì lực đàn hồi sẽ kéo các điểm biên này quay trở lại. Do đó, trong thời gian tới ta có thể mong rằng các điểm biên sẽ luôn luôn nằm sát với đường biên mong muốn của chúng. IV. Kết quả mô phỏng số Để xác minh tính phù hợp của phương pháp biên nhúng cho bài toán có biên di chuyển, ta sẽ tiến hành mô phỏng dòng chảy lưu chất nhớt không nén được qua trụ tròn dao động theo phương ngang. Với kết quả mô phỏng số đạt được sẽ so sánh với kết quả từ các phương pháp số khác. Dao động của trụ tròn được cho dưới dạng là một dao động điều hòa bởi phương trình : X(t)= A (sin2πft) (4.1) Trong đó : A- biên độ dao động ; f - tần số dao động Với miền lưu chất Ω = [0,3] x [0,1.6], lưới tính toán là 512x216, đường kính trụ tròn D=0.1, khối lượng riêng của lưu chất 1.0, vận tốc lưu chất trên toàn miền u 1.0 , hệ số Re=200, biên độ dao động A=0.15, tần số dao động f=0.2 [4]. Các đại lượng này đều là những đại lượng không thứ nguyên. Trong trường hợp này phương trình (4.1) sẽ được áp đặt trên toàn bộ các điểm điều khiển mà cụ thể hơn là các điểm trên biên nhúng được chia trên lưới Lagrangian, và sự chuyển động của bề mặt trụ tròn theo phương x đều không tồn tại. a. Z.C.Zheng & N.Zhang,2007 b. Hiện tại
7. Hình 4.1 Đường bao xoáy của trụ tròn dao dộng sau 40 chu kỳ dao động tại Re=200, (a: Zhang & Zheng; b: Hiện tại) Hình 4.1 là đường bao xoáy sau 40 chu kỳ dao động của trụ tròn tại Re=200 và được so sánh với kết quả mô phỏng ở [4]. Trong hình này ta thấy dòng xoáy Von Karman đã phát triển gần như là đầy đủ và có thể dễ dàng quan sát được. Với kết quả này cho ta thấy phương pháp biên nhúng cho ra kết quả mô phỏng khá phù hợp với kết quả của Z.C.Zheng & N.Zhang, 2007 [4]. Re=200: HỆ SỐ CẢN TRƯỜNG Re=200: HỆ SỐ CẢN TRƯỜNG HỢP TRỤ TRÒN CỐ ĐỊNH HỢP TRỤ TRÒN DAO ĐỘNG Hình 4.2 Hệ số cản CD tại Re=200 đối với trụ tròn cố định và trụ tròn dao động Giá trị hệ số cản CDđạt được CD=1.29 ±0.002 nhỏ hơn so với trường hợp trụ tròn cố định CD=1.32 ±0.060 khoảng 2,27% và thời gian đạt trạng thái dao động ổn định ở giây 70 nhanh hơn so với trụ tròn cố định đạt tại giây 100. Re=200: HỆ SỐ NÂNG TRƯỜNG Re=200: HỆ SỐ NÂNG TRƯỜNG HỢP TRỤ TRÒN CỐ ĐỊNH HỢP TRỤ TRÒN DAO ĐỘNG Hình 4.3 Hệ số nâng CL tại Re=200 đối với trụ tròn cố định và trụ tròn dao động Giá trị hệ số nâng CLđạt được CL= ±0.543 lớn hơn so với trường hợp trụ tròn cố định CL=±0.495 khoảng 8,8% và thời gian đạt trạng thái dao động ổn định ở giây 65 nhanh hơn so với trụ tròn cố định đạt tại giây 115.
8. V. Kết luận Bài toán mô phỏng dòng chảy lưu chất nhớt không nén được qua trụ tròn dao động tại hệ số Re=200 đã được tính toán và mô phỏng bằng phương pháp biên nhúng. Ở bài toán này vấn đề chính là biên di chuyển mà không cần tái tạo lại lưới tính toán đã cho ra kết quả mô phỏng phù hợp với các kết quả mô phỏng từ các phương pháp số khác từ các tài liệu có sẵn. Với các kết quả từ mô phỏng số được so sánh, đối chứng với các kết quả từ các phương pháp số khác và kết quả từ thực nghiệm đã chứng tỏ tính chính xác, tính hợp lý của phương pháp biên nhúng. TÀI LIỆU THAM KHẢO [1] C. S. Peskin, The immersed boundary method, Acta Numer.11 (2) 479–517 (2002). [2] C. S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 220– 252 (1977). [3] M. Griebel, T. Dornseifer, T. Neunhoeffer, Numerical simulation in fluid dynamics: A practical introduction. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, (1998). [4] Zhang, N., Zheng, Z.C.,. An improved direct-forcing immersed-boundary method for finite difference applications. Journal of Computational Physics 221, (2007) 250–268.
9. BÀI BÁO KHOA HỌC THỰC HIỆN CÔNG BỐ THEO QUY CHẾ ĐÀO TẠO THẠC SỸ Bài báo khoa học của học viên có xác nhận và đề xuất cho đăng của Giảng viên hướng dẫn Bản tiếng Việt ©, TRƯỜNG ĐẠI HỌC SƯ PHẠM KỸ THUẬT TP. HỒ CHÍ MINH và TÁC GIẢ Bản quyền tác phẩm đã được bảo hộ bởi Luật xuất bản và Luật Sở hữu trí tuệ Việt Nam. Nghiêm cấm mọi hình thức xuất bản, sao chụp, phát tán nội dung khi chưa có sự đồng ý của tác giả và Trường Đại học Sư phạm Kỹ thuật TP. Hồ Chí Minh. ĐỂ CÓ BÀI BÁO KHOA HỌC TỐT, CẦN CHUNG TAY BẢO VỆ TÁC QUYỀN! Thực hiện theo MTCL & KHTHMTCL Năm học 2016-2017 của Thư viện Trường Đại học Sư phạm Kỹ thuật Tp. Hồ Chí Minh.