CN105843977A - 一种固定外边界的结构动网格生成方法 - Google Patents

Abstract

本发明属于CFD网格生成技术领域，具体涉及一种可用于强迫振荡、六自由度自由飞CFD计算的结构动网格生成方法。该方法包括如下步骤：(1)生成初始网格，并将初始网格分为内场与外场两个区间；(2)求解初始网格的源项；(3)旋转内场网格，给定旋转轴和转角，刚性旋转内部网格，得到新状态的内部网格；(4)求解椭圆形方程光顺外场网格。本发明所生成的网格保证了外边界的固定不动，同时保证了壁面处网格的正交性以及法向高度等特性，外场网格通过求解椭圆形方程光顺质量较好，在大角度转动后仍能保证较好的网格质量。

Description

一种固定外边界的结构动网格生成方法

技术领域

本发明属于CFD网格生成技术领域，具体涉及一种可用于强迫振荡、六自由度自由飞CFD计算的结构动网格生成方法。

背景技术

近年来，随着技术的进步和飞行器性能需求的升级，基于非定常流动的飞行器气动特性研究越来越成为关注的热点，而强迫运动和自由飞是研究飞行器非定常气动特性的重要途径。随着计算机技术和数值计算方法的快速发展，CFD已经成为解决上述非定常问题的有效途径，采用CFD计算非定常动态特性必须首先解决网格运动的问题。

现有结构动网格生成方法主要有两种：刚性旋转法和超限插值法。前者网格与物面刚性固连，随物体一起运动，网格形状不发生改变，该方法简单，适用于任意外形、任意运动，网格质量好，完全符合几何守恒率，但由于网格随体运动，使远场边界的位移和速度很大，导致运动远场边界处理不易；后者将插值法和刚性旋转法加权混合的动网格生成方法，接近外边界时，插值法网格的权重，靠近物面时，刚性旋转法网格权重，该方法保持外边界的静止，对复杂外形具有较好的适应能力，但是对于大幅度的运动会出现网格扭曲及负体积等现象，如大迎角俯仰运动。

发明内容

本发明的目的在于提供一种固定外边界的结构动网格生成方法，对初始网格固定外边界，划分内、外场，同时通过对靠近壁面的内场刚性旋转、对外场求解椭圆型方程进行光顺，获得高质量的结构动网格。

为达到上述目的，本发明所采取的技术方案为：

一种固定外边界的结构动网格生成方法，包括如下步骤：

(1)生成初始网格，并将初始网格分为内场与外场两个区间，s1为内场，s2为外场；

(2)求解初始网格的源项P、Q，

ξ_xx+ξ_yy＝(x_ξ ²+y_ξ ²)P(ξ,η)

η_xx+η_yy＝(x_η ²+y_η ²)Q(ξ，η) (1)

式中ξ、η为网格点在曲线坐标系中的坐标，x、y为网格点在物理坐标系中的坐标；

(3)旋转内场网格，给定旋转轴和转角，刚性旋转内部网格，得到新状态的内部网格；

(4)求解椭圆形方程光顺外场网格。

所述步骤(1)初始网格由Gridgen/Pointwise、ICEM工程软件生成；初始网格应保证良好的正交性及光滑过渡特性，同时保证外场的大小满足CFD计算要求，并可适当增大外场。

所述步骤(2)初始网格源项求解过程如下：

第一步：将式(1)由物理平面转化到计算平面，即将自变量和因变量对调，得到式(2)：

式中：其中为单位正交向量，x,y为网格点在物理平面内坐标值； $g_{11}={\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{r}}_{\mathrm{\ξ}}=x_{\mathrm{\ξ}}^2+y_{\mathrm{\ξ}}^2;g_{22}={\stackrel{\→}{r}}_{\mathrm{\η}}\·{\stackrel{\→}{r}}_{\mathrm{\η}}=x_{\mathrm{\η}}^2+y_{\mathrm{\η}}^2;g_{12}={\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{r}}_{\mathrm{\η}}=x_{\mathrm{\ξ}}{x}_{\mathrm{\η}}+y_{\mathrm{\ξ}}{y}_{\mathrm{\η}};$

初始网格在边界处满足正交性条件，见式(3)：

将式(3)代入式(2)得式(4)：

分别用r_ξ和r_η点积上式，解出该边界处源项的如下两个表达式：

式中，r_η，r_ηη，r_ξ，r_ξξ根据初始网格得坐标值差分得到；

第二步：得到边界源项后，内部源项采用以边界源项为基础的指数函数内插方式获得，见式(7)：

