1.一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，具体包括如下步骤：S1，采集目标层位不同岩性岩石，进行分类标注；按岩石类型将岩石试块制备成长方体；分析裂隙面力学成因，采用轴向劈裂、剪切试验方法制备粗糙裂隙，制备粗糙裂隙上、下壁面；

S2，采用三维激光扫描技术对粗糙裂隙上、下壁面进行扫描，获取粗糙裂隙上、下壁面形貌特征数据集；

S3，对粗糙裂隙上、下壁面形貌特征数据集进行预处理，建立裂隙上、下壁面高度数据集Ztop与Zbot；

S4，根据步骤S3获取的裂隙上、下壁面高度数据集，进行裂隙开度分布特征分析；

S5，根据步骤S3和步骤S4获取的数据集，进行粗糙裂隙上、下壁面间匹配程度分析；

S6，根据S3获取的高度数据集合，进行裂隙粗糙度分析；

S7，基于步骤S4～S6分析并建立研究区不同岩石类型的粗糙裂隙集合Rockname，用于确定生成的裂隙模型是否满足在裂隙开度、粗糙度特征上的要求；

S8，基于步骤S4～S6分析，生成与实际粗糙裂隙特征、属性相仿的粗糙裂隙数据集；

S9，基于matlab with comsol平台，建立粗糙裂隙数值模型，开展粗糙裂隙非线性渗流求解；

S10，循环步骤S8～S9，获取研究区域不同岩石类型在不同粗糙裂隙内的渗流特征，建立研究区针对不同岩石类型及裂隙粗糙程度的裂隙岩体非线性渗流指标。

2.根据权利要求1所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S3具体为：利用Surfer软件数据插值功能，通过设置x、y方向固定的节点数量M、N，将裂隙上、下壁面各划分为M×N的正方形网格，建立裂隙上、下壁面高度数据集Ztop与Zbot：式中， 表示裂隙上壁面的网格节点(M，N)位置处的高度值， 表示裂隙下壁面的网格节点(M，N)位置处的高度值。

3.根据权利要求2所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S4具体包括如下步骤：S4.1，基于Ztop与Zbot，将上、下裂隙面位置坐标z值对应相减，得到自然状态下裂隙开度数据集Aperture：式中，ApM，N表示裂隙的网格节点(M,N)位置处裂隙上壁面与下壁面间的裂隙开度；

S4.2，对裂隙开度数据集Aperture中的裂隙开度数据进行分布特征分析，绘制裂隙开度频率直方图，分析裂隙开度分布符合的概率密度函数，选取概率密度函数进行描述，并确定分布均值Mean和标准差SD。

4.根据权利要求3所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S5具体包括如下步骤：S5.1，将裂隙上、下壁面高度数据集和裂隙开度数据集导入Matlab软件，采用Matlab中的快速傅里叶变换函数，计算裂隙上、下壁面高度数据集以及裂隙开度数据集的功率谱密度PSD；快速傅里叶变换函数所需计算参数包括：傅里叶分量中的波数k、采样频率v、波长λ，其中k＝1/λ，v＝2π/λ，波长λ等于裂隙剖面的长度；

S5.2，基于求得的各数据集PSD曲线，利用公式(4)计算相关系数PSDR，采用双对数坐标轴形式绘制PSDR与波长λ的关系曲线，进行裂隙匹配程度分析；

式中，PSD裂隙开度为裂隙开度数据集功率谱密度曲线，PSD上壁面为裂隙上壁面功率谱密度曲线，PSD下壁面为裂隙下壁面功率谱密度曲线；

S5.3，根据相关系数PSDR与波长λ的关系曲线，确定裂隙上、下壁面间的最小匹配系数MinMF、最大匹配系数MaxMF、失配长度ML以及过渡长度TL。

5.根据权利要求4所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S6具体包括如下步骤：S6.1，基于步骤S3中的裂隙上、下壁面高度数据集合，采用Z2s描述法进行裂隙粗糙度评价，Z2s表示为：式中，Δm和Δn为分别为x、y方向阶段间距；zm，n代表高度值；

当裂隙上、下壁面不同时，用上、下裂隙面的 平均值表示，公式为：式中 表示裂隙平均Z2s粗糙度， 为裂隙上、下壁面Z2s粗糙度；

S6.2，采用迂曲度Rs表征法，从高度分布特征角度，对裂隙面的粗糙程度进行表征，Rs为裂隙面实际面积Sreal与裂隙面投影面积S之比，裂隙面投影面积S’可以通过测量获得，裂隙面实际面积Sreal用式(7)表示：当裂隙上下表面不同时，利用上、下裂隙壁面的平均值 表示裂隙面的粗糙程度，公式为：式中Rstop与Rsbot分别为裂隙面上、下表面的Rs值，Sreal‑top表示裂隙上壁面实际面积，Sreal‑bot表示裂隙下壁面实际面积，Stop表示裂隙上壁面投影面积，Sbot表示裂隙下壁面投影面积。

