Sử dụng phương pháp biên nhúng (IBM) xây dựng mô hình 3D tính toán khí động cho cánh máy bay

Nghiên cứu khoa học công nghệ 
Tạp chí Nghiên cứu KH&CN quân sự, Số Kỷ niệm 55 năm Viện KHCNQS, 10 - 2015 101 
SỬ DỤNG PHƯƠNG PHÁP BIÊN NHÚNG (IBM) 
XÂY DỰNG MÔ HÌNH 3D TÍNH TOÁN KHÍ ĐỘNG 
CHO CÁNH MÁY BAY 
Mai Xuân Cảnh1*, Lã Hải Dũng2, Nguyễn Trang Minh3, Phạm Minh Vương4 
Tóm tắt: Trên cơ sở tính toán khí động cho cánh máy bay với mô hình 2D bằng 
phương pháp biên nhúng (IBM) đã được công bố [1], nhóm tác giả nghiên cứu xây 
dựng mô hình tính toán cho cánh với mô hình 3D. Trong bài báo này, phương pháp 
biên nhúng được tiếp tục cải tiến về mặt phương pháp tính, xây dựng mô hình 3D và 
thuật giải tính toán song song, để có thể tính được bài toán phức tạp với tốc độ tính 
nhanh hơn. Chương trình tính đã được thử nghiệm trên nhiều mô hình 3D mà đặc 
biệt là cánh máy bay mà các tác giả trình bày dưới đây. Kết quả tính toán cho thấy 
mô hình tính toán xác định được các hệ số khí động, cũng như mô phỏng được cấu 
trúc dòng chảy trên cánh với độ chính xác cho phép. 
Từ khóa: Tính toán khí động, Hệ số khí động, Cấu trúc dòng chảy, IBM 
1. ĐẶT VẤN ĐỀ 
Cơ sở lý thuyết của phương pháp biên nhúng và việc áp dụng phương pháp này để 
nghiên cứu tác động của dòng chảy lên vật thể hai chiều đã được các tác giả công bố [1]. 
Trong bài nghiên cứu trước, các tác giả đã trình bày các kết quả về vận tốc, đường dòng và 
phân bố áp suất ở nhiều góc tấn khác nhau. Các kết quả đã khẳng định độ chính xác của 
phương pháp biên nhúng dùng cho trường hợp dòng chảy qua cánh 2D. Trong phần nghiên 
cứu này, các tác giả tiếp tục ứng dụng phương pháp tính biên nhúng cho mô hình cánh 3D 
để có thể tính toán được dòng xoáy, phân bố áp suất cũng như áp suất trên cánh. Phần 
trình bày cụ thể về phương pháp và các kết quả tính sẽ được giới thiệu dưới đây. 
2. GIỚI THIỆU PHƯƠNG PHÁP TÍNH 
Phương pháp biên nhúng sử dụng cho nghiên cứu này được phát triển từ phần nghiên 
cứu tính toán khí động với mô hình cánh 2D mà các tác giả đã nghiên cứu trước đó [1]. 
Điều khác biệt là trong phần nghiên cứu này các tác giả xây dựng mô hình cho bài tính 
toán 3 chiều. Về mặt lý thuyết, việc áp dụng phương pháp biên nhúng để tính toán khí 
động cho cánh máy bay nhờ các điểm lực và điểm ảo. Đây là các điểm ở hai phía gần bề 
mặt. Chia miền tính toán thành hai miền: miền vật lý và miền ảo. Miền vật lý là miền của 
dòng chảy. Miền ảo nằm bên trong vật thể. Các điểm lực là các điểm của dòng chảy, 
nhưng không thể có được các giá trị tại các điểm này vì chúng ở gần bề mặt và sự rời rạc 
là không khả thi. Vì vậy, phương pháp IBM sẽ tính toán các giá trị tại các điểm này mà 
không rời rạc. Còn các điểm ảo là các điểm trong vật thể và chỉ nhằm mục đích phục vụ 
tính toán chứ không có nghĩa vật lý. Giá trị tính toán trên các điểm ảo được ngoại suy từ 
các giá trị trên biên và trên các điểm nằm trên miền lưu chất. 
Các phương trình tổng quát chú trọng đến sự bảo toàn mô men động lượng và khối 
lượng. Một thành phần lực trên biên fi được thêm vào phương trình mô men động lượng 
[4], [6]: 
 0,j
j
u
x


 (1) 
iji
i
j
Fu
f
t x

 
 (2) 
