Nguyên cứu đặc trưng biến dạng bất đẳng hướng của khối đá nứt nẻ tại Km119+000 quốc lộ 3B, tỉnh Bắc Kạn theo phương pháp EFC

ĐỊA KỸ THUẬT SỐ 4 - 2020 43 NGUYÊN CỨU ĐẶC TRƯNG BIẾN DẠNG BẤT ĐẲNG HƯỚNG CỦA KHỐI ĐÁ NỨT NẺ TẠI KM119+000 QUỐC LỘ 3B, TỈNH BẮC KẠN THEO PHƯƠNG PHÁP EFC ĐẶNG HÔNG LAM*, PHÍ HỒNG THỊNH* A study of anisotropic mechanical properties of fractured rock mass at Km119+000 in 3B highway, Backan province by EFC method Abstract: The stability of rock slope is almost studied based on the assumption of homogenous and isotropic material properties of fractured rock mass by equivalent equilibrium method. However, the complicated distribution of fractures give significantly the anisotropy of mechanical properties of fractured rock as results in the literature. The article presents the anisotropic mechanical properties of fractured rock mass at Km119+000 in 3B highway, Backan province by EFC method in which the mechanical properties such as Young modulus and Poisson ratio are obtained from mechanical properties of intact rock and fractures. By the site survey and re-generation of fractures distribution at the Km119+000 in 3B highway, the anisotropy of fractured rock mass is highlighted. 1. GIỚI THIỆU CHUNG* Nghiên cứu trượt lở ở Việt Nam được tiến hành từ đầu những năm 1985 và tập trung chủ yếu ở Viện Địa chất-Viện Hàn lâm Khoa học và Công nghệ Việt Nam, Viện Khoa học Địa chất và Khoáng sản, Viện Khoa học và Công nghệ Giao thông Vận tải. Phần lớn các nghiên cứu đều được tiến hành dựa vào xử lý số liệu ảnh vệ tinh, địa hình, địa mạo... để xây dựng bản đồ phân vùng và dự báo nguy cơ trượt lở. Các nghiên cứu đã kết luận trong những năm gần đây tại các tỉnh miền núi phía Bắc, dạng tai biến này diễn ra ngày càng nguy hiểm hơn với 12/16 tỉnh nằm trong vùng có nguy cơ trượt – lở cao. Theo hướng này, phải kể đến các công trình: Nghiên cứu đánh giá tai biến trượt lở khu vực các tỉnh miền núi phía Bắc và các giải pháp phòng tránh (Lê Thị Nghinh và nnk, 2003). Nghiên cứu đánh giá trượt lở, lũ quét- lũ bùn đá tại vùng trọng điểm tỉnh Lào Cai (H. Bát Xát, H. Sapa và TP. Lào Cai - tỉnh Lào Cai) và kiến nghị các giải pháp phòng tránh giảm nhẹ thiệt hại (Nguyễn * Khoa Công trình-Đại học Giao thông Vận tải 03 Cầu Giấy-Láng Thượng-Đống Đa-Hà Nội Email: [email protected] Trọng Yên, 2006). Có thể phân biệt hai nhóm phương pháp mô hình nghiên cứu trượt lở: - Nhóm phương pháp vật lý dựa trên các phương trình toán lý mô phỏng bản chất vật lý của quá trình trượt. Trong đó phương pháp phân tích cân bằng giới hạn được sử dụng phổ biến. Một số tác giả khác của Việt Nam đã phát triển Lý thuyết khối của Goodman và Shi (1985) để phân tích xác định khối nguy cơ trượt lở trên bề mặt mái dốc đá dựa vào số đo thế nằm khe nứt và thế nằm của bề mặt mái dốc. - Nhóm phương pháp thống kê dựa trên quan hệ thống kê giữa các điểm trượt lở và các yếu tố được cho là nguyên nhân gây nên trượt lở. Phí Trường Thành (2018) đã nhận dạng các kiểu trượt và mô hình hóa 3 chiều (3D) các khối nguy cơ trượt trên một số mái dốc đá điển hình trong mối liên quan với hoạt động tân kiến tạo: Áp dụng cho tuyến Quốc lộ 4B, đoạn Tiên Yên - Lạng Sơn. Trần Mạnh Liểu (2013) đã dự báo nguy cơ và cường độ phát triển trượt lỡ khu vực Thị xã Bắc Kạn. Tuy nhiên, các nghiên cứu ở trên vẫn chưa làm nổi bật sự ảnh hưởng của vết nứt đến tính ổn định của bờ dốc đá nứt nẻ. Cụ thể là do sự phân bố của vết nứt làm cho tính chất biến dạng ĐỊA KỸ THUẬT SỐ 4 - 2020 44 của khối đá theo các phương là khác nhau. Một số nghiên cứu của Dang (2018) và Dang và nnk (2019) cho khối đá ở trên thế giới như Sellafield chỉ ra rằng ứng xử của khối đá nứt nẻ là không đẳng hướng. Do đó mục đích của bài báo là tìm hiểu sự ảnh hưởng khe nứt đến tính chất biến dạng của khối đá ở quốc lộ 3B thông qua kết quả khảo sát hiện trường về tình trạng nứt nẻ trên Quốc lộ 3B (QL3B) tại Km119. Cấu trúc của bài báo như sau: tiếp theo phần giới thiệu là lý thuyết xác định các thông số biến dạng của đá nứt nẻ. Tiếp theo nhóm tác giả tiến hành khảo sát hiện trường QL3B tại Km119 để xác định thế nằm và các vết nứt nẻ của đá. Trên cơ sở kết quả hiện trường, nhóm tác giả tiến hành xây dựng lại hệ thống vết nứt và các thông số đặc trưng của vết nứt. Từ đó, nhóm tác dụng áp dụng phương pháp EFC (Dang, 2018; Dang và nnk, 2019) để xác định đặc trưng biến dạng của khối đá. Kết quả chỉ ra tính chất bất đẳng hướng rất lớn của khối đá thực tế theo các phương ở tuyến QL3B. Cuối cùng là một số kết luận chung của bài báo. 2. PHƢƠNG PHÁP EFC TRONG VIỆC XÁC ĐỊNH TÍNH CHẤT CƠ HỌC TƢƠNG ĐƢƠNG CỦA ĐÁ NỨT NẺ Mô đun đàn hồi tương đương của đá nứt nẻ Ý tưởng để xác định mô đun đàn hồi tương đương của đá nứt nẻ ở đây được đề xuất trong tài liệu (Dang và nnk, 2019) trong đó phương pháp nhúng vết nứt vào môi trường liên tục (EFC) được sử dụng để mô hình hóa mạng lưới vết nứt dưới ứng xử cơ học. Điểm cơ bản của phương pháp này, so với các phương pháp đã đề cập trước đây, là việc sử dụng khái niệm phần tử tương đương. Theo đó, mỗi vết nứt được mô hình hóa bởi môi trường gồm một nhóm các ô lưới (phần tử vết nứt) giao nhau (xem Hình 1) (Dang và nnk, 2019). Để thực hiện chi tiết, chúng tôi phân biệt ô lưới thành 2 loại sau: (1) ô lưới có chứa vết nứt (fracture cell) và (2) ô lưới không chứa vết nứt (matrix cell) (Dang và nnk, 2019). Cũng lưu ý rằng trong phương pháp EFC, chúng tôi phân chia loại ô lưới có chứa vết nứt thành 2 loại nhỏ: loại I đại diện cho tế bào được giao nhau bởi chỉ một vết nứt trong khi đó loại II là tập hợp tất cả ô lưới có chứa nhiều vết nứt (nhiều hơn một vết nứt) như được chỉ ra trong Hình 1. Hình 1: Sơ họa khái niệm ô lưới có chứa vết nứt và ô lưới không chứa vết nứt trong phương pháp EFC (a), loại I biểu thị các ô lưới chỉ chứa một vết nứt (b) và loại II biểu thị ô lưới chứa nhiều hơn một vết nứt (c). Một trong những nghiên cứu đầu tiên trong phương pháp này được thực hiện bởi Figueiredo và nnk (2015) trong đó mô đun đàn hồi của đá nứt nẻ là đẳng hướng được tính từ mô đun đàn hồi của đá gốc (đá không nứt nẻ), độ cứng của khe nứt và khoảng cách của khe nứt (công thức 1). Hơn nữa, hệ số Poisson của đá nứt nẻ được đề xuất bằng hệ số Poisson của đá gốc 1 1 1 ; fc iso m n fc iso m E E k h      (1) Trong đó Em, vm, kn và h tương ứng là mô đun đàn hồi, hệ số Poisson của đá gốc, độ cứng pháp tuyến và kích thước của ô lưới. Nghiên ĐỊA KỸ THUẬT SỐ 4 - 2020 45 cứu của Dang và nnk (2019) đã đề xuất công thức xác định mô đun đàn hồi tương đương của bề mặt đá có nứt nẻ theo công thức (2)(3): Đối vớiphần tử loại I có một vết nứt, mô đun đàn hồi theo các phương của phần tử được tính theo công thức (2): 1 2 2 41 cos .sin cos ; ; . . fc fc iso iso m m t n E E k d k d                (2) Trong đó:  là góc nghiêng của vết nứt minh họa theo hình 1b. kn và kt lần lượt là độ cứng vuông góc (chống nén) và tiếp tuyến (chống cắt) của vết nứt. d là khoảng cách các vết nứt. Tương tự, đối với phần tử loại II có nhiều hơn một vết nứt, mô đun đàn hồi theo các phương cũng được xác định theo công thức (3): 1 2 2 4cos .sin cos1 ( ) ; ; . . fc fci i i iso iso mi i im t i n i E E k d k d                 (3) Trong đó chỉ số i là chỉ số ứng với vết nứt thứ i. Ý nghĩa của các đại lượng tương tự như ở công thức (2) và được minh họa ở hình 1c. 3. ÁP DỤNG PHƢƠNG PHÁP EFC ĐỂ XÁC ĐỊNH MÔ ĐUN ĐÀN HỒI TƢƠNG ĐƢƠNG CỦA ĐÁ NỨT NẺ Ở QUỐC LỘ 3B TỈNH BẮC KẠN Phương pháp xác định tính mô đun đàn hồi tương đương Mô đun đàn hồi tương đương của các khối đá nứt nẻ được tính toán trong bài báo này dựa trên giả thuyết biến dạng mặt phẳng 2D. Cách xác định ma trận biến dạng đã được đề xuất trong (Min và Jing, 2003) (Dang và nnk, 2019) được áp dụng ở đây. Với giả thuyết biến dạng phẳng, mô đun đàn hồi của vật liệu tổng quát được viết dưới dạng ma trận như công thức (4) sau: 11 12 13 14 21 22 23 24 31 32 33 34 41 42 43 44 zz zz xx xx yy yy xy xy S S S S S S S S S S S S S S S S                                          (4) Lưu ý ở đây là tọa độ Oxyz thể hiện tọa độ tổng thể, tọa độ phần tử được liên hệ với tọa độ tổng thể như sau z = 1, x = 2, y = 3 Giả sử rằng mô đun đàn hồi và hệ số Poisson theo hướng z bằng với các giá trị của đá gốc ( zx zy m    , z mE E ) và tuân theo 12 21 13 31 / ,m mS S S S E     11 1/ mS E . Nghiên cứu của (Min và Jing, 2003) cũng chỉ ra rằng các giá trị S14,S24,S41,S42 rất nhỏ so với các thành phần khác trong khi S34,S43bằng 0 do thực tế là ứng suất cắt σxy không ảnh hưởng đến biến dạng hướng z. Do đó, ma trận biến dạng trong phương trình (4) có thể được viết lại dưới dạng như công thức (5) sau: 22 23 32 33 44 1 0 0 02 0 0 0 m m m m m zz zz m xx xx m yy yy m xy xy m E E E S S E S S E S                                                         (5) (a) (b) (c) Hình 2: Ba bài toán: nén theo hướng x (a), nén theo hướng y (b), cắt thuần túy (c) để xác định năm thành phần S22,S23,S32,S33,S44 của khối đá nứt nẻ ĐỊA KỸ THUẬT SỐ 4 - 2020 46 Do đó, việc tính toán các đặc trưng tính chất đàn hồi tương đương của các khối đá bị nứt nẻchỉ còn làviệc xác định năm thành phần S22, S23, S32, S33, S44. Các thông số này được thực hiện thông qua ba bài toán minh họa trong Hình 2. Bài toán đầu tiên (Hình 2a) thể hiện cho bài toán nén đơn phương với ứng suất pháp tuyến f a x được áp dụng theo phương x. Hai thành phần S22 và S32được tính toán từ bài toán này theo công thức (6)trong đó có xét đến điều kiện ( ) 0 a zz  (do đó a a zz m xf  ): ( ) ( )2 2 11 21( ) ( ) ; a a xx yym m a a x m x m S S f E f E        ( ) ( ) ( ) ( ) 22 32( ) ( ) ; a a a am m xx m x yy m x m m a a x x f f E E S S f f          (6) Tương tự, hai thành phần S23 và S33 được tính từ bài toán thứ hai theo công thức (7) (Hình 2b): ( ) ( )2 2 21 22( ) ( ) ; b b xx yym m b b y m y m S S f E f E        (b) (b) (b) (b) 23 33(b) (b) ; m m xx m y yy m y m m y y f f E E S S f f          (7) Thành phần cuối cùng S44 được xác định từ bài toán thứ ba theo công thức (8) (Hình 2c): ( ) ( ) 44 2 c cxy xy S       (8) Lưu ý rằng ( ) ( ),a axx yy  và (b) (b),xx yy  thể hiện biến dạng trung bình dưới ứng suất theo hướng x và y tương ứng trong khi (c)xy chỉ ra biến dạng cắt trung bình theo tải trọng cắt. Từ các thành phần này của ma trận biến dạng, chúng ta suy ra mô đun đàn hồi, hệ số Poisson và mô đun chống cắt của vật liệu tương đương (đá nứt nẻ) được xác định theo công thức (9): 23 32 44 22 33 1 1 ; ; ; ; 1/x y xy y yx x xyE E S E S E G S S S       (9) Kết quả mô đun đàn hồi của khối đá nứt nẻ trên QL3B tỉnh Bắc Kạn Tiến hành đo vẽ vết nứt tại Km 119 xác định được dạng vết nứt được sơ họa theo hình 3. Tiến hành đo góc vết nứt và độ nghiêng của vết nứt theo minh họa ở hình thu được thông số vết nứt theo như bảng1. Bằng phép thống kê số lượng vết nứt cho mặt bờ dốc có kích thước 5m2, ta được thống kê số lượng vết nứt như bảng 2. Hình 3. Sơ họa mặt nứt của đá Trong đó: • Đường phương (a) là đường thẳng nằm ngang trên mặt lớp đá (chỉ phương kéo dàicủa lớp đá) • Đường hướng dốc (b‟) là đường thẳng nằm ngang, vuông góc với đường phương vàcắm theo hướng dốc của lớp đá (chỉ hướng cắm của lớp đá). Góc của đường hướng dốc so với phương bắc được ký hiệu là () • Góc dốc (α) là góc tạo bởi giữa mặt lớp đá với mặt phẳng nằm ngang Bảng 1: Thế nằm của các lớp đá trên bờ dốc đá tại Km119 Lớp Góc phương vị đường hướng dốc () (độ) Góc dốc (θ) (độ) Khoảng cách (m) 1 125 60 3 2 345 40 10 3 95 70 0,2 ĐỊA KỸ THUẬT SỐ 4 - 2020 47 Bảng 2: Sự phân bố chiều dài vết nứt tại Km 119 Chiều dài vết nứt (m) Số lượng vết nứt (cái) Góc phương vị đường hướng dốc (θ) (độ) Min. Max. 0,5 0,6 42 95 0,6 0,7 24 95 0,7 0,8 15 95 0,8 1 17 95 1 1,5 16 95 1,5 2 5 95 2 3 4 95 3 4 1 125 4 5 1 125 5 10 1 345 Khi đó thông số đặc trưng cho mật độ vết nứt được xác định như sau: P20 = 5.04 vết nứt / m 2 P21 = 4.78 m dài / m 2 Giống với sự phân bố vết nứt của đá được nghiên cứu bởi Min và Jing (2003), sự phân bố chiều dài của khe nứt được phân bố theo quy luật power với công thức (10) trong đó trong đó NF là số lượng khe nứt (với chiều dài gãy lớn hơn chiều dài L) trên một đơn vị diện tích; C, D là các hệ số đến mật độ phân bố của vết nứt . D FN C L  (1 0) Đối với đá ởquốc lộ 3B, C=1.1 và D=2.2 được lựa chọn. Khi đó sự phân bố theo lý thuyết và thực tế của chiều dài vết nứt được biểu thị ở Hình 4. Dựa trên các cuộc thảo luận này, trong công trình này, kích thước đặc trưng của khối đá 5m được chọn để tính toán mô đun đàn hồi của khối đá nứt nẻ này bằng cách sử dụng phương pháp EFC. Trong đó sự phân bố thực tế của khối đá 5m được chụp lại ở hình 5. Hình 6 trình bày một mẫu hình vuông (L = Lx = Ly = 5m) của đá nứt nẻ ở QL3B được tạo ra tương ứng số lượng vết nứt là 118 theo phương pháp đã trình bày bởi Dang và Phi (2020) 0 20 40 60 80 100 120 0 1 2 3 4 5 S ố lƣ ợ n g v ết n ứ t Chiều dài (m) Lý thuyết Thực tế Hình 4: So sánh sự phân bố chiều dài vết nứt theo lý thuyết và thực tế Hình 5: Sự phân bố vết nứt của đá lại Km119 ĐỊA KỸ THUẬT SỐ 4 - 2020 48 Hình 6: Kết quả tạo lại mẫu đá nứt nẻ có cạnh 5m Vấn đề tiếp theo trong việc xác định các tính chất biến dạng tương đương của khối nứt ở quốc lộ 3B đó là các đặc trưng biến dạng cần được xác định.Nghiên cứu của Ahmed Alnuaim và nnk (2019) đã chỉ ra giá trị cường độ kháng nén nở hông (qu) của đá vôi thường dao động trong khoảng từ 100MPa đến 200MPa đồng thời kiến nghị giá trị mô đun đàn hồi Emtrong khoảng từ 47.0 GPa đến 58.0 GPa, hệ số poisson từ 0.19 đến 0.28. Do đó, nhóm tác giả chọn giá trị Em = 50.0 GPa và hệ số poisson bằng 0.24 cho đá vôi ở quốc lộ 3B trong nghiên cứu này được thống kê ở bảng 3. Về độ cứng (pháp tuyến và tiếp tuyến) của khe nứt, đây là vấn đề khó xác định trong điều kiện Việt Nam khi chưa thể tiến hành thí nghiệm cho khối đá nguyên dạng (khối đá bao gồm cả khe nứt). Trong điều kiện của bài báo này, nhóm tác giả xác định độ cứng của khe nứt đá ở QL3B dựa theo các nghiên cứu đã có trên thế giới, cụ thể như sau. Bandis và nnk (1983) đã đề xuất công thức(11) tính độ cứng pháp tuyến của vết nứt đá theo công thức sau: (11) Trong đó: JRC là hệ số độ nhám của vết nứt. Barton and Choubey (1977) để xuất thang hệ số JRC ứng với 10 hình dạng vết nứt khác. Theo nghiên cứu của Du Shigui và nnk (2011) đã chỉ ra cho đá vôi thường có giá trị JRC trong khoảng từ 8.43 đến 9.20, do đó trong nghiên cứu này nhóm tác giả đề xuất hệ số JRC = 9.0 JCS là cường độ nén của vết nứt (MPa), theo kết quả nghiên cứu của Du Shigui và nnk (2011) đã chỉ ra rằng JCS của đá vôi ở QL3B có giá trị trung bình khoảng 70 MPa. Giá trị JCS = 70 MPa này sẽ được lựa chọn trong tính toán ở bài báo này. a0 là bề rộng ban đầu của vết nứt (mm). Kết quả khảo sát hiện trường cho thầy vết nứt đã được lấp đầy bằng vật liệu lớp nhét. Nghiên cứu của Benmadi và nnk (2018) đã chỉ ra rằng với đá vôi thì các vết nứt ở đây thước có kích thước từ 0.05 đến 0.33 mm trong đó giá trị nhiều nhất ở 0.1mm. Trong trường hợp tại quốc lộ 3B, quan sát và đo trực tiếp ở hiện trường rút ra được bề rộng vết nứt trung bình ở đây là 0.1mm và giá trị này sẽ được sử dụng trong tính toán ở bài báo này. Kết quả tính toán theo công thức (11) trên cho được giá trị kn=22.6 MPa/mm = 22.6GPa/m Về giá trị độ cứng chống cắt (ks), đã có nhiều nghiên cứu về giá trị này cho đá như Barton và Chubey (1977), Małkowski, Piotr. (2015) đã đề xuất công thức (12)tính độ cứng chống chắt (ks) như sau: (12) Trong đó: Lx là khoảng cách giữa các vết nứt. σn là ứng suất pháp. JRC là hệ số độ nhám của vết nứt. φr là góc ma sát trong của vết nứt. Phần mềm UDEC (2000) đề xuất giá trị ks trong khoảng từ 10 đến 100 MPa/m cho trường hợp khe nứt được lấp phủ bởi đất sét trong khi đó giá trị này được tăng lên rất nhiều cho đá Magma là 100GPa/m. Giá trị mặc định ở phần mềm Phase 2 (2013) là kn = 100 GPa/m và ks = ĐỊA KỸ THUẬT SỐ 4 - 2020 49 10 GPa/m. Nghiên cứu của Baohua Guo và Hangyu Dong (2018) về cường độ kháng cắt của đá đưa gia giá trị trong khoảng từ 100MPa/m đến 800 MPa/m. Với các nghiên cứu trên, nhóm tác giả đề xuất giá trị dộ cứng chống chắt của khe nứt đá ở QL3B là giá trị tương ứng với 1% độ ứng chống nén, tức là ks = 0.226 GPa/m Bảng 3: Giá trị biến dạng của đá gốc và vết nứt ở Quốc lộ 3B Đá gốc Mô đun đàn hồi Em [GPa] 50,0 Hệ số poisson νm 0,24 Vết nứt Độ cứng pháp tuyến kn [GPa/m] 22,6 Độ cứng tiếp tuyến kt [GPa/m] 0,226 Tiến hành xác định mô đun và hệ số poisson tương đương của mẫu đá có vết nứt với kích thước vuông 5m theo các công thức đã nêu ra ở trên, kết quả chuyển vị của mẫu đá ứng với các bài toán (a), (b), (c) minh họa ởhình 7 như sau. Từ đó xác định các đặc trưng biến dạng tương đương của mẫu đá theo công thức (11) cho kết quả thông kê ở bảng 4. Bảng 4: Mô đun đàn hồi tƣơng đƣơng (Ex, Ey, υyz, υxy, Gxy) xE (GPa) yE (GPa) yx xy xyG (GPa) 40,45 1,26 0,26 0,23 16,45 Hình 7: Kết quả chuyển vị của mẫu theo các bài toán nêu ra ở hình 2 4. KẾT LUẬN Các thông số biến dạng của đá nứt nẻ (mô đun đàn hồi và hệ số poisson) ở quốc lộ 3B địa phận Tỉnh Bắc Kạn đã được nghiên cứu thông qua phương pháp EFC. Kết quả nghiên cứu đã chỉ ra rằng mô đun đàn hồi của khối đá nứt nẻ này chỉ còn lại khoảng 80% và 2,5% so với đá gốc khi xét đến sự ảnh hưởng của khe nứt trong khối đá nứt nẻ theo các phương khác nhau. Điều này chứng tỏ tính bất đẳng hướng rất lớn theo 2 phương khác nhau ở khối đá thực tế đang xảy ra ở Km119 thuộc quốc lộ 3B tỉnh Bắc Kạn.Kết quả nghiên cứu này chỉ ra sự cần thiết phải đưa mô hình vật liệu bất đẳng hưởng trong xem xét sự ổn định bờ dốc đá trên tuyến quốc lộ 3B cho các nghiên cứu tiếp theo sau này. LỜI CẢM ƠN Nghiên cứu này được tài trợ bởi Trường đại học Giao thông Vận tải (ĐH GTVT) trong đề tài mã số T2020-CT-024. ĐỊA KỸ THUẬT SỐ 4 - 2020 50 TÀI LIỆU THAM KHẢO 1. Alnuaim, Ahmed & Hamid, Wagdi & Alshenawy, Abdulhafiz. (2019). Unconfined Compressive Strength and Young's Modulus of Riyadh Limestone. Electronic Journal of Geotechnical Engineering. 24. 707-717. 2. Bandis, Stavros & Lumsden, A.C. & Barton, Nick. (1983). Fundamental of joint deformation. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts. 20. 249-268. 10.1016/0148- 9062(83)90595-8. 3. Barton, Nick & Choubey, Vishnu. (1977). The shear strength of rock joints in theory and practice. Rock Mechanics Felsmechanik Mecanique des Roches. 10. 1-54. 10.1007/BF01261801. 4. Dang Hong-Lam (2018), A hydro- mechanical modeling of double porosity and double permeability. PhD thesis, University of Orleans (France), ftp://ftp.univ- orleans.fr/theses/honglam-dang_3747.pdf 5. Dang, Hong-Lam & Do, Duc Phi & Hoxha, Dashnor. (2019). Effective Elastic and Hydraulic Properties of Fractured Rock Masses with High Contrast of Permeability: Numerical Calculation by an Embedded Fracture Continuum Approach. Advances in Civil Engineering. 2019. 1-21. 10.1155/2019/7560724. 6. Dang, Hong-Lam and Thinh, Phi-Hong (2020), A methodology of re-generating a representative element volume of fractured rock mass. Transp