6.根据权利要求5所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S7具体为：式中，Ap‑Meanmax为不同粗糙裂隙模型裂隙开度分布均值最大值,Ap‑Meanmin为不同粗糙裂隙模型裂隙开度分布均值最小值，Ap‑SDmax为不同粗糙裂隙模型裂隙开度分布标准差最大值，Ap‑SDmin为不同粗糙裂隙模型裂隙开度分布标准差最小值，Frac‑Z2smax为不同粗糙裂隙模型裂隙粗糙度指标Z2s最大值，Frac‑Z2smin为不同粗糙裂隙模型裂隙粗糙度指标Z2s最小值，Frac‑Rsmax为不同粗糙裂隙模型裂隙粗糙度指标Rs最大值，Frac‑Rsmin为不同粗糙裂隙模型裂隙粗糙度指标Rs最小值。

7.根据权利要求6所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S8具体为：S8.1，采用基于空间频率分布的反傅里叶变换建模法进行粗糙裂隙的建模，其中反傅里叶变换建模法需要的参数包括：序列集合F(m，n)与节点数量M、N，其中序列集合F(m,n)具体表示为：F(m，n)＝{Rad cos(Phase)，Rad sin(Phase)}；

式中，Rad、Phase分别为序列集合F(m，n)中具体参数；

S8.2，重复步骤S8.1，明确反傅里叶变换函数中的两个计算参数Rad与Phase；利用matlab中的反傅里叶函数function IFFT2D(F(m，n)，M，N)，自动生成粗糙裂隙上、下壁面高度数据集；

S8.3，将生成的粗糙裂隙上、下壁面高度数据集重复步骤S2～S6，计算裂隙开度分布特征与粗糙度特征，与步骤S7中构建的粗糙裂隙集合Rockname进行对比，当裂隙参数满足标准裂隙开度分布、粗糙度特征集合标准时，判定为有效粗糙裂隙模型，并对粗糙裂隙集合数据文件进行保存；当不满足标准集合要求时，判定为无效粗糙裂隙模型，清除集合数据，并重复步骤S8。

8.根据权利要求7所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S8.1具体包括：S8.1.1，输入粗糙裂隙上、下壁面间匹配系数、失配长度和过渡长度参数；

S8.1.2，令i＝1～M，j＝1～N，循环迭代执行步骤S8.1.2～S8.1.4，求解粗糙裂隙区域范围内的反傅里叶函数参数集合F(m,n)；

S8.1.3，循环迭代计算中，参数Rad为衰减函数，即粗糙裂隙壁面功率谱，按照函数A(m,m)分布逐渐衰减，表示为：式中，为频谱指数，表示衰减的速度， D为粗糙裂隙分形维数；

S8.1.4，裂隙壁面参数phase根据粗糙裂隙上、下壁面间匹配程度，通过比较节点位置

2 2 1/2

长度A(m,m)＝(m+n) 与失配长度ML大小进行确定：当节点位置长度A(m,m)＜ML/2时，粗糙裂隙壁面参数phase均等于phase＝2πRan1；

当节点位置长度A(m，m)≥ML/2时，粗糙裂隙上、下壁面参数phase上和phase下采用随机数Ran1、Ran2和Ran3进行计算，分别为：phase上＝2πRan1；

phase下＝2πRan3；

Ran1与Ran3间关系具体判断如下：

Ran3＝γRan1+(1‑γ)Ran2；

式中，kc为失配波长的波数，kc＝2π/λc，λc为失配波长，k为傅里叶分量中的波数。

9.根据权利要求8所述的一种基于真实粗糙裂隙的模型构建与非达西渗流模拟方法，其特征在于，步骤S9具体包括：S9.1，基于matlab with comsol平台，加载步骤S8生成的粗糙裂隙集合数据文件，利用参数化扫描功能，实现上、下粗糙裂隙壁面构建；随后，建立长方体，利用布泽尔运算法，将长方体与上、下粗糙裂隙壁面求并集，建立粗糙裂隙数值模型；

S9.2，利用matlab with comsol平台，调用软件平台的控制方程，其中调用纳维‑斯托克斯方程作为流体流动的控制方程，在流体密度恒定且为不可压缩流体时，N‑S方程表示为：式中，u为流体速度，p为流体压力，ρ为流体密度，μ为流体动力黏度，F为体积力，即作用于流体上的外力，为拉普拉斯算子；

调用表示质量守恒的连续性方程，与公式(12)形成耦合方程组，连续性方程表示为：同时，采用湍流Standradk‑ε模型，来模拟粗糙裂隙内流体的非线性渗流过程，其中湍流动能kT和湍流耗散ε的约束性方程分别为：其中，Pk为湍动能生成项，表示为：

2

式中，湍流黏度μT＝ρCμk/ε，湍流模型参数采取默认参数：Cε1＝1.44；Cε2＝1.44；Cμ＝

0.09；

S9.3，为粗糙裂隙数值设定边界条件，以x方向作为流体流动方向，在模型入口与出口边界分别设置为压力边界条件，并将y方向壁面作为不透水边界；

S9.4，利用comsol模拟软件进行模型网格划分，调用求解器对流场和压力场进行求解；

求解结束后，计算粗糙裂隙模型入口处渗流速度，通过在入口处对流速积分，获取不同裂隙开度分布情况下的稳态体积流速Q，具体为：Q＝Au         (17)；

式中：A＝∫dh，A为垂直于流体流动方向的横截面积，u为流体通过横截面积的速度；

S9.5，采用立方定律，利用模型入口压力P与体积流速Q，计算不同体积流速下的等效裂隙开度，作为研究区非线性渗流参考指标，具体计算公式为：式中，he为等效裂隙开度，W为裂隙模型宽度。