基于逐步剖分法求解渗流自由面的方法技术

本发明专利技术涉及一种基于逐步剖分法求解渗流自由面的方法，其属于坝工建设技术领域。它首先对渗流场做有限单元划分，依据横向能量损失率最小计算渗流溢出点；然后由溢出点向上游汇入点逐层单元推进，每推进一层单元，将完成计算的一层单元的渗流虚域切掉，基于能量损失率最小计算此层的自由面点，逐层推进，直到汇入点；最后将每一层的自由面点连线得到完整的自由面曲线。本发明专利技术的有益效果是：对有电模拟试验解的矩形坝、有甘油模型试验解的矩形坝和有解析解的梯形坝的计算，具有很高的计算精度，避免了现有自由面曲线求解方法中需要迭代计算的缺陷。

A Method for Solving Free Surface of Seepage Flow Based on Stepwise Partition Method

【技术实现步骤摘要】
基于逐步剖分法求解渗流自由面的方法
本专利技术涉及一种基于逐步剖分法求解渗流自由面的方法，其属于坝工建设

技术介绍
坝体渗流自由面位置的确定在坝工建设中非常重要，坝体渗流自由面也被称为坝体的生命线，只有计算出自由面的位置，才能得到准确的渗流区范围。采用有限元法计算坝体稳定性，通常以坝体自由面为分割线，自由面以下视为饱和区，自由面以上为非饱和区。在稳定渗流求解中，自由面应同时满足水头边界条件和流量边界条件，而自由面位置是事先未知的。国内外学者解决该问题时，多根据水头边界条件或流量边界条件逐步迭代逼近计算渗流自由面，如自由面适应网格法是一种同时逼近水头和流量边界条件的迭代计算；虚单元法、丢单元法则是通过逼近水头边界条件的迭代计算；初流量法、改进初流量法、改进截止负压法是逼近流量边界条件的迭代计算。采用以上方法计算渗流自由面时，为减少误差，使计算结果更加接近实际的自由面，需提高迭代次数，计算量大。
技术实现思路
本专利技术提供一种基于逐步剖分法求解渗流自由面的方法，基于渗流域能量损失率最小求解稳定渗流场自由面的逐步剖分法，将渗流自由面求解转化为求解渗流域能量损失率最小值。本专利技术解决上述技术问题的技术方案如下：一种基于逐步剖分法求解渗流自由面的方法，根据达西定律和地下水运动的连续性条件，不考虑土体和水体的压缩性，均质各向异性土体的二维稳定渗流满足以下控制微分方程式中h是水头函数；kx、kx分别为x方向、z方向的渗透系数。对稳定渗流场，已知的边界条件有：①上、下游的水头边界条件，②自由渗出段的水头边界条件和流量边界条件，③底部不透水边界的流量边界条件，④渗流自由面的水头边界条件和流量边界条件，该边界是未知的。采用有限单元法计算渗流问题的关键是运用变分原理确定渗流微分方程的泛函。因此，根据给定的边界条件由(3)式求解渗流水头函数等价于求渗流域Ω内泛函的极值问题：在(2)式中，分别是x、z方向的渗流速度，分别是x、z方向的水力坡降，与渗透力成正比；相当于单位时间渗透力做的功(即功率)，数值上等于渗流的能量损失率。因此，式(2)的积分I(h)与整个渗流域的能量损失率成正比，后文以I(h)表示能量损失率。对式(2)求极小值，将所得线性微分方程组汇总成如下的矩阵形式[K]{h}＝{f}(1)式(1)中，[K]为总渗透矩阵，[h]是已知常数项列阵；{f}为节点水头列阵。渗流自由面是渗流域的上部边界，位于其上的部分实际不存在渗流，所以I(h)应不包含这个区域。基于此将整个渗流域分为两部分：渗流实域，即渗流自由面以下的部分BDAF；和渗流虚域，即渗流自由面以上的部分BDC，如图1所示。在求解渗流自由面时，式(2)的泛函只考虑渗流实域，即由上游汇入边界AB、下游流出边界EF、自由渗出边界DE，和渗流自由面BD围成的区域。依据最小功原理，在已确定上、下游及溢出面边界的渗流场中，当渗流达到稳定状态时，域内的能量损失率最小，即满足δI＝0。以图1所示的矩形坝AFCB为例。设上游水位为h1，下游水位为h2，AF为不透水边界，D为溢出点。以A为坐标原点，水平向右为x轴的正方向，垂直向上为z轴正方向建立直角坐标系。计算渗流场的能量损失率时，以下游水平面上一点o为坐标原点，水平向右为正方向，表示自由面节点位置不同时能量损失率的大小；垂直向上为正方向，表示不同自由面节点位置，建立如图1所示的直角坐标系。当渗流场达到稳定状态时，设B3D为真实的渗流自由面，渗流实域B3DFA内的能量损失率为I3；若其他条件不变，将渗流自由面的3点抬高到1点，自由面为B1D时，渗流实域B1DFA的能量损失率为I1，基于以上分析此时应有I1＞I3。若其他条件不变，将渗流自由面的3点压低到2点，自由面为B2D时，计算出渗流实域B2DFA的能量损失率为I2，此时也应有I2＞I3。因此，在上、下游及溢出边界确定的情况下，稳定渗流场中自由面总是位于使得域内能量损失率最小的位置。逐步剖分法原理：由渗流溢出点向上游入渗点逐层单元逆推，每逆推一层单元，将假定自由面以上单元丢弃，依据实域能量损失率最小确定待求渗流自由面节点的位置，直至逆推到入渗点位置，将每层求得的自由面节点连成光滑的曲线即为真实的自由面曲线。下面以图2a所示的均质矩形坝为例，介绍逐步剖分法求解渗流自由面的具体步骤。首先，依据横向能量损失率最小计算出渗流溢出点位置，如图2b中节点1所示。然后，以节点1(即溢出点)为起点，向上游方向逆推一层单元，设节点1左上方的相邻的节点2为此层单元的待求自由面节点。假定线段1-2为局部渗流自由面，赋予其自由面的边界条件，并将其以上单元丢弃，形成新的渗流计算模型，如图2b所示。通过依次调整假定自由面节点2的高程，基于能量损失率最小求解自由面点原理计算出实域能量损失率最小所对应的节点位置，即为该层单元真实的自由面点。如图3所示，当自由面节点分别位于2a、2b、2c位置时，实域能量损失率分别为I2a(h)、I2b(h)、I2c(h)。因为图3中节点位于2a处的实域能量损失率最小，所以2a即为所求自由面位置。依据上述步骤，依次逆推求解出各层单元渗流自由面节点的位置。将求解出的每层自由面节点、溢出点以及入渗点连成一条光滑曲线，即为真实的渗流自由面。逐步剖分法的网格划分要求，选取六节点三角形单元计算二维稳定渗流问题，对渗流域进行横向与竖向的规则剖分，所有网格均划分为四边形；然后将每个四边形单元沿对角线划分为两个三角形单元。为方便计算过程中自由面节点的调整，避免自由面节点出现跨单元调整造成单元畸形，竖向单元划分的尺寸取式(4)中，h1为上游水位，h3为溢出点的高程，L为上游汇入点与溢出点之间的水平距离，L0为单元横向尺寸。如图2所示，横向网格剖分时，距离溢出点1m的初始段内，应至少划分两层单元网格。因为以溢出点为起始点进行逐步剖分推进计算，在初始段，计算实域能量损失时渗流虚域的面积存在较大误差，对渗流自由面的计算结果影响较大，当此段内横向网格仅划分一个单元时，渗流虚域的面积误差对渗流自由面的结果集中在此段内，计算误差过大，也会增加后续计算中自由面节点的误差积累；当该段横向网格划分加密，在推进求解时，渗流虚域在数学模型中所占比例的逐渐减少会消减一部分误差，该段内自由面节点的计算误差会减小，后续计算中的积累误差也会减小。因此加密初始段横向单元划分，可提高整体自由面的计算精度。根据上述原理，采用Fortran语言编制了求解平面稳定渗流自由面的有限元程序，计算步骤如下：①选取六节点三角形单元对计算模型进行有限元划分，设在x轴方向共划分了m层单元。依据横向能量损失率最小确定溢出点位置；②由节点1(即溢出点)位置往上游逆推一层网格，如图2b所示，假定自由面节点的位置为高于溢出点位置的节点2，连接1-2并丢弃其上的单元：引入单元函数zb(i)以区分1-2自由面以上或以下的单元，zb(i)＝0表示以上的单元，将此类单元丢弃；zb(i)＝1表示以下的单元，保留此类单元；③将已知的上、下游和溢出面各节点水头值以及假定的自由面1-2的水头值代入到式(3)中，求解计算域内所有节点的水头值。引入单元分类函数pb(i)，若单元完全位于渗流虚区，即单元内三个角节点的水头值与高程的差值均小于零，则pb(i)＝0；若单元本文档来自技高网...

