一种预报超疏水表面曲面边界层转捩位置的方法

摘要

本发明涉及一种预报超疏水表面曲面边界层转捩位置的方法，包括下列步骤：划分计算域网格，其中，在曲面几何体壁面的近壁区，绘制一层或者多层垂直曲面的网格；模拟曲面几何体壁面的的超疏水表面，列出壁面滑移速度边界条件；在计算域网格上求解NS方程，计算超疏水表面曲面边界层的流场，并插值到曲面几何体表面沿法向的正交化网格下，获得各处的壁面切向速度分量剖面；基于线性稳定性理论，分析超疏水表面曲面边界层流场的流动稳定性，获得各频率扰动波幅值沿流向的增长率；计算其最大增长倍数，预报超疏水表面曲面边界层的自然转捩位置。

权利要求

1.一种预报超疏水表面曲面边界层转捩位置的方法，包括下列步骤：第一步，划分计算域网格，其中，在曲面几何体壁面的近壁区，绘制一层或者多层垂直曲面的网格；

第二步，根据第一步划分的计算域网格，模拟曲面几何体壁面的超疏水表面，列出壁面滑移速度边界条件，方法如下：使用第一步中曲面几何体壁面网格点的数据，通过计算得到曲面几何体任意位置表面的切线方向，设定此方向为超疏水表面当地的滑移速度方向；

在每一个垂直曲面的近壁区网格单元内计算当地的壁面法向的速度梯度；使用线性滑移速度边界条件假设，如(1)式所示，列出超疏水表面曲面几何体的壁面滑移速度边界条件：* *

式中 为超疏水表面曲面几何体沿壁面切向的滑移速度，U为切向速度，y 为壁面法*向，λ为滑移长度，上标“*”表示有量纲的物理量；

第三步，根据第二步设定的壁面滑移速度边界条件，在第一步中的计算域网格上求解NS方程，计算超疏水表面曲面边界层的流场，并插值到曲面几何体表面沿法向的正交化网格下，获得各处的壁面切向速度分量剖面；

第四步，基于线性稳定性理论，分析第三步中超疏水表面曲面边界层流场的流动稳定性，获得各频率扰动波幅值沿流向的增长率，方法如下：从不可压缩流动的无量纲NS方程出发，将流动的瞬时量分解成基本流与扰动量之和，T即 其中φ＝(v1,v2,v3,p) 为瞬时量矢量， 为基本流矢量，φ'＝T

(v'1,v'2,v'3,p') 为扰动量矢量，v1,v2,v3分别为流向、法向和周向的速度分量，p为压力，上标“‑”代表基本流，上角标“'”代表扰动量；通过适用于水介质的流向和周向带曲率的曲面边界层流动的线性稳定性分析，得到扰动的线性稳定性方程，即：式中 ， 代 表 将扰 动 φ '写 成 行 进 波形 式 后的 特 征 函 数矢 量 ，即其中t为时间，q1,q2和q3分别为正交曲线坐标系下当地的流向、法向和周向空间坐标，ω为圆频率，α为流向波数，β为周向波数，c.c.为复共轭；在空间模式中，α为复数，表示为α＝αr+iαi,其中，‑αi为扰动的空间增长率，αr为实数形式的流向波数； 为系数矩阵，包括基本流、壁面的流向曲率drx和周向曲率drz、扰动圆频率、波数和雷诺数参数；

扰动在超疏水表面曲面边界层传播时，其流向速度分量v'1和周向速度分量v'3在壁面处适用第二步中列出的滑移速度边界条件，线性稳定性方程(2)的在q2＝0处的壁面边界条件为：基于第三步中得到的超疏水表面曲面边界层流场，在每一个流向位置处求解线性稳定性方程(2)，得到超疏水表面边界层内各频率扰动波幅值沿流向的空间增长率‑αi；

