低雷诺数翼型层流分离泡长度的快速确定方法与流程

本发明涉及低雷诺数工具装置技术领域，具体涉及一种低雷诺数翼型层流分离泡长度的快速确定方法。

背景技术：

在新能源领域对低雷诺数翼型的实际使用中，层流分离泡的出现实际上增厚了翼型边界层，导致翼型的型阻增加等不良效应。学术界通常使用计算流体力学方法(cfd)来预测低雷诺数翼型上的层流分离泡尺寸，但是要想计算准确，需要花费相当长的网格绘制时间以及流场计算时间，且需要正确选择湍流模型来保证结果的准确度。

技术实现要素：

(一)要解决的技术问题

本发明要解决的技术问题是：如何设计一种低雷诺数翼型层流分离泡长度的快速确定方法。

(二)技术方案

为了解决上述技术问题，本发明提供了一种低雷诺数翼型层流分离泡长度的快速确定方法，包括以下步骤：

第一步、使用薄翼型理论预测翼型表面的无粘流动速度

1.3.确定弯度和厚度方程

设翼型的弦线是ab，其中a是前缘点，b是尾缘点，ab的中点o被作为原点，x轴的方向被定义为沿着弦线，以o为圆心，ab为直径作圆aob。直线l1和l2分别与圆相切于a点和b点，翼型完全处于直线l1和l2之间，直线pp’与x轴正交，与翼型的上下表面分别相交于pu和pl两点，pu和pl的偏心角分别定义为pob＝θ和p’ob＝-θ；

设弦长为c，翼型的形状用以下方程表示：

f(θ)表示以θ为自变量的函数，此处θ的变化区间是0到π，因此翼型的上下表面的曲线表示为：

yu＝c·f(θ),yl＝c·f(-θ)

弯度方程定义为：

厚度方程定义为：

1.4.预测二阶精度翼型表面无粘流动速度

考虑一个置于不可压直匀流中的翼型，以u表示直匀流的速度，q表示翼型表面的当地流动速度，为无粘流动速度，α表示翼型的攻角；

对于表面当地流动速度q的一阶精度近似解q1表示为：

±分别表示翼型的上、下表面，u表示自由流速度，与厚度相关的水平速度小扰动ut1表示为：

与弯度和攻角相关的水平速度小扰动uc1表示为：

yt′和y′c分别表示厚度方程和弯度方程的一阶导数；

rle表示翼型的前缘半径，表示q1的修正值。

在二阶近似解中，二阶弯度和厚度方程分别定义为：

相应的，对应于二阶弯度和厚度的二阶近似水平速度小扰动分别表示为：

因此，表面速度的二阶精度近似解q2表示为：

改进后的修正方程为：

λ表示翼型弯度线的初始角，表示q2的修正值；

第二步、使用有粘边界层方法求解层流分离泡尺寸

2.1预测层流分离点

thwaites方法的线性方程为：

q表示翼型表面的无粘流动速度，通过上一步骤中的计算公式得到，θ代表边界层的动量厚度，ν表示空气的运动粘度，a、b为常数，thwaites方法的线性方程被积分为：

下标i代表初始条件，在从前缘点开始计算时为0，一旦θ从以上方程中求解得到，层流分离点就从thwaites的层流分离参数m中求解得到，m的方程如下：

2.2预测层流-湍流转捩点

根据当地边界层的动量厚度，使用michel经验公式来计算层流-湍流转捩点的位置：

当基于动量厚度的雷诺数满足方程(1)时，转捩发生，这里，l是翼型表面当地位置与驻点间的曲线长度，称为特征长度；

reθ是一个函数，rex是另一个函数，而reθtr、rextr分别是一个值，当这两个函数同时满足式(1)时，reθ＝reθtr，rex＝rextr，下角标tr表示转捩的位置，当reθ和rex同时满足以上方程时，得到的层流-湍流转捩点的位置，是当公式(1)成立时，通过得到的rextr的值求出的l值；

2.3预测湍流再附着点

湍流再附着点x′使用horton模型来求解得到，改进的horton模型表示为：

λr是湍流再附着参数，根据湍流再附着参数计算得到湍流再附着点，两等号之间的r也表示下标，θ′指的是湍流边界层中的动量厚度；

二维湍流边界层中，thwaites方法用来求解湍流边界层中的动量厚度θ′：

xtr是之前求解出来的转捩点的位置，初始条件为x′＝xtr；

根据该初始条件，计算得到const，然后利用const计算得到动量厚度θ′；根据得到的层流分离泡的湍流再附着点x′和层流分离点x的坐标，得到的两者之差为层流分离泡的长度。

