计算水力学06.ppt

计算水力学 第六章二维水流模拟 计算水力学 李光炽 第1节分步方法 李光炽 计算水力学 1基本方程描述二维水流运动的基本方程是浅水波方程 李光炽 计算水力学 其中x y t为平面坐标和时间 Z u v为水位和流速 为水深 为河床高程 分别为沿x和y方向的单宽流量 为风应力系数 为风速的大小 f为柯氏力系数 g为重力加速度 初始条件 李光炽 计算水力学 边界条件 开边界 陆边界 李光炽 计算水力学 第一分步 李光炽 计算水力学 第二分步 李光炽 计算水力学 对于第一分步 其等价方程可写成 即流进控制体的动量之和等于控制体中的动量增量 可视作在没有外力作用下的控制体中的动量守恒原理 李光炽 计算水力学 对于第二分步 其等价方程可写成 采用数值方法 对方程线性化可得 李光炽 计算水力学 即 和其中A为形成控制体四周的面积 为控制体自由水面面积 方程表示的是考虑水为不可压缩 流进控制体的水量之和等于控制体的体积变化 即控制体的质量守恒原理 2数值方法 李光炽 计算水力学 图1计算网格图Fig 1ComputationnetworksSketch p w e n s 李光炽 计算水力学 以W边为例 通过该控制面的动量通量为 式中函数表达使取决于所采用的差分格式 李光炽 计算水力学 若对流项采用逆风格式 扩散项采用中心差分 于是有 李光炽 计算水力学 同理有 李光炽 计算水力学 离散得u分量的方程 因此 对每一个控制面计算动量通量 向对应的单元 控制体 叠加 形成以单元变量为基本未知量的代数方程 每一个单元都可以建立这样的方程 由此形成一代数方程组 可唯一求解流速u 李光炽 计算水力学 同理 对流速v亦有 求解线性方程组 可得 对于方程可用类似的方法离散为 李光炽 计算水力学 李光炽 计算水力学 对每一个单元 都可建立这样的方程 对整个计算区域 形成以单元中心水位为基本未知量的线性代数方程组 可唯一地求解 由动量方程的线性表达式可求得 在建立方程式的过程中 可采用计算每一边的通量 向对应的单元叠加的方法形成 这样可使计算单元编码任意 无序化 方便灵活 易于扩充 以流量形式表达 与一维计算采用水位和流量作为基本变量的表达形式相同 这样使得二维计算易于同一维计算耦合联解 李光炽 计算水力学 本文的求解方法 最终归结为流速u的方程 流速v的方程和水位z的方程的求解 它们是高阶线性代数方程组 对于宽广的计算区域 单元众多 方程组的阶数很高 是一个巨大的矩阵 但是 该矩阵是高稀疏矩阵 方程组的求解效率决定着方法的成败 对于这一类方程组 可用高效的矩阵标识法求解 3矩阵标识法 第2节正交边界拟合坐标变换 李光炽 计算水力学 1正交边界拟合坐标由于河道的计算区域边界复杂 且长 宽尺度相差悬殊 因而在直角坐标系下对二维浅水问题进行求解存在着复杂边界不易拟合 网格多等困难 为此引进正交边界拟合坐标变换 将复杂的计算区域变换成规则的求解区域进行求解 在变换过程中可以根据需要布置网格的疏密 李光炽 计算水力学 边界拟合坐标是美国Mississippi州立大学J F Thompson等人提出的 其特点是可以把物理平面上任意形状的区域变换成计算平面上的规则区域 该方法的主导思路是寻找一种适当的变换 1 将物理平面上的复杂边界变换成计算平面上的规则边界 李光炽 计算水力学 假定计算平面的坐标系统与原物理平面的直角坐标系统之间满足如下的泊松方程 2 通过方程式 2 的变换 可以把 坐标平面上复杂的计算域转换成 平面上的矩形域 x y平面上的曲线网格变成 平面上间距为 的正方形网格 李光炽 计算水力学 为收缩因子 对变换结果有显著影响 其作用有二 适当选择 可使网格疏密根据需要分布 或使曲线网格正交 考察恒定 有压有势的二维水流运动方程 3 式中u v为垂线平均流速 h为水深 为流体密度 r为摩阻系数 为压力势函数 李光炽 计算水力学 对于二维有势运动 存在流函数 x y 流函数 与流速u v的关系为 4 由流函数与势函数的性质可知 曲线与曲线正交 李光炽 计算水力学 李光炽 计算水力学 从方程 3 可得 从方程 4 可得 李光炽 计算水力学 所得到的变换是正交变换 可以把物理平面上的复杂区域的正交曲线网格变换到计算平面上的规则均匀网格 2求解方法 李光炽 计算水力学 为了实现正交边界拟合坐标变换 必须对方程 2 进行数值求解 由于物理平面上的区域边界复杂 直接差分求解存在一定困难 因此 选择在计算平面上求解 在计算平面上 自变量为 因变量为x y 利用如下的微分关系 李光炽 计算水力学 把方程 2 转换为 13 李光炽 计算水力学 李光炽 计算水力学 李光炽 计算水力学 李光炽 计算水力学 李光炽 计算水力学 采用线Gauss Seidel迭代法 求解步骤如下 根据边界条件和工程计算的要求 事先假定区域的初始条件 即给定初始计算网格 在计算平面上对方程 13 进行离散 由于式 13 为一非线性的偏微分方程组 采用线性化的方法将其离散 其离散过程中 主要是将用初始计算网格进行计算 而其它项则采用隐式中心差商进行离散 李光炽 计算水力学 3 对于离散后的差分方程 用线Gauss Seidel迭代法求解新的节点坐标值 4 根据正交条件 滑动边界节点坐标值 以便边界处满足正交条件 5 重复2 3 4的步骤 直到计算的曲线网格全部满足正交条件为止 李光炽 计算水力学 第3节河道二维水流模型 李光炽 计算水力学 1 基本方程描述平面二维水流运动的基本方程组为 1 李光炽 计算水力学 4 李光炽 计算水力学 式中u v分别为沿 和 方向的流速 分别为曲线网格的长度和宽度 为曲线网格的面积 李光炽 计算水力学 与u v之间的变换关系为 通过正交变换 把原来在x y坐标系统中利用方程 1 求解变量Z u v变为在 坐标系统中利用 4 求解Z u v 下面为了书写方便 略去下标 2 边界条件 李光炽 计算水力学 1 固壁边界 采用不穿透条件 V n 0 n为固体边界的法向矢量 2 自由边界上边界 Z上为已知或 u为已知或 为已知及假定水位无横比降下边界 Z下为已知或 u为已知或 为已知及假定水位无横比降 3河道二维水流计算模型 李光炽 计算水力学 节点布置 图2物理平面上的节点示意图 李光炽 计算水力学 图3变换平面上的节点示意图 李光炽 计算水力学 为了便于边界条件的处理 变量采用交错布置 即在网格中心布置水位变量 在网格四边布置相应的流速变量 图 2 图 3 为物理平面及变换平面上的节点变量示意图 从图 3 中可见 流速u有N M个节点 流速v有 N 1 M 1 个节点 水位Z有 N 1 M个节点 李光炽 计算水力学 方程离散 连续方程式的差分在点 2k 1 2j 对方程 4 的连续方程 采用如下格式差分 非线性项 g hu和g hv的线性化采用如下公式 李光炽 计算水力学 采用中心差分 将上述差分近似和边界条件代入到连续方程可得如下线性差分方程组 上式写成矢量形式为 A1k Z2k 1 B1k Z2k 1 C1k Z2k 3 D1k V2k 1 E1k U2k F1k U2k 2 H1k k 1 2 N 式中 A1 B1k C1 E1 F1 为M M矩阵 D1k为M M 1 矩阵 H1 为M维矢量 动量方程的差分 李光炽 计算水力学 李光炽 计算水力学 将上述的差分近似及以边界条件代入到 方向的动量方程中得如下线性方程组 在点 2k 1 2j 1 对 方向动量方程类似于 方向动量方程的处理可得如下线性方程式 李光炽 计算水力学 A2k Z2k 1 B2k Z2k 1 C2k U2k 2 D2k U2k E2k U2k 2 F2k V2k 1 G2k V2k 1 H2k k 1 2 N 9 A3k Z2k 1 B3k U2k C3k U2k 2 D3k V2k 1 E3k V2k 1 F3k V2k 3 H3k k 1 2 N 式中 A2 B2 C2 D2 E2 为M M 矩阵 F2 G2 为M M 矩阵 A3 B3 C3 为 M 1 M矩阵 D3 E3 F3 为 M 1 M 1 矩阵 H2 为M维矢量 H3 为 M 1 维矢量 李光炽 计算水力学 方程求解 矩阵追赶法是一种针对二维河道水流计算特性所提所出的一种求解离散后的线性差分方程组的一种算法 V1 0上边界水位无横比降无横向流速 11 A21 Z1 B21 Z3 D21 U2 E21 U4G21 V3 H21 12 Z1 Z1 t 已知上游边界水位条件 10 李光炽 计算水力学 A1k Z2k 1 B1k Z2k 1 C1k Z2k 3 D1k V2k 1 E1k U2k F1k U2k 2 H1k 13 A3k Z2k 1 B3k U2k C3k U2k 2 D3k V2k 1 E3k V2k 1 F3k V2k 3 H3k 14 A2k 1 Z2k 1 B2k 1 Z2k 3 C2k 1 U2k D2k 1 U2k 2 E2k 1 U2k 4 F2k 1 V2k 1 G2k 1 V2k 3 H2k 1 15 k 1 2 N 1 李光炽 计算水力学 V2N 1 0下边界水面无横比降无横向流速 16 Z2N 1 Z2N 1 t 已知下边界水位过程 17 将边界条件 10 11 代入至式 12 求解得 U2 UZ1 Z3 UU1 U4 UV1 V3 UF1 18 将式 10 11 及 18 代入到式 14 可求得 V2k 1 VZk Z2k 1 VUk U2k 2 VVk V2k 3 VFk k 1 2 N 1 19 将式 10 18 及 19 代入到式 13 可求得 Z2k 1 ZZk Z2k 3 ZUk U2k 2 ZVk V2k 3 ZFk k 1 2 N 1 20 李光炽 计算水力学 将式 18 19 及 20 代入到式 15 可求得 U2k 2 Zk 1 Z2k 3 UUk 1 U2k 4 UVk V2k 3 UFk 1 k 1 2 N 1 21 按上述步骤逐步递推追赶直到k N 1可以得 U2N UZN Z2N 1 UUN U2N 2 UVN V2N 1 UFN 由下边界条件 16 17 式及对称假定U2N 2 U2N代入式 21 式求得U2N后逐步回代可求得Z2N 1 V2N 1 U2N 2 U2 完成该步长的求解 李光炽 计算水力学 其中0为零矢量 UZk UUk ZZk ZUk为M M矩阵 UVk ZVk为M M 1 矩阵 VZk VUk为 M 1 M矩阵 VVk为 M 1 M 1 矩阵 UFk ZFk为M维矢量 VFk为 M 1 维矢量 1 通过分析可以发现矩阵A1k B1k C1k D1k E1k F1k A2k B2k C2k D2k E2k F2k G2k A3k B3k C3k D3k E3k F3k均为三对角矩阵 且其中有许多为对角阵 因此进行运算时充分利用这些特性 其实际计算的工作量并不是很大 李光炽 计算水力学 2 矩阵追赶法主要是针对河道一类问题提出 而河道一类问题来说 河宽方向上的网格数据一般要比沿纵向要小得多 即M取值比N要小得多 所以前面导出的矩阵运算是小尺度的矩阵运算 3 由于本模型采用的是高精度离散格式 且其精度高 稳定性好 因而计算的时间步长可以取得很大 这样相对地减少了计算工作量 根据在各地的算例 Courant数可以达到400后 仍能满足精度要求 李光炽 计算水力学 李光炽 计算水力学 李光炽 计算水力学 李光炽 计算水力学 李光炽 计算水力学 第4节ADI方法 李光炽 计算水力学 在数值计算中应用比较广泛的一种格式为交替方向隐格式 这种格式为Peaceman与Rachford 1955 和Douglas 1955 同时提出 以简单的热传导方程为例 考查这种格式的优越性 1 李光炽 计算水力学 构造差分格式 这个格式第一个时间步为隐式逼近 而为显式逼近 这个格式第二