Mô phỏng khuếch tán chất trong môi trường rỗng kép không bão hòa bằng phương pháp đồng nhất hóa

 68 SCIENCE & TECHNOLOGY DEVELOPMENT JOURNAL, Vol 20, No.K7- 2017 
 Mô phỏng khuếch tán chất trong 
 môi trường rỗng kép không bão hòa 
 bằng phương pháp đồng nhất hóa 
 Trần Ngọc Tiến Dũng, Phạm Minh Quang, Nguyễn Thống, 
 Bùi Phạm Phương Thanh, Trần Văn Tiếng 
 Từ khóa— khuếch tán, truyền chất, môi trường 
 Tóm tắt— Khuếch tán chất là một quá trình then rỗng kép, đồng nhất hóa. 
chốt được chú trọng nghiên cứu trong nhiều lĩnh vực 
như khoa học vật liệu hay kỹ thuật môi trường. 
Khuếch tán của một chất hóa học trong một môi 1 MỞ ĐẦU 
trường rỗng thường được mô phỏng một cách truyền 
 ơ chế khuếch tán (diffusion) được Fick (1855) 
thống bằng định luật truyền chất Fick. Tuy nhiên, 
môi trường địa chất tự nhiên thường không đống C đề xuất, diễn ra khi có chênh lệch nồng độ và 
nhất và vì vậy cần có những mô hình tiên tiến cho được áp dụng rộng rãi trong nhiều lĩnh vực. 
phép mô tả ứng xử truyền chất trong môi trường Tuy nhiên, nó được quan tâm nhiều trong kỹ thuật 
như vậy. Quan niệm môi trường “rỗng kép” có thể địa chất môi trường từ những năm 1970, để giải 
đại diện cho một lớp môi trường rỗng không đồng quyết các bài toán xây dựng các công trình chôn 
nhất mà ở đó tồn tại hai miền rỗng lớn và bé với đặc giữ chất thải rắn sinh hoạt và nguy hại [1]. Hệ số 
trưng thủy lực rất khác biệt nhau, ví dụ đá nứt nẻ, thấm bé của môi trường địa chất (đất, đá nức nẻ) 
đất kết cục,... Trong bài báo này, chúng tôi trình bày chỉ là điều kiện cần nhưng chưa đủ để đảm bảo an 
cách phát triển một mô hình “vĩ mô” khuếch tán 
 toàn tránh rò rỉ. Khi đó chất ô nhiễm lan truyền 
chất trong môi trường đất rỗng kép không bão hòa 
nước bằng phương pháp đa tỷ lệ. Mô hình vĩ mô trong môi trường bị chi phối mạnh mẽ bởi cơ chế 
nhận được sẽ là một hệ phương trình phối hợp ở cả khuếch tán, hơn là những cơ chế lan truyền khác. 
hai thang tỷ lệ vi mô - vĩ mô và một ten-xơ khuếch Trong môi trường rỗng có hệ số thấm bé như vậy, 
tán hiệu dụng đại diện cho cả môi trường rỗng không rò rỉ của chất ô nhiễm sẽ nhanh hơn bởi sự khuếch 
đồng nhất. Mô hình phát triển được kiểm tra khi so tán chất, trong khi đó cơ chế đối lưu thuần túy 
sánh với lời giải tham chiếu thông qua một ví dụ số (convection), đối lưu với khuếch tán và phân tán 
3D của bài toán địa chất thủy văn. cơ học (dispersion) cho kết quả tính toán sai lệch 
 [2]. 
 Truyền chất khuếch tán trong môi trường rỗng 
 đồng nhất (homogeneous) chưa bão hòa cũng được 
 mô phỏng bằng định luật Fick [3], ở đó hệ số 
 Bản thảo nhận ngày 02 tháng 4 năm 2017, hoàn chỉnh sửa khuếch tán phụ thuộc vào độ khúc khuỷu của môi 
chữa ngày 16 tháng 10 năm 2017 trường rỗng mà nó có thể ước tính bằng công thức 
 Nghiên cứu này được tài trợ bởi Quỹ Phát triển khoa học và 
công nghệ Quốc gia (NAFOSTED) trong đề tài mã số 107.01- thực nghiệm Millington và Quirk [4], thông qua độ 
2015.25. ẩm của môi trường. Tuy nhiên, tồn tại những 
 Trần Ngọc Tiến Dũng, Phạm Minh Quang – CARE/Khoa trường hợp, đặc biệt có môi trường không đồng 
