一种基于虚拟单元法的有自由面渗流问题的分析方法

本发明涉及岩土工程、水利工程渗流计算分析，具体地指一种基于虚拟单元法的有自由面渗流问题的分析方法。

背景技术：

1、自由面渗流问题是岩土工程、水利工程等领域常遇到的一个重要研究课题，其分析难点在于自由面位置事先是未知的，一般均需要迭代求解。求解该问题的有限元法可分为调整网格法和固定网格法两类。调整网格法将自由面作为分析区域的可移动边界，迭代过程中不断调整网格形状，修改自由面位置直至自由面位置稳定。调整网格法存在计算量大、当假定自由面位置与实际自由面位置相差较大时易造成网格畸变或交错重叠等重大缺陷，因此有被固定网格法全面取代的趋势。虽然如此，由于调整网格法简单直观，仍有不少学者针对其进行改进。

2、国内外学者发展了多种格式的固定网格法，郑宏等将其分为直觉化方法和变分不等式法两类。直觉化方法有调整流量、调整单元渗透矩阵、局部修改网格三种路径。调整流量类似于应力分析中的初应力法，执行时必须计算或考虑自由面。调整单元渗透矩阵类似于切线刚度法，使用heaviside函数将仅可适用于湿区的darcy定律延拓至全域，对于自由面经过的单元，常引入人为参数对heaviside函数进行微调为避免求解振荡。局部修改网格，如虚单元法、子单元法和丢单元法等的做法，形式上保持固定网格不变，在迭代过程中将自由面附近的网格进行修改，自由面以下为分析问题的有效区域。但是，三种路径的直觉化方法在求解过程中均无法绕开确定自由面的问题。

3、固定网格的变分不等式法，一般在固定区域内构造一个新问题，使得自由面及其上的条件不被显式包含在问题中，但该问题一旦被解出，便可使用某些函数条件确定出自由面，如截止负压法和signorini变分不等式提法等。值得注意的是，郑宏等将潜在出渗面上的边界条件提为signorini型条件，从理论上消除了出渗点的奇异性。变分不等式法迭代求解过程中无须关注自由面位置，数学基础更加严密，但理论相对复杂。

4、此外，li等人开展了用无网格法求解的尝试，但最大的障碍来自于本质边界条件的准确施加。zheng等通过无网格与流形单元法结合，破解此瓶颈后成功求解了有自由面渗流问题，对于无网格类方法在此问题上的应用具有重要意义。

5、从工程应用角度来看，传统的调整网格法虽仍有广泛应用但存在诸多缺陷；固定网格法的变分不等式法由于理论复杂尚未被工程师广泛应用；无网格法求解的研究尚待丰富深入。直觉化的固定网格法三种路径中，局部修改网格无需引入额外物理概念或数学函数，思路更为直观，但局部修改网格时由于常规有限元对单元形状的限制，往往需要对自由面切割的单元进行慎重处理。

6、适用于一般多边形网格的虚拟单元法(virtual element method,vem)，其通过投影构造单元方程中的矩阵和向量，与多边形有限元相比计算方便快捷,其最大优势在于其对单元形状的包容性，执行局部修改网格的固定网格法时被自由面切割后形成的多边形单元，vem无须进行额外处理可同样执行计算分析。因此，若能将vem与局部网格修改法结合求解渗流问题对实际工程意义重大。

技术实现思路

1、本发明的目的在于克服上述不足，提供一种基于虚拟单元法的有自由面渗流问题的分析方法，克服传统局部修改网格的局限性，在迭代过程中实现自由面上节点的自动加密，效率较高且能保证很好的数值精度。

2、本发明为解决上述技术问题，所采用的技术方案是：一种基于虚拟单元法的有自由面渗流问题的分析方法，它包括以下步骤：

3、步骤一：采集可观测数据，以其视为基本信息，建立有自由面的土石坝渗流模型；

4、步骤二：对渗流区域进行区域离散化，生成固定网格；

5、步骤三：计算单元渗透矩阵并集合为整体渗透矩阵，同时计算整体流量向量；

6、步骤四：施加边界条件、求解渗流自由面；

7、步骤五：定义误差阈值，若计算结果与实际结果的差值小于阈值，则认为计算结果即为真实自由面位置，执行步骤七；若误差检验未通过，则执行步骤六；

8、步骤六：以生成的自由面为基准，局部修改网格，更新整体渗透矩阵和流量向量，并再次执行步骤四；

9、步骤七：输出渗透域内参数。

10、优选地，所述步骤一具体包括：

11、s1.1：对拟研究对象进行观测，确定外形尺寸；

12、s1.2：确定上下游水位、坝体宽度并定义为h1、h2、h3；

13、s1.3：确定岩层分布以及渗透系数，若所研究问题为均质问题，取相同渗透系数；若问题为非均质问题，取不同的渗透系数；

14、s1.4：以上述信息为基础，建立有自由面的土石坝渗流模型。

15、优选地，所述步骤二具体包括：