Tên lửa & Thiết bị bay 
M.X.Cảnh, L.H.Dũng, ,“Sử dụng phương pháp biên nhúng cánh máy bay.” 102 
Ở đây: 
- ui,j (i, j = 1, 2, 3) là thành phần vận tốc theo các trục; 
- xj (j=1, 2, 3) là tọa độ không gian 3 chiều. 
- Fij = uiuj + Pδij - ν( j j
j i
u u
x x
 
 
) (3) 
Trong đó: - P là áp suất dòng chảy có mật độ dòng ρ, độ nhớt động học ν = μ/ ρo; 
- δij là tensor đường chéo đơn vị (Kronecker delta). 
Việc giải phương trình mô men động lượng tại mỗi thời điểm phụ thuộc vào điều kiện 
biên. Lực là một hàm của thời gian phụ thuộc vào vị trí và vận tốc của dòng chảy. Thực tế 
xi không trùng khớp với các điểm nút, nhưng lực phải được ngoại suy từ các điểm nút này. 
Lực đặt trên biên nhúng fi bằng 0 trong vùng dòng chảy và khác không ở trong vùng điểm 
ảo [2], [3], [5]. Lực trên bề mặt để xác định biên vật thể được đặt vào trong miền tính được 
tính theo biểu thức dưới đây [1], [6]: 
1n n
i
V u
f R SH
t
 (4)
Trong đó: - Vn+1 là vận tốc tại điểm lực được nội suy từ vận tốc trên bề mặt của biên; 
- RSH là thành phần bên phải của phương trình Navier - Stockes trong đó có chứa 
thành phần đối lưu và lực nhớt; 
- un là vận tốc tại thời điểm n; Δt là bước thời gian tính toán. 
Trình tự giải bài toán bằng phương pháp biên nhúng được tiến hành theo các bước sau: 
Xác định các điểm lực; nội suy để tìm ra các giá trị lực cần thiết để áp đặt điều kiện biên; 
tính toán vận tốc ở bước dự đoán; giải phương trình Poisson để đảm bảo các điều kiện liên 
tục; cập nhật trường vận tốc cho các bước tiếp theo. 
3. SƠ ĐỒ THUẬT TOÁN 
Phương pháp giải bài toán theo phương pháp của Mac-Cormack, gồm hai bước tính 
toán vận tốc dự đoán (n) và tính toán vận tốc chính xác (n+1). Sơ đồ thuật toán tổng quát 
được thể hiện như trên hình 1. 
 Hình 1. Sơ đồ thuật toán tổng quát. 
 Sau khi xây dựng mô hình tính toán, các điểm tính được phân bố trên các máy tính để 
tính toán. Hai bước tính toán chính được thực hiện. Bước phỏng đoán dựa trên việc tính 
toán vận tốc dự đoán từ phương trình động lượng, ước tính lực trên các điểm lực và giải 
Tính vận 
tốc dự 
đoán 
Ước tính 
lực trên 
điểm lực 
Giải phương 
trình áp suất 
(Poisson) 
Giải ph.trình 
Poisson (dựa 
trên v.tốc c.xác) 
Tính lực trên điểm 
lực (dựa trên vận 
tốc chính xác) 
Xây dựng 
mô hình 
3D 
Phân bố 
điểm trên 
các CPU 
Hiệu chỉnh 
vận tốc 
chính xác 
K.tra đ/k 
ổn định 
Kết 
thúc 
Bắt 
đầu 
Sai 
Đúng 
Bước phỏng đoán 
Trao đổi dữ 
liệu song song 
Bước hiệu chỉnh 
Nghiên cứu khoa học công nghệ 
Tạp chí Nghiên cứu KH&CN quân sự, Số Kỷ niệm 55 năm Viện KHCNQS, 10 - 2015 103 
phương trình áp suất Poisson. Bước hiệu chỉnh được thực hiện nhờ việc trao đổi dữ liệu 
với bước phỏng đoán để tính toán lại vận tốc chính xác theo sơ đồ Runge - Kutta bậc 2, 
sau đó tính toán lực và giải phương trình Poisson dựa trên vận tốc chính xác này. Điều 
kiện ổn định được kiểm tra để kết thúc quá trình tính toán. 
4. MÔ HÌNH VÀ KẾT QUẢ TÍNH TOÁN 
4.1. Mô hình tính toán 
Tiến hành xây dựng mô hình số để tính toán khí động cho cánh máy bay có biên dạng 
kiểu NaCa 0012, với các thông số thiết kế: 
+ Vận tốc, v (m/s): 53 + Độ dài dây cung gốc cánh, bo (m): 2,15 
+ Độ cao, H (m): 500 + Độ dài dây cung mút cánh, bk (m): 1,08 
+Chiều dài nửa sải cánh, l (m): 4,65 + Độ dài dây cung khí động trung 
 bình, ba, (m): 