第五步，使用第四步中所获得的不同频率扰动波幅值的增长率，计算其最大增长倍数，预报超疏水表面曲面边界层的自然转捩位置，方法如下：对第四步所获得的波幅值沿流向的空间增长率‑αi沿流向积分，得到某一流向位置q1处扰动波幅值的放大因子N，即：式中，f为某一扰动的频率大小，q10为积分的起始位置，取该频率扰动开始失稳的流向位置；

计算不同频率的扰动波在每一个流向位置处的最大增长倍数 其中Nmax为所有频率扰动在当地q1处N值的最大值，即满足当最大增长倍数 达到转捩判据给出的临界值 时，认为层流转捩为湍流，此时所对应的流向位置q1就是发生转捩的位置。

2.根据权利要求1所述的预报超疏水表面曲面边界层转捩位置的方法，其特征在于，第一步中的近壁区垂直网格划分用于计算每一个位置处壁面法向的速度梯度。

3.根据权利要求1所述的预报超疏水表面曲面边界层转捩位置的方法，其特征在于，第五步中的转捩判据取

4.根据权利要求1所述的预报超疏水表面曲面边界层转捩位置的方法，其特征在于，将第四步中的流向曲率和周向曲率均设定为零，预报超疏水表面平面边界层转捩位置。

说明书

一种预报超疏水表面曲面边界层转捩位置的方法

技术领域

[0001] 本发明属于水动力学研究技术领域，具体地，涉及一种预报超疏水表面曲面边界层转捩位置的方法。

背景技术

[0002] 进入本世纪以来，超疏水表面的滑移特性引起了人们的广泛关注和研究。通过在水下航行器壳体材料表面制备超疏水涂层，可以直接影响壁面剪切应力状态，使得壁面速度大于零，最终导致湍流转捩点向下游移动，这也使其在减阻降噪等方面表现出了极为诱人的应用前景。而回转体/水翼作为典型的水下曲面航行器，与潜艇在水中的运动、舰船在水面的航行以及鱼雷在水下的发射等实际问题密切相关，那么如何可行且合理地预报超疏水表面曲面边界层的转捩位置，对于水下曲面航行器超疏水涂层技术的发展和应用显得尤为关键。

[0003] 由于曲面边界层流动比平面边界层流动更复杂，因此相关的转捩预报研究相对较少，而超疏水表面的曲面边界层又比普通表面的曲面边界层更为复杂，目前尚未见到相关的转捩预报研究工作。究其原因，针对这一问题的转捩预报技术主要存在两个难点：第一，如何模拟曲面上的超疏水表面，并计算出曲面边界层的流场；第二，如何分析超疏水表面曲面边界层的流动稳定性，预报转捩位置。近些年，已经发展出一些相关的技术和方法，为解决第一个难点奠定了一定的基础。超疏水表面平板边界层的计算相对简单，已经有了相应的计算方法(Liu,B.and Zhang,Y..A numerical study on the natural transition locations  in the flat‑plate boundary  layers on  superhydrophobic surfaces.Phys.Fluids  32,124103,2020；Liu,B.and  Zhang,Y..Numerical 
investigation of the natural transition in flat‑plate boundary layers on superhydrophobic surfaces considering the influence of the leading edge region.AIP Advances 12,035140,2022)。有关第二个难点的研究，需要同时考虑壁面滑移速度条件和壁面曲率，目前尚未见到这样的工作发表。但是，单独考虑壁面滑移速度条件的工作，已经有人做了，针对的是超疏水表面平板边界层(Liu,B.and Zhang,Y..A numerical study on the natural transition locations in the flat‑plate boundary layers on superhydrophobic surfaces.Phys.Fluids 32,124103,2020；Liu,B.and Zhang,Y..Numerical investigation of the natural transition in flat‑plate boundary layers on superhydrophobic surfaces considering the influence of the leading edge region.AIP Advances 12,035140,2022)；单独考虑壁面曲率的工作也已经发表，针对的是无滑移条件的普通表面曲面边界层，曲面外形比较复杂，壁面的流向和周向曲率都已考虑在内(Liu,J.,Chu,X.,and Zhang,Y..Numerical investigation of natural transitions  of  bow  boundary layers  over  underwater axisymmetric 
bodies.Phys.Fluids 33,074101,2021)。虽然这些工作没有同时考虑超疏水表面滑移速度条件和壁面曲率，但是为解决第二个难点奠定了基础。