Môi trường và Tài nguyên, Trường Đại học Bách Khoa – 
ĐHQG-HCM, Việt Nam (e-mail: tntdung@hcmut.edu.vn). nhất (heterogeneous), mà chúng ta không thể áp 
 Nguyễn Thống - Khoa Kỹ thuật Xây dựng, Trường Đại học dụng định luật Fick “truyền thống” [5]. Thật vậy, 
Bách Khoa – ĐHQG-HCM. môi trường địa chất tự nhiên thường có cấu trúc 
 Bùi Phạm Phương Thanh - Khoa Khoa học Quản lý, Trường 
Đại học Thủ Dầu Một, Bình Dương. không đồng nhất do những tác động lý-hóa-sinh 
 Trần Văn Tiếng - Khoa Xây dựng, Trường Đại học Sư phạm chất chồng. Đối với một tập hợp môi trường địa 
Kỹ thuật TP. HCM 
TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K7-2017 69 
chất không đồng nhất mà được đặc trưng bởi 2 = x (x1, x2, x3) và y = y (y1, y2, y3) [-]: biến vi mô 
miền rỗng (rỗng lớn và rỗng bé) có tính chất thủy (y = X/ ℓ, biến nhanh) và biến vĩ mô (x = X/ L, 
lực và truyền chất khác biệt nhau như đất kết cục biến chậm), với X = X (X1, X2, X3) [L] là biến 
(aggregated soil) hay đá nức nẻ (fractured rock), không gian. Vì vậy toán tử vi phân: 
chúng ta có thể dùng quan niệm “môi trường rỗng ∇ y → ∇ y + ε ∇ x . (2) 
kép” [6] để đại diện chúng. Một mô hình toán học 
 Trong môi trường có cấu trúc vi mô (REV) 
khuếch tán trong môi trường composite (tương tự 
 mang tính chu kỳ, những biến đại lượng vật lý của 
môi trường rỗng kép) dưới điều kiện bão hòa đã bài toán được biễu diễn dưới dạng khai triển tiệm 
được giới thiệu trong [5]. Mô hình này nhận được cận: 
dưới dạng một phương trình, bằng phương pháp * * (0) * (1) * 2 (2) *
đồng nhất hóa tiệm cận [7], nhưng có khó khăn ψ (x, y,t ) =ψ (x, y,t ) + εψ (x, y,t ) + ε ψ (x, y,t ) + ...
chính trong việc xác định biến đổi Laplace trong (3) 
mô hình này. ở đây “*” chỉ định những đại lượng phi thứ 
 Mục đích của bài báo này là trình bày mô hình nguyên. Để đạt được mô hình ở tỷ lệ vĩ mô, quá 
hóa truyền chất khuếch tán trong một môi trường trình đồng nhất hóa trải qua những bước sau: 
rỗng kép không bão hòa nước bằng phương pháp - Bước 1: Mô tả bài toán vật lý ở tỷ lệ vi mô; 
đồng nhất hóa tiệm cận và thực thi số để kiểm tra - Bước 2: Chuẩn hóa các phương trình chủ đạo; 
mô hình nhận được. - Bước 3: Đánh giá số phi thứ nguyên đặc 
 trưng theo lũy thừa ε; 
 2 PHƯƠNG PHÁP ĐỒNG NHẤT HÓA VÀ - Bước 4: Khai triển tiệm cận cho mọi đại 
 THÀNH LẬP BÀI TOÁN. lượng vật lý phi thứ nguyên; 
 - Bước 5: Giải các bài toán giá trị biên cho các 
2.1 Phương pháp lũy thừa ε khác nhau; 
 Phương pháp đồng nhất hóa (Homogenization) - Bước 6: Quay lại mô tả bài toán vĩ mô với 
cho phép mô tả ứng xử của một môi trường không đầy đủ thứ nguyên. 
đồng nhất bằng cách thay thế nó bởi một môi 
 2.2 Giả thiết bài toán 
trường đồng nhất tương đương. Quá trình này 
được gọi là sự chuyển đổi thang tỷ lệ và nguyên lý Chúng ta xem xét hiện tượng khuếch tán của 
của phương pháp dựa trên những khai triển tiệm chất tan có nồng độ bé diễn ra trong một môi 
cận. Từ phương pháp toán học này [8,9], Auriault trường rỗng kép (double-porosity medium) có kích 
[7] đã đề xuất một cách tiếp cận gọi là đồng nhất thức L, mang tính chu kỳ, không bão hòa và không 
hóa vật lý (physical homogenization), được thực bị biến dạng (rắn). Môi trường này cấu thành từ 
hiện trên những số phi thứ nguyên đặc trưng cho những REV giống nhau có kích thước ℓ, chứa 2 
những quá trình vật lý liên quan. Phương pháp của loại vật liệu rỗng và thỏa mãn phương trình (1) 
GS. Auriault được sử dụng trong bài báo này và (Hình 1). Hai loại vật liệu rỗng này có đặc trưng 
giới thiệu tổng quát ở đây (chi tiêt có thể tìm đọc ở thủy lực và truyền chất rất khác biệt nhau: 
 2
[10]). D2 / D1 = O(ε ) (4) 
 2 -1
 Đồng nhất hóa bắt đầu từ việc mô tả những hiện với D1 và D2 [L T ] là hệ số khuếch tán của 