$P(\mathrm{\ξ},\mathrm{\η})=p\left(\mathrm{\ξ}\right)e^{-\mathrm{a\η}}+r\left(\mathrm{\ξ}\right)e^{-c({\mathrm{\η}}_{\max }-\mathrm{\η})}$

$Q(\mathrm{\ξ},\mathrm{\η})=q\left(\mathrm{\ξ}\right)e^{-\mathrm{b\η}}+s\left(\mathrm{\ξ}\right)e^{-d({\mathrm{\η}}_{\max }-\mathrm{\η})}---\left(7\right)$

式中，p(ξ)和r(ξ)分别是由式(5)决定的P(ξ,η_min)和P(ξ,η_max)，q(ξ)和s(ξ)分别是由式(6)决定的Q(ξ,η_min)和Q(ξ,η_max)；a，b，c，d都是正常数，表征由边界向流场内部源项的衰减速率，较小的常数表示较慢的衰减。

所述步骤(4)，外场网格的内外边界已经给定，外场各网格点初值为初始网格，通过迭代求解方程来获得高质量的外场网格，求解方法可以采用点松弛迭代方法、线松弛迭代方法。

所述点松弛迭代方法的步骤如下：

求解方程式(1)，外边界固定不动，内边界通过旋转求出，则η＝η_max,η_min上边界条件给定，外场采用O型网格，在ξ＝const边界上采用周期性边界条件；设I_max,J_max分别为ξ,η方向上的网格节点总数，两个方向上的步长均取为1，方程中涉及到的偏导数均采用中心差分格式，即：

逐点超松弛方法步骤如下：

第一步：先采用Gauss－Seidel改进迭代，计算的循环安排是j＝2,3,…,J_max-1为外循环，i＝1,2,3,…I_max为内循环，则方程(1)在(i,j)点尽量采用最新值的Gauss-Seidel改进迭代的第n场的差分方程为：

式中，g₁₁,g₂₂,g₁₂也尽量采用最新值算得，则可求得将之作为中间结果

第二步：引入超松弛迭代因子，将与前一步解进行加权平均，得到逐点超松弛迭代公式：

w的取值范围为[0,2]。

本发明所取得的有益效果为：

本发明将初始网格分为靠近物面的内场以及带外边界的外场两部分，由于壁面网格对法向第一层高度、正交性等要求较高，因此对内场采取刚性旋转的方法生成动网格，保证了初始网格在壁面处的网格质量，同时外边界固定不动，外场通过求解椭圆型方程进行光顺的方式得到，其中方程中的源项采用初始网格的源项，在方程迭代求解中减少了一层内迭代，缩短了网格生成时间。所生成的网格保证了外边界的固定不动，同时保证了壁面处网格的正交性以及法向高度等特性，外场网格通过求解椭圆形方程光顺质量较好，在大角度转动后仍能保证较好的网格质量。

附图说明

图1为网格内外场划分示意图；

图2为物理平面上的网格示意图；

图3为计算平面上的网格示意图；

图4为采用Gridgen软件生成的初始网格示意图；

图5为旋转30度示意图；

图6为旋转180度示意图；

图7为三维OREX飞船外形效果示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

本发明所述固定外边界的结构动网格生成方法包括如下步骤：

1、生成初始网格，并将初始网格分为内场与外场两个区间，见图1，其中s1为内场，s2为外场，ξ、η分别为二维曲线坐标系的两个坐标符。由于网格生成软件的普及，初始网格一般由Gridgen/Pointwise、ICEM等工程软件生成，以NACA0012翼型为例，采用Gridgen软件生成的初始网格见图4。

初始网格应保证良好的正交性及光滑过渡特性，同时保证外场的大小满足CFD计算要求，并可适当增大外场。

2、求解初始网格的源项，以二维网格为例，源项即式(1)中的P、Q，为后续外场求解椭圆型方程提供源项值。

ξ_xx+ξ_yy＝(x_ξ ²+y_ξ ²)P(ξ,η)

η_xx+η_yy＝(x_η ²+y_η ²)Q(ξ,η) (1)

式中ξ、η为网格点在曲线坐标系中的坐标，x、y为网格点在物理坐标系中的坐标，在求解椭圆型方程(式(1))光顺网格的过程中，源项起着控制网格的正交性以及过渡特性，且需再加入一层内迭代来求解，由于初始网格具有较好质量，为了缩短计算时间，后续源项可直接采用初始网格的源项而不需内迭代。初始网格源项求解过程如下：