[0004] 综上所述，针对超疏水表面曲面边界层转捩位置，探索一种可行且合理的预报技术和方法具有重要指导意义。

发明内容

[0005] 本发明的目的在于提供一种可行且合理的预报超疏水表面曲面边界层转捩位置的方法。该方法不仅为回转体/水翼等水下航行体曲面边界层转捩的准确预报提供了一种有效的方法，而且还为超疏水表面曲面滑移速度边界条件的实现提供了一种可行的建模方式。为实现上述目的，本发明所采用下述技术方案：

[0006] 一种预报超疏水表面曲面边界层转捩位置的方法，包括下列步骤：

[0007] 第一步，划分计算域网格，其中，在曲面几何体壁面的近壁区，绘制一层或者多层垂直曲面的网格。

[0008] 第二步，根据第一步划分的计算域网格，模拟曲面几何体壁面的的超疏水表面，列出壁面滑移速度边界条件，方法如下：

[0009] 使用第一步中曲面几何体壁面网格点的数据，通过计算得到曲面几何体任意位置表面的切线方向，设定此方向为超疏水表面当地的滑移速度方向；

[0010] 在每一个垂直曲面的近壁区网格单元内计算当地的壁面法向的速度梯度；使用线性滑移速度边界条件假设，如(1)式所示，列出超疏水表面曲面几何体的壁面滑移速度边界条件：

[0011]* *

[0012] 式中 为超疏水表面曲面几何体沿壁面切向的滑移速度，U 为切向速度，y为壁*面法向，λ为滑移长度，上标“*”表示有量纲的物理量。

[0013] 第三步，根据第二步设定的壁面滑移速度边界条件，在第一步中的计算域网格上求解NS方程，计算超疏水表面曲面边界层的流场，并插值到曲面几何体表面沿法向的正交化网格下，获得各处的壁面切向速度分量剖面。

[0014] 第四步，基于线性稳定性理论，分析第三步中超疏水表面曲面边界层流场的流动稳定性，获得各频率扰动波幅值沿流向的增长率，方法如下：

[0015] 从不可压缩流动的无量纲NS方程出发，将流动的瞬时量分解成基本流与扰动量之和，即 其中 为瞬时量矢量， 为基本流矢量，φ'T
＝(v'1,v'2,v'3,p') 为扰动量矢量，v1,v2,v3分别为流向、法向和周向的速度分量，p为压力，上标“‑”代表基本流，上角标“'”代表扰动量。通过适用于水介质的流向和周向带曲率的曲面边界层流动的线性稳定性分析，得到扰动的线性稳定性方程，即：

[0016]

[0017] 式中，代表将扰动φ'写成行进波形式后的特征函数矢量，即其中t为时间，q1,q2和q3分别为正交曲线坐标系
下当地的流向、法向和周向空间坐标，ω为圆频率，α为流向波数，β为周向波数，c.c.为复共轭；在空间模式中，α为复数，表示为α＝αr+iαi,其中，‑αi为扰动的空间增长率，αr为实数形式的流向波数； 为系数矩阵，包括基本流、壁面的流向曲率drx和周向曲率drz、扰动圆频率、波数和雷诺数参数；

[0018] 扰动在超疏水表面曲面边界层传播时，其流向速度分量v'1和周向速度分量v'3在壁面处适用第二步中列出的滑移速度边界条件，线性稳定性方程(2)的在q2＝0处的壁面边界条件为：