tượng vật lý ở thang tỷ lệ vật chất không đồng miền rỗng 1 và miền rỗng 2. 
nhất, mà thang tỷ lệ này rất bé so với thang tỷ lệ 
công trình mà chúng ta thực hiện. Giả thiết nền 
tảng của phương pháp đồng nhất là thỏa mãn thông 
số tách tỷ lệ (scale separation) của những đại Ω1
 Γ
lượng hình học và vật lý: Ω2
 ε =  / L << 1 (1) 
 ở đây ℓ và L [L] lần lượt là chiều dài đặc trưng ℓ
ở thang tỷ lệ “vi mô” và “vĩ mô”. Điều kiện (1) 
 L 
tương đương với việc tồn tại một thể tích phân tố Hình 1. Môi trường rỗng kép chu kỳ có 2 miền: miền rỗng 
đại diện (representative elementary volume, REV). Ω1 có tính thấm lớn và liên tục, miền rỗng Ω2 tính thấm kém 
Trong môi trường hiện diện 2 thang tỷ lệ, chúng ta và giao diện Γ. 
có thể đưa vào 2 biến không gian phi thứ nguyên x 
 70 SCIENCE & TECHNOLOGY DEVELOPMENT JOURNAL, Vol 20, No.K7- 2017 
 Trong mỗi vật liệu rỗng này, định luật Fick 2 2
 L θ1c  θ1c 2
được áp dụng để tính thông lượng trung bình. Bỏ tc = . Vì vậy số P1l = = O(ε ) và theo 
 D D t
qua dòng chảy trong môi trường (số Peclet rất bé). 1c 1c c
 2
Hiện tượng hấp phụ, cũng như là bất kỳ phản ứng  θ2c θ2c
 (4), P2l = = O(1) với = O(1) (1) 
sinh hóa nào không được xem xét ở đây. D2ctc θ1c
2.3 Mô tả khuếch tán ở thang tỷ lệ địa phương Những phương trình (10) – (13) được viết lại với 
 Phương trình chủ đạo của truyền chất khuếch những xem xét (14): 
tán trong mỗi miền rỗng của môi trường rỗng kép *
 2 ∂(θ1 C1 ) * *
không bão hòa như sau: ε = ∇ y ⋅ (D1 ⋅∇ yC1 ) trong miền Ω1, (2) 
 ∂t*
 ∂(θ1C1 ) *
 = ∇ X ⋅(D1(θ1)⋅∇ X C1 ) trong miền Ω1, (5) ∂(θ2C2 ) * *
 ∂t = ∇ y ⋅ (D2 ⋅∇ yC2 ) trong miền Ω2, (3) 
 ∂t*
 ∂(θ2C2 )
 = ∇ ⋅(D (θ ) ⋅∇ C ) trong miền Ω2, (6) * * 2 * *
 ∂t X 2 2 X 2 (D1 ⋅∇ yC1 )⋅ N = ε ⋅ (D2 ⋅∇ yC2 )⋅ N trên Γ, (4) 
cùng với điều kiện liên tục tại giao diện: * *
 C1 = C2 trên Γ, (5) 
 (D (θ )⋅∇ C )⋅ N = (D (θ )⋅∇ C )⋅ N 
 1 1 X 1 2 2 X 2 
 trên Γ, (7) 
 2.6 Phân tích đồng nhất hóa 
 C1 = C2 trên Γ, (8) 
 ở đây C [ML-3] là nồng độ chất tan; D [L2T-1] là Mục đích của việc phân tích là đưa ra phương trình 
ten-xơ khuếch tán; [L3L-3] là độ ẩm môi trường; t chủ đạo chứa thông số có hiệu mô tả truyền chất 
[T] là thời gian; và N [-] là vec-tơ pháp tuyến đơn khuếch tán cho một môi trường liên tục từ môi 
 trường thực tế rỗng kép. 
vị với chiều dương hướng ra ngoài Ω . “1” và “2” 
 1 Thay thế khai triển tiệm cận (3) vào (15) – (18) và 
