Báo cáo Khoa học Ứng dụng kỹ thuật tạo lưới trong bài toán mô phỏng dòng phun rối xoáy hai pha không đẳng nhiệt

Bỏo cỏo khoa học: Ứng dụng kỹ thuật tạo lưới trong bài toỏn mụ phỏng dũng phun rối xoỏy hai pha khụng đẳng nhiệt Tạp chí KHKT Nông nghiệp 2006 Tập IV, số 6:111-115 Đại học Nông nghiệp I ứng dụng kỹ thuật tạo l−ới trong bài toán mô phỏng dòng phun rối xoáy hai pha không đẳng nhiệt An application of grid generation technology to simulate Two-phase Non-isothermal swirling turbulent Flow Nguyễn Thanh Nam1, Nguyễn Thanh Hào2, Hoàng Đức Liên3 SUMMARY The numerical solution of partial differential equations requires some discretization of the field into a collection of points. The differential equations are approximated by a set of algebraic equations on this collection. This system of algebraic equations is then solved to produce a set of discrete values which approximate the solution of the partial differential system over the field. The discretization of the field requires some organization for the solution thereon to be efficient, it must be possible to readily identify the points neighboring the computation site. Furthermore, the discretization must conform to the boundaries of the region in such a way that boundary conditions can be accurately represented. The boundaries of the flame are not straight lines, how to determine grid points inside physical region which be used to solves equations of two- phase non-isothermal swirling turbulent flow in industrial combustion chamber is presented in this paper. Generalizing from the above consideration, the computational region may be treated as follows: ξ(x,y) and η(x,y) on the boundaries of the flame in manner: Set, η = constant, ξ = monotonically varying along the boundary of the physical region, and set, ξ = constant, η = monotonically varying along the boundary of the physical region. The grid points inside computational region will be deformed into rectangle to form the transformed region. After that, the inside points of physical region will be determined as follows: x(ξ,η) and y(ξ,η)... Here we consider the situations in which Cartesian coordinates are used both in the physical and computational regions. Key words: grid generation technology, two-phase, non-isothermal swirling turbulent flow 1. ĐặT VấN Đề Sử dụng ph−ơng pháp sai phân hữu hạn có thể dễ dàng giải các bài toán có miền khảo sát dạng hình chữ nhật. Tuy nhiên, khi giải bài toán mô phỏng dòng phun rối xoáy hai pha không đẳng nhiệt bằng ph−ơng pháp này lại gặp rất nhiều khó khăn vì biên dạng của dòng phun là đ−ờng cong đối xứng qua trục toạ độ (Hình 3.1a). Xét ph−ơng trình vi phân tổng quát của dòng phun rối xoáy hai pha không đẳng nhiệt có dạng: ( ) ( ) 0. . .. . ..... =+      ∂ ∂ ∂ ∂ −      ∂ ∂ ∂ ∂ −            ∂ ∂ ∂ ∂ −      ∂ ∂ ∂ ∂ ϕ ϕ ϕ ϕ ϕϕ ϕϕψϕψϕ dr r c rb rz c rb zzrrz a (1.1) Trong đó: aϕ, bϕ, cϕ, dϕ là các hệ số, ψ là hàm của đ−ờng dòng và ϕ là biến số. 1 Khoa Cơ khí, Đại học Bách khoa TP HCM 2 Khoa Cơ khí, Đại học Công nghiệp TP HCM 3 Khoa Cơ - Điện, Đại học Nông nghiệp I Theo Ph−ơng pháp sai phân hữu hạn, ph−ơng trình tổng quát (1.1) sẽ đ−ợc rời rạc hoá bằng cách thay thế các biểu thức vi phân bằng các tỷ sai phân t−ơng ứng, để chuyển ph−ơng trình vi phân (1.1) về dạng ph−ơng trình đại số. Sau đó sử dụng thuật toán nội suy kết hợp với ph−ơng pháp lặp giữa các điểm nằm trên biên với các điểm bên trong miền vật lý để giải hệ ph−ơng trình đại số. Nh−ng việc nội suy này không thực hiện đ−ợc hoặc phải chấp nhận sai số rất lớn, do biên dạng của ph−ơng trình tổng quát (1.1) là một đ−ờng cong, nên khi tiến hành chia l−ới miền khảo sát dạng l−ới hình chữ nhật sẽ có một số nút l−ới không nằm trên biên của miền vật lý mà chúng chỉ nằm gần biên hoặc rơi ra khỏi miền vật lý. Vấn đề này sẽ đ−ợc giải quyết khi ta chuyển hệ trục toạ độ từ miền vật lý có biên dạng là một đ−ơng cong về miền tính toán có biên dạng là đ−ờng thẳng (hình 3.1), trong đó khoảng cách giữa các nút l−ới theo ph−ơng x là đều nhau, còn khoảng cách giữa các nút l−ới theo ph−ơng y là không đều. 2. MÔ HìNH TíNH TOáN CHUYểN ĐổI LƯớI Trình tự chuyển đổi hệ trục toạ độ từ miền vật lý (x,y) sang miền tính toán (ξη) bao gồm các b−ớc cơ bản nh− sau (M.Necati Ozisik, 2000): - Xác định mối quan hệ giữa hệ trục toạ độ từ miền vật lý (x,y) và hệ trục toạ độ tính toán (ξη) bới các ph−ơng trình vi phân Laplas hoặc ph−ơng trình Poison của elliptic. - Chuyển đổi toạ độ từ miền vật lý (x,y) sang miền tính toán (ξ,η) trong hệ trục toạ độ Đề các. - Chuyển đổi các ph−ơng trình vi phân trong miền vật lý thành các ph−ơng trình vi phân trong miền tính toán. - Giải các ph−ơng trình trong miền tính toán, sau đó chuyển kết quả tìm đ−ợc trong miền tính toán thành miền vật lý ban đầu. Xét một ph−ơng trình vi phân riêng phần có các biến độc lập x,y trong miền vật lý. Phép biến đổi từ các biến x,y sang các biến ξ,η có thể đ−ợc biểu diễn nh− sau: ξ ≡ ξ(x,y); η ≡ η(x,y) (2.1a,b) Và phép biến đổi ng−ợc là x ≡ ξ,η; y ≡ ξ,η (2.2a,b) Phép biến đổi Jacobi J đ−ợc đ−a ra nh− sau (Courant, 1956): 0 , , ≠−= − − =      = ξηηξ ηη ξξ ηξ yxyxyx yxyxJJ (2.3) Trong đó: ηξ ηξ ∂ ∂ = ∂ ∂ = yyxx , (2.4a,b) Theo định luật Cramer, ta có: ηη ξξ xJyJ yx 1 , 1 −== (2.5a,b) ξξ ηη xJ y J yx 1 , 1 =−= (2.5c,d) Quan hệ giữa các nút l−ới trên biên của miền vật lý trong hệ toạ độ Đề các (x,y) và các nút l−ới trên biên của miền tính toán trong hệ toạ độ Đề các (ξ,η) (Joe F.Thompson, Z.U.A.Warsi, C.Wayne Mastin, 1985) là: η = hằng số, ξ = thay đổi tuyến tính dọc theo biên của miền vật lý. ξ = hằng số, η = thay đổi tuyến tính dọc theo biên của miền vật lý. B−ớc tiếp theo là chuyển đổi các toạ độ l−ới bên trong miền vật lý sang miền tính toán với điều kiện các đ−ờng toạ độ có xu h−ớng cách đều nhau ở trong miền và các giá trị ξ,η thoả msn ph−ơng trình Poisson (Joe F.Thompson, Z.U.A.Warsi, C.Wayne Mastin, 1985): 02 2 2 2 = ∂ ∂ + ∂ ∂ yx ξξ (2.6a) 02 2 2 2 = ∂ ∂ + ∂ ∂ yx ηη (2.6b) Tuy nhiên, khi tiến hành giải các ph−ơng trình (2.6) bằng ph−ơng pháp sai phân hữu hạn để xác định tính chất của dòng phun, ta phải tiến hành giải bài toán ng−ợc đó là xác định giá trị các toạ độ x,y t−ơng ứng với các giá trị các toạ độ ξ,η đs biết trong miền tính toán. Khi đó ph−ơng trình (2.6) trở thành: 02 2 22 2 2 = ∂ ∂ + ∂∂ ∂ − ∂ ∂ η γ ηξβξα yyy (2.7a) 02 2 22 2 2 = ∂ ∂ + ∂∂ ∂ − ∂ ∂ η γ ηξβξα xxx (2.7b) Trong đó các hệ số hình học a, b, g và ma trận Jacobi đ−ợc xác định nh− sau: 22       ∂ ∂ +      ∂ ∂ = ηη α yx (2.8a) ηξηξβ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ = yyxx (2.8b) 22       ∂ ∂ +      ∂ ∂ = ξξγ yx (2.8c) ξηηξ ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ = yxyxJ (2.8d) Giả sử các b−ớc l−ới Dξ = Dh = 1. áp dụng khai triển chuỗi Taylor, ta có: )( 2 1)( ,1,1, jijiji fff −+ −=ξ (2.9a) )( 2 1)( 1,1,, −+ −= jijiji fffη (2.9b) )2()( ,1,,1, jijijiji ffff −+ +−=ξξ (2.9c) )2()( 1,,1,, −+ +−= jijijiji ffffηη (2.9d) )( 4 1)( 1,11,11,11,1, −−−++−++ +−−= jijijijiji fffffξη (2.9e) Trong đó f ≡ x hoặc y và các chỉ số i và j t−ơng ứng liên hệ với ξ và h. Các biểu thức sai phân hữu hạn cho bởi ph−ơng trình (2.9) đ−ợc thay vào ph−ơng trình (2.7), ta có: )()(5.0)([5.0 1,1,1,11,11,11,1,1,1, −+−−−++−++−+ +++−−−++= jijijijijijijijiji fffffffff γβαγα (2.10) Trong đó đại l−ợng a, b, g và ma trận J đ−ợc coi là các hệ số và đ−ợc tính bằng sai phân hữu hạn sau khi làm trễ một b−ớc lặp: ( )21,1,, 4 1 −+ −= jijiji ffα (2.11) ( )( )1,1,,1,1, 4 1 −+−+ −−= jijijijiji ffffβ (2.12) ( )21,1,, 4 11 −+ −+= jijiji ffγ (2.13) ( )1,1,, 2 1 −+ −= jijiji ffJ (2.14) 3. KếT QUả TíNH TOáN Việc tính toán đ−ợc thực hiện trên máy tính, ch−ơng trình tạo l−ới tự động đ−ợc xây dựng bằng phần mềm Matlab 6.5 trong đó toạ độ các nút l−ới trong miền vật lý (x, y) hoàn toàn đ−ợc xác định t−ơng ứng với l−ới hình chữ nhật trong miền tính toán (ξ, η). Với giả thiết khoảng cách giữa các nút l−ới theo ph−ơng x là đều nhau, còn khoảng cách các nút l−ới theo ph−ơng y là không đều cho ta ma trận điểm và đồ thị chia l−ới sau khi chạy ch−ơng trình. Số l−ợng các nút l−ới theo ph−ơng x, y và dung sai cho phép (độ hội tụ) đ−ợc nhập vào theo yêu cầu của ng−ời sử dụng. Ma trận kết quả của biến y = f(x) a) Chia l−ới miền vật lý b) Chia l−ới miền tính toán Hình 3-1. Kết quả phân bố l−ới của dòng phun rối xoáy hai pha không đẳng nhiệt sau khi chạy ch−ơng trình 4. KếT LUậN ứng dụng kỹ thuật tạo l−ới cho phép ta giải bài toán mô phỏng dòng phun rối xoáy hai pha không đẳng nhiệt trong buồng đốt công nghiệp một cách rất dễ dàng và chính xác bằng ph−ơng pháp sai phân hữu hạn. Nghiệm nhận đ−ợc từ ch−ơng trình đạt độ chính xác mong muốn vì dung sai cho phép và số nút l−ới trên biên của miền vật lý là do ng−ời sử dụng ch−ơng trình trực tiếp nhập vào. Ch−ơng trình còn có thể ứng dụng trong các bài toán dẫn nhiệt trong mặt hình học không đều, đối l−u tự nhiên trong hình bao không đều... Ngoài ra kỹ thuật tạo l−ới còn có thể ứng dụng trong việc chia l−ới các miền vật lý có hình dạng phức tạp khác trong tự nhiên cũng nh− trong kỹ thuật. Tài liệu tham khảo Joe F.Thompson, Z.U.A.Warsi, C.Wayne Mastin (1985). Numerical Grid Generation Foundations and Applications, Elsevier Science Publishing Co. - Inc. pp. 7. P.D.Thomas, J.F.Middlecoff (1979). Direct Control of the Grid Point Distribution in Meshes Generated by Elliptic Equations. AIAA Journal Vol.18 - No.6. pp. 1462. M.Necati Ozisik (2000). Finite Difference Methods in Heat Transfer. CRC Press. Pp. 307 ữ 333. Courant, R. (1956). Differential and Integral Calculus. Blackie & Son, Ltd., London. Pp. 133. Nguyễn Hoài Sơn (chủ biên), Đỗ Thanh Việt, Bùi Xuân Lâm (2002). ứng dụng Matlab trong tính toán kỹ thuật. Nhà xuất bản Đại học Quốc gia Tp. HCM. Trang 13 ữ 49. Công trình đ−ợc sự hỗ trợ quí báu của ch−ơng trình nghiên cứu cơ bản trong khoa học tự nhiên. Các tác giả xin chân thành cảm ơn!