[0019]

[0020] 基于第三步中得到的超疏水表面曲面边界层流场，在每一个流向位置处求解线性稳定性方程(2)，得到超疏水表面边界层内各频率扰动波幅值沿流向的空间增长率‑αi。

[0021] 第五步，使用第四步中获得的不同频率扰动波幅值的增长率，计算其最大增长倍数，预报超疏水表面曲面边界层的自然转捩位置，方法如下：

[0022] 对第四步所获得的波幅值沿流向的空间增长率‑αi沿流向积分，得到某一流向位置q1处扰动波幅值的放大因子N，即：

[0023]

[0024] 式中，f为某一扰动的频率大小，q10为积分的起始位置，取该频率扰动开始失稳的流向位置。

[0025] 计算不同频率的扰动波在每一个流向位置处的最大增长倍数 其中Nmax为所有频率扰动在当地q1处N值的最大值，即满足

[0026]

[0027] 当最大增长倍数 达到转捩判据给出的临界值 时，认为层流转捩为湍流，此时 所对应的流向位置q1就是发生转捩的位置。

[0028] 进一步地，第一步中的近壁区垂直网格划分用于计算每一个位置处壁面法向的速度梯度。

[0029] 进一步地，第五步中的转捩判据取

[0030] 进一步地，将第四步中的流向曲率和周向曲率均设定为零，预报超疏水表面平面边界层转捩位置。

[0031] 与现有技术相比，本发明至少具有以下有益效果：

[0032] 本发明提供了一种模拟曲面壁的超疏水表面的方法，可以设定复杂曲面几何体表面的滑移速度方向，实现壁面滑移速度边界条件。

[0033] 本发明提出了一种计算超疏水表面曲面边界层流场的方法。

[0034] 本发明建立了超疏水表面曲面边界层流动的线性稳定性分析方法。

[0035] 本发明建立的转捩位置预报方法，不局限于预报超疏水表面曲面边界层，而且还可以预报超疏水表面平面边界层和普通表面曲面边界层的转捩位置。

[0036] 本发明预报的超疏水表面曲面边界层转捩位置，为进一步发展水下曲面航行器的减阻降噪等技术提供了理论基础和支撑。

[0037] 综上所述，本发明提供了一种可行且合理的预报超疏水表面曲面边界层转捩位置的方法。

[0038] 本发明的其他特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

[0039] 附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

[0040] 图1是根据本发明的一种预报超疏水表面曲面边界层转捩位置方法的流程图；

[0041] 图2是根据本发明的一种实施方式的超疏水表面SUBOFF模型近壁区网格示意图；

[0042] 图3是根据本发明的一种实施方式的超疏水表面SUBOFF模型壁面处的滑移速度边界条件示意图；

[0043] 图4是根据本发明的一种实施方式的超疏水表面SUBOFF模型壁面处的滑移速度；

[0044] 图5是根据本发明的一种实施方式的超疏水表面SUBOFF模型距离前缘0.1m处壁面切向速度分量剖面；

[0045] 图6是根据本发明的一种实施方式的超疏水表面SUBOFF模型边界层内扰动幅值的增长率分布等值线图；

[0046] 图7是根据本发明的一种实施方式的超疏水表面SUBOFF模型边界层内扰动幅值的N值分布等值线图；

[0047] 图8是根据本发明的一种实施方式的超疏水表面SUBOFF模型边界层内扰动幅值的最大增长倍数曲线。

具体实施方式

[0048] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0049] 下面将结合附图来说明本发明所提供的一种预报超疏水表面曲面边界层转捩位置的方法的应用实例，如图1所示，包括以下步骤：

[0050] 本实施例模型为无附体的SUBOFF模型(Groves N.,Huang  T.,and Chang M..Geometric characteristics of DARPA Suboff models.DTRC/SHD‑1298‑01.David Taylor Research Center–Ship Hydromechanics Department,Department of the Navy.1989)，是研究水下曲面航行器绕流的标准模型之一，轴对称的几何结构主要由钝型的首部、圆柱形的中间体组成；工况为零攻角，航行速度10m/s。