lần lượt là những chỉ số cho miền rỗng lớn và bé. 
 nhóm những số hạng có cùng bậc ε, dẫn đến: 
2.4 Chuẩn hóa phương trình chủ đạo ∂ θ *C(0)
 ε 2 ( 1 1 ) =
 Bài toán vật lý truyền chất khuếch ở thang tỷ lệ ∂t*
địa phương trên đây sẽ được biễu diễn theo những * (0) * (1) * (0)
 = ∇ y ⋅[D1 ⋅∇ yC1 + ε (D1 ⋅∇ yC1 + D1 ⋅∇xC1 )]+ 
đại lượng không chiều: 2 * (2) * (1) * (0)
 + ∇ y ⋅[ε (D1 ⋅∇ yC1 + D1 ⋅∇xC1 )]ε ∇x ⋅ D1 ⋅∇ yC1 +
 * C1 * C2 * D1 * D2
C = ; C = ; D = ; D = ; 2 * (1) (0)
 1 2 1 2 + ε ∇ ⋅ D ⋅ (∇ C + ∇ C )
 Cc Cc D1c D2c x 1 y 1 x 1
 θ θ t trong miền Ω1, (6) 
θ * = 1 ; θ * = 2 và t* = , (9) 
 1 2 ∂ θ * C(0) + ε C(1) + ε 2 C(2)
 θc θc tc [ 2 ( 2 2 2 )]
 * =
 Chỉ số “c” biểu thị cho các đại lượng đặc trưng ∂t
 * (0) (1) 2 (2)
mà chúng ta quan sát được. Từ (9) những phương = ∇ y ⋅ [D2 ⋅∇ y (C2 + ε C2 + ε C2 )]+
trình (5) – (8) được chuẩn hóa đối với chiều dài vi * (0) 2 (1)
 + ∇ y ⋅ [D2 ⋅∇x (ε C2 + ε C2 )]+
mô ℓ: 
 + ε ∇ D* ⋅∇ C(0) + ε 2∇ D* ⋅∇ C(1) + ε 2∇ D* ⋅∇ C(0)
 2 * x 2 y 2 x 2 y 2 x 2 x 2
  θ1c ∂(θ C ) * *
 1 1 = ∇ ⋅ D ⋅∇ C trong miền Ω1, (10) 
 * y ( 1 y 1 ) trong miền Ω2, (20) 
 D1c tc ∂t
  * (0) (1) 2 (2) 
 2 * D1 ⋅∇ y (C1 + ε C1 + ε C1 )+
  θ2c ∂(θ C ) * *
 2 2 = ∇ ⋅ (D ⋅∇ C ) trong miền Ω2, (11)   ⋅ N =
 D t * y 2 y 2  * (0) 2 (1)  
 2c c ∂t + D1 ⋅∇ x (ε C1 + ε C1 ) 
 * * D2c * * 2 * (0)
 (D ⋅∇ C )⋅ N = ⋅ (D ⋅∇ C )⋅ N trên Γ, (12) = [ε D2 ⋅∇ y C2 ]⋅ N
 1 y 1 D 2 y 2
 1c trên Γ , (21) 
C* = C* trên Γ, (13) 
 1 2 (C (0) + ε C (1) + ε 2 C (2) )= (C (0) + ε C (1) + ε 2 C (2) ) 
 1 1 1 2 2 2
 trên Γ , (22) 
2.5 Đánh giá số phi thứ nguyên Thực hiện đồng nhất hóa (19) – (22) theo bậc ε sẽ 
 Chúng ta xem xét thời gian đặc trưng: nhận được các bài toán cho ra những kết quả quan 
TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K7-2017 71 
tâm. Chi tiết quá trình biến đổi của các bài toán 3.3 Bài toán vĩ mô – Mô hình khuếch tán trong 
được trình bày trong [11]. môi trường rỗng kép 
 Để xác định mô hình thang tỷ lệ vĩ mô, chúng ta 
 3 KẾT QUẢ VÀ THẢO LUẬN. viết phương trình (19) và (21) ở bậc ε2: 