1,64 + Diện tích nửa cánh, S (m2): 7,5 
Các điều kiện biên tính toán được trình bày trên hình 1 trong đó vận tốc đầu vào là 53 
m/s. Biên trên và biên dưới của miền tính được sử dụng điều kiện biên đối xứng. Ở đầu ra, 
điều kiện biên áp suất ra được sử dụng. Miền tính được chia lưới đều với số phần tử lưới 
là: 105 106 phần tử; kích thước lưới dx = 0,0204 m, dy = 0,0125 m, dz= 0,027 m (trong mô 
hình c = 1,62 m). Bài toán được chạy trên hệ thống máy tính song song với 64 CPU và 
thời gian tính cho một bài toán là 7 ngày để có được kết quả tính toán ổn định. 
Hình 2. Mô hình tính toán. 
Hình 3. Hình lưới trên cánh 3D. 
4.2. Kết quả tính toán 
Tiến hành tính toán cho hai trường hợp với 2 góc tấn khác nhau là 2o và 12o. Phương 
pháp tính toán dòng không dừng đã được sử dụng với hai bước tính khác nhau. Trước tiên 
dòng khí được nghiên cứu trong một khoảng thời gian nhất định, gọi là thời gian chuyển 
Tên lửa & Thiết bị bay 
M.X.Cảnh, L.H.Dũng, ,“Sử dụng phương pháp biên nhúng cánh máy bay.” 104 
tiếp để đảm bảo dòng chảy vào được ổn định. Sau đó, bước kế tiếp là chạy để tính giá trị 
trung bình về lực nâng và lực cản. Phương pháp mô phỏng dòng rối bằng cách tính các 
cấu trúc lớn dùng mô hình rối LES đã được áp dụng để có thể bắt được các cấu trúc thực 
của dòng chảy. Phương pháp này cho phép bắt được hầu hết các cấu trúc của dòng và cho 
kết quả tốt nhất. Tuy nhiên, để có thể bắt được các cấu trúc dòng và tùy từng trường hợp 
bài toán, chúng tôi đã đạt đến giá trị y+ trên thành nhỏ hơn 10. Với phương pháp mô phỏng 
rối LES, thuật giải tính toán không dừng là bắt buộc để có thể bắt được các cấu trúc dòng 
như đã trình bày ở trên. 
Các kết quả về mô phỏng dòng quanh cánh, vận tốc, áp suất cũng như cường độ xoáy 
được trình bày dưới đây: 
Hình 4. Vận tốc ở các góc tấn khác nhau. 
Hình 5. Phân bố cường độ xoáy. 
Hình 6. Phân bố áp suất xung quanh biên dạng cánh. 
Hình 4 thể hiện kết quả phân bố vận tốc trên một mặt cắt 2D qua giữa sải cánh bên 
phải. Tương tự hình 5 thể hiện phân bố cường độ xoáy và hình 6 thể hiện phân bố áp suất. 
Kết quả tính cho thấy rằng dòng chảy sau cánh là dòng chảy rối và các cấu trúc dòng sau 
cánh đã được xé nhỏ thành các xoáy nhỏ khác nhau. Khi góc tới càng lớn thì vùng xoáy 
Nghiên cứu khoa học công nghệ 
Tạp chí Nghiên cứu KH&CN quân sự, Số Kỷ niệm 55 năm Viện KHCNQS, 10 - 2015 105 
sau cánh càng rộng như trong trường hợp góc tới là 12o. Phân bố áp suất trong trường hợp 
hình 6 cũng thể hiện được rằng áp suất phía dưới cánh và mũi cánh lớn hơn áp suất phía 
trên cánh. Đây cũng là nguyên lý cơ bản để tạo ra lực nâng và lực cản trên cánh. 
Hình 7 thể hiện phân bố hệ số áp suất Cp trên toàn bộ bề mặt cánh ở 2 góc tấn khác 
nhau là 2o và 12o. Ở đây: 
21
2
P
P
C
V 
 Hình 7. Phân bố hệ số áp suất trên cánh. 