【技术保护点】
1.一种基于逐步剖分法求解渗流自由面的方法，其特征在于，采用Fortran语言编制求解平面稳定渗流自由面的有限元程序，其计算步骤如下：①选取六节点三角形单元对计算模型进行有限元划分，设在x轴方向共划分了m层单元，确定溢出点位置；②由溢出点即节点1位置往上游逆推一层网格，假定自由面节点的位置为高于节点1位置的节点2，连接节点1和节点2得到自由面1‑2，并丢弃其上的单元：引入单元函数zb(i)以区分节点1‑2面以上或以下的单元，zb(i)＝0表示以上的单元，将此类单元丢弃，zb(i)＝1表示以下的单元，保留此类单元；③将已知的上、下游和溢出面各节点水头值以及假定的自由面1‑2的水头值代入到由线性微分方程组汇总的矩阵形式[K]{h}＝{f}   (1)中，求解计算域内所有节点的水头值；在式(1)中，[K]为总渗透矩阵、[h]是已知常数项列阵、{f}为节点水头列阵；引入单元分类函数pb(i)，若单元完全位于渗流虚区，即单元内三个角节点的水头值与高程的差值均小于零，则pb(i)＝0；若单元完全位于渗流实区，即单元内三个角节点的水头值与高程的差值均大于零，则pb(i)＝1；④依据式