3.1 Bài toán biến nồng độ thang vĩ mô và khuếch * (0)
 ∂(θ1 C1 ) * (2) * (1) * (1) (0)
 tán thang vi mô = ∇ y ⋅ (D1 ⋅∇ yC1 + D1 ⋅∇ xC1 )+ ∇ x ⋅ D1 ⋅ (∇ yC1 +∇ xC1 ) 
 ∂t*
 Ở bậc ε0, phương trình (19) – (22) trở thành: 
 trong miền Ω1, (13) 
 * (0) * (2) * (1) * (0)
 ∇ y ⋅ (D1 ⋅∇ yC1 )= 0 trong miền Ω1, (23) 
 (D1 ⋅∇ yC1 + D1 ⋅∇ xC1 )⋅ N = (D2 ⋅∇ yC2 )⋅ N 
 * (0)
 ∂(θ2C2 ) * * trên Γ , (14) 
 = ∇ y ⋅ (D2 ⋅∇ yC2 ) trong miền Ω2, (7) 
 ∂t* Thực hiện hàng loạt các biến đổi cho (34) [11], 
 * (0) chúng ta nhận được: 
 (D1 ⋅∇ yC1 )⋅ N = 0 trên Γ, (8) 
 * (0) * (0)
 (0) (0) ∂(θ1 C1 ) * (0) 1 ∂(θ 2C2 )
 C = C trên Γ, (9) φ1 = ∇ x ⋅ (Deff ⋅∇ xC1 )− d Ω 
 1 2 ∂t* Ω ∫ ∂t*
 Ω
 Lời giải của bài toán (23) và (25) có dạng [12]: 2
 (15) 
 (0) (0) *
 C1 = C1 (x,t ) , (10) với ten-xơ khuếch tán có hiệu của môi trường rỗng 
 (0) kép được định nghĩa như sau: 
 nghĩa là ở bậc đầu tiên C1 đã là biến vĩ mô vì 
 1
không phụ thuộc vào y, và trở thành C (0) . Phương D* = D* ⋅(∇ χ + I )d Ω (16) 
 eff Ω ∫ 1 y
trình (24) mô tả bài toán khuếch tán ở tỷ lệ địa Ω1
 (0)
phương (vi mô) và C2 không là hằng số trong Ω1
 và φ1 = . (17) 
miền Ω2. Như vậy, ta có 2 biến nồng độ vĩ mô Ω
 (0) (0)
C và vi mô C2 , tương tác với nhau qua điều Cuối cùng quay lại bài toán vật lý thứ nguyên, 
kiện (26). chúng ta đạt được mô hình khuếch tán vĩ mô 
3.2 Bài toán giá trị biên địa phương không cân bằng địa phương với C1 và C2 với một 
 thông số có hiệu Deff của môi trường: 
 Phương trình (19) và (21) được phân tích ở bậc 
 1
ε khi kể đến lời giải (27): ∂(θ1C1 ) ∂ θ 2C2
 φ1 = ∇ X ⋅ (D ⋅∇ xC1 )− , (18) 
 * (1) * (0) ∂t eff ∂t
∇ y ⋅(D1 ⋅∇ yC1 + D1 ⋅∇ xC1 )= 0 trong miền Ω1, 
 ∂(θ2C2 )
 (11) = ∇ X ⋅(D2 ⋅∇ X C2 ) , (40) 
 * (1) * (0) ∂t
(D1 ⋅∇ yC1 + D1 ⋅∇ xC1 )⋅ N = 0 trên Γ , (12) 
 C1 = C2 trên Γ , (41) 
Hệ phương trình (28) và (29) có nghiệm dưới 1
 với D (θ ) = D (θ )⋅(∇ χ + I )d Ω . (42) 
dạng: eff Ω ∫ 1
 (1) (0) (1) * Ω1
C1 = χ( y)∇ xC1 + C1 (x,t ) , (30) 
 Phương pháp đồng nhất hóa đưa đến kết quả 
ở đây χ( y) là hàm vectơ chu kỳ và là một 
 nhiều bài toán khác nhau. Nó cho phép giải thích 
 (1) *
C1 (x,t ) hàm bất kỳ. Vectơ χ( y) nhận được từ mối quan hệ giữa các tỷ lệ, tùy thuộc vào độ lớn 
việc giải bài toán giá trị biên địa phương sau: của thông số tách tỷ lệ ε. Mô hình nhận được (39) 
 * – (42) khác với những mô hình giới thiệu trước 
 ∇ y ⋅[D1 ⋅(∇ y χ + I )]= 0 trong miền Ω1, (31) 
 đây (ví dụ [5]) khi cho phép mô phỏng bài toán 
 *
 [D1 ⋅(∇ y χ + I )]⋅ N = 0 trên Γ , (32) môi trường chưa bão hòa. Để kiểm tra tính dự 
 1 đoán mô hình này, chúng ta sẽ so sánh với lời giải 
 χ = χ d Ω = 0 ( Ω = Ω + Ω ) (33) 
 Ω ∫ 1 2 tham chiếu (reference solution) của mô hình không 
 Ω1 đồng nhất hóa (mô hình 1 phương trình khuếch tán 
 truyền thống) để đối chứng thông qua một ví dụ 
 số. 
 Môi trường rỗng kép của ví dụ số 3D bao gồm 
 72 SCIENCE & TECHNOLOGY DEVELOPMENT JOURNAL, Vol 20, No.K7- 2017 