第一步，将式(1)由物理平面转化到计算平面，即将自变量和因变量对调，得到式(2)，物理平面与计算平面上得网格转换见图2、图3。

式中：其中为单位正交向量，x,y为网格点在物理平面内坐标值

$g_{11}={\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{r}}_{\mathrm{\ξ}}=x_{\mathrm{\ξ}}^2+y_{\mathrm{\ξ}}^2$

$g_{22}={\stackrel{\→}{r}}_{\mathrm{\η}}\·{\stackrel{\→}{r}}_{\mathrm{\η}}=x_{\mathrm{\η}}^2+y_{\mathrm{\η}}^2$

$g_{12}={\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{r}}_{\mathrm{\η}}=x_{\mathrm{\ξ}}{x}_{\mathrm{\η}}+y_{\mathrm{\ξ}}{y}_{\mathrm{\η}}$

由于初始网格在边界处满足正交性条件，见式(3)，

将式(3)代入式(2)得式(4)：

分别用r_ξ和r_η点积上式，可解出该边界处源项的如下两个表达式：

式中，r_η，r_ηη，r_ξ，r_ξξ可根据初始网格得坐标值差分得到。

第二步，得到边界源项后，内部源项采用以边界源项为基础的指数函数内插方式获得，见式(7)：

$P(\mathrm{\ξ},\mathrm{\η})=p\left(\mathrm{\ξ}\right)e^{-\mathrm{a\η}}+r\left(\mathrm{\ξ}\right)e^{-c({\mathrm{\η}}_{\max }-\mathrm{\η})}$

$Q(\mathrm{\ξ},\mathrm{\η})=q\left(\mathrm{\ξ}\right)e^{-\mathrm{b\η}}+s\left(\mathrm{\ξ}\right)e^{-d({\mathrm{\η}}_{\max }-\mathrm{\η})}---\left(7\right)$

式中，p(ξ)和r(ξ)分别是由式(5)决定的P(ξ,η_min)和P(ξ,η_max)，q(ξ)和s(ξ)分别是由式(6)决定的Q(ξ,η_min)和Q(ξ,η_max)。a，b，c，d都是正常数，表征由边界向流场内部源项的衰减速率，较小的常数表示较慢的衰减。

类似的，可以求出另一边界(ξ＝const)上的源项分布。

3、旋转内场网格。给定旋转轴和转角，刚性旋转内部网格，得到新状态的内部网格，由于外场外边界固定不动，此时已经固定了外场网格的内外边界(η＝η_max,η_min)。

4、求解椭圆形方程光顺外场网格。外场网格的内外边界已经给定，外场各网格点初值为初始网格，通过迭代求解方程来获得高质量的外场网格，求解方法可以采用点松弛迭代方法、线松弛迭代方法等，为了使迭代稳定或加速收敛可以采用松弛法进行，方程中各项采用中心差分形式，以二维NACA0012翼型网格为例，采用逐点超松弛迭代方法的步骤推导如下：

求解方程式(1)，由于外边界固定不动，内边界通过旋转求出，则η＝η_max,η_min上边界条件给定，由于外场采用O型网格，在ξ＝const边界上采用周期性边界条件，而不必特别指定边界条件。设I_max,J_max分别为ξ,η方向上的网格节点总数，两个方向上的步长均取为1，方程中涉及到的偏导数均采用中心差分格式，即：

逐点超松弛(SOR)方法步骤如下：

第一步，先采用Gauss－Seidel改进迭代，尽可能利用最新的值进行计算，计算的循环安排是j＝2,3,…,J_max-1为外循环，i＝1,2,3,…I_max为内循环，则方程(1)在(i,j)点尽量采用最新值的Gauss-Seidel改进迭代的第n场的差分方程为：

式中，g₁₁,g₂₂,g₁₂也尽量采用最新值算得，则可求得将之作为中间结果

第二步，引入超松弛迭代因子，将与前一步解进行加权平均，得到逐点超松弛迭代公式：

w的取值范围为[0,2]。当w＝1时就是Gauss-Seidel改进迭代。

三维动网格生成方法可依此类推。

NACA0012翼型旋转30度后生成的网格见图5，旋转180度后生成的网格见图6。

三维OREX飞船外形旋转80度后生成的网格见图7。

Claims (5)

1.一种固定外边界的结构动网格生成方法，其特征在于：包括如下步骤：

(1)生成初始网格，并将初始网格分为内场与外场两个区间，s1为内场，s2为外场；