[0051] 1、第一步，划分计算域网格，其中在复杂曲面几何体的壁面附近，绘制一层或者多层垂直曲面的网格。远离壁面处的网格，则没有这一特殊要求，可以根据具体情况进行划分。

[0052] 首先设定超疏水表面SUBOFF模型的径向计算域为5倍模型长度大小，然后使用FORTRAN编写程序划分计算域网格，其中近壁区的网格需要垂直曲面，如图2所示。这样的近壁区垂直网格是为了在第二步中方便计算每一个位置处的壁面法向的速度梯度，实现超疏水表面曲面滑移速度边界条件。

[0053] 2、第二步，根据第一步划分的计算域网格，模拟曲面壁的超疏水表面，实现壁面滑移速度边界条件。

[0054] 首先使用第一步中超疏水表面SUBOFF模型壁面网格点的数据，通过FORTRAN编写程序计算得到SUBOFF模型任意位置表面的切线方向，设定这个切线方向为超疏水表面当地的滑移速度方向。然后再使用壁面网格单元的数据，可以获得单元中心距离壁面的垂直距离，使用ANSYS Fluent软件中的用户自定义函数在每一个垂直曲面的近壁区网格单元内计算当地的壁面法向的速度梯度。最后使用1823年Navier提出的线性滑移速度边界条件假设(该假设是已有的公共知识，具体可见参考文献(Navier,C.L.M.H..Mémoire sur les lois du mouvement des fluids.Mém.Acad.R.Sci.Inst.6,389–440,1823))，如(1)式所示，实现如图3所示的超疏水表面曲面几何体的壁面滑移速度边界条件：

[0055]* *

[0056] 式中 为超疏水表面曲面几何体沿壁面切向的滑移速度，U为切向速度，y为壁*面法向，λ为滑移长度，上标“*”表示有量纲的物理量。本应用实例中的滑移长度取自实验中的测量值，且为常数，大小为21微米(Ou J.,Perot B.,and Rothstein J P..Laminar drag reduction in microchannels using ultrahydrophobic surfaces.Phys.Fluids 
16,4635,2004)。

[0057] 3、第三步，根据第二步设定的壁面滑移速度边界条件，在第一步中的计算域网格上求解NS方程，计算超疏水表面曲面边界层的流场，获得各处的壁面切向速度分量的剖面。

[0058] 首先，将第二步中的超疏水表面SUBOFF模型壁面滑移速度边界条件在ANSYS Fluent软件中完成配置。然后，使用ANSYS Fluent软件中的NS方程求解器在第一步中的计算域网格上求解NS方程，得到超疏水表面曲面边界层的流场，其中超疏水表面SUBOFF模型壁面处的滑移速度如图4中黑色实线所示。为了验证超疏水表面曲面滑移速度边界条件得到实现，图4中黑色圆点画出了使用滑移长度与当地壁面法向的速度梯度乘积所得的滑移速度，可以看出，两者完全一致。最后，将边界层的流场插值到曲面几何体表面沿法向的正交化网格下，最终获得各处的壁面切向速度分量剖面，图5为距离模型前缘0.1m处壁面切向速度分量剖面。

[0059] 4、第四步，基于线性稳定性理论，分析第三步中超疏水表面曲面边界层流场的流动稳定性，获得各频率扰动波幅值沿流向的增长率。