những viên sét hình tròn có bán kính R sắp xếp trường rỗng kép là Deff = 9.05 x 10-10 m2/s, được 
đều đặn trong đất cát (Hình 2). Điều kiện biên và tính toán từ bài toán giá trị biên địa phương (31) và 
ban đầu được tóm tắt trên Hình 2. Thông số hình (32). (Không trình bày chi tiết ở đây). 
học và vật lý của 2 loại vật liệu được trình bày 
trong Bảng 1. Hệ số khuếch tán có hiệu của môi 
 t > 0 t > 0
 X =0 X =L
 C = 1 g/l C = 0
 t = 0; C = 0
 L = 0.02 m
 Hình 2. Kích thước hình học và tình huống truyền chất của ví dụ số. 
 X = 0
 X = L
 Thang 
 t = 1 × 103 s t = 2× 104 s t = 1 × 104 s t = 1 × 106 s
 nồng độ 
 Hình 3. Nồng độ chất trong môi trường ở thời gian khác nhau tính toán từ mô hình đối chứng. (Xem bản số của bài báo để phân 
biệt thang nồng độ với mã màu khác nhau.) 
TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K7-2017 73 
 1.00 1.00
 Mô hình đối chứng
 Mô hình đối chứng
 Mô hình đồng nhất hóa_miền vĩ mô
 0.75 Mô hình đồng nhất hóa_miền vĩ mô 0.75
 t = 1e3 = t s
 0.50 =2e4t s 0.50
 ở
 ở
 C [g/l] C [g/l] 0.25
 C [g/l] C [g/l] 0.25
 0.00
 0.000 0.005 0.010 0.015 0.020 0.00
 0.000 0.005 0.010 0.015 0.020
 x [m] x [m]
 1.00 1.00
 Mô hình đối chứng Mô hình đối chứng
 Mô hình đồng nhất hóa_miền vĩ mô Mô hình đồng nhất hóa_miền vĩ mô
 0.75 Mô hình đồng nhất hóa_miền vi mô 0.75
 t = 1e6 = t s
 t = 8e4 = t s 0.50 0.50
 ở
 ở
 C [g/l] C [g/l]
 C [g/l] C [g/l] 0.25 0.25
 0.00 0.00
 0.000 0.005 0.010 0.015 0.020 0.000 0.005 0.010 0.015 0.020
 x [m] x [m]
 Hình 4. Nồng độ chất ở vị trí trung tâm môi trường ở thời gian khác nhau tính toán từ mô hình đối chứng và nồng độ vĩ mô từ mô 
