含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法

1.含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法，其特征在于，所述方法包括以下步骤：S1、构造与边界面平行的相似面，定义比例边界有限元SBFEM坐标系：环向坐标 η 、 ζ 平行边界面，径向坐标 ξ 垂直边界面，给出问题域中任意点的坐标表达式，引入形函数对边界面和相似面上的点进行离散，求得面单元内问题域中任意点在SBFEM坐标系中的映射，通过坐标变换的雅可比矩阵实现笛卡尔坐标系与SBFEM坐标系的转换，进一步求得雅可比矩阵的逆矩阵并推导得到瞬态渗流问题的微分算子表达式；S2、以刻画系统能量状态的总能量泛函为变分基础，推导瞬态渗流的能量泛函表达式，在SBFEM坐标系内离散化场变量，获取问题域内任意点水头和梯度的插值形式，将以上插值代入总能量泛函，并在环向坐标 上进行积分，将总能量泛函表达式转换为仅依赖于径向坐标的离散形式，通过应用第一变分原理，令总能量泛函的驻点的一阶导数为零，并采用分部积分法从弱形式方程中提取微分方程，得到了关于水头的SBFEM坐标系下的控制方程；S3、基于虚功原理与加权残值法获取内部等效节点通量，根据内外部等效节点通量的相反性，描述节点水头与外部等效节点通量的关系，并进一步推导得到瞬态渗流矩阵的控制方程，引入无量纲频率和无量纲渗流矩阵，对无量纲渗流矩阵分别求关于径向坐标 和频率 的偏导，通过变量替换将原方程转化为频率域内的瞬态渗流矩阵控制方程；S4、应用连续分式法将 趋于无穷远的瞬态渗流矩阵分解为幂级数形式，通过求解一般特征值问题获取系数矩阵 ，将高阶项系数设为零，通过求解广义特征值问题和Lyapunov方程得到关于归一化频率 的待求解函数表达式，再通过幂级数展开剩余的低阶项，递归计算后续系数矩阵，设定终止条件截断展开式，得到瞬态渗流矩阵的表达式；S5、令频率 ，得到稳态渗流问题的控制方程，将控制方程表述为特征值问题，通过Schur分解求解得到稳态渗流矩阵，描述问题域边界面上通量与水头之间的关系，并应用傅里叶逆变换得到时间域的水头表达式，采用改进精细时程积分法MPTSIM，采用矩阵指数计算得到任意时间步的水头值，通过执行递归循环，获得边界面上水头各时间步的水头值，定义水头关于径向坐标的函数，采用四阶Runge-Kutta法求解半无限域水头分布。 2.如权利要求1所述的含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法，其特征在于，在所述步骤S1中，构造与边界面平行的相似面，定义比例边界有限元SBFEM坐标系：环向坐标 η 、 ζ 平行边界面，径向坐标 ξ 垂直边界面，给出问题域中任意点的坐标表达式，具体如下：其中， , 和 表示问题域中点的坐标， , 和 表示相似面上点的坐标， , 和 表示边界面上点的坐标， 表示SBFEM坐标系的径向坐标；引入形函数对边界面和相似面上的点进行离散：其中， 表示面单元的形函数， ， ， 分别为面单元内相似面上点的节点坐标向量； ， ， 表示面单元内边界面上点的节点坐标向量；求得面单元内问题域中任意点在SBFEM坐标系中的映射：其中， ， ， 表示面单元内边界面上点的节点坐标向量；通过坐标变换的雅可比矩阵实现笛卡尔坐标系（ x , y , z ）与SBFEM坐标系（ ξ , η , ζ ）的转换，雅可比矩阵的表达式如下：其中， 表示雅克比矩阵，下标 ， ， 分别表示对径向坐标 和环向坐标 ， 求一阶导数；求得雅可比矩阵的逆矩阵：其中，上标-1表示矩阵的逆， 是雅可比矩阵的行列式， ， ， ， ， ， ， ， ， 为雅克比矩阵逆矩阵中的所有元素；推导得到瞬态渗流问题的微分算子表达式：其中， 表示微分算子； ， ， 分别表示对 x ， y ， z 求偏导； ， ， 分别表示对 ξ ， η ， ζ 求偏导； ， ， 表示子系数向量； 。 3.如权利要求2所述的含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法，其特征在于，在所述步骤S2中，以刻画系统能量状态的总能量泛函为变分基础，推导瞬态渗流的能量泛函表达式，具体表达式为：其中， 表示总能量泛函， 表示耗散能量泛函， 表示储存能量泛函， 表示水头， 表示贮水系数， k 表示渗透系数， 表示水力梯度， 表示频率， 表示整个计算域， dV 表示对体积分， i 表示虚数，上标 T 表示矩阵的转置；在SBFEM坐标系内离散化场变量，获取问题域内任意点水头和梯度的插值形式：其中， 表示SBFEM坐标系下的水头插值形式； 表示面单元的水头形函数矩阵； ， 表示坐标变换矩阵，其定义如下： ， 将以上插值代入总能量泛函，并在表面坐标 上进行积分，将总能量泛函表达式转换为仅依赖于径向坐标的离散形式，获得以下关系：其中， 表示总能量泛函的一阶变分； ， 表示径向坐标的起点值和终点值； 表示连接相似面单元节点与边界面单元节点的射线方向水头向量； 表示渗透系数矩阵； 表示对径向坐标积分； ， ， 表示瞬态渗流系数矩阵； 表示瞬态渗流质量矩阵，表达式定义如下：其中， ， 表示对环向坐标 ， 积分， 表示在整个域求和；通过应用第一变分原理，令总能量泛函的驻点的一阶导数为零，具体表达式为：其中， 为变分符号；采用分部积分法从弱形式方程中提取微分方程，得到了关于水头的SBFEM坐标系下的控制方程；其中， 表示对径向坐标 求二阶导数。 4.如权利要求3所述的含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法，其特征在于，在所述步骤S3中，基于虚功原理与加权残值法获取内部等效节点通量，具体表达式为：其中， 为内部等效节点通量；根据内外部等效节点通量的相反性，描述节点水头与外部等效节点通量的关系：其中， 为外部等效节点通量， 为瞬态渗流矩阵；推导得到瞬态渗流矩阵的控制方程：引入无量纲频率：其中， 是无量纲频率， r 是环向边界面的面积，上标2表示平方；引入无量纲渗流矩阵，无量纲渗流矩阵与瞬态渗流矩阵存在如下关系：其中， 为无量纲渗流矩阵；对无量纲渗流矩阵分别求关于 和 的偏导，得到如下表达式：其中， 为 的简称， 为对频率 求偏导；通过变量替换将原方程转化为频率域内的瞬态渗流矩阵控制方程：其中， 为 趋于无穷远的瞬态渗流系数矩阵和质量矩阵。 5.如权利要求4所述的含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法，其特征在于，在所述步骤S4中，应用连续分式法将渗流矩阵分解为幂级数形式：其中， 表示 趋于无穷远的瞬态渗流矩阵； 表示关于 的 o 阶函数； ， 为幂级数展开中的系数矩阵； 表示归一化频率， 与频率 之间存在如下表达式：通过求解一般特征值问题获取系数矩阵 ：其中， 表示特征向量矩阵； 表示特征值矩阵；将高阶项系数设为零，通过求解广义特征值问题和Lyapunov方程得到关于 的待求解函数的表达式：其中， 表示 o 阶函数的初始值； 表示 o 阶函数的递归计算值； 表示关于 的 o +1阶函数；通过幂级数展开剩余的低阶项，并递归计算后续系数矩阵，设定终止条件截断展开式，得到瞬态渗流矩阵的表达式为：其中， 表示1阶函数的初始值； 表示1阶函数的递归计算值； 表示2阶函数的初始值； 表示2阶函数的递归计算值。 6.如权利要求5所述的含三维非均质半无限域堤坝瞬态渗流的时域积分求解方法，其特征在于，在所述步骤S5中，令 ，得到稳态渗流问题的控制方程：其中 是稳态渗流矩阵；将控制方程表述为特征值问题：其中， 为特征向量矩阵， 为特征值矩阵， 为哈密顿矩阵；上式显式表示为：其中， 为单位矩阵； ， ， ， 为特征向量子矩阵； 为特征值子矩阵；通过Schur分解求解得到稳态渗流矩阵：描述问题域边界面上通量与水头之间的关系，具体表达式为：其中， 为 的简写； 为 的简写； 为边界面的质量矩阵；并将傅里叶逆变换应用于上述方程得到时间域的水头表达式：其中， 为时间域的水头， 为水头的一般形式；采用改进精细时程积分法MPTSIM，采用矩阵指数精确计算得到任意时间步的水头值，表达式如下：其中， t 为任意历史时间点， ， 分别为第 k 和 k +1的时间步对应的时间点； ， 分别为第 k 和 k +1的时间步对应的时间点的水头； 为初始水头； 为状态矩阵； 为外部源项； dt 为对时间积分；exp为指数函数；通过执行递归循环：其中， 为迭代次数； N 为总迭代次数； 为泰勒展开指数矩阵；获得边界面上水头各时间步的水头值：其中， 为指数矩阵； 为时间步长； n 为泰勒展开级数； 为2 n- 1次的高阶可忽略项； 和 分别为高斯积分的权重和高斯点； 为从1到 l 进行求和， l 表示积分点个数；定义水头关于径向坐标的函数为：其中， 表示对应于点 处水头关于径向坐标的函数； 表示对应于点 处的瞬态渗流矩阵； 表示对应于点 处的外部等效节点通量；采用四阶Runge-Kutta法求解半无限域中的水头分布，具体表达式如下：其中， 表示第 j 步求得的 处水头； 表示第 j -1步求得的 处水头； 为迭代步长； , , , 分别对应于点 , , , 的梯度，具体表达式为：其中， 表示第 j 步时 处水头关于径向坐标的函数； , , 表示第 j -1步时对应于点 , , 处的瞬态渗流矩阵； , , 表示第 j -1步时对应于点 , , 处的外部等效节点通量。