一种壁面加热/冷却的水下平板边界层转捩位置的预测方法

1.本发明涉及水动力学研究技术领域，具体涉及一种壁面加热/冷却的水下平板边界层转捩位置的预测方法。


背景技术：

2.流体在物体表面流动过程中，表面边界层内的流动存在层流和湍流两种流动状态。层流状态和湍流状态的性质有着显著的区别，如阻力、噪声、热流等这些对实际工程问题有重要影响的物理量。流动状态从层流转变为湍流的过程称为转捩。控制边界层的转捩对航海、航空、航天等工程领域中的实际问题有着重要的影响。
3.前人对平板边界层有较多的的实验、计算以及理论研究成果，所以可以从平板边界层入手，研究边界层转捩控制。控制平板边界层转捩的方法有很多，加热/冷却壁面是其中一种重要的控制手段，而预测壁面加热/冷却的水下平板边界层转捩位置对这一控制手段的具体实施是非常关键的。
4.目前对于壁面加热/冷却的水下平板边界层的转捩预测，现有可能的方法有三种。一是采用实验办法，但是对实验环境和设备的要求较高，实验花费成本较高；二是采用直接数值模拟方法，但是计算量太大，当前的计算水平无法用于实际问题；三是采用经验或半经验的方法，其中半经验的en方法是目前看来可用的一种转捩预测方法。


技术实现要素：