16、s2.1：对拟分析对象进行固定网格划分，并存储固定网格信息、节点信息；

17、s2.2：根据经验假设一条自由渗流面的位置，将其设定为初始自由面。

18、优选地，所述步骤三具体包括：

19、s3.1：计算单元渗透矩阵ke，单元渗透矩阵的具体表达形式为：

20、

21、s3.2：矩阵元素kij的计算格式为：

22、式中，d为渗透张量，分别表示i、j方向上的基函数，且其自身无显示表达，为哈密顿算子；

23、s3.3：按照一一对应的原则将各个单元渗透矩阵填充到整体渗透矩阵k内；

24、s3.4：计算整体流量向量fe，具体为：

25、式中，为规范化的单元基函数向量：

26、qn为流量；所以的投影具体表达为：

27、

28、此处的li和li+1为与节点i关联的两条邻边的长度，ni和ni+1为两条邻边的单位外法线向量，ni＝(nix,niy)t，a为多变形面积，xi、yi为节点的横纵坐标值；n等于多边形的边数；

29、s3.5：建立基函数在多项式函数空间p1(ωe)的投影后，将基函数分解为：

30、

31、由于且

32、式中，p1(ωe)表示多项式函数空间，v1(ωe)表示测试函数空间；

33、由于投影残差则有：

34、

35、根据构造投影时所采用的正交条件，可知上式子右侧的第二项和第三项为零，于是进一步简化为：

36、

37、s3.6：确定整体渗透矩阵形式，上式右侧的第一项是可准确计算的，虚拟单元法称其为kij的一致项；右侧的第二项则称为kij的稳定项，使用投影残差在单元节点处的值构建稳定项：

38、

39、上式中的λ为用于保证稳定项数量级与一致项具有可比性的正近似乘子，直接取1，为合理的考虑渗透张量d的影响，鉴于元素kij的计算由一致项和稳定项组成，将单元渗透矩阵也分为相应的两部分：

40、ke＝kc+λks          (10)

41、s3.7：推导乘子表达形式，此处kc的元素为kij一致项，ks的元素为不使用乘子λ情况下按式(9)计算的kij稳定项；式(10)中的乘子λ计算思路为：kc代表了基函数的投影对单元渗透矩阵的贡献，λks则代表了基函数的投影残差对单元渗透矩阵的贡献；既然式(9)中用投影残差在单元节点处值估计其所做的贡献，那么投影的贡献可采用同样方式估计，为：

42、

43、s3.8：确定乘子，引入参考矩阵krefer，其元素为不使用乘子λ情况下按式(11)计算的kij的参考项，那么乘子λ可按下面的方式确定：trace表示求取矩阵的迹

44、

45、s3.9：存储固定网格的整体渗透矩阵k和流量向量信息fe。

46、优选地，所述步骤四具体包括：

47、s4.1：根据此刻网格分布施加边界条件；

48、s4.2：求解方程组：kehe＝fe(13)，确定自由面的位置；

49、s4.3：采用二次曲线相交法确定出渗点位置。

50、优选地，所述步骤五具体包括：根据实际工程情况以及预期结果精度定义误差阈值ε，对所求自由面位置进行校核，自由面应满足h-y<ε；

51、式中h为自由面上任一点的水头值，y为该点纵坐标。

52、优选地，所述步骤六具体包括：

53、s6.1：读取固定网格信息，使用新自由面对其进行切割；

54、s6.2：将单元分为三类：删除单元(全部节点在自由面之上)、原生单元(全部节点在自由面之下)、贯穿单元(部分节点在自由面之下，部分节点在自由面之上)；

55、s6.3：定义自由面与贯穿单元的交点为新生节点，新生节点与贯穿单元位于自由面以下的节点构成新生单元；

56、s6.4：淘汰删除单元和贯穿单元，添加新生单元，保留原生单元；

57、s6.5：仿照单元形式对节点进行同样处理，最终形成新的渗流区域网格的单元节点信息；

58、s6.6：根据新生成的渗流域网格的单元节点信息，更新整体渗透矩阵和流量向量。

59、优选地，所述步骤七具体包括：

60、s7.1：输出渗流自由面的位置；

61、s7.2：输出各节点的水头值；

62、s7.3：计算出渗点的位置以及迭代误差。

63、本发明的有益效果：本发明克服了传统局部修改网格的局限性，在迭代过程中实现自由面上节点的自动加密，效率较高且能保证很好的数值精度；本发明是基于虚拟单元法原理求解自由面的位置的数值计算方法，避免了网格畸变对计算精度的影响；同时迭代过程中采用局部网格修改的原则，降低了存储要求；大大的提高了求解效率；其主要应用于岩土工程、水利工程等领域的渗流分析计算中；在设计土坝时，需要通过渗流计算来确定渗漏损失和合理的防渗排渗措施，而坝型和坝的断面尺寸也经常需要借助渗流计算来比较选定，故而采用本方法准确计算自由面的位置对实际工程意义重大。