优选地，常数a等于0.45。

优选地，常数b等于0.6。

优选地，在混合边界层中，如果层流分离发生在层流-湍流转捩之前，并且湍流再附着也随之发生，此时满足层流分离泡形成的条件，在有些情况下，层流分离泡不能形成，因为湍流没有再附着于翼型表面，这种情况包括翼型在低攻角下工作，或者在低雷诺数下工作。

优选地，当m达到0.082时层流开始分离，或者，当m达到0.09时层流开始分离。

(三)有益效果

本发明通过进行无粘薄翼型理论计算，输出包括翼型表面速度在内的流场信息，再通过相关边界层模型的计算，快速得到有粘边界层内的流场信息，有效输出层流分离泡的分离、转捩和在附着状态。该方法采用fortran90自编译实现，用以研究低雷诺数下主导翼型气动性能的层流分离泡。

附图说明

图1是薄翼型理论中的翼型几何形状描述；

图2是直匀流中等攻角下的翼型示意图；

图3是零度攻角下naca0012翼型的表面流动速度无粘预测图；

图4是雷诺数100000时e387翼型的表面流动速度无粘预测及实验数据对比图；

图5是本发明的方法流程图；

图6是雷诺数100000(左)和200000(右)下，本发明方法预测的层流分离泡位置，及其与cfd预测结果和实验结果的对比。

具体实施方式

为使本发明的目的、内容、和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。

为快速预测低雷诺数翼型的层流分离泡尺寸，本发明提出一种基于薄翼型和边界层理论的有粘无粘交互的工程分析方法，该方法能够快速确定低雷诺数翼型层流分离泡长度，包括以下步骤：

2.使用薄翼型理论预测翼型表面的无粘流动速度

2.1.确定弯度和厚度方程

考虑一个典型的薄翼型，如图1所示。设翼型的弦线是ab，其中a是前缘点，b是尾缘点。ab的中点o被作为原点，x轴的方向被定义为沿着弦线。以o为圆心，ab为直径作圆aob。直线l1和l2分别与圆相切于a点和b点。翼型完全处于直线l1和l2之间。直线pp’与x轴正交，与翼型的上下表面分别相交于pu和pl两点。pu和pl的偏心角分别定义为pob＝θ和p’ob＝-θ。

设弦长为c，根据上文描述，翼型的形状可以用以下方程表示：

f(θ)表示以θ为自变量的函数，翼型不同，f(θ)不同，此处θ的变化区间是0到π。因此翼型的上下表面的曲线可以表示为：

yu＝c·f(θ),yl＝c·f(-θ)

弯度方程可以定义为：

从上式可知，弯度曲线函数为偶函数。

厚度方程可以定义为：

2.2.预测二阶精度翼型表面无粘流动速度

考虑一个置于不可压直匀流中的翼型，如图2所示。图中的u表示直匀流的速度，q表示翼型表面的当地流动速度(为无粘流动速度)，α表示翼型的攻角。

对于表面当地流动速度q的一阶精度近似解q1可表示为：

这里的±符号分别表示翼型的上、下表面，u表示自由流速度，下标1表示一阶精度的近似解。与厚度相关的水平速度小扰动ut1可以被表示为：

与弯度和攻角相关的水平速度小扰动uc1可以被表示为：

方程中，yt′和y′c分别表示厚度方程和弯度方程的一阶导数。

从以上方程得知，薄翼型理论的一阶近似解会导致前缘点a处出现奇点。lighthill学者提出了一种通用修正方法来解决这个问题，方程如下：

rle表示翼型的前缘半径，表示q1的修正值。

在精度更高的二阶近似解中，二阶弯度和厚度方程被分别定义为：

相应的，对应于二阶弯度和厚度的二阶近似水平速度小扰动分别可以表示为：

因此，表面速度的二阶精度近似解q2可以表示为：

相似于一阶近似解，lighthill为二阶近似解也存在的前缘奇点问题提出了修正。这个修正后来被vandyke加以改进，给出了翼型前缘部分表面速度更精确的解。改进后的修正方程为：

这里λ表示翼型弯度线的初始角，表示q2的修正值。

无粘流动速度的预测效果：

前面详细介绍了本发明的计算方法中的无粘部分，该部分可用来快速计算翼型表面的无粘流场的速度分布。这里用两个算例来验证无粘速度的预测效果。

使用本发明中的薄翼型理论的无粘预测方法，图三给出了naca0012翼型在0度攻角下的表面流动速度分布。