(2)求解初始网格的源项P、Q，

ξ_xx+ξ_yy＝(x_ξ ²+y_ξ ²)P(ξ,η)

(1)

η_xx+η_yy＝(x_η ²+y_η ²)Q(ξ,η)

式中ξ、η为网格点在曲线坐标系中的坐标，x、y为网格点在物理坐标系中的坐标；

(3)旋转内场网格，给定旋转轴和转角，刚性旋转内部网格，得到新状态的内部网格；

(4)求解椭圆形方程光顺外场网格。

2.根据权利要求1所述的固定外边界的结构动网格生成方法，其特征在于：所述步骤(1)初始网格由Gridgen/Pointwise、ICEM工程软件生成；初始网格应保证良好的正交性及光滑过渡特性，同时保证外场的大小满足CFD计算要求，并可适当增大外场。

3.根据权利要求1所述的固定外边界的结构动网格生成方法，其特征在于：所述步骤(2)初始网格源项求解过程如下：

第一步：将式(1)由物理平面转化到计算平面，即将自变量和因变量对调，得到式(2)：

式中：其中为单位正交向量，x,y为网格点在物理平面内坐标值； $g_{11}={\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{r}}_{\mathrm{\ξ}}=x_{\mathrm{\ξ}}^2+y_{\mathrm{\ξ}}^2;g_{22}={\stackrel{\→}{r}}_{\mathrm{\η}}\·{\stackrel{\→}{r}}_{\mathrm{\η}}=x_{\mathrm{\η}}^2+y_{\mathrm{\η}}^2;g_{12}={\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{r}}_{\mathrm{\η}}=x_{\mathrm{\ξ}}{x}_{\mathrm{\η}}+y_{\mathrm{\ξ}}{y}_{\mathrm{\η}};$

初始网格在边界处满足正交性条件，见式(3)：

将式(3)代入式(2)得式(4)：

分别用r_ξ和r_η点积上式，解出该边界处源项的如下两个表达式：

式中，r_η，r_ηη，r_ξ，r_ξξ根据初始网格得坐标值差分得到；

第二步：得到边界源项后，内部源项采用以边界源项为基础的指数函数内插方式获得，见式(7)：

$P(\mathrm{\ξ},\mathrm{\η})=p\left(\mathrm{\ξ}\right)e^{-\mathrm{a\η}}+r\left(\mathrm{\ξ}\right)e^{-c({\mathrm{\η}}_{\max }-\mathrm{\η})}$

(7)

$Q(\mathrm{\ξ},\mathrm{\η})=q\left(\mathrm{\ξ}\right)e^{-\mathrm{b\η}}+s\left(\mathrm{\ξ}\right)e^{-d({\mathrm{\η}}_{\max }-\mathrm{\η})}$

式中，p(ξ)和r(ξ)分别是由式(5)决定的P(ξ,η_min)和P(ξ,η_max)，q(ξ)和s(ξ)分别是由式(6)决定的Q(ξ,η_min)和Q(ξ,η_max)；a，b，c，d都是正常数，表征由边界向流场内部源项的衰减速率，较小的常数表示较慢的衰减。

4.根据权利要求1所述的固定外边界的结构动网格生成方法，其特征在于：所述步骤(4)，外场网格的内外边界已经给定，外场各网格点初值为初始网格，通过迭代求解方程来获得高质量的外场网格，求解方法可以采用点松弛迭代方法、线松弛迭代方法。

5.根据权利要求4所述的固定外边界的结构动网格生成方法，其特征在于：所述点松弛迭代方法的步骤如下：

求解方程式(1)，外边界固定不动，内边界通过旋转求出，则η＝η_max,η_min上边界条件给定，外场采用O型网格，在ξ＝const边界上采用周期性边界条件；设I_max,J_max分别为ξ,η方向上的网格节点总数，两个方向上的步长均取为1，方程中涉及到的偏导数均采用中心差分格式，即：

逐点超松弛方法步骤如下：

第一步：先采用Gauss－Seidel改进迭代，计算的循环安排是j＝2,3,…,J_max-1为外循环，i＝1,2,3,…I_max为内循环，则方程(1)在(i,j)点尽量采用最新值的Gauss-Seidel改进迭代的第n场的差分方程为：

式中，g₁₁,g₂₂,g₁₂也尽量采用最新值算得，则可求得将之作为中间结果

第二步：引入超松弛迭代因子，将与前一步解进行加权平均，得到逐点超松弛迭代公式：

w的取值范围为[0,2]。