一种基于Spectral/hp法的方腔顶盖驱动流的计算方法

摘要

本发明涉及一种基于Spectral/hp法的方腔顶盖驱动流的计算方法，属于计算流体力学技术领域。构建方腔顶盖驱动流的计算模型，该步骤包括：对无障碍物和有方形障碍物的方腔顶盖驱动流模型进行几何模型的建立；进行网格划分；流入速度与边界条件的施加；对计算模型进行高阶离散化，以划分的网格作为求解域，在该求解域上构造高阶多项式插值函数；对计算模型求解Navier‑Stokes方程，该步骤包括：求解弱压泊松问题，求解离散点处的速度。本发明能有效减少前处理时间，提供更高效的计算效率，得到具有较高精度的计算结果。

说明书

技术领域

本发明涉及一种基于Spectral/hp法的方腔顶盖驱动流的计算方法，属于计算流体 力学技术领域。

背景技术

方腔驱动流是由一个或多个壁面运动产生的内部再循环流动。方腔驱动流问题的求解不仅在技术上很重要，而且具有重要的科学意义，因为它在最简单的几何设置中 显示了许多重要的的流体力学现象，并且可以在相同的封闭几何结构中进行研究。

在h型有限元方法中，每个单元使用固定阶多项式作为插值基函数，并通过减小单元尺寸的大小来实现收敛。这就是经典的h型有限元法，其中h表示单元的尺寸， 该方法具有极高的几何灵活性。在p型有限元法中，采用固定网格，通过在每个单元 中增加插值多项式的阶数来提高数值解的收敛性，这就是p型有限元法，其中p表示 单元中插值多项式的阶数。高阶插值多项式有助于数值解的快速收敛。如果将整个求 解域作为一个单元来处理，那么p型有限元法就变成了谱方法。而hp型有限元法则是 通过同时减小单元尺寸的大小和提高插值多项式的阶数来加速数值解的收敛性。

与传统的低阶有限元法相比，Spectral/hp单元法理论上可以在解足够光滑的假设下 提供任意阶的空间逼近精度。而且，Spectral/hp单元法的紧凑性非常适合利用现代多核 计算硬件。所有上述特性都将Spectral/hp方法定位为一种用于以相对较低的计算成本 获得高精度数值解的计算方法。Spectral/hp离散化是连续和不连续Galerkin公式的基础 近似，可以在一维、二维和三维空间中构建。一般来说，计算网格可以由不同形状的 混合单元组成，这些单元可以是二维的三角形和四边形，也可以是三维的四面体、棱 锥、棱柱和六面体。在每个单元中使用多项式插值，对于连续伽辽金展开，这些单元 的插值函数通常可分为点模式、边模式和内部模式插值函数。

Spectral/hp单元法将经典h型有限元法的几何灵活性与谱方法和hp有限元法的高 精度数值特性相结合，在粗糙的有限元网格上采用高阶分段多项式基函数，空间近似基于正交多项式，例如勒让德或切比雪夫多项式。在计算和理论上，通过细分网格h 和增加多项式阶数p，可以获得高精度数值解并具有快速收敛的特性，特别是在某些条 件和假设下，数值解和精确解之间的近似误差可以实现指数级减小。

发明内容

针对现有技术存在的问题及不足，本发明提供一种基于Spectral/hp法的方腔顶盖 驱动流的计算方法。本发明通过以下技术方案实现。

一种基于Spectral/hp法的方腔顶盖驱动流的计算方法，其按以下步骤进行：

步骤1、构建方腔顶盖驱动流的计算模型，该步骤包括：对无障碍物和有方形障碍物的方腔顶盖驱动流模型进行几何模型的建立；进行网格划分；流入速度与边界条件 的施加；

步骤2、对计算模型进行高阶离散化，以步骤1划分的网格作为求解域，在该求解域上构造高阶多项式插值函数，其具体步骤如下：

一个求解域Ω，将它划分为一个包含N

对于域Ω＝{x|0＜x＜l}，一维情况下用点x

第e个单元定义为：

Ω

引入一维标准单元：

Ω

根据局部坐标ξ定义一个在Ω

标准单元Ω

χ

通过将标准单元Ω

在标准单元Ω

二维情况下的多项式插值类似于一维的情况，所有多项式插值都将在标准区域Ω

式中

以二维情况下的四边形单元为例，如图7所示：

四个顶点A、B、C、D的点模式基函数构造如下：

对应于四边形四个边AB、CD、AC、BD的边模式定义如下：

