Nghiên cứu dự báo sa bồi luồng sông Hậu theo phương pháp của Van Rijn's Method
Tóm tắt
Trong công tác chỉnh trị cửa sông ven biển phục vụ cho phát triển kinh tế ngành cảng, đường
thủy ở Việt Nam, việc xác định được mức độ sa bồi tại khu vực bể cảng bao gồm cửa vào
và khu nước là yếu tố then chốt để đánh giá được hiệu quả đầu tư. Trong những năm gần
đây nghiên cứu đánh giá sa bồi đã có những chuyển biến tích cực giải quyết được phần lớn
các yêu cầu trên cơ sở nền tảng cơ sở khoa học cao. Các nghiên cứu đã phát triển dần từ
việc đánh giá sa bồi dựa trên phương pháp truyền thống (chập bình đồ đo đạc thực tế) đến
việc kết hợp với sử dụng mô hình toán với các phần mềm tiên tiến. Tuy nhiên mỗi phương
pháp đều có các ưu, nhược điểm riêng. Phương pháp truyền thống chỉ có thể đánh giá sa
bồi trong điều kiện hiện trạng mà không đáp ứng được việc đánh giá sa bồi đối với các
phương án đầu tư trong tương lai. Phương pháp nghiên cứu trên mô hình toán đã dần được
nghiên cứu và đưa vào sử dụng để đáp ứng được cả hai yêu cầu trên nhưng đối với phương
pháp này đòi hỏi bộ số liệu đầu vào phải được đo đạc đầy đủ và đồng bộ (bao gồm các số
liệu địa hình, địa chất, thủy hải văn). Điều này dẫn tới chi phí dành cho công tác điều tra cơ
bản là rất tốn kém và kéo dài thời gian nghiên cứu. Đây là nhược điểm rất lớn cho những
trường hợp phát sinh đòi hỏi thời gian nghiên cứu nhanh và kinh phí không nhiều, thực tế
này đã xảy ra trong quá trình thực hiện một số dự án lớn của Việt Nam như dự án “Đầu tư
xây dựng luồng cho tầu biển trọng tải lớn vào sông Hậu”.
Xuất phát từ lý do trên nhóm nghiên cứu đã tiếp cận phương pháp lý thuyết bán kinh nghiệm
để giải quyết nhanh bài toán sa bồi (schematized box method). Kết quả nghiên cứu của
phương pháp này được so sánh với kết quả thực đo của các công ty tư vấn thực hiện và kết
quả so sánh đã cho thấy sự phù hợp của cơ sở lý thuyết với số liệu thực đo tại các khu vực
cảng nghiên cứu
Tóm tắt nội dung tài liệu: Nghiên cứu dự báo sa bồi luồng sông Hậu theo phương pháp của Van Rijn's Method
CHÚC MỪNG NĂM MỚI 2019 40 Tạp chí khoa học Công nghệ Hàng hải Số 57 - 01/2019 NGHIÊN CỨU DỰ BÁO SA BỒI LUỒNG SÔNG HẬU THEO PHƯƠNG PHÁP CỦA VAN RIJN'S METHOD RESEARCH ON THE SILTATION OF HAU RIVER'S CHANNEL ACCORDING TO VAN RIJN'S METHOD LÊ THỊ HƯƠNG GIANG, ĐỖ THỊ MINH TRANG, LÊ THỊ LỆ Khoa Công trình, Trường Đại học Hàng hải Việt Nam Email liên hệ: gianglh.ctt@vimaru.edu.vn Tóm tắt Trong công tác chỉnh trị cửa sông ven biển phục vụ cho phát triển kinh tế ngành cảng, đường thủy ở Việt Nam, việc xác định được mức độ sa bồi tại khu vực bể cảng bao gồm cửa vào và khu nước là yếu tố then chốt để đánh giá được hiệu quả đầu tư. Trong những năm gần đây nghiên cứu đánh giá sa bồi đã có những chuyển biến tích cực giải quyết được phần lớn các yêu cầu trên cơ sở nền tảng cơ sở khoa học cao. Các nghiên cứu đã phát triển dần từ việc đánh giá sa bồi dựa trên phương pháp truyền thống (chập bình đồ đo đạc thực tế) đến việc kết hợp với sử dụng mô hình toán với các phần mềm tiên tiến. Tuy nhiên mỗi phương pháp đều có các ưu, nhược điểm riêng. Phương pháp truyền thống chỉ có thể đánh giá sa bồi trong điều kiện hiện trạng mà không đáp ứng được việc đánh giá sa bồi đối với các phương án đầu tư trong tương lai. Phương pháp nghiên cứu trên mô hình toán đã dần được nghiên cứu và đưa vào sử dụng để đáp ứng được cả hai yêu cầu trên nhưng đối với phương pháp này đòi hỏi bộ số liệu đầu vào phải được đo đạc đầy đủ và đồng bộ (bao gồm các số liệu địa hình, địa chất, thủy hải văn). Điều này dẫn tới chi phí dành cho công tác điều tra cơ bản là rất tốn kém và kéo dài thời gian nghiên cứu. Đây là nhược điểm rất lớn cho những trường hợp phát sinh đòi hỏi thời gian nghiên cứu nhanh và kinh phí không nhiều, thực tế này đã xảy ra trong quá trình thực hiện một số dự án lớn của Việt Nam như dự án “Đầu tư xây dựng luồng cho tầu biển trọng tải lớn vào sông Hậu”. Xuất phát từ lý do trên nhóm nghiên cứu đã tiếp cận phương pháp lý thuyết bán kinh nghiệm để giải quyết nhanh bài toán sa bồi (schematized box method). Kết quả nghiên cứu của phương pháp này được so sánh với kết quả thực đo của các công ty tư vấn thực hiện và kết quả so sánh đã cho thấy sự phù hợp của cơ sở lý thuyết với số liệu thực đo tại các khu vực cảng nghiên cứu. Từ khóa: Bồi lắng, luồng tàu, cảng, bê cảng, phương pháp hình hộp. Abstract In training estuaries for port development and naviation purposes in Vietnam, the estimation of sediment deposition in harbor basins, including harbor entrances, is crucial for evaluating the effectiveness of investment. In recent years, studies on sediment deposition has remarkable advances, meeting most requirements with scientifically proven methods. The researches has evolved from traditional method (overlaying surveyed bathymetries) to mathematical modelling with the applications of advanced computer programs. However each method has its own advantage and disadvantages. The traditional method can only predict sediment deposition in present situation and not suitable for future investment projects. The methods based on mathematical modelling is recently more concerned and applied to meet the above two requirements, but this method requires complete and synchronized surveyed data (including bathymetry, geology, and oceanography). This has lead to huge cost of basic survey, and lengthen the research time. This is a major disadvantage for ad-hoc situations where time and money budgets are limited. This type of situation happens during the implementation of certain large projects in Vietnam e.g. the investment project on excavation of a navigation channel for high capacity vessels into Hau river. For the above reason the author adopted a semi-empirical method (aka schematized box method) to quickly solve the sediment deposition problem. The calculated result is compared with surveyed data carried out by consultant companies, and showed an agreement between the result with measured data in the basin. Keywords: Siltation, channel, port, basin, schematized box method. 1. Giới thiệu về dự án nạo vét luồng sông Hậu [1] Luồng cho tầu biển trọng tải lớn vào sông Hậu được triển khai xây dựng vào năm 2009 với hạng mục nạo vét và thi công kè bờ (hai bên bờ kênh) với chiều dài đoạn kênh khoảng 4300 m. Quá trình thực hiện đã gặp phải khó khăn về vốn nên đã phải dừng thi công trong thời gian dài, từ 2012 đến cuối năm 2014 mới triển khai tiếp. Khi đó việc thi công nạo vét luồng tầu đã hoàn thành được 90%. Đến tháng 7 năm 2014 khi có đủ điều kiện về vốn, gói thầu được tiếp tục thực hiện. Tuy nhiên CHÚC MỪNG NĂM MỚI 2019 Tạp chí khoa học Công nghệ Hàng hải Số 57 - 01/2019 41 đã gặp phải vấn đề lớn về sa bồi. Khối lượng bồi lại rất lớn với khối lượng khoảng 1,16 triệu m3. Việc xác định nguyên nhân sa bồi trở lại mang ý nghĩa rất lớn về kinh tế của Nhà thầu và Chủ đầu tư, kinh phí để nạo vét lại khối lượng sa bồi này được xác định như thế nào, Chủ đầu tư hay Nhà thầu chịu trách nhiệm. Để xác định được trách nhiệm chi trả kinh phí phải xác định đượng khối lượng bồi lắng tự nhiên trong quá trình dừng giãn thi công. Với những phương pháp xác định sa bồi đang được áp dụng, khối lượng sa bồi tự nhiên này có thể xác định được bằng phương pháp mô phỏng trên mô hình toán. Tuy nhiên để thực hiện được cần phải có các số liệu đo đạc tự nhiên đồng bộ (địa hình, thủy văn mẫu đất mặt....), kinh phí ước tính khoảng 3 tỷ đến 4 tỷ đồng. Ngoài ra phải tốn thời gian khảo sát, đo đạc 5 đến 6 tháng. Trong trường hợp này cả thời gian và kinh phí đều không khả thi. Trước tình hình này nhóm nghiên cứu đã tiếp cận phương pháp bán kinh nghiệm của Van Rijn để xác định sa bồi trong thủy vực cảng. 2. Mô hình nghiên cứu của Vanrjin Căn cứ vào điều kiện dòng chảy - sóng, vận chuyển bùn cát và các vấn đề sa bồi có thể lựa chọn một phương pháp tổng hợp để tính toán sa bồi. Với bài toán thuỷ vực một cửa vào: - Nếu muốn tính tốc độ (hay bề dày) bồi lắng trung bình trong bể cảng thì sử dụng phương pháp mô hình bán kinh nghiệm đơn giản; - Nếu muốn tính phân bố sa bồi trong bể cảng thì cần dùng phương pháp mô hình toán 2 chiều về thuỷ động lực kết hợp với sóng. Sau đó sử dụng kết quả trường dòng chảy này để phân tích diễn biến vận chuyển bùn cát cũng như bồi xói bể cảng; - Nếu muốn mô phỏng với độ chính xác cao hơn nữa thì cần tiến hành thí nghiệm mô hình vật lý, tuy nhiên chi phí tương đối đắt; - Các mô hình toán 3 chiều mặc dù phản ánh được bản chất vật lí của hiện tượng rối động nhưng mô hình chạy lâu và không thể lập luận độ chính xác vì khó tiến hành thu thập đủ và đúng số liệu tại hiện trường. Phương pháp của Van Rijn (2016) dựa trên cơ sở giản hóa địa hình bể cảng thành một hình hộp (‘Bể cảng 1’ trong Hình 1). Ngoài ra, có thể bổ sung một hộp phụ nữa (‘Bể cảng 2’) trong trường hợp thủy vực kéo dài và co hẹp một chút. Trong thời gian duy trì dòng triều dâng giữa các mực nước thấp và mực nước cao, nước tiến vào thủy vực. Còn khi dòng triều rút giữa các mực nước cao và mực nước thấp, nước rút khỏi thủy vực. Khi đó, cân bằng khối lượng trong bể cảng ‘Basin 1’ có thể được viết như sau: - Đối với dòng triều dâng: A×d(hc)/dt = (Qtide1 + Qhor + Qvert1 + Qvert2)co − (Qhor + Qvert1 + Qvert2)c – Acws,eff − Qtide2c + pQriver (criver2 – c) = Qtide1 co + Qhor(co − c) + (Qvert1 + Qvert2)(co − c) – Acws,eff − Qtide2c + pQriver (criver2 – c) (1) Trong đó các thành phần (số hạng) ở vế phải được giải thích như sau: - Trầm tích có nồng độ co vào thủy vực theo dòng triều, do trao đổi nước theo phương ngang và thẳng đứng; - Trầm tích có nồng độ c thoát ra khỏi cửa thủy vực, do trao đổi nước theo phương ngang và đứng, và dòng chảy sông; - Trầm tích có nồng độ c đọng xuống đáy thủy vực; - Trầm tích có nồng độ c thoát sang bể phụ 2; - Trầm tích có nồng độ criver2 vào bể cảng theo dòng chảy sông. - Đối với dòng triều rút: d(hc)/dt = (Qhor + Qvert1 + Qvert2)co − (Qtide1 + Qhor + Qvert1 + Qvert2)c − Acws,eff − Qtide2c2 + pQriver (criver2 – c) = –Qtide1 c + Qhor(co − c) + (Qvert1 + Qvert2)(co − c) − Acws,eff − Qtide2c + pQriver (criver2 – c) (2) Trong đó các thành phần (số hạng) ở vế phải được giải thích như sau: - Trầm tích có nồng độ co vào thủy vực do trao đổi nước theo phương ngang và đứng; - Trầm tích có nồng độ c thoát ra khỏi cửa thủy vực, do thoát nước, trao đổi nước theo phương ngang và đứng, và dòng chảy sông; - Trầm tích có nồng độ c2 vào thủy vực từ bể phụ (c2 = αco với α = 0,1 đến 0,5); - Trầm tích có nồng độ c đọng xuống đáy thủy vực; - Trầm tích có nồng độ criver2 vào bể cảng theo dòng chảy sông. CHÚC MỪNG NĂM MỚI 2019 42 Tạp chí khoa học Công nghệ Hàng hải Số 57 - 01/2019 Hình 1. Bể cảng (mặt bằng và mặt cắt ngang) giản hóa dùng cho mô phỏng [5] Lưu ý: Thủy vực đang xét là Bể cảng 1 (Basin 1) nơi có bồi lắng, khu chứa phụ là Bể cảng 2 (Basin 2), ngoài biển là ‘Outside’ nơi có dao động mực nước triều; Inflow of fresh water and mud: dòng nước ngọt và bùn chảy vào (từ sông 2). Các giá trị lưu lượng nước có thể được tính bởi: Triều dâng: Qtide1 = (Abasin + Abasin2) dηo /dt Qtide2 = Abasin2 dηo /dt Qhor = f1,flood b h uo Qvert1 = 0,5 f3 b h [(Δρo /ρo)gh]0,5 Qvert2 = p (fmix - 1) Qriver2 Triều rút: Qtide1 = (Abasin + Abasin2) dηo /dt Qtide2 = Abasin2 dηo /dt Qhor = f1,ebb b h uo (thường lấy f1,ebb = 0 Qvert1 = 0,5 f3 b h [(Δρo /ρo)gh]0,5 Qvert2 = p (fmix - 1) Qriver2 Kì triều dâng được xác định là thời gian có giá trị dηo /dt > 0 (hay mực nước tăng dần bắt đầu từ t = 0). Theo đó, các tham số triều cơ bản có thể được biểu diễn dưới dạng: ηo = - ηo,max cos(ωt); uo = ur - uo,max cos(ωt); Δρo = - Δρo,max cos(ωt); co = co1 - co2 cos(ωt). Còn tốc độ chìm lắng hiệu quả của trầm tích được cho bởi: ws,eff = αws ws,o Ý nghĩa của các thông số và biến số đầu vào như sau: Abasin1 = diện tích bể cảng; Abasin2 = diện tích khu chứa phụ; ηo = mực nước triều so với mực biển trung bình; ηo,max = biên độ mực nước triều; uo = vận tốc trung bình độ sâu của dòng chảy sông+triều ngoài cửa cảng; uo,max = biên độ của vận tốc triều bên ngoài và song song với cửa vào; ur = vận tốc dòng chảy sông trung bình độ sâu bên ngoài (nếu có) và song song với cửa vào; b = bề rộng cửa vào; bbasin = bề rộng bể cảng; ho = độ sâu nước dưới MN biển trung bình (hằng số); co = mật độ trầm tích trung bình độ sâu ngoài bể cảng; c2 = mật độ trầm tích trung bình độ sâu trong bể phụ; ws,o = tốc độ chìm hạt trầm tích trong nước tĩnh; ucr,d = lưu tốc phân giới cho bồi lắng (≈ 0,3 m/s); Hs/h = chiều cao sóng tương đối trong bể cảng (≈ 0,05 đến 0,15); f1,flood = hệ số trao đổi phương ngang kì triều dâng (≈ 0,025); f1,ebb = hệ số trao đổi phương ngang kì triều rút (≈ 0) ; f3 = hệ số trao đổi phương đứng (= 0,3); Δρo = chênh lệch mật độ nước ngoài - trong cảng; Δρo,max = biên độ của chênh lệch nêu trên; ρo = mật độ nước bên ngoài; p = phần trăm thời gian có dòng chảy sông; fmixing = hệ số hoà trộn; CHÚC MỪNG NĂM MỚI 2019 Tạp chí khoa học Công nghệ Hàng hải Số 57 - 01/2019 43 Qriver2 = lưu lượng dòng chảy sông vào thuỷ vực; ρd = mật độ khô của bùn lắng đọng; Tmud,lat = lượng vận chuyển bùn dọc bờ thuỷ vực (nếu có bờ là đất bùn); Lm = chiều dài bờ đất bùn ở đó có dòng bùn chảy vào. Các biến lượng trong quá trình tính toán bao gồm: h = ho + ηo = độ sâu nước (cả trong và ngoài bể cảng); Qtide1 = lưu lượng dòng triều qua cửa vào; Qtide2 = lưu lượng dòng triều qua cửa ra (exit), tức là cửa nối với bể phụ; Qhor = lưu lượng trao đổi do dòng chảy qua cửa vào; Qvert1 = lưu lượng trao đổi do dòng mật độ (giá trị tuyệt đối); Qvert2 = lưu lượng trao đổi do dòng chảy qua cửa vào; ubasin = Qtide1/(h bbasin) = lưu tốc trong bể cảng; c01 = hằng số nồng độ trầm tích ngoài cửa cảng; c02 = biên độ biến thiên nồng độ bùn ngoài cửa vào; criver2 = mật độ bùn trung bình của sông 2. Biến số kết quả là: c = nồng độ trầm tích trung bình độ sâu trong bể cảng. Tổng lượng trao đổi nước (chảy vào và ra) qua cửa vào bể cảng do những xoáy cỡ lớn và do chênh lệch mật độ nước thì bằng 0. Thể tích nước chảy vào luôn bằng thể tích nước chảy ra khi không có dòng chảy trao đổi. Vì nồng độ bùn bên trong bể cảng (c) nhỏ hơn ở bên ngoài (co) nên luôn có dòng bùn liên tục chảy vào bể cảng, cả khi triều dâng lẫn triều rút. Nếu nồng độ bùn trong bể lớn hơn ở ngoài bể do nạo vét (có khuấy động đáy biển), thì sẽ có lượng tịnh dòng chảy bùn hướng ra do các dòng chảy triều và trao đổi. Việc nạo vét có thể được mô phỏng bởi một nồng độ bùn ban đầu cao ở trong bể cảng. Thể tích bùn cát hằng năm (m3/năm) chuyển vào qua cửa vào và cửa nối với bể phụ được cho bởi: Vs,in = [Ntides/ρd] [floodΣ{Qtide1 co Δt} + ebbΣ{Qtide2 c2 Δt} + + tideΣ{(Qhor + Qvert) (co - c) Δt} + tideΣ{Tmud,lat Lm Δt} + + tideΣ{p Qriver Criver Δt}] (3) Thể tích bùn cát hằng năm (m2/năm) chuyển ra qua hai cửa này được cho bởi: Vs,out = [Ntides/ρd] [ebbΣ {Qtide1 c Δt} + floodΣ {Qtide2 c Δt} + tideΣ{p Qriver c Δt}] (4) Thể tích sa bồi tịnh hằng năm trong bể cảng có thể được xác định bằng hiệu số các lưu lượng vận chuyển qua cửa vào (1) và cửa phụ (2) của bể cảng, từ đó ta có: Vs,bồi = Vs,in - Vs,out (5) Bề dày lớp bồi lắng hằng năm được cho bởi: Δzb = Vs,bồi /A (6) Phương pháp ‘hộp’ này đặc biệt phù hợp với các bể cảng nhỏ có mặt bằng gần như chữ nhật, khi đó lớp sa bồi sẽ gần như đồng đều trong bể cảng. Các phương trình (1 đến 6) có thể được giải bằng cách số trị nếu cho trước điều kiện biên, và đã được Van Rijn lập bảng tính mẫu bằng Excel. Đạo hàm thời gian của nồng độ được xác định tại t, lấy giá trị c từ thời gian trước t - Δt. Kết quả tính độ nhạy cho thấy tốc độ bồi lắng tăng gần như là tuyến tính với nồng độ bùn (co) bên ngoài bể cảng. Với luồng sông Hậu, các số liệu đầu vào được chọn như sau: - Diện tích của bể cảng được tính xấp xỉ: A = 2,4 km × 2 km = 4,8 km2. - Biên độ triều η0 = ½ độ lớn triều = 1,75 m. - Chu kì triều T = 44700 s (bán nhật triều, 12h25’). Để tiện việc trình bày trong bảng tính, lấy giá trị 12 h tức là T = 43200 s. - Vận tốc cực đại của dòng triều ngoài bể cảng u0m = 0,4 m/s (lấy giá trị của hải lưu, mặc dù đây chỉ là lựa chọn sơ bộ). - Nồng độ trầm tích cực đại ngoài bể cảng c0 = 0,35 g/l. - Biên độ nồng độ trầm tích ngoài bể cảng c02 = 0,1 g/l. - Chênh lệch cực đại nồng độ trong và ngoài bể cảng (Δρm): được lựa chọn trong quá trình tính toán để thu được kết quả hợp lý. - Chọn các giá trị hệ số trao đổi động lượng như mặc định. - Dung trọng nước ρ = 1.025 kg/m3; - Tốc độ chìm đều ws được tính theo đường kính hạt 0,06 mm; được ws = 0,0031 m/s; - Mật độ khối của trầm tích chìm lắng (có kể đến độ xốp rỗng), ρdry, có thể được ước tính từ đường kính hạt (Wu và Wang, 2006): ρdry =1,0 1,5 tấn/m3; chọn ρdry = 1.300 kg/m3. CHÚC MỪNG NĂM MỚI 2019 44 Tạp chí khoa học Công nghệ Hàng hải Số 57 - 01/2019 3. Kết quả tính toán Với sự trợ giúp của phần mềm bảng tính Excel ‘Harbour sedimentation’ được Van Rijn thiết lập dựa trên cơ sở lý thuyết đã trình bày ở trên. Nhóm nghiên cứu thu được các đặc trưng bồi lắng sau 1 năm như sau: - Khối lượng bồi lắng m = 1,58 × 109 kg/năm; - Thể tích bồi lắng V = 1,21 × 106 m3/năm; - Bề dày bồi lắng z = 0,252 m/năm. Kết quả nghiên cứu của phương pháp này được đối chiếu, so sánh với kết quả thực đo là z = 0,244 m/năm. Kết quả này được Tư vấn TEDIPORT tính toán, Tư vấn Portcoast thẩm tra xác nhận.Kết quả so sánh đã cho thấy sự phù hợp của cơ sở lý thuyết với số liệu thực đo tại khu vực cảng nghiên cứu. 4. Kết luận Một số kết luận quan trọng sau đây được rút ra từ các kết quả nghiên cứu nêu trên: - Đã giải quyết được vấn đề về khối lượng bồi lắng trở lại, đánh giá định lượng được khối lượng bồi lắng theo nguyên nhân khách quan (sa bồi theo chế độ dòng chảy tự nhiên) làm cơ sở để các nhà quản lý đánh giá trách nhiệm về tài chính đối với khối lượng bồi lắng trong quá trình dừng thi công; - Giải pháp đã giúp giải quyết nhanh và tiết kiệm được 3 đến 4 tỷ đồng khảo sát số liệu điều kiện tự nhiên (địa hình, thủy hải văn, bùn cát lơ lửng, nếu tính toán bằng phương pháp mô phỏng trên mô hình toán. TÀI LIỆU THAM KHẢO [1] Báo cáo nghiên cứu phân tích, đánh giá nguyên nhân và xác định khối lượng sa bồi sau dừng giãn thuộc Dự án “Đầu tư xây dựng công trình luồng cho tàu biển trọng tải lớn vào sông Hậu”, Công ty Cổ phần Tư vấn Đầu tư Phát triển hạ tầng Quảng Ninh. [2] Thanh V.Q., Reyns J., Wackerman C., Eidam E.F., Roelvink D. Modelling suspended sediment dynamics on the subaqueous delta of the Mekong River, Continental Shelf Research, 2017. [3] Van Rijn L.C. Principles of sediment transport in rivers, estuaries and coastal seas. Aqua Publication, 1993. [4] Van Rijn L.C. Principles of Sedimentation and Erosion Engineering in Rivers, Estuaries and Coastal Seas. Aqua Publications, 2005. [5] Van Rijn L.C. Harbour siltation and control measures. Technical Note, 2016. Ngày nhận bài: 19/11/2018 Ngày nhận bản sửa: 07/01/2019 Ngày duyệt đăng: 11/01/2019
File đính kèm:
- nghien_cuu_du_bao_sa_boi_luong_song_hau_theo_phuong_phap_cua.pdf