hình đồng nhất hóa. (Nồng độ chất ở miền vi mô chỉ trình bày ờ t = 8e4 s, ở t = 1e3s hay 2e4 s rất bé (~ 0), và ở t = 1e6 s nồng độ ở 
2 miền xấp xỉ nhau (trạng thái cân bằng).) 
 5.E-08
 Mô hình đối chứng
 1]
 - 4.E-08 Mô hình đồng nhất hóa
 2 s
 -
 [g [g m 3.E-08
 2.E-08
 x = 0.02 m 0.02 m x =
 1.E-08
 Thông lượng tại khuếch tán Thông
 0.E+00
 1.E-01 1.E+01 1.E+03 1.E+05 1.E+07
 t [s] 
 Hình 5. Thông lượng khuếch tán chất quan sát ở lối ra của môi trường mô phỏng từ mô hình đối chứng và mô hình đồng nhất hóa. 
 Bảng 1: Thông số hình học và vật lý của 2 vật liệu trong môi trường rỗng kép ví dụ. 
 Thông số L ℓ R φ1 φ2 θ1 θ2 D1 D2 
 Đơn vị m m m - - m3/m3 m3/m3 m2/s m2/s 
 Giá trị 0.02 0.002 0.0005 0.935 0.065 0.5 0.3 1×10-13 1×10-9 
 Thực thi số (numerical implementation) của mô phần tử hữu hạn phát triển với Mathlab. Thời gian 
hình đồng nhất hóa có thể tham khảo ở [13]. Miền giải là < 1 s cho mô hình đồng nhất và 16 s cho mô 
vĩ mô được phân chia thành 50 phẩn tử và 1200 hình không đồng nhất với bộ xử lý của máy tính 
phần tử tam giác cho miền vi mô, trong khi đó mô 3.30 GHz và RAM 4.0 GB. 
hình không đồng nhất được giải ở tỷ lệ mịn (fine Kết quả nồng độ chất ở một số thời gian được 
scale) bằng việc phân chia miền bài toán thành chọn giới thiệu, Hình 3. Hình 4 cho chúng ta thấy 
~13400 phần tử tứ diện (tetrahedral). Cả hai bài nồng độ chất mô phỏng bởi mô hình đối chứng dao 
toán của mô hình đồng nhất và không đồng nhất động vì sự khuếch tán lần lượt qua miền vật liệu 1 
được thực thi bằng một phần mềm thương mại (cát) và miền vật liệu 2 (sét). Kết quả mô hình 
 74 SCIENCE & TECHNOLOGY DEVELOPMENT JOURNAL, Vol 20, No.K7- 2017 
đồng nhất hóa cho ra đường cong nồng độ vĩ mô – 
vi mô trùng với kết quả mô hình đối chứng lần TÀI LIỆU THAM KHẢO 
lượt ở miền 1 và 2. Trung bình nồng độ ở một mặt 
cắt tại những vị trí X sẽ rất giống nhau từ 2 mô [1] Shackelford C.D., “The ISSMGE Kerry Rowe Lecture: 
hình, Hình 5. Như vậy mô hình phát triển trong bài The role of diffusion in environmental geotechnics”, 
báo này không chỉ được kiểm chứng, mà cách thực Can. Geotech. J., vol. 51, no. 11, pp. 1219–1242, 2014. 
thi số cho mô hình đồng nhất ở đây chứng tỏ hiệu [2] C.D. Shackelford and P.L. Redmond, “Solute 
quả cho phép tính toán nổng độ ở cả hai tỷ lệ vĩ breakthrough curves for processed kaolin at low flow 
 rates”, Journal of Geotechnical Engineering, vol. 121, 
mô và vi mô. no. 1, pp. 17-32, 1995. 
 Thông qua ví dụ số trên đây, chúng ta nhận thấy [3] Rowe R.K., Geotechnical and Geoenvironmental 
mô hình đồng nhất cho phép biết được sự khuếch Engineering Handbook, Springer, New York, 2001. 
tán chất ở những thang tỷ lệ hiện diện trong môi [4] R.J. Millington and J.M. Quirk, “Permeability of porous 
trường. Ngoài ra, mô hình đồng nhất vĩ mô giúp solids”, Transaction of Faraday Soc., vol. 57, pp. 1200–
tiết kiệm tài nguyên, đặc biệt trong trường hợp mô 1207, 1961. 
phỏng cho bài toán thực tế. [5] J.L. Auriault and J. Lewandowska, “Non-gaussian 
 diffusion modeling in composite porous medium by 
 homogenization: tail effect”, Transp. Porous Media, vol. 
 4 KẾT LUẬN. 21, pp. 47–74, 1995. 
 Mô hình truyền chất khuếch tán trong môi [6] G. Barenblatt, I. Zheltov and I. Kochina, “Basic concepts 
trường rỗng kép không bão hòa đã được phát triển in the theory of seepage of homogeneous liquids in the 
 fissured rocks”, J. Appl. Math, vol. 24, no. 5, pp. 1286–
bằng cách áp dụng phương pháp đa thang tỷ lệ 1303, (1960). 
đồng nhất hóa. Mô hình nhận được là một hệ [7] J.L. Auriault, “Heterogeneous medium. Is an equivalent 
phương trình cho phép mô tả khuếch tán ở cả hai macroscopic description possible?”, Int. J. Engg. Sci., 
thang tỷ lệ vi mô và vĩ mô mà nồng độ của chúng vol. 29, no. 7, pp. 785–795, 1991. 
tương tác với nhau thông qua biên tiếp xúc giữa [8] A. Bensoussan, J.L. Lions and G. Papanicolaou, 
chúng. Sự truyền chất qua lại giữa hai miền trong Asymptotic Analysis for Periodic Structures, North-
 Holland, Amsterdam, 1987. 
môi trường và hệ số khuếch tán có hiệu đại diện 
cho toàn bộ môi trường hiện diện trong mô hình [9] Sanchez-Palencia E., “Non-homogeneous Media and 
 Vibration Theory”, Lecture Notes in Phys. Springer-
một cách chính xác mà không chịu bất kỳ điều Verlag, Berlin, 1980. 
kiện áp đặt nào. Mô hình đề xuất được xác nhận có [10] J.L. Auriault, C. Boutin and C. Geindreau, 
khả năng mô phỏng truyền chất khuếch tán trong “Homogenization of Coupled Phenomena in 
môi trường rỗng kép không đồng nhất, thông qua Heterogenous Media”, ISBN: 978-1-84821-161-2, 476 
việc so sánh kết quả nồng độ của mô hình đồng pages, Wiley-ISTE, 2009. 
nhất hóa và mô hình truyền thống. [11] N. T.D. Tran, Transport de solutés dans un milieu à 
 double-porosité non saturé. Modélisation par 
 Cuối cùng, phương pháp mô hình hóa giới thiệu homogénéisation et application, Thèse doctorale, 
trong bài báo này có thể mở rộng cho những Université Joseph Fourier-Grenoble, 185 p (in French), 
trường hợp truyền chất khác tính đến sự biến đổi 2008. 
hóa sinh của chất ô nhiễm, cũng như là áp dụng [12] Auriault J.L., Upscaling by multiscale asymptotic 
cho những bài toán đa vật lý trong môi trường địa expansions, CISM Lecture 480 Applied micromechanics 
 of porous materials, L. Dormieux and F.-J. Ulm Eds 
chất nhiều pha và nhiều thành phần. Springer, 2005. 
 [13] N.T.D. Tran, J. Lewandowska, M. Vauclin and H. 
 Bertin, “Two-scale modeling of solute dispersion in 
 unsaturated double-porosity media: Homogenization and 
 experimental validation”, Int. J. Numer. Anal. Meth. 
 Geomech., vol. 35, no. 14, pp. 1536-1559, 2011. 
TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K7-2017 75 
Trần Ngọc Tiến Dũng tốt nghiệp kỹ sư xây dựng Nguyễn Thống tốt nghiệp tiến sĩ Viện ĐH Bách 
Trường Đại học Bách Khoa, ĐHQG-HCM 2001, khoa Grenoble (Pháp) 1991 và đạt học hàm PGS 
thạc sĩ kỹ thuật Môi trường ĐH INSA Lyon (Pháp) 2005. Hiện là chủ nhiệm bộ môn Kỹ thuật Tài 
2004 và tiến sĩ Khoa học Trái đất và Môi trường nguyên Nước, Khoa Kỹ thuật Xây dựng, Trường 
ĐH Grenoble (Pháp) 2008. Là nghiên cứu viên- Đại học Bách Khoa – ĐHQG-HCM. 
giảng viên Trung tâm Châu Á Nghiên cứu về 
Nước (CARE), Trường Đại học Bách Khoa, Bùi Phạm Phương Thanh tốt nghiệp thạc sĩ 
ĐHQG-HCM. TS. Dũng quan tâm thực hiện các 2013, là giảng viên Trường Đại học Thủ Dầu Một, 
nghiên cứu liên ngành về Địa chất Thủy văn và Bình Dương từ năm 2010. 
Thủy văn Đô thị. 
 Trần Văn Tiếng tốt nghiệp tiến sĩ Đại học 
Phạm Minh Quang nghiên cứu viên Trung tâm Grenoble (Pháp) 2011, là giảng viên Trường Đại 
CARE, Trường Đại học Bách Khoa, ĐHQG-HCM học Sư phạm Kỹ thuật TP. HCM từ 2012. Hướng 
đến 2017. Nay ThS. Quang công tác tại Công Ty nghiên cứu chính của TS. Tiếng là phương pháp 
CP Nhựa Tân Phú. phần tử rời rạc để mô phỏng ứng xử của địa vật 
 liệu bão hòa và chưa bão hòa. 
 A model of solute diffusion in unsaturated 
 double-porosity medium by 
 homogenization 
 Tran Ngoc Tien Dung, Pham Minh Quang, Nguyen Thong, 
 Bui Pham Phuong Thanh, Tran Van Tieng 
Abstract— Solute diffusion is a key process in many fields like for example material science or 
environmental engineering. Diffusion mechanism in porous media is often described by Fick’s law. 
However, we could not use this law for nonstandard diffusion behaviors occurring in cases of 
heterogeneous media. The conception of double-porosity medium can be applied to a class of such media. 
The double-porosity medium is characterized by two distinct pore sizes: macro-porosity domain and 
micro-porosity domain, respectively, having the contrasted hydraulic properties. This paper presents the 
development of a macroscopic model for the solute diffusion in unsaturated double-porosity medium, by 
using homogenization method. The obtained macroscopic model is a system of two equations coupling on 
the interface of the macro- and micro-porosity domain for diffusion. This model contains the effective 
diffusion tensor representing for the entire medium. The developed model is verified by comparing with 
the reference solution of the fine scale model through a 3D numerical example of hydrogeology problem. 
Index Terms— diffusion, mass transport, double-porosity medium, homogenization.