Hình 8 thể hiện phân bố hệ số lực nâng (CL) và lực cản (CD) trên cánh theo thời gian ở 
hai góc tới khác nhau. Hệ số lực cản trung bình là 0,041 cho trường hợp góc tấn 2° và 0,16 
cho trường hợp góc tấn 12°. Hệ số lực nâng trung bình là 0.07 cho trường hợp góc tấn 2° 
và 0,6 cho trường hợp góc tấn 12°. Các kết quả tính toán được so sánh với một số kết quả 
nghiên cứu trước đó của của các tác giả trên các hình dạng khác nhau đã khẳng định được 
tính chính xác của chương trình tính [1], [4] với sai số nhỏ hơn 8,5%. 
α = 2o α = 12o 
Hình 8. Đồ thị hệ số lực nâng và lực cản theo thời gian. 
5. KẾT LUẬN 
Nội dung nghiên cứu trong bài báo này đã xây dựng được mô hình tính toán khí động 
trên cánh máy bay với mô hình 3D bằng phương pháp biên nhúng (IBM). Các bước lập 
trình cũng như chương trình tính đã được hoàn thiện. Kết quả cho thấy mô hình tính toán 
cho phép tính được các hệ số khí động, cũng như cấu trúc dòng chảy trên cánh. 
Các kết quả tính toán sẽ tiếp tục được kiểm chứng bằng mô hình thí nghiệm mà các tác 
giả sẽ tiến hành trong thời gian tới. 
Phương pháp tính toán dòng chảy bằng phương pháp biên nhúng, dùng thuật giải tính 
toán không dừng, trong trường hợp cánh 3D có thể không phải là một trong các phương 
pháp hiệu quả nhất. Tuy nhiên, phương pháp này có một ưu điểm mà các phương pháp 
tính truyền thống chưa thực hiện được, đó là với phương pháp biên nhúng, khi cácnh dịch 
chuyển thì không cần phải chia lại lưới và điều này đảm bảo độ chính xác cũng như hiệu 
Tên lửa & Thiết bị bay 
M.X.Cảnh, L.H.Dũng, ,“Sử dụng phương pháp biên nhúng cánh máy bay.” 106 
quả tính toán. Phương pháp này sẽ được nghiên cứu phát triển để tính cho cánh có cánh lái 
và tính phân bố áp suất trên cánh cho các bài toán đàn hồi khí động. 
TÀI LIỆU THAM KHẢO 
[1]. Lã Hải Dũng, Nguyễn Trang Minh, Mai Xuân Cảnh, Đào Duy Trường, Phạm Minh 
Vương (2013), "Xây dựng mô hình số tính toán khí động trên cánh máy bay bằng 
phương pháp biên nhúng (IBM)", Tuyển tập Hội nghị Khoa học Cơ học thủy khí, 
Quảng Bình (2013). 
[2]. Iaccarino, G. and Verzicco, R., "Immersed boundary technique for turbulent flow 
simulations", Appl. Mech. Rev., Vol. 56, pp. 331–347, 2003. 
[3]. Minh Vuong Pham, Frédéric Plourde, and Son Doan Kim, "Large-eddy simulation of 
a pure thermal plume under rotating conditions", Journal of Physics of Fluids, 
(2007). 
[4]. M.V. Pham, F. Plourde *, S.K. Doan (2008) "Turbulent heat and mass transfer in 
sinusoidal wavy channels", International Journal of Heat and Fluid Flow. 
[5]. Truong D. D, V. M. Pham, N. T. T. Le, "New modeling in the fluid-structure 
interaction flutter phenomenon analysis", Tuyển tập Hội nghị Cơ học toàn quốc lần 
thứ IX , Hà Nội, (2012). 
[6]. Yu-Heng Tseng , Joel H. Ferziger(2003), "A ghost-cell immersed boundary method 
for flow in complex geometry", Journal of Computational Physics, pp. 593–623. 
[7]. Nguyễn Chí Công, Phạm Minh Vương, Lê Thị Minh Nghĩa "Mô phỏng dòng qua 
hình trụ vuông bằng phương pháp LES-Mô hình SMAGORINSKY", Tuyển tập Hội 
nghị Khoa học Cơ học thủy khí, Hà Tiên (2004). 
ABSTRACT 
USING THE IMMERSED BOUNDARY METHOD (IBM) IN ORDER TO BUILD THE 
3D CACULATION MODEL FOR THE AERODYNAMIC OF AIRCRAFT WINGS 
Based on the aerodynamic calculations for wing with 2D models by Immersed 
Boundary Method (IBM) had been published [1], the authors study to build 
computational models for the wings with 3D models. In this paper, the immersed 
boundary method is further improved in terms of methodology, building 3D models 
and parallel computing algorithms to calculate the complex problems with faster 
calculated speed. The program has been tested on multiple 3D models which are 
particularly aircraft wing that the authors present below. The results show that the 
computational model had identified the aerodynamic coefficients, as well as 
simulating the flow structure on the wing with allowing precision. 
Keywords: Aaerodynamic calculations, Aerodynamic coefficients, Flow structure, IBM. 
Nhận bài ngày 15 tháng 07 năm 2015 
Hoàn thiện ngày 05 tháng 9 năm 2015 
Chấp nhận đăng ngày 10 tháng 9 năm 2015 
Địa chỉ: 
1Cục Kỹ thuật Quân chủng PK-KQ; 
 2 Viện Nghiên cứu phát triển Viettel; 
 3Viện Khoa học và Công nghệ quân sự; 
 4Công ty DFM-ENGINEERING. 
 *Email: canhmaixuan75@gmail.com