【技术特征摘要】
1.一种基于逐步剖分法求解渗流自由面的方法，其特征在于，采用Fortran语言编制求解平面稳定渗流自由面的有限元程序，其计算步骤如下：①选取六节点三角形单元对计算模型进行有限元划分，设在x轴方向共划分了m层单元，确定溢出点位置；②由溢出点即节点1位置往上游逆推一层网格，假定自由面节点的位置为高于节点1位置的节点2，连接节点1和节点2得到自由面1-2，并丢弃其上的单元：引入单元函数zb(i)以区分节点1-2面以上或以下的单元，zb(i)＝0表示以上的单元，将此类单元丢弃，zb(i)＝1表示以下的单元，保留此类单元；③将已知的上、下游和溢出面各节点水头值以及假定的自由面1-2的水头值代入到由线性微分方程组汇总的矩阵形式[K]{h}＝{f}(1)中，求解计算域内所有节点的水头值；在式(1)中，[K]为总渗透矩阵、[h]是已知常数项列阵、{f}为节点水头列阵；引入单元分类函数pb(i)，若单元完全位于渗流虚区，即单元内三个角节点的水头值与高程的差值均小于零，则pb(i)＝0；若单元完全位于渗流实区，即单元内三个角节点的水头值与高程的差值均大于零，则pb(i)＝1；④依据式计算渗流实域的能量损...

【专利技术属性】
技术研发人员：侯兴民，郑珊珊，
申请(专利权)人：烟台大学，
类型：发明
国别省市：山东,37