四边形单元插值函数的内部模式定义为：

步骤3、对计算模型求解Navier-Stokes方程，该步骤包括：求解弱压泊松问题， 求解离散点处的速度，步骤3具体步骤包括：

在求解域内取粘性不可压N-S方程内积建立压力泊松方程；

首先写出如下形式的黏性不可压Navier-Stokes方程：

其中u为速度，p为压力，v为运动粘度，f为外力，t为时间，式中将密度简化为ρ＝1并合到压力p中。

步骤3.1、通过取方程在解域Ω上与

其中

通过将方程最后一项设置为零来使散度为零。将方程中的第1、2和最后一项使用分部积分，得到弱压方程：

其中

其中

最后，引入了非线性项和旋涡项的一致外推形式：

最终得到弱压力近似公式：

这个公式可以重新转换为压力泊松方程的等效强形式：

一致的Neumann边界条件规定为：

式中

步骤3.2、在时间为n+1处，对离散方程使用上一步中n+1时间处压力来求解速 度u

这就变成了Helmholtz问题：

式中

求解这个Helmholtz方程，就得到了速度场，完成了对时间n+1时的速度u

所述步骤3.2求解得到离散点处的速度后，通过下列公式求得主涡、二次涡中心位置的坐标、涡量和流函数；

涡量ω的定义：

流函数

通过带入不同位置的x方向上的速度u，y方向上的速度v，从而得到不同位置的 涡量与流函数。

流函数的等值线是流线，可以直观的表达流场的流动特性，从而判断涡产生的位置坐标，得到主涡，二次涡中心位置的坐标。

上述各标号为本领域技术人员公知的标号含义。

本发明的有益效果是：

(1)本发明基于Spectral/hp单元法求解方腔顶盖驱动流，提供了一种新的计算思路；

(2)本发明能有效减少前处理时间，提供更高效的计算效率，得到具有较高精度的计算结果。

附图说明

图1是本发明实施例1的方腔顶盖驱动流几何参数示意图；

图2是本发明实施例1的方腔顶盖驱动流涡结构命名法和边界条件示意图；

图3是本发明实施例1的无障碍物Re＝100时不同网格数下中心垂直线上速度分量u图；

图4是本发明实施例1的无障碍物不同Re稳定状态下的流函数等值线；

图5是本发明实施例2的有方形障碍物的方腔顶盖驱动流几何参数示意图；

图6是本发明有方形障碍物的l＝0.2时不同雷诺数下流线图；

图7是二维四边形单元示意图。

具体实施方式

下面结合附图和具体实施方式，对本发明作进一步说明。

实施例1

该基于Spectral/hp法的方腔顶盖驱动流的计算方法，按以下步骤进行：

步骤1、构建方腔顶盖驱动流的计算模型，该步骤包括：对无障碍物的方腔顶盖驱动流模型进行几何模型的建立；进行网格划分；流入速度与边界条件的施加；

方腔顶盖驱动流几何参数如图1所示，高为H＝1、宽W＝1，高宽比B＝H/W＝1， 在腔顶部边界上速度恒定u＝1、v＝0，边界条件设置为：u＝1，v＝0，

步骤2、对计算模型进行高阶离散化，以步骤1划分的网格作为求解域，在该求解域上构造高阶多项式插值函数；

四边形单元插值函数的内部模式定义为：

在Re＝100情况下，分别采用2×2、4×4、8×8、16×16、32×32五种不同数量的网格进行计算，每种网格情况下分别计算得到了p从1到8的稳定解；

步骤3、对计算模型求解Navier-Stokes方程，该步骤包括：求解弱压泊松问题， 求解离散点处的速度；

所述步骤3具体步骤包括：

步骤3.1、在求解域内取粘性不可压N-S方程内积建立压力泊松方程；

建立的压力泊松方程为

该方程一致的Neumann边界条件规定为

式中

步骤3.2、利用步骤3.1中的得到的压力泊松方程求解得到压力，然后再通过压力求解离散点处的速度；

将步骤3.1的压力泊松方程使用时间倒数的近似表示，获得新的方程如下：

式中

通过新方程得到了速度场，完成了对时间n+1时的速度u

并将方腔中心垂直线(x＝0.5)上的水平速度分量u的计算值，如图3所示，

所述步骤3.2求解得到离散点处的速度后，通过下列公式求得主涡、二次涡中心位置的坐标、涡量和流函数；

涡量ω的定义：

流函数