[0060] 首先，从不可压缩流动的无量纲NS方程出发，将流动的瞬时量分解成基本流与扰T动量之和，即 其中矢量φ＝(v1,v2,v3,p)为瞬时量矢量， 为基本
T
流矢量，φ'＝(v'1,v'2,v'3,p') 为扰动量矢量，v1,v2,v3分别为流向、法向和周向的速度分量，p为压力，上标“‑”代表基本流，上角标“'”代表扰动量。然后使用Liu等(Liu,J.,Chu,X.,and Zhang,Y..Numerical investigation of natural transitions of bow boundary layers over underwater axisymmetric bodies.Phys.Fluids 33,074101,2021)建立的适用于水介质的流向和周向带曲率的曲面边界层流动的线性稳定性分析方法，可以得到扰动的线性稳定性方程，即：

[0061]

[0062] 式 中 是 将 扰 动 φ '写 成 行 进 波形 式 后 的 特 征函 数 矢 量 ，即其中t为时间，q1,q2和q3分别为流向、法向和周向空间坐标，ω为圆频率，α为流向波数，β为周向波数，c.c.为复共轭，在空间模式中，α为复数，即α＝αr+iαi,其中‑αi为扰动的空间增长率； 为系数矩阵，其中包含基本流、壁面的流向曲率drx和周向曲率drz、扰动圆频率、波数和雷诺数等参数。

[0063] 然后要注意的是，扰动在超疏水表面曲面边界层传播时，其流向速度分量v'1和周向速度分量v'3在壁面处也要使用滑移速度边界条件，其中滑移长度的大小也为21微米，那么线性稳定性方程(2)的在q2＝0处的壁面边界条件为：

[0064]

[0065] 最后，基于第三步中得到的超疏水表面曲面边界层流场，采用数值方法在每一个流向位置处求解线性稳定性方程(2)，得到超疏水表面SUBOFF模型边界层内各频率扰动波幅值沿流向的有量纲增长率分布等值线图，如图6所示。其中数值方法为本领域公知的手段，具体见《流动稳定性》，周恒和赵耕夫，国防工业出版社。

[0066] 5、第五步，使用第四步中获得的不同频率扰动波幅值的增长率，计算其最大增长倍数，预报超疏水表面曲面边界层的自然转捩位置。

[0067] 首先，认为来流湍流度很低的边界层内不同频率的扰动波可以持续增长到湍流转捩点，那么对其沿流向积分第四步中获得的增长率‑αi，就可以得到扰动波幅值的放大因子N值分布等值线图，如图7所示，N的大小由(4)式给出：

[0068]

[0069] 式中，f为某一扰动的频率大小，q10为扰动波失稳的位置。

[0070] 然后，对不同频率下的放大因子N值取包络，这样就可以得到不同频率的扰动波在每一个流向位置处的最大放大因子Nmax，即：

[0071]

[0072] 最后，可以得到超疏水表面SUBOFF模型边界层扰动幅值的最大增长倍数eNmax曲Nmax Ntr线，如图8所示。当最大增长倍数e 首次等于转捩判据给定的e 时，认为层流转捩为湍流，从图中可以看出，此时所对应的流向位置为0.17m，那么超疏水表面SUBOFF模型边界层转捩位置则为0.17m。为了清晰地看出超疏水表面对SUBOFF模型边界层转捩位置的影响，我们将第二步和第四步中的滑移长度均设定为零，则得到普通表面SUBOFF模型边界层转捩位置为0.12m。转捩位置结果显示本发明可以很好地考虑超疏水表面对SUBOFF模型边界层转捩位置的影响。

[0073] 本发明提出的超疏水表面曲面边界层转捩位置的预报方法，不仅为回转体/水翼等水下航行体曲面边界层转捩的准确预报提供了一种有效的方法，而且还为超疏水表面曲面滑移速度边界条件的实现提供了一种可行的建模方式。

[0074] 本领域的技术人员应当理解，上述实施方式仅仅是为了清楚地说明本发明，而并非是对本发明的范围进行限定。对于所属领域的技术人员而言，在上述发明的基础上还可以做出其它变化或变型，并且这些变化或变型仍处于本发明的范围内。

[0075] 本发明未详细说明部分为本领域技术人员公知技术。