总体来看二阶近似解比一阶近似解更加精确，因为考虑了二阶水平速度小扰动的缘故。由于考虑了前缘半径的修正，附带vandyke修正的二阶近似解给出了比二阶近似解更为准确的预测，也是四种方法中最准确的预测。因此在下一个算例中，选用附带vandyke修正的二阶近似解来预测e387翼型的表面速度分布。

如图4所示，在雷诺数100000时，除了翼型吸力面0.6至0.9倍弦长区间的层流分离泡区域，无粘薄翼型理论方法预测的e387翼型表面速度分布与nasa的实验数据很好的吻合。

无粘的假设无法预测出包括层流分离泡尺寸在内的边界层粘性特性，这将在本发明方法的有粘部分详细求解。

3.使用有粘边界层方法求解层流分离泡尺寸

预测层流分离泡的尺寸，必须解决三个关键位置参数：层流分离点，层流湍流转捩点和湍流再附着点。而在计算其长度之前，必须确认层流分离泡存在。在数值计算中，仅当层流分离点发生在层流湍流转捩点之前，层流分离泡才会存在。相反，如果转捩过程发生在层流分离之前，翼型上发生自然转捩，则不会出现层流分离泡。

3.1.预测层流分离点

thwaites方法被广泛采用来求解层流分离点x的位置，其近似于一种线性方程如下所示：

这里q表示翼型表面的无粘流动速度，通过上一步骤中的计算公式得到。θ代表边界层的动量厚度，ν表示空气的运动粘度。thwaites建议方程中的常数a等于0.45，b等于0.6，这些都基于实验数据得到，因此方程可以被积分为：

积分时，总会从初始条件积分到最终条件，这里下标i代表初始条件，在从前缘点开始计算时为0。一旦θ从以上方程中求解得到，层流分离点就可以从thwaites的层流分离参数m中求解得到，m的方程如下：

thwaites建议当m达到0.082时层流开始分离，但这个参数后来被curle修改为0.09。

3.2.预测层流-湍流转捩点

根据当地边界层的动量厚度，使用michel经验公式来计算层流-湍流转捩点的位置：

当基于动量厚度的雷诺数满足以上方程时，转捩发生。这里,l是翼型表面当地位置与驻点间的曲线长度，称为特征长度；

reθ是一个函数，rex是另一个函数，而reθtr、rextr分别是一个值，当这两个函数同时满足上式(1)时，reθ＝reθtr，rex＝rextr。下角标tr表示转捩的位置。当reθ和rex同时满足以上方程时，得到的层流-湍流转捩点的位置，是当公式(1)成立时，通过得到的rextr的值求出的l值。需要注意的是，该方程仅能用于层流的情况。

在混合边界层中，如果层流分离发生在层流-湍流转捩之前，并且湍流再附着也随之发生，此时满足层流分离泡形成的条件。在有些情况下，层流分离泡不能形成，因为湍流没有再附着于翼型表面。这种情况包括翼型在低攻角下工作，或者在低雷诺数下工作。

3.3.预测湍流再附着点

湍流再附着点x′可以使用horton模型来求解得到。改进的horton模型可以表示为：

这里λr是湍流再附着参数，该参数由大量实验数据总结归纳而成，根据湍流再附着参数可以计算得到湍流再附着点，两等号之间的r也表示下标。这里的θ′指的是湍流边界层中的动量厚度，而非层流边界层，因此之前层流边界层中求解θ的方法这里不再适用。

二维湍流边界层中，thwaites方法可以用来求解湍流边界层中的动量厚度θ′：

这里xtr是之前求解出来的转捩点的位置，在本算例中初始条件即为x′＝xtr；

根据该初始条件，计算得到const，然后可利用const计算得到动量厚度θ′。

本发明的预测方法，通过结合上面阐述的翼型表面的无粘速度预测和本步骤阐述的有粘边界层方法，完成对层流分离泡的尺寸预测。预测方法的整体步骤如图5所示。通过图中的步骤，得到层流分离泡的湍流再附着点x′和层流分离点x的坐标，两者之差可以近似为层流分离泡的长度。

算例结果和讨论：

本发明选用e387翼型作为低雷诺数翼型的代表，使用本发明中介绍的计算方法来计算其层流分离泡的位置和尺寸。为了对比计算结果，选用美国nasa低湍流度风洞在雷诺数100000和200000的实验数据作比较验证。cfd计算也被应用于相同算例，用以进一步验证和对比。

在本发明的计算方法中，使用图5中的层流分离位置和层流-湍流转捩位置，来求解湍流再附着点的位置，以此预测层流分离泡的尺寸。而在cfd计算中，使用流场中的流动分离流线以及翼型表面的几何坐标来圈定层流分离泡的尺寸。层流分离点的坐标和湍流再附着点的坐标，即由此得到。