通过带入不同位置的x方向上的速度u，y方向上的速度v，从而得到不同位置的 涡量与流函数。

流函数的等值线是流线，可以直观的表达流场的流动特性，从而判断涡产生的位置坐标，得到主涡，二次涡中心位置的坐标。

根据速度将计算所得到的主涡、二次涡中心位置的坐标、涡量、流函数的结果如表1所示。

表1无障碍物主涡和二次涡的特性：流函数、涡量与位置

由图3可以看出，对于所有的网格情况计算值曲线都随着多项式阶数p的提高曲线更加平滑，说明数据的有效性。在16×16和32×32的网格下，多项式阶p＝1时的计 算结果就平滑，说明小步长h下更容易收敛。而在2×2网格下，多项式阶直到p＝8时 计算结果误差依然很大，对于4×4的网格也是如此，过度粗糙的网格划分未能在多项 式阶p≤8的情况下得到有效的解，需要继续提高插值多项式阶数来减小误差。在8×8 的网格下，随着p的增加，水平速度分量u的快速收敛，p≥4时的计算结果就与其他 有效结果吻合。因此，对Re＝100的情况，8×8的网格和p＝4的多项式阶是最为高效 的选择。

由图4可以看出，随着雷诺数的增大，流动趋于复杂，在Re＝1000时左下角和右 下角的二次涡BL1、BR1更加明显；Re＝2500时，左上角出现了一个二次涡TL1；Re＝5000 时，左下角和右下角二次涡增加至两个，出现了BL2与BR2；Re＝7500时，BL2与BR2 更加明显。

实施例2

该基于Spectral/hp法的方腔顶盖驱动流的计算方法，按以下步骤进行：

步骤1、构建方腔顶盖驱动流的计算模型，该步骤包括：对加入一个正方形障碍物的方腔顶盖驱动流模型进行几何模型的建立；进行网格划分；流入速度与边界条件的 施加；

加入一个正方形障碍物的方腔顶盖驱动流几何参数如图5所示，方腔为边长W＝1的正方形，腔内部中心位置处有一个边长为l＝0.2的正方形障碍物，在腔顶部边界上 速度恒定u＝1、v＝0，边界条件设置为：u＝1，v＝0，

步骤2、对计算模型进行高阶离散化，以步骤1划分的网格作为求解域，在该求解域上构造高阶多项式插值函数；

四边形单元插值函数的内部模式定义为：

在Re＝100、500、1000、5000情况下，采用16×16的网格进行计算，计算得到了 p为8的稳定解；

步骤3、对计算模型求解Navier-Stokes方程，该步骤包括：求解弱压泊松问题， 求解离散点处的速度；

所述步骤3具体步骤包括：

步骤3.1、在求解域内取粘性不可压N-S方程内积建立压力泊松方程；

建立的压力泊松方程为

该方程一致的Neumann边界条件规定为

式中

步骤3.2、利用步骤3.1中的得到的压力泊松方程求解得到压力，然后再通过压力求解离散点处的速度；

将步骤3.1的压力泊松方程使用时间倒数的近似表示，获得新的方程如下：

式中

通过新方程得到了速度场，完成了对时间n+1时的速度u

所述步骤3.2求解得到离散点处的速度后，通过下列公式求得主涡、二次涡中心位置的坐标、涡量和流函数；

涡量ω的定义：

流函数

通过带入不同位置的x方向上的速度u，y方向上的速度v，从而得到不同位置的 涡量与流函数。

流函数的等值线是流线，可以直观的表达流场的流动特性，从而判断涡产生的位置坐标，得到主涡，二次涡中心位置的坐标。

根据速度将计算所得到的主涡、二次涡中心位置的坐标、涡量、流函数的结果如表2所示。

表2有障碍物主涡和二次涡的特性：流函数、涡量与位置

如图6所示，在低雷诺数(Re＝100)时，腔内有一个明显的主涡，且随着雷诺数 增大逐渐靠近障碍物的上边界。在Re＝5000时，不仅左下角和右下角的角涡范围变大， 还在左上角附近产生了一个诱导涡，在方柱的的左壁和下壁分别出现了一个黏性效应 导致的次级涡。而主涡随着雷诺数的增大而逐渐向着腔内方柱的右上角收缩。

综合以上两个实施例可以得出如下结论：本发明的Spectral/hp单元法求解方腔顶 盖驱动流问题是有效正确的，进行计算时所需网格数较少(前处理工作量大大减少)，且收敛速度较快。

以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实 施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。