PHƯƠNG PHÁP TẬP MỨC KHÔNG LƯỚI: CƠ SỞ TOÁN HỌC VÀ KHẢ NĂNG ỨNG DỤNG TRONG NGÀNH DẦU KHÍ
Độ tin cậy và thiết thực của việc mô phỏng một quá trình vật lý không những phụ thuộc vào mô hình toán học mô tả quá trình, thường ở dạng những phương trình vi phân, mà còn phụ thuộc vào độ chính xác và tính hiệu quả của phương pháp số dùng để giải các phương trình vi phân đó.
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
PHƯƠNG PHÁP TẬP MỨC KHÔNG LƯỚI: CƠ SỞ TOÁN HỌC VÀ KHẢ
NĂNG ỨNG DỤNG TRONG NGÀNH KỸ THUẬT DẦU KHÍ
MESHLESS LEVEL SET METHOD: MATHEMATICAL
FUNDAMENTALS AND POTENTIAL APPLICATIONS IN
PETROLEUM ENGINEERING
Mai Cao Lân*, Trần Công Thành**
* Khoa Kỹ thuật Địa chất & Dầu khí, Đại học Bách khoa Tp. Hồ Chí Minh, Việt Nam
** University of Southern Queensland, Australia
---------------------------------------------------------------------------------------------------------------------------
TÓM TẮT
Độ tin cậy và tính thiết thực của việc mô phỏng một quá trình vật lý không những phụ thuộc vào
mô hình toán học mô tả quá trình, thường ở dạng những phương trình vi phân, mà còn phụ thuộc vào
độ chính xác và tính hiệu quả của phương pháp số dùng để giải các phương trình vi phân đó. Bài báo
này trình bày cơ sở lý thuyết một phương pháp số mới mang tên phương pháp Tập mức Không lưới
(Meshless Level set method) trong đó những tính năng ưu việt của 2 nhóm phương pháp không lưới
(meshless) và tập mức (level set) được tích hợp để giải các bài toán biên di động. Một số bài toán mẫu
giải bằng phương pháp này được trình bày trong bài báo để minh họa cho độ chính xác và tính hiệu
quả của nó cũng như khả năng ứng dụng của phương pháp trong ngành kỹ thuật dầu khí.
ABSTRACT
The reliability and usefulness of the numerical simulations of a physical process depend not only
on the mathematical model representing the process, normally in forms of differential equations, but
also on the accuracy and efficiency of the numerical methods for solving such equations. This paper
presents the theoretical basics of a new numerical method, namely Meshless Level Set method, in
which the advantageous features of meshless methods and level set methods are integrated to solve
moving boundary problems. Some benchmark problems solved by the method are presented to
demonstrate the accuracy and efficiency of the method as well as its potential applications in
petroleum engineering.
thông thường không thể dùng phương pháp giải
1. GIỚI THIỆU
tích được. Thay vào đó, người ta sử dụng các
Đa số mô hình toán mô tả một quá trình vật phương pháp số để tìm nghiệm gần đúng của
lý thường ở dạng các phương trình vi phân. Đối chúng. Hiện nay các phương pháp số được sử
với bài toán đa biến, ta có các phương trình vi dụng phổ biến gồm có phương pháp sai phân
phân riêng phần. Việc tìm nghiệm của những hữu hạn (finite difference method - FDM), phần
phương trình này nói chung là phức tạp nên tử hữu hạn (finite element method - FEM), khối
175
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
hữu hạn (finite volume method - FVM), v.v…. động (moving boundaries). Có hai nhóm
Xin xem [Tannehill et al. (1997), Chung (2002)] phương pháp số được sử dụng cho các bài toán
để biết thêm chi tiết. Các phương pháp này được dạng này: Nhóm phương pháp dựa trên lưới di
gọi chung là phương pháp rời rạc hóa theo động và nhóm phương pháp sử dụng lưới cố
không gian. Đối với các bài toán phụ thuộc thời định. Phương pháp Tập mức (level set method)
gian, ta cần thêm công cụ số để rời rạc hóa thuộc nhóm phương pháp thứ hai, do Osher and
phương trình vi phân theo biến thời gian. Xin Sethian (1988) đề xuất. Phương pháp này ban
xem [Quarteroni and Valli (1994), Quarteroni et đầu được thiết lập để sử dụng với nhóm các
al. (2000)] để biết thêm chi tiết về các phương phương pháp dựa vào lưới như FDM, FEM,
pháp này. FVM [Sethian (1999), Osher and Fedkiw
(2003)]. Trong bài báo này, phương pháp tập
Nếu như các phương pháp FDM, FEM,
mức được triển khai trên nền tảng của phương
FVM, v.v… rời rạc hóa phương trình vi phân
pháp không lưới IRBFN.
trên cơ sở chia nhỏ miền tính toán thành một
lưới (mesh) gồm những phần tử ràng buộc lẫn
nhau trên lưói theo những nguyên tắc xác định 2. CƠ SỞ TOÁN HỌC
(ta gọi chung các phương pháp này là nhóm
2.1. Phương pháp Tập mức
phương pháp dựa vào lưới) thì đối với các
phương pháp Không lưới, miền tính toán được Trong phương pháp Tập mức, biên di động
chia thành một tập hữu hạn các điểm rời rạc, có Γ(t) của miền Ω ⊂ ℜ2 được xem là tập mức
thể bố trí tùy ý (unstructured) và không có bất không (zero) của một hàm φ(x,t), gọi là hàm tập
kỳ mối ràng buộc nào về vị trí tương đối giữa mức, trong không gian ℜ3
chúng trong quá trình tính toán. Kết quả là các
phương pháp không lưới rất thích hợp cho các Γ(t ) = {x ∈ ℜ 2 | φ ( x, t ) = 0} (1)
bài toán có biến dạng lớn (như trong cơ học rạn Hàm φ(x,t) có thể chọn tùy ý với điều kiện
nứt) hoặc các bài toán có biên di động (như dự phải là hàm trơn. Trong [Sethian (1999), Osher
đoán quá trình điền khuôn đúc hoặc mô phỏng
and Fedkiw (2003)] , φ(x,t) được chọn là hàm
mặt tiến dầu-nước/khí-dầu trong quá trình bơm
khoảng cách sao cho
ép/thu hồi tăng cường dầu) trong khi đối với các
+
phương pháp dựa vào lưới, việc giải các bài toán ⎧+ d ( x, t ), x ∈ Ω
⎪
này sẽ rất phức tạp (đôi khi làm giảm độ chính φ ( x, t ) = ⎨ 0 x∈Γ (2)
xác của lời giải) do phải thường xuyên điều ⎪− d ( x, t ), x ∈ Ω −
⎩
chỉnh lưới bị biến dạng trầm trọng. Có nhiều
phương pháp không lưới [Kansa (1990a,b), Trong đó d(x,t) là khoảng cách từ điểm x đến
Aluri (2002)], trong đó có phương pháp Indirect biên di động; Ω+ và Ω- là miền bên ngoài và bên
Radial Basis Function Networks (IRBFN) [Mai- trong biên tương ứng. Như vậy, trong phương
Duy and Tran-Cong (2001,2003)] dùng để giải pháp tập mức, đối tượng nghiên cứu là hàm tập
các phương trình vi phân không lệ thuộc thời mức φ(x,t) chuyển động với vận tốc “mở rộng”
gian. Phương pháp này gần đây đã được mở (extended velocity) V thay vì là biên Γ(t) di
rộng để giải các bài toán phụ thuộc thời gian chuyển với tốc độ F [Osher and Sethian (1988)].
[Mai-Cao and Tran-Cong (2003,2004,2005)]. Phương trình chuyển động của hàm tập mức
Các phương pháp số để giải bài toán biên di tương ứng với dịch chuyển của biên trong
động đã và đang được các nhà nghiên cứu quan trường vận tốc V của môi trường xung quanh
tâm vì tính phức tạp của bản thân các biên di như sau:
176
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
∂φ w(t ) = G −1U (t ) vào phương trình (4), ta có
+ V ⋅ ∇φ = 0 (3)
∂t công thức xấp xỉ hàm
Ở một thời điểm bất kỳ, thông tin về biên di uˆ ( x, t ) = g T ( x)G −1U (t ) (6)
động (vị trí, hình dáng, độ cong, v.v…) có thể
được tái tạo từ hàm tập mức φ(x,t) bằng cách Trong đó G là ma trận được xác định bằng
xác định tập hợp các đoạn trên Γ(t) sao cho cách áp dụng (5) tại M điểm trong miền tính
φ(x,t) triệt tiêu. toán với tập các hàm cơ sở đã cho g(x). Trong
phương trình (6), giá trị hàm tại các điểm nút
Do phương trình (3) được giải bằng phương
U(t) là biến cần tìm. Đạo hàm bậc nhất và bậc
pháp số nên chỉ sau một bước thời gian φ(x,t) sẽ hai của hàm u(x,t) được xấp xỉ bằng cách lấy
không còn là hàm khoảng cách. Vì vậy việc tái đạo hàm phương trình (6) tương ứng:
thiết lập hàm tập mức thỏa điều kiện (2) là một
T
bước cần thiết và được thực hiện bằng cách tìm uˆ , j ( x, t ) = g , j ( x)G −1U (t ) (7a)
lời giải dừng (steady) cho bài toán sau [Sussman
T
et al. (1994)] uˆ , j ..l ( x, t ) = g , j ...l ( x)G −1U (t ) ) (7b)
∂φ Nếu như trong phương pháp Kansa (1990a),
= S ε (φ )(1− | ∇φ |)
∂t (4a) hàm cơ sở g(x) trong phương trình (5) được
φ ( x , t = 0) = φ ( x ) chọn là hàm multiquadrics (MQ) thì trong
IRBFN, g(x) là đạo hàm bậc k của hàm
ở đó Sε là một hàm trơn sao cho multiquadrics hoặc thin plate splines (TPS). Đa
φ số các bài toán trong lĩnh vực Cơ học Chất lỏng
S ε (φ ) = (4b) Tính toán (Computational Fluid Dynamics -
φ 2 +ε2 CFD) được giải với k=2. Chi tiết về cở sở lý
với ε là khoảng cách ngắn nhất giữa một thuyết của phương pháp IRBFN cũng như ứng
điểm bất kỳ với các điểm khác trong miền tính dụng của nó để giải các toán phụ thuộc thời gian
toán. Cơ sở lý thuyết cũng như các ứng dụng đã được trình bày trong [Mai-Cao and Tran-
tiêu biểu của phương pháp này được trình bày Cong (2005)].
chi tiết trong [Sethian (1999), Osher and Fedkiw
3. PHƯƠNG PHÁP TẬP-MỨC KHÔNG-
(2003)].
LƯỚI VÀ CÁC BÀI TOÁN MINH HỌA
2.2. Phương pháp Không lưới IRBFN
Quy trình giải một bài toán biên di chuyển bị
Xấp xỉ uˆ ( x, t ) của hàm u(x,t) có thể được động dưới tác dụng của trường vận tốc không
viết ở dạng tổ hợp tuyến tính của N hàm cơ sở đổi bằng phương pháp Tập mức không lưới bao
gồm các bước sau:
N
u(x, t) ≈ uˆ(x, t) = ∑wi (t)gi (x) = gT (x)w(t) (5) Bước 1: Xây dựng hàm tập mức ban đầu là
i=1
hàm khoảng cách thỏa phương trình (2);
Trong đó g(x)=[g1(x),g2(x),…,gN(x)]T là tập
Bước 2: Thực hiện dịch chuyển hàm tập mức
các hàm cơ sở cho trước; w(t)=[w1(t),…,wN(t)]T
trong một bước thời gian bằng cách giải phương
là tập N trọng số cần tìm. Với một tập hợp M
trình (3);
điểm trong miền tính toán và giá trị hàm tương
ứng tại các điểm đó tại thời điểm t, U(t)=[U1(t), Bước 3: Tái thiết lập hàm tập mức thỏa
U2(t),…,UM(t)], bằng cách thay phương trình (2) bằng cách tìm lời giải dừng cho
phương trình (4a,b). Thông tin về biên di động ở
177
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
thời điểm đang xét có thể tái tạo bằng giải thuật u=-sin(πx) cos(πy)
lấy đường đồng mức zero của hàm φ(x,t);
v=-cos(πx) sin(πy)
Bước 4: Quay lại bước 2 cho bước thời gian
Kết quả mô phỏng ở nhiều thời điểm khác
kế tiếp hoặc kết thúc quá trình khi thời gian mô
nhau được trình bày trong hình 1. Ở mỗi thời
phỏng đạt đến giá trị giới hạn cho trước.
điểm, giải thuật trích đường đồng mức zero của
Trong bài báo này, các phương trình vi phân hàm tập mức cho ta biên dạng bọt có dạng đa
ở bước 2 và 3 được giải bằng các lược đồ số dựa giác khép kín. Diện tích của hình đa giác này
trên phương pháp IRBFN mô tả trong [Mai-Cao chính là diện tích của bọt ở thời điểm tương
and Tran-Cong (2005)]. ứng. Kết quả tính toán cho thấy tỉ lệ phần trăm
3.1. Bài toán bọt xoay tròn thay đổi về diện tích của hình tròn trong suốt
quá trình mô phỏng không vượt quá 2% với mật
Xét một bọt hình tròn bán kính r=0.15 ban độ điểm trong miền tính toán là 32 x 32.
đầu được đặt tại vị trí (0.5,0.7) trong miền chữ
nhật [0,1] x [0,1] có trường vận tốc xoáy (u,v)
được xác định như sau:
Hình 1: Bài toán bọt xoay tròn
178
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
Hình 2: Bài toán bọt xoay tròn (tiếp theo)
giữa các bọt được mô phỏng hoàn toàn theo quy
3.2. Bài toán 4 bọt di động trong dòng chảy
trình 4-bước tổng quát mô tả ở trên mà không
xoáy
cần phải xử lý cho từng trường hợp riêng biệt
Các bọt ban đầu được bố trí ngẫu nhiên như như trong các phương pháp truyền thống khác.
ở hình 2 trong một trường vận tốc xoáy giới hạn
trong miền [-1,1] x [-1,1]. Các hình bên trái của 4. KẾT LUẬN & HƯỚNG PHÁT TRIỂN
hình 2 thể hiện biên di động là đường đồng mức CỦA ĐỀ TÀI
zero (màu xanh dương, trong cùng) ở các thời
Phương pháp Tập mức Không lưới được xây
điểm ban đầu và t=4.1333. Bên phải là hàm tập
dựng trên cơ sở triển khai phương pháp Tập
mức ở các thời điểm tương ứng trên đó biên
mức trên nền không lưới của phương pháp
dạng của 4 bọt di động được gắn vào. Như vậy
IRBFN. Qua các bài toán mẫu trình bày trong
thay vì theo dõi sự chuyển động và biến dạng
bài báo, phương pháp mới cho thấy độ chính xác
của bản thân 4 bọt di động, ta quan sát hàm tập
và tính hiệu quả cao của nó khi giải các bài toán
mức di chuyển theo quy luật (3) và trích đường
phụ thuộc thời gian, trong đó mô hình các bọt di
đồng mức zero của nó để có biên dạng của các
động có thể được mở rộng để mô phỏng chế độ
bọt ở thời điểm cần quan tâm. Với phương pháp
dòng chảy của hỗn hợp dung dịch trong ống
Tập mức Không lưới, sự kết dính và tách rời
khai thác.
179
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
Hình 3: Bài toán 4 bọt di động trong dòng chảy xoáy
Đây chính là hướng phát triển của đề tài 3. Kansa, E.J. Multiquadrics - A Scattered
trong đó có xét tới sự phân chia thành những bọt Data Approximation Scheme with
thứ cấp cũng như sự kết hợp giữa các bọt khí Applications to Computational Fluid-
trong quá trình đi từ đáy giếng lên bề mặt. Ngoài Dynamics--I. Surface Approximations and
ra, các mô hình chất lỏng phi Newton (non- Partial Derivative Estimates, Computers and
Newtonian fluid) cũng sẽ được xem xét để mô tả Mathematics with Applications 19 (1990a),
ứng xử phức tạp của hỗn hợp dung dịch khai pp. 27-145.
thác.
4. Kansa, E.JMultiquadrics - A Scattered Data
TÀI LIỆU THAM KHẢO Approximation Scheme with Applications to
Computational Fluid-Dynamics--II.
1. Atluri, S.N. and Shen, S. The Meshless
Solutions to Parabolic, Hyperbolic and
Local Petrov-Galerkin (MLPG) Method,
Elliptic Partial Differential Equations,
Tech Science Press, Encino, USA (2002).
Computers and Mathematics with
2. Chung, T.J. Computational Fluid Dynamics, Applications 19 (1990b), pp. 147-161
Cambridge University Press, UK (2002).
180
Hội nghị khoa học và công nghệ lần thứ 9, Trường Đại học Bách khoa Tp. HCM, 11/10/2005
5. Mai-Cao, L. and Tran-Cong, T. Solving Quarteroni, A. and Sacco, R. and Saleri, F.:
Time-Dependent PDEs with a Meshless Numerical Mathematics, Vol. 37 of Texts in
IRBFN-based Method. In: Alves, C.J.S. and Applied Mathematics, Springer-Verlag,
Chen, C.S. and Leitao, V. (eds): New York (2000).
International Workshop on MeshFree
13. Sethian, J.A. Level Set Methods and Fast
Methods, July 21-23, Lisbon, Portugal
Marching Methods: Evolving Interfaces in
(2003)
Computational Geometry, Fluid Mechanics,
6. Mai-Cao, L. and Tran-Cong, T. Element- Computer Vision, and Materials Science,
Free Simulation for non-Newtonian Flows. Cambridge University Press, New York
In: Atluri, S.N. and Beskos, D.E. and (1999).
Polyzos, D. (eds): International Conference
14. Sussman, M. and Smereka, P. and Osher,
on Computational & Experimental
S.J. A Level Set Approach for Computing
Engineering & Sciences, ICCES, July 26-
Solutions to Incompressible Two-Phase
29, Madeira, Portugal (2004).
Flow, Journal of Computational Physics 114
7. Mai-Cao, L. and Tran-Cong, T. Meshless (1994), pp.146-159.
IRBFN-Based Method for Transient
15. Tannehill, J.C. and Anderson, D.A. and
Problems, Computer Modeling in
Pletcher, R.H. Computational Fluid
Engineering & Sciences 7 (2005), pp. 149-
Mechanics and Heat Transfer, Taylor &
171.
Francis,USA (1997).
8. Mai-Duy, N. and Tran-Cong, T. Numerical
Solution of Differential Equations Using
Multiquadric Radial Basis Function
Networks, Neural Networks 14 (2001),
pp.185-199.
9. Mai-Duy, N. and Tran-Cong, T.
Approximation of Function and its
Derivatives Using Radial Basis Function
Networks, Applied Mathematical Modelling
27 (2003), pp. 197-220.
10. Osher, S. and Fedkiw, R.: Level Set
Methods and Dynamic Implicit Surfaces,
Springer, New York (2003).
11. Osher, S. and Sethian, J.A. Fronts
Propagating with Curvature-Dependent
Speed: Algorithms Based on Hamilton-
Jacobi Formulations, Journal of
Computational Physics 79 (1988), pp. 12-
49.
12. Quarteroni, A. and Valli, A.: Numerical
Approximation of Partial Differential
Equations, Springer-Verlag, New York
181