5.本发明的目的在于提供一种壁面加热/冷却的水下平板边界层转捩位置的预测方法。本发明的技术方案如下：
6.一种壁面加热/冷却的水下平板边界层转捩位置的预测方法，包括以下步骤：
7.步骤一，确定计算工况，包括确定来流温度、来流速度以及需要计算的水下平板壁面加热/冷却时的壁面温度；
8.步骤二，针对壁面加热/冷却的水下平板边界层，考虑能量方程和温度变化，考虑温度和速度的耦合，且考虑粘性系数和热传导系数随温度的变化，而密度和比热则视为常数，建立层流基本流二维定常有量纲的边界层方程，再进一步推导出考虑粘度、热传导系数随温度变化的blasius方程，从而获得步骤一中所确定工况下的壁面加热/冷却的水下平板边界层的层流基本流流场，方法如下：
9.1)建立层流基本流二维定常有量纲的边界层方程(1)：
[0010][0011]
公式中，带“*”的量为有量纲量，不带“*”的量为无量纲量，有量纲的边界层方程(1)中第一式为连续性方程，第二式为流向动量方程，第三式为法向动量方程，第四式为能量方程；x
*
、y
*
分别为直角坐标系中流向坐标和法向坐标；u
*
、v
*
、t
*
、p
*
分别为流向速度、法向速度、温度、压强；ρ
*
、c
*
、μ
*
、k
*
分别为密度、比热、动力学粘性系数、热传导系数；
[0012]
2)引入变量f、g和坐标η，且变量f和g均为坐标η的函数，f对η一阶导数为无量纲的速度，为无量纲的温度，为壁面处的温度变量g；带下角标“e”的量表示来流中的量，带下角标“w”的量表示壁面处的量，为来流速度，为来流中的动力学粘性系数，为壁面温度，为来流温度；将f和g的表达式代入有量纲的边界层方程(1)，考虑流向压力梯度为零的情况，即得到考虑粘度、热传导系数随温度变化的blasius方程(2)：
[0013][0014]
考虑粘度、热传导系数随温度变化的blasius方程(2)的系数分别为
[0015]
设壁面处使用无滑移条件和等温条件，远离壁面处使用来流速度和来流温度条件，则考虑粘度、热传导系数随温度变化的blasius方程(2)的边界条件(3)：
[0016][0017]
3)使用打靶法和四阶的runge-kutta法，求解具有边界条件(3)的考虑粘度、热传导系数随温度变化的blasius方程，得到壁面加热/冷却的水下平板边界层的层流基本流流场；
[0018]
步骤三，基于线性稳定性理论，考虑能量方程和温度变化，考虑温度和速度的耦合，且考虑粘性系数和热传导系数随温度的变化，而密度和比热则视为常数，并且考虑温度的扰动，对壁面加热/冷却的水下平板边界层进行稳定性分析，得到各频率小扰动增长率沿
着流向的分布即增长率云图；
[0019]
步骤四，计算壁面加热/冷却的水下平板边界层小扰动沿流向最大增长倍数，方法如下：
[0020]
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数，由式(7)给出：
[0021][0022]
其中a表示小扰动在位置x处的幅值，a0表示中性点的小扰动幅值，x0为中性点对应的流向位置，n为幅值放大因子；
[0023]
2)通过从中性曲线开始，对幅值增长率沿水下平板壁面积分，获得各频率小扰动的幅值放大因子n沿流向坐标的分布，用一条曲线将不同频率波的n值沿流向坐标分布的最大值包络，包络线的各点取值代表各种频率扰动在此处能达到的最大幅值放大因子n
max
，获得n
max
曲线，最大增长倍数的大小为
[0024]
步骤五，预测壁面加热/冷却的水下平板边界层的转捩位置；
[0025]
1)认为只要某频率的小扰动波幅在沿流向发展的过程中，幅值达到初始幅值的倍时层流即发展为湍流，作为转捩判据的n
tr
值的大小由实验标定；在步骤四中获得的n
max
曲线上找到n
tr
对应流向位置，即为预测出来的壁温加热/冷却的水下平板的转捩位置。
[0026]
2.根据权利要求1所述的壁面加热/冷却的水下平板边界层转捩位置的预测方法，其特征在于，步骤五中，转捩判据取为n
tr
＝7。
[0027]
进一步地，步骤三的方法如下：
[0028]
1)建立无量纲的三维流动的线性扰动方程(4)：
[0029][0030]
其中，x、y、z分别为直角坐标系中的流向坐标、法向坐标、展向坐标，t为时间，φ'＝(u',v',w',t',p')
t
为扰动矢量，且φ'为关于t、x、y、z的函数，u'为流向速度扰动量、v'为法向速度扰动量、w'为展向速度扰动量、t'温度扰动量、p'为压强扰动量，γ'、a'、b'、c'、d'、v
11
'、v
22
'、v
33
'、v
12
'、v
13
'、v
23
'为系数矩阵，系数矩阵中包含基本流的信息，基本流通过步骤二中求解考虑粘度、热传导系数随温度变化的blasius方程(2)得到；系数矩阵也包含了无量纲参数雷诺数re、埃克特数ec、普朗特数pr；
[0031]
2)在边界层流向上使用近似平行流假设，得到扰动展开式(5)：
[0032]
φ'＝φ
·ei(αx iβ-ωt)
c.c.
ꢀꢀꢀꢀꢀꢀꢀ
(5)
[0033]
其中，为扰动量的特征函数，且φ仅为关于y的函数，u为流向速度扰动量的特征函数、为法向速度扰动量的特征函数、w为展向速度扰动量的特征函数、t温度扰动量的特征函数、p为压强扰动量的特征函数；ω为无量纲圆频率，对应的有量纲频率为f
*
，β为展向波数，c.c.表示复共轭；空间模式下α＝αr iαi为复数形式波数,αr为实数形式波数，-αi为空间增长率；
[0034]
3)将扰动展开式(5)代入线性扰动方程(4)，得到线性稳定性方程(6)：
[0035][0036]
其中，v
22
、b、d为系数矩阵，表达式为：
[0037]v22
＝v
22
'，
[0038]
b＝b' iαv
12
' iβv
23
'，
[0039]
d＝-iωγ' iαa' iβc' d'-α
2v11
'-β
2v33
'-αβv
13
'；
[0040]
4)系数矩阵v
22
、b、d中包含基本流的信息和参数雷诺数re、复波数α、展向波数β、圆频率ω、埃克特数ec、普朗特数pr；采用数值方法计算线性稳定性方程(6)，获得基本流不同流向站位的特征值，即扰动展开式(5)中的波数α、圆频率ω与基本流流场雷诺数re的参数关系；通过特征值的计算得到各层流基本流的中性曲线，并得到各频率小扰动增长率-αi沿着流向的分布，即增长率云图。
[0041]
进一步地，步骤四的方法如下：
[0042]
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数，由式(7)给出：
[0043][0044]
其中a表示小扰动在位置x处的幅值，a0表示中性点的小扰动幅值，x0为中性点对应的流向位置，n为幅值放大因子；
[0045]
2)通过从中性曲线开始，对幅值增长率沿水下平板壁面积分，获得各频率小扰动的幅值放大因子n沿流向坐标的分布，用一条曲线将不同频率波的n值沿流向坐标分布的最大值包络，包络线的各点取值代表各种频率扰动在此处能达到的最大幅值放大因子n
max
，从而获得n
max
曲线，最大增长倍数的大小为
[0046]
本发明与现有技术相比的有益效果：
[0047]
(1)本发明基于线性稳定性分析理论，建立了考虑壁面加热/冷却的水下平板边界层的转捩位置预测方法；
[0048]
(2)本发明可用于通过加热/冷却壁面来控制水下平板边界层的转捩位置。
附图说明
[0049]
图1为本发明步骤框图
[0050]
图2为所计算工况的速度剖面
[0051]
图3为所计算工况的温度剖面
[0052]
图4为所计算工况的中性曲线
[0053]
图5为加热壁面的增长率云图
[0054]
图6为加热壁面的n值分布云图
[0055]
图7为所计算工况的n
max
曲线
具体实施方式
[0056]
下面结合具体实施方式及附图对本发明进行详细说明。在下面的描述中，出于解
释而非限制性的目的，阐述了具体细节，以帮助全面地理解本发明。然而，对本领域技术人员来说显而易见的是，也可以在脱离了这些具体细节的其它实施例中实践本发明。
[0057]
在此需要说明的是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的设备结构和/或处理步骤，而省略了与本发明关系不大的其他细节。
[0058]
本发明实施例提供了一种壁面加热/冷却的水下平板边界层转捩位置的预测方法，如图1所示，包括下述步骤：
[0059]
步骤一，根据实际情况，确定计算工况，即确定来流温度、来流速度以及要计算的水下平板壁面加热/冷却时的壁面温度。
[0060]
该步骤中，取计算工况为：零攻角，来流速度来流温度平板壁面温度取为4.4℃(壁面冷却)、15.6℃(壁面温度等于来流温度)、32.2℃(壁面加热)。
[0061]
步骤二，计算步骤一中所确定工况下的水下平板的层流基本流，获得层流基本流流场。
[0062]
1)针对壁面加热/冷却的水下平板边界层，考虑能量方程和温度变化，考虑温度和速度的耦合，且考虑粘性系数和热传导系数随温度的变化，而密度和比热则视为常数。从二维定常有量纲的n-s方程出发，推导出有量纲的边界层方程，再进一步推导出blasius方程，获得层流基本流流场。方法如下：
[0063]
2)层流基本流满足二维定常有量纲n-s方程(n-s方程的形式为本领域公知技术，具体见《流体力学》，周光炯等，高等教育出版社)。设定两个条件，一是边界层厚度在量级上小于到平板前缘距离，二是边界层内的粘性力和惯性力的量级相同。对二维定常有量纲n-s方程进行量级分析，忽略高阶小量，得到如下的有量纲的边界层方程(1)(这一步推导过程为本领域公知技术，具体见《流体力学》，周光炯等，高等教育出版社)：
[0064][0065]
公式中，带“*”的量为有量纲量，不带“*”的量为无量纲量，方程(1)中第一式为连续性方程，第二式为流向动量方程，第三式为法向动量方程，第四式为能量方程。x
*
、y
*
分别为直角坐标系中流向坐标和法向坐标；u
*
、v
*
、t
*
、p
*
分别为流向速度、法向速度、温度、压强；ρ
*
、c
*
、μ
*
、k
*
分别为密度、比热、动力学粘性系数、热传导系数。
[0066]
3)引入变量f、g和坐标η，且变量f和g均为坐标η的函数，f对η一阶导
数为无量纲的速度，为无量纲的温度，为壁面处的温度变量g；带下角标“e”的量表示来流中的量，带下角标“w”的量表示壁面处的量，为来流速度，为来流中的动力学粘性系数，为壁面温度，为来流温度；将f和g的表达式代入有量纲的边界层方程(1)，考虑流向压力梯度为零的情况，即得到考虑粘度、热传导系数随温度变化的blasius方程(2)(这一步推导过程为本领域公知技术，具体见《高超声速和高温气体动力学》，john d.anderson jr等，航空工业出版社)：
[0067][0068]
方程(2)的系数分别为
[0069]
设壁面处使用无滑移条件和等温条件，远离壁面处使用来流速度和来流温度条件，则blasius方程(2)的边界条件为：
[0070][0071]
4)使用打靶法和四阶的runge-kutta法(这一方法为本领域公知技术，具体见《高超声速和高温气体动力学》，john d.anderson jr等，航空工业出版社)，求解具有边界条件(3)的blasius方程(2)，求解过程中使用适用于0.1mpa、0-100℃的水的动力学粘性系数μ
*
和热传导系数k
*
随温度t
*
的变化关系；再使用关系式和得到壁面加热/冷却的水下平板边界层的层流基本流流场。
[0072]
5)在该步骤中，适用于0.1mpa、0-100℃的水的动力学粘性系数μ
*
和热传导系数k
*
随温度t
*
的变化关系使用关系式(4)和(5)(具体见文献huber m l,perkins r a,laesecke a,et al.new international formulation for the viscosity of h2o[j].journal of physical&chemical reference data,2009,38(2):101-125.和huber m l,perkins r a,friend d g,et al.new international formulation for the thermal conductivity of h2o[j].journal of physical&chemical reference data,2012,41(3):387-440.)：
[0073]
[0074][0075]
6)该步骤中，层流基本流的边界层内需要保证足够的网格数量，以保证后续稳定性分析的需求。
[0076]
7)该步骤可得到所计算工况的基本流剖面，即速度剖面如图2所示和温度剖面如图3所示。
[0077]
步骤三，基于线性稳定性理论，对壁面加热/冷却的水下平板边界层进行稳定性分析，得出小扰动的增长率。
[0078]
1)针对壁面加热/冷却的水下平板边界层，考虑能量方程和温度变化，考虑温度和速度的耦合，且考虑粘性系数和热传导系数随温度的变化，而密度和比热则视为常数，此外，还要考虑温度的扰动。从无量纲的三维流动的n-s方程(n-s方程的形式为本领域公知技术，具体见《流体力学》，周光炯等，高等教育出版社)出发，使用小扰动假设推导出线性扰动方程，再使用近似平行流假设推导出线性稳定性方程，从而进行稳定性分析。方法如下：
[0079]
2)考虑小扰动的情况，将瞬时量表示为层流基本流和小扰动之和，即
[0080]
(us,vs,ws,ts,ps,μs,ks)
t
＝(u,v,w,t,p,μ,k)
t
(u',v',w',t',p',μ',k')
t
。
[0081]
其中，us、vs、ws、ts、ps、μs、ks为瞬时量，且ws为展向速度；u、v、w、t、p、μ、k为层流基本流的量，u'、v'、w'、t'、p'、μ'、k'为扰动量；由粘性系数和热传导系数随温度的变化关系，可得到粘性系数和热传导系数的扰动量随温度扰动量的变化关系式，表示为μ'＝τt'，k'＝σt'，其中τ和σ为基本流温度t的函数。
[0082]
3)瞬时量和基本流都满足无量纲的三维流动的n-s方程，用瞬时量的方程减去基本流的方程，再略去非线性的高阶小量，得到线性扰动方程(6)，用矢量形式表示为：
[0083][0084]
其中，x、y、z分别为直角坐标系中的流向坐标、法向坐标、展向坐标，t为时间，φ'＝(u',v',w',t',p')
t
为扰动矢量，γ'、a'、b'、c'、d'、v
11
'、v
22
'、v
33
'、v
12
'、v
13
'、v
23
'为系数矩阵，系数矩阵中包含基本流的信息，基本流通过步骤二中求解blasius方程(2)得到；系数矩阵也包含了无量纲参数雷诺数re、埃克特数ec、普朗特数pr。
[0085]
4)考虑到边界层在流向上变化缓慢，在流向上使用近似平行流假设，将扰动描述为行进波的形式，得到扰动展开式(7)：
[0086]
φ'(t,x,y,z)＝φ(y)
·ei(αx iβ-ωt)
c.c.
ꢀꢀꢀꢀꢀꢀꢀ
(7)
[0087]
其中，为扰动的特征函数；ω为无量纲圆频率，对应的有量纲频率为f
*
(单位是hz)，β为展向波数，c.c.表示复共轭。空间模式下α＝αr iαi为复数形式波数,αr为实数形式波数，-αi为空间增长率。将扰动展开式(7)代入线性扰动方程(6)，得到线性稳
定性方程(8)：
[0088][0089]
其中，v
22
、b、d为系数矩阵，表达式为：
[0090]v22
＝v
22
'，
[0091]
b＝b' iαv
12
' iβv
23
'，
[0092]
d＝-iωγ' iαa' iβc' d'-α
2v11
'-β
2v33
'-αβv
13
'。
[0093]
5)系数矩阵v
22
、b、d中包含基本流的信息和参数雷诺数re、复波数α、展向波数β、圆频率ω、埃克特数ec、普朗特数pr。采用数值方法计算线性稳定性方程(8)，获得基本流不同流向站位的特征值，即扰动展开式中的波数α、圆频率ω与基本流流场雷诺数re的参数关系。通过特征值的计算可以得到各层流基本流的中性曲线；也可以得到各频率小扰动增长率-αi沿着流向的分布，即增长率云图。
[0094]
6)该步骤中，稳定性方程特征值分析方法可采用本领域公知的手段，具体见《流动稳定性》，周恒等，国防工业出版社。
[0095]
7)该步骤中得到的有量纲的中性曲线如图4所示，时的增长率云图如图5所示。
[0096]
步骤四，计算壁面加热/冷却的水下平板边界层小扰动沿流向最大增长倍数。
[0097]
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数，由式(7)给出：
[0098][0099]
其中a表示小扰动在位置x处的幅值，a0表示中性点的小扰动幅值，x0为中性点对应的流向位置，n为幅值放大因子。
[0100]
2)该步骤中，通过从中性曲线开始，对幅值增长率沿水下平板壁面积分，获得各频率小扰动的幅值放大因子n沿流向坐标的分布，工况的n值分布云图如图6所示。
[0101]
3)该步骤中，用一条曲线将不同频率波的n值沿流向坐标分布的最大值包络，包络线代表各种频率扰动在此处能达到的最大幅值放大因子，获得n
max
曲线，所计算工况的n
max
曲线如图7所示。n
max
曲线表示出不同频率的小扰动可以达到的最大增长倍数，这个倍数的大小为
[0102]
步骤五，预测壁面加热/冷却的水下平板边界层的转捩位置。
[0103]
1)该步骤中，根据对边界层自然转捩的认识，可以认为只要某频率的小扰动波幅在沿流向发展的过程中，幅值达到初始幅值的倍时层流即发展为湍流。
[0104]
2)该步骤中，作为转捩判据的n
tr
值的大小由实验标定出来，转捩判据取为n
tr
＝7。
[0105]
3)该步骤中，在步骤四中获得的n
max
曲线上找到n
tr
对应流向位置，即为预测出来的壁温加热/冷却的水下平板的转捩位置。
[0106]
4)该步骤中，可得到每个工况的转捩位置以及转捩延迟距离，将结果列在表1中：
[0107]
表1不同壁温时的转捩预测结果(来流温度为15.6℃)
[0108][0109][0110]
表1中，是壁温为时的转捩位置，是壁温为15.6℃时的转捩距离，为转捩延迟距离，其表达式为
[0111]
如上针对一种实施例描述和/或示出的特征可以以相同或类似的方式在一个或更多个其它实施例中使用，和/或与其它实施例中的特征相结合或替代其它实施例中的特征使用。
[0112]
应该强调，术语“包括/包含”在本文使用时指特征、整件、步骤或组件的存在，但并不排除一个或更多个其它特征、整件、步骤、组件或其组合的存在或附加。
[0113]
本发明以上的装置和方法可以由硬件实现，也可以由硬件结合软件实现。本发明涉及这样的计算机可读程序，当该程序被逻辑部件所执行时，能够使该逻辑部件实现上文所述的装置或构成部件，或使该逻辑部件实现上文所述的各种方法或步骤。本发明还涉及用于存储以上程序的存储介质，如硬盘、磁盘、光盘、dvd、flash存储器等。
[0114]
这些实施例的许多特征和优点根据该详细描述是清楚的，因此所附权利要求旨在覆盖这些实施例的落入其真实精神和范围内的所有这些特征和优点。此外，由于本领域的技术人员容易想到很多修改和改变，因此不是要将本发明的实施例限于所例示和描述的精确结构和操作，而是可以涵盖落入其范围内的所有合适修改和等同物。
[0115]
本发明未详细说明部分为本领域技术人员公知技术。