由本发明计算方法计算得到的层流分离泡的详细信息如图6所示，图中总结了雷诺数100000和200000条件下层流分离泡的层流分离点和湍流再附着点的位置。

预测结果的合理性分析：

在图6展示的两个雷诺数下，通过与实验数据相比较发现，本发明所使用的计算方法，均对e387翼型上的层流分离点的位置提供了很好的预测结果。同时对于湍流再附着点的预测，本发明方法的预测结果与实验数据趋势相同，但位置较实验数据略偏下游。在雷诺数为100000时，湍流再附着点的预测误差约10％。在雷诺数为200000时，湍流再附着点的预测结果，在较低攻角下略偏下游约8％。但当攻角大于三度时，本发明方法预测结果与实验结果吻合。

湍流再附着点的预测误差来源于再附着点的判定标准。该判定标准由经验公式得出。判定标准的提出者horton指出，湍流中的大速度梯度可以快速改变湍流再附着区域，这使得实验测量非常困难。因此，准确的湍流再附着位置非常难以确定。判定标准的改进者也指出，另一个湍流再附着位置预测误差的原因是不能准确的预测湍流边界层中的转捩和动量厚度。但值得一提的是，作为一种快速响应的计算工具，本发明计算方法能够提供误差10％以内的层流分离泡位置快速预测。误差必然存在是因为该计算方法是一种通用方法，而非为了某一种特定翼型而专门研发的计算方法。

为了进一步验证本发明的计算方法，同时使用cfd方法计算层流分离泡的位置并总结于图6中。粗略看来cfd计算结果与实验结果的吻合度总体要高于本发明计算方法的吻合度，这是因为cfd计算中选择了较为合适的湍流模型，而且花费了大量的网格生成和计算时间。

在雷诺数为100000时，与实验相比，本发明的计算方法对层流分离点的预测比cfd方法预测更为准确。cfd的计算结果在攻角大于6度时开始快速向前缘方向移动。这种过早的层流分离预测直接导致了更大的层流分离泡的形成，最终造成cfd预测的升力系数偏小，阻力系数偏大。

在雷诺数为200000时，如图6的右图所示，本发明的计算结果和cfd计算结果在与实验结果对比时均提供了合理的层流分离点预测。在较低攻角下，本发明的计算结果稍偏下游。两种计算结果随攻角增加的趋势相同。

图6中两种方法预测的层流分离泡长度被总结在表1中，与实验对比的计算误差与计算时间对比也包含在表一中。在两个雷诺数下，除了7度攻角，分离泡长度在可接受范围内被有效预测出。总的来看，在雷诺数为100000时本发明的预测方法更为精确。

表1

预测结果的准确性分析：

如表1所示，两种方法对于e387翼型在七度攻角下的层流分离泡长度的预测结果，与实验结果相比较，均达到20％以上的误差。对于这种现象一个可能的解释是层流分离泡的敏感性。因为七度攻角已经非常接近于该翼型的自然转捩攻角，因此此时层流分离泡变得非常敏感且不稳定。该翼型在七度攻角左右层流分离泡出线的敏感性问题，被众多研究者所发现。比如lin和pauley发现，风洞实验的自由流湍流会推迟层流的分离，因此使得层流分离泡更短。他们也发现，cfd计算会过高的估计层流分离泡在七度攻角时的长度。

另外，本发明中的计算方法的不准确性主要来自于湍流再附着区域预测的不准确。对于湍流再附着预测不准确的一个合理解释是，经典的再附着判定标准并没有将湍流度考虑在内，而湍流度本身又对层流分离泡的形成起着至关重要的作用，尤其是湍流再附着点的位置。况且，预测湍流再附着点位置的经验方程，是基于其他翼型的实验结果，并没有将e387翼型的自然转捩攻角考虑在内。因此两种不同的计算方法都过长地估算了层流分类泡的长度。

预测结果的高效性分析：

两种计算方法对层流分离泡长度的计算误差和计算时长也被总结于表一中。与实验数据相比，两种方法都给出了可以接受的长度预测，cfd方法的预测结果稍微更精确一些。但是，本发明的计算方法再保证计算误差可接受的情况下，可以提供及其快速的响应计算时间。在同时使用英特尔i76700处理器的情况下，本发明计算方法仅需单核单线程就可以在五秒内完成单次运算，即使cfd方法采用了四核八线程并行运算，本发明方法依然比cfd方法快了百倍以上，这其中还不包含网格生成和网格独立性调查的时间。这就使得本发明的计算方法成为一种快速便捷响应的估算低雷诺数翼型层流分离泡长度的计算工具。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。