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

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

Info

Publication number
CN114528665B
CN114528665B CN202210174746.9A CN202210174746A CN114528665B CN 114528665 B CN114528665 B CN 114528665B CN 202210174746 A CN202210174746 A CN 202210174746A CN 114528665 B CN114528665 B CN 114528665B
Authority
CN
China
Prior art keywords
temperature
equation
boundary layer
flow
disturbance
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210174746.9A
Other languages
English (en)
Other versions
CN114528665A (zh
Inventor
张永明
吕军波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN202210174746.9A priority Critical patent/CN114528665B/zh
Publication of CN114528665A publication Critical patent/CN114528665A/zh
Application granted granted Critical
Publication of CN114528665B publication Critical patent/CN114528665B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • G06F30/17Mechanical parametric or variational design
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08Thermal analysis or thermal optimisation
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

本发明涉及一种壁面加热/冷却的水下平板边界层转捩位置的预测方法,包括以下步骤:确定计算工况;针对壁面加热/冷却的水下平板边界层,建立层流基本流二维定常有量纲的边界层方程,再进一步推导出考虑粘度、热传导系数随温度变化的Blasius方程,从而获得步骤一中所确定工况下的壁面加热/冷却的水下平板边界层的层流基本流流场;对壁面加热/冷却的水下平板边界层进行稳定性分析,得到各频率小扰动增长率沿着流向的分布即增长率云图;计算壁面加热/冷却的水下平板边界层小扰动沿流向最大增长倍数;预测壁面加热/冷却的水下平板边界层的转捩位置。

Description

一种壁面加热/冷却的水下平板边界层转捩位置的预测方法
技术领域
本发明涉及水动力学研究技术领域,具体涉及一种壁面加热/冷却的水下平板边界层转捩位置的预测方法。
背景技术
流体在物体表面流动过程中,表面边界层内的流动存在层流和湍流两种流动状态。层流状态和湍流状态的性质有着显著的区别,如阻力、噪声、热流等这些对实际工程问题有重要影响的物理量。流动状态从层流转变为湍流的过程称为转捩。控制边界层的转捩对航海、航空、航天等工程领域中的实际问题有着重要的影响。
前人对平板边界层有较多的实验、计算以及理论研究成果,所以可以从平板边界层入手,研究边界层转捩控制。控制平板边界层转捩的方法有很多,加热/冷却壁面是其中一种重要的控制手段,而预测壁面加热/冷却的水下平板边界层转捩位置对这一控制手段的具体实施是非常关键的。
目前对于壁面加热/冷却的水下平板边界层的转捩预测,现有可能的方法有三种。一是采用实验办法,但是对实验环境和设备的要求较高,实验花费成本较高;二是采用直接数值模拟方法,但是计算量太大,当前的计算水平无法用于实际问题;三是采用经验或半经验的方法,其中半经验的eN方法是目前看来可用的一种转捩预测方法。
发明内容
本发明的目的在于提供一种壁面加热/冷却的水下平板边界层转捩位置的预测方法。本发明的技术方案如下:
一种壁面加热/冷却的水下平板边界层转捩位置的预测方法,包括以下步骤:
步骤一,确定计算工况,包括确定来流温度、来流速度以及需要计算的水下平板壁面加热/冷却时的壁面温度;
步骤二,针对壁面加热/冷却的水下平板边界层,考虑能量方程和温度变化,考虑温度和速度的耦合,且考虑粘性系数和热传导系数随温度的变化,而密度和比热则视为常数,建立层流基本流二维定常有量纲的边界层方程,再进一步推导出考虑粘度、热传导系数随温度变化的Blasius方程,从而获得步骤一中所确定工况下的壁面加热/冷却的水下平板边界层的层流基本流流场,方法如下:
1)建立层流基本流二维定常有量纲的边界层方程(1):
Figure GDA0003890316690000021
公式中,带“*”的量为有量纲量,不带“*”的量为无量纲量,有量纲的边界层方程(1)中第一式为连续性方程,第二式为流向动量方程,第三式为法向动量方程,第四式为能量方程;x*、y*分别为直角坐标系中流向坐标和法向坐标;u*、v*、T*、p*分别为流向速度、法向速度、温度、压强;ρ*、c*、μ*、k*分别为密度、比热、动力学粘性系数、热传导系数;
2)引入变量f、g和坐标η,且
Figure GDA0003890316690000022
变量f和g均为坐标η的函数,f对η一阶导数
Figure GDA0003890316690000023
为无量纲的速度,
Figure GDA0003890316690000024
为无量纲的温度,
Figure GDA0003890316690000025
为壁面处的温度变量g;带下角标“e”的量表示来流中的量,带下角标“w”的量表示壁面处的量,
Figure GDA0003890316690000026
为来流速度,
Figure GDA0003890316690000027
为来流中的动力学粘性系数,
Figure GDA0003890316690000028
为壁面温度,
Figure GDA0003890316690000029
为来流温度;将f和g的表达式代入有量纲的边界层方程(1),考虑流向压力梯度为零的情况,即得到考虑粘度、热传导系数随温度变化的Blasius方程(2):
Figure GDA00038903166900000210
考虑粘度、热传导系数随温度变化的Blasius方程(2)的系数分别为
Figure GDA00038903166900000211
设壁面处使用无滑移条件和等温条件,远离壁面处使用来流速度和来流温度条件,则考虑粘度、热传导系数随温度变化的Blasius方程(2)的边界条件(3):
Figure GDA0003890316690000031
3)使用打靶法和四阶的Runge-Kutta法,求解具有边界条件(3)的考虑粘度、热传导系数随温度变化的Blasius方程,得到壁面加热/冷却的水下平板边界层的层流基本流流场;
步骤三,基于线性稳定性理论,考虑能量方程和温度变化,考虑温度和速度的耦合,且考虑粘性系数和热传导系数随温度的变化,而密度和比热则视为常数,并且考虑温度的扰动,对壁面加热/冷却的水下平板边界层进行稳定性分析,得到各频率小扰动增长率沿着流向的分布即增长率云图;
步骤四,计算壁面加热/冷却的水下平板边界层小扰动沿流向最大增长倍数,方法如下:
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数,由式(7)给出:
Figure GDA0003890316690000032
其中A表示小扰动在位置x处的幅值,A0表示中性点的小扰动幅值,x0为中性点对应的流向位置,N为幅值放大因子;
2)通过从中性曲线开始,对幅值增长率沿水下平板壁面积分,获得各频率小扰动的幅值放大因子N沿流向坐标的分布,用一条曲线将不同频率波的N值沿流向坐标分布的最大值包络,包络线的各点取值代表各种频率扰动在此处能达到的最大幅值放大因子Nmax,获得Nmax曲线,最大增长倍数的大小为
Figure GDA0003890316690000035
步骤五,预测壁面加热/冷却的水下平板边界层的转捩位置;
1)认为只要某频率的小扰动波幅在沿流向发展的过程中,幅值达到初始幅值的
Figure GDA0003890316690000034
倍时层流即发展为湍流,作为转捩判据的Ntr值的大小由实验标定;在步骤四中获得的Nmax曲线上找到Ntr对应流向位置,即为预测出来的壁温加热/冷却的水下平板的转捩位置。
进一步地,步骤五中,转捩判据取为Ntr=7。
进一步地,步骤三的方法如下:
1)建立无量纲的三维流动的线性扰动方程(4):
Figure GDA0003890316690000033
其中,x、y、z分别为直角坐标系中的流向坐标、法向坐标、展向坐标,t为时间,φ'=(u',v',w',T',p')T为扰动矢量,且φ'为关于t、x、y、z的函数,u'为流向速度扰动量、v'为法向速度扰动量、w'为展向速度扰动量、T'温度扰动量、p'为压强扰动量,Γ'、A'、B'、C'、D'、V11'、V22'、V33'、V12'、V13'、V23'为系数矩阵,系数矩阵中包含基本流的信息,基本流通过步骤二中求解考虑粘度、热传导系数随温度变化的Blasius方程(2)得到;系数矩阵也包含了无量纲参数雷诺数Re、埃克特数Ec、普朗特数Pr;
2)在边界层流向上使用近似平行流假设,得到扰动展开式(5):
φ'=φ·ei(αx+iβ-ωt)+c.c. (5)
其中,
Figure GDA0003890316690000041
为扰动量的特征函数,且φ仅为关于y的函数,u为流向速度扰动量的特征函数、
Figure GDA0003890316690000042
为法向速度扰动量的特征函数、w为展向速度扰动量的特征函数、T温度扰动量的特征函数、p为压强扰动量的特征函数;ω为无量纲圆频率,对应的有量纲频率为f*,β为展向波数,c.c.表示复共轭;空间模式下α=αr+iαi为复数形式波数,αr为实数形式波数,-αi为空间增长率;
3)将扰动展开式(5)代入线性扰动方程(4),得到线性稳定性方程(6):
Figure GDA0003890316690000043
其中,V22、B、D为系数矩阵,表达式为:
V22=V22',
B=B'+iαV12'+iβV23',
D=-iωΓ'+iαA'+iβC'+D'-α2V11'-β2V33'-αβV13';
4)系数矩阵V22、B、D中包含基本流的信息和参数雷诺数Re、复波数α、展向波数β、圆频率ω、埃克特数Ec、普朗特数Pr;采用数值方法计算线性稳定性方程(6),获得基本流不同流向站位的特征值,即扰动展开式(5)中的波数α、圆频率ω与基本流流场雷诺数Re的参数关系;通过特征值的计算得到各层流基本流的中性曲线,并得到各频率小扰动增长率-αi沿着流向的分布,即增长率云图。
进一步地,步骤四的方法如下:
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数,由式(7)给出:
Figure GDA0003890316690000051
其中A表示小扰动在位置x处的幅值,A0表示中性点的小扰动幅值,x0为中性点对应的流向位置,N为幅值放大因子;
2)通过从中性曲线开始,对幅值增长率沿水下平板壁面积分,获得各频率小扰动的幅值放大因子N沿流向坐标的分布,用一条曲线将不同频率波的N值沿流向坐标分布的最大值包络,包络线的各点取值代表各种频率扰动在此处能达到的最大幅值放大因子Nmax,从而获得Nmax曲线,最大增长倍数的大小为
Figure GDA0003890316690000052
本发明与现有技术相比的有益效果:
(1)本发明基于线性稳定性分析理论,建立了考虑壁面加热/冷却的水下平板边界层的转捩位置预测方法;
(2)本发明可用于通过加热/冷却壁面来控制水下平板边界层的转捩位置。
附图说明
图1为本发明步骤框图
图2为所计算工况的速度剖面
图3为所计算工况的温度剖面
图4为所计算工况的中性曲线
图5为加热壁面的增长率云图
图6为加热壁面的N值分布云图
图7为所计算工况的Nmax曲线
具体实施方式
下面结合具体实施方式及附图对本发明进行详细说明。在下面的描述中,出于解释而非限制性的目的,阐述了具体细节,以帮助全面地理解本发明。然而,对本领域技术人员来说显而易见的是,也可以在脱离了这些具体细节的其它实施例中实践本发明。
在此需要说明的是,为了避免因不必要的细节而模糊了本发明,在附图中仅仅示出了与根据本发明的方案密切相关的设备结构和/或处理步骤,而省略了与本发明关系不大的其他细节。
本发明实施例提供了一种壁面加热/冷却的水下平板边界层转捩位置的预测方法,如图1所示,包括下述步骤:
步骤一,根据实际情况,确定计算工况,即确定来流温度、来流速度以及要计算的水下平板壁面加热/冷却时的壁面温度。
该步骤中,取计算工况为:零攻角,来流速度
Figure GDA0003890316690000061
来流温度
Figure GDA0003890316690000062
平板壁面温度
Figure GDA0003890316690000063
取为4.4℃(壁面冷却)、15.6℃(壁面温度等于来流温度)、32.2℃(壁面加热)。
步骤二,计算步骤一中所确定工况下的水下平板的层流基本流,获得层流基本流流场。
1)针对壁面加热/冷却的水下平板边界层,考虑能量方程和温度变化,考虑温度和速度的耦合,且考虑粘性系数和热传导系数随温度的变化,而密度和比热则视为常数。从二维定常有量纲的N-S方程出发,推导出有量纲的边界层方程,再进一步推导出Blasius方程,获得层流基本流流场。方法如下:
2)层流基本流满足二维定常有量纲N-S方程(N-S方程的形式为本领域公知技术,具体见《流体力学》,周光炯等,高等教育出版社)。设定两个条件,一是边界层厚度在量级上小于到平板前缘距离,二是边界层内的粘性力和惯性力的量级相同。对二维定常有量纲N-S方程进行量级分析,忽略高阶小量,得到如下的有量纲的边界层方程(1)(这一步推导过程为本领域公知技术,具体见《流体力学》,周光炯等,高等教育出版社):
Figure GDA0003890316690000064
公式中,带“*”的量为有量纲量,不带“*”的量为无量纲量,方程(1)中第一式为连续性方程,第二式为流向动量方程,第三式为法向动量方程,第四式为能量方程。x*、y*分别为直角坐标系中流向坐标和法向坐标;u*、v*、T*、p*分别为流向速度、法向速度、温度、压强;ρ*、c*、μ*、k*分别为密度、比热、动力学粘性系数、热传导系数。
3)引入变量f、g和坐标η,且
Figure GDA0003890316690000065
变量f和g均为坐标η的函数,f对η一阶导数
Figure GDA0003890316690000066
为无量纲的速度,
Figure GDA0003890316690000067
为无量纲的温度,
Figure GDA0003890316690000068
为壁面处的温度变量g;带下角标“e”的量表示来流中的量,带下角标“w”的量表示壁面处的量,
Figure GDA0003890316690000071
为来流速度,
Figure GDA0003890316690000072
为来流中的动力学粘性系数,
Figure GDA0003890316690000073
为壁面温度,
Figure GDA0003890316690000074
为来流温度;将f和g的表达式代入有量纲的边界层方程(1),考虑流向压力梯度为零的情况,即得到考虑粘度、热传导系数随温度变化的Blasius方程(2)(这一步推导过程为本领域公知技术,具体见《高超声速和高温气体动力学》,John D.Anderson Jr等,航空工业出版社):
Figure GDA0003890316690000075
方程(2)的系数分别为
Figure GDA0003890316690000076
设壁面处使用无滑移条件和等温条件,远离壁面处使用来流速度和来流温度条件,则Blasius方程(2)的边界条件为:
Figure GDA0003890316690000077
4)使用打靶法和四阶的Runge-Kutta法(这一方法为本领域公知技术,具体见《高超声速和高温气体动力学》,John D.Anderson Jr等,航空工业出版社),求解具有边界条件(3)的Blasius方程(2),求解过程中使用适用于0.1MPa、0-100℃的水的动力学粘性系数μ*和热传导系数k*随温度T*的变化关系;再使用关系式
Figure GDA0003890316690000078
Figure GDA0003890316690000079
得到壁面加热/冷却的水下平板边界层的层流基本流流场。
5)在该步骤中,适用于0.1MPa、0-100℃的水的动力学粘性系数μ*和热传导系数k*随温度T*的变化关系使用关系式(4)和(5)(具体见文献Huber M L,Perkins R A,LaeseckeA,et al.New International Formulation for the Viscosity of H2O[J].Journal ofPhysical&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 Conductivityof H2O[J].Journal of Physical&Chemical Reference Data,2012,41(3):387-440.):
Figure GDA0003890316690000081
Figure GDA0003890316690000082
6)该步骤中,层流基本流的边界层内需要保证足够的网格数量,以保证后续稳定性分析的需求。
7)该步骤可得到所计算工况的基本流剖面,即速度剖面如图2所示和温度剖面如图3所示。
步骤三,基于线性稳定性理论,对壁面加热/冷却的水下平板边界层进行稳定性分析,得出小扰动的增长率。
1)针对壁面加热/冷却的水下平板边界层,考虑能量方程和温度变化,考虑温度和速度的耦合,且考虑粘性系数和热传导系数随温度的变化,而密度和比热则视为常数,此外,还要考虑温度的扰动。从无量纲的三维流动的N-S方程(N-S方程的形式为本领域公知技术,具体见《流体力学》,周光炯等,高等教育出版社)出发,使用小扰动假设推导出线性扰动方程,再使用近似平行流假设推导出线性稳定性方程,从而进行稳定性分析。方法如下:
2)考虑小扰动的情况,将瞬时量表示为层流基本流和小扰动之和,即
(us,vs,ws,Ts,pss,ks)T=(u,v,w,T,p,μ,k)T+(u',v',w',T',p',μ',k')T
其中,us、vs、ws、Ts、ps、μs、ks为瞬时量,且ws为展向速度;u、v、w、T、p、μ、k为层流基本流的量,u'、v'、w'、T'、p'、μ'、k'为扰动量;由粘性系数和热传导系数随温度的变化关系,可得到粘性系数和热传导系数的扰动量随温度扰动量的变化关系式,表示为μ'=τT',k'=σT',其中τ和σ为基本流温度T的函数。
3)瞬时量和基本流都满足无量纲的三维流动的N-S方程,用瞬时量的方程减去基本流的方程,再略去非线性的高阶小量,得到线性扰动方程(6),用矢量形式表示为:
Figure GDA0003890316690000091
其中,x、y、z分别为直角坐标系中的流向坐标、法向坐标、展向坐标,t为时间,φ'=(u',v',w',T',p')T为扰动矢量,Γ'、A'、B'、C'、D'、V11'、V22'、V33'、V12'、V13'、V23'为系数矩阵,系数矩阵中包含基本流的信息,基本流通过步骤二中求解Blasius方程(2)得到;系数矩阵也包含了无量纲参数雷诺数Re、埃克特数Ec、普朗特数Pr。
4)考虑到边界层在流向上变化缓慢,在流向上使用近似平行流假设,将扰动描述为行进波的形式,得到扰动展开式(7):
φ′(t,x,y,z)=φ(y)·ei(αx+iβ-ωt)+c.c. (7)
其中,
Figure GDA0003890316690000093
为扰动的特征函数;ω为无量纲圆频率,对应的有量纲频率为f*(单位是Hz),β为展向波数,c.c.表示复共轭。空间模式下α=αr+iαi为复数形式波数,αr为实数形式波数,-αi为空间增长率。将扰动展开式(7)代入线性扰动方程(6),得到线性稳定性方程(8):
Figure GDA0003890316690000094
其中,V22、B、D为系数矩阵,表达式为:
V22=V22',
B=B'+iαV12'+iβV23',
D=-iωΓ'+iαA'+iβC'+D'-α2V11'-β2V33'-αβV13'。
5)系数矩阵V22、B、D中包含基本流的信息和参数雷诺数Re、复波数α、展向波数β、圆频率ω、埃克特数Ec、普朗特数Pr。采用数值方法计算线性稳定性方程(8),获得基本流不同流向站位的特征值,即扰动展开式中的波数α、圆频率ω与基本流流场雷诺数Re的参数关系。通过特征值的计算可以得到各层流基本流的中性曲线;也可以得到各频率小扰动增长率-αi沿着流向的分布,即增长率云图。
6)该步骤中,稳定性方程特征值分析方法可采用本领域公知的手段,具体见《流动稳定性》,周恒等,国防工业出版社。
7)该步骤中得到的有量纲的中性曲线如图4所示,
Figure GDA0003890316690000101
时的增长率云图如图5所示。
步骤四,计算壁面加热/冷却的水下平板边界层小扰动沿流向最大增长倍数。
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数,由式(7)给出:
Figure GDA0003890316690000102
其中A表示小扰动在位置x处的幅值,A0表示中性点的小扰动幅值,x0为中性点对应的流向位置,N为幅值放大因子。
2)该步骤中,通过从中性曲线开始,对幅值增长率沿水下平板壁面积分,获得各频率小扰动的幅值放大因子N沿流向坐标的分布,工况
Figure GDA0003890316690000103
的N值分布云图如图6所示。
3)该步骤中,用一条曲线将不同频率波的N值沿流向坐标分布的最大值包络,包络线代表各种频率扰动在此处能达到的最大幅值放大因子,获得Nmax曲线,所计算工况的Nmax曲线如图7所示。Nmax曲线表示出不同频率的小扰动可以达到的最大增长倍数,这个倍数的大小为
Figure GDA0003890316690000104
步骤五,预测壁面加热/冷却的水下平板边界层的转捩位置。
1)该步骤中,根据对边界层自然转捩的认识,可以认为只要某频率的小扰动波幅在沿流向发展的过程中,幅值达到初始幅值的
Figure GDA0003890316690000105
倍时层流即发展为湍流。
2)该步骤中,作为转捩判据的Ntr值的大小由实验标定出来,转捩判据取为Ntr=7。
3)该步骤中,在步骤四中获得的Nmax曲线上找到Ntr对应流向位置,即为预测出来的壁温加热/冷却的水下平板的转捩位置。
4)该步骤中,可得到每个工况的转捩位置以及转捩延迟距离,将结果列在表1中:
表1不同壁温时的转捩预测结果(来流温度为15.6℃)
Figure GDA0003890316690000106
Figure GDA0003890316690000111
表1中,
Figure GDA0003890316690000112
是壁温为
Figure GDA0003890316690000113
时的转捩位置,
Figure GDA0003890316690000114
是壁温为15.6℃时的转捩距离,
Figure GDA0003890316690000115
为转捩延迟距离,其表达式为
Figure GDA0003890316690000116
如上针对一种实施例描述和/或示出的特征可以以相同或类似的方式在一个或更多个其它实施例中使用,和/或与其它实施例中的特征相结合或替代其它实施例中的特征使用。
应该强调,术语“包括/包含”在本文使用时指特征、整件、步骤或组件的存在,但并不排除一个或更多个其它特征、整件、步骤、组件或其组合的存在或附加。
本发明以上的装置和方法可以由硬件实现,也可以由硬件结合软件实现。本发明涉及这样的计算机可读程序,当该程序被逻辑部件所执行时,能够使该逻辑部件实现上文所述的装置或构成部件,或使该逻辑部件实现上文所述的各种方法或步骤。本发明还涉及用于存储以上程序的存储介质,如硬盘、磁盘、光盘、DVD、flash存储器等。
这些实施例的许多特征和优点根据该详细描述是清楚的,因此所附权利要求旨在覆盖这些实施例的落入其真实精神和范围内的所有这些特征和优点。此外,由于本领域的技术人员容易想到很多修改和改变,因此不是要将本发明的实施例限于所例示和描述的精确结构和操作,而是可以涵盖落入其范围内的所有合适修改和等同物。
本发明未详细说明部分为本领域技术人员公知技术。

Claims (3)

1.一种壁面加热/冷却的水下平板边界层转捩位置的预测方法,包括以下步骤:
步骤一,确定计算工况,包括确定来流温度、来流速度以及需要计算的水下平板壁面加热/冷却时的壁面温度;
步骤二,针对壁面加热/冷却的水下平板边界层,考虑能量方程和温度变化,考虑温度和速度的耦合,且考虑粘性系数和热传导系数随温度的变化,而密度和比热则视为常数,建立层流基本流二维定常有量纲的边界层方程,再进一步推导出考虑粘度、热传导系数随温度变化的Blasius方程,从而获得步骤一中所确定工况下的壁面加热/冷却的水下平板边界层的层流基本流流场,方法如下:
1)建立层流基本流二维定常有量纲的边界层方程(1):
Figure FDA0003901733600000011
公式中,带“*”的量为有量纲量,不带“*”的量为无量纲量,有量纲的边界层方程(1)中第一式为连续性方程,第二式为流向动量方程,第三式为法向动量方程,第四式为能量方程;x*、y*分别为直角坐标系中流向坐标和法向坐标;u*、v*、T*、p*分别为流向速度、法向速度、温度、压强;ρ*、c*、μ*、k*分别为密度、比热、动力学粘性系数、热传导系数;
2)引入变量f、g和坐标η,且
Figure FDA0003901733600000012
变量f和g均为坐标η的函数,f对η一阶导数
Figure FDA0003901733600000013
为无量纲的速度,
Figure FDA0003901733600000014
为无量纲的温度,
Figure FDA0003901733600000015
为壁面处的温度变量g;带下角标“e”的量表示来流中的量,带下角标“w”的量表示壁面处的量,
Figure FDA0003901733600000016
为来流速度,
Figure FDA0003901733600000017
为来流中的动力学粘性系数,
Figure FDA0003901733600000018
为壁面温度,
Figure FDA0003901733600000019
为来流温度;将f和g的表达式代入有量纲的边界层方程(1),考虑流向压力梯度为零的情况,即得到考虑粘度、热传导系数随温度变化的Blasius方程(2):
Figure FDA0003901733600000021
考虑粘度、热传导系数随温度变化的Blasius方程(2)的系数分别为
Figure FDA0003901733600000022
Figure FDA0003901733600000023
设壁面处使用无滑移条件和等温条件,远离壁面处使用来流速度和来流温度条件,则考虑粘度、热传导系数随温度变化的Blasius方程(2)的边界条件(3):
Figure FDA0003901733600000024
3)使用打靶法和四阶的Runge-Kutta法,求解具有边界条件(3)的考虑粘度、热传导系数随温度变化的Blasius方程,得到壁面加热/冷却的水下平板边界层的层流基本流流场;
步骤三,基于线性稳定性理论,考虑能量方程和温度变化,考虑温度和速度的耦合,且考虑粘性系数和热传导系数随温度的变化,而密度和比热则视为常数,并且考虑温度的扰动,对壁面加热/冷却的水下平板边界层进行稳定性分析,得到各频率小扰动增长率沿着流向的分布即增长率云图,方法如下:
1)建立无量纲的三维流动的线性扰动方程(4):
Figure FDA0003901733600000025
其中,x、y、z分别为直角坐标系中的流向坐标、法向坐标、展向坐标,t为时间,φ'=(u',v',w',T',p')T为扰动矢量,且φ'为关于t、x、y、z的函数,u'为流向速度扰动量、v'为法向速度扰动量、w'为展向速度扰动量、T'温度扰动量、p'为压强扰动量,Γ'、A'、B'、C'、D'、V11'、V22'、V33'、V12'、V13'、V23'为系数矩阵,系数矩阵中包含基本流的信息,基本流通过步骤二中求解考虑粘度、热传导系数随温度变化的Blasius方程(2)得到;系数矩阵也包含了无量纲参数雷诺数Re、埃克特数Ec、普朗特数Pr;
2)在边界层流向上使用近似平行流假设,得到扰动展开式(5):
φ'=φ·ei(αx+iβ-ωt)+c.c. (5)
其中,
Figure FDA0003901733600000031
为扰动量的特征函数,且φ仅为关于y的函数,u为流向速度扰动量的特征函数、
Figure FDA0003901733600000032
为法向速度扰动量的特征函数、w为展向速度扰动量的特征函数、T温度扰动量的特征函数、p为压强扰动量的特征函数;ω为无量纲圆频率,对应的有量纲频率为f*,β为展向波数,c.c.表示复共轭;空间模式下α=αr+iαi为复数形式波数,αr为实数形式波数,-αi为空间增长率;
3)将扰动展开式(5)代入线性扰动方程(4),得到线性稳定性方程(6):
Figure FDA0003901733600000033
其中,V22、B、D为系数矩阵,表达式为:
V22=V22',
B=B'+iαV12'+iβV23',
D=-iωΓ'+iαA'+iβC'+D'-α2V11'-β2V33'-αβV13';
4)系数矩阵V22、B、D中包含基本流的信息和参数雷诺数Re、复波数α、展向波数β、圆频率ω、埃克特数Ec、普朗特数Pr;采用数值方法计算线性稳定性方程(6),获得基本流不同流向站位的特征值,即扰动展开式(5)中的波数α、圆频率ω与基本流流场雷诺数Re的参数关系;通过特征值的计算得到各层流基本流的中性曲线,并得到各频率小扰动增长率-αi沿着流向的分布,即增长率云图;
步骤四,计算壁面加热/冷却的水下平板边界层小扰动沿流向最大增长倍数
Figure FDA0003901733600000041
步骤五,预测壁面加热/冷却的水下平板边界层的转捩位置,方法为:根据只要某频率的小扰动波幅在沿流向发展的过程中,幅值达到初始幅值的
Figure FDA0003901733600000042
倍时层流即发展为湍流的认定,设定作为转捩判据的Ntr值;在步骤四中获得的Nmax曲线上找到Ntr对应流向位置,即为预测出来的壁温加热/冷却的水下平板的转捩位置。
2.根据权利要求1所述的壁面加热/冷却的水下平板边界层转捩位置的预测方法,其特征在于,步骤五中,转捩判据取为Ntr=7。
3.根据权利要求1所述的壁面加热/冷却的水下平板边界层转捩位置的预测方法,其特征在于,步骤四的方法如下:
1)给定频率小扰动的幅值从中性点至流向某位置的放大倍数,由式(7)给出:
Figure FDA0003901733600000043
其中A表示小扰动在位置x处的幅值,A0表示中性点的小扰动幅值,x0为中性点对应的流向位置,N为幅值放大因子;
2)通过从中性曲线开始,对幅值增长率沿水下平板壁面积分,获得各频率小扰动的幅值放大因子N沿流向坐标的分布,用一条曲线将不同频率波的N值沿流向坐标分布的最大值包络,包络线的各点取值代表各种频率扰动在此处能达到的最大幅值放大因子Nmax,从而获得Nmax曲线,最大增长倍数的大小为
Figure FDA0003901733600000044
CN202210174746.9A 2022-02-24 2022-02-24 一种壁面加热/冷却的水下平板边界层转捩位置的预测方法 Active CN114528665B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210174746.9A CN114528665B (zh) 2022-02-24 2022-02-24 一种壁面加热/冷却的水下平板边界层转捩位置的预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210174746.9A CN114528665B (zh) 2022-02-24 2022-02-24 一种壁面加热/冷却的水下平板边界层转捩位置的预测方法

Publications (2)

Publication Number Publication Date
CN114528665A CN114528665A (zh) 2022-05-24
CN114528665B true CN114528665B (zh) 2022-12-09

Family

ID=81624902

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210174746.9A Active CN114528665B (zh) 2022-02-24 2022-02-24 一种壁面加热/冷却的水下平板边界层转捩位置的预测方法

Country Status (1)

Country Link
CN (1) CN114528665B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115146383B (zh) * 2022-07-01 2023-01-24 天津大学 一种预报超疏水表面曲面边界层转捩位置的方法
CN117407635B (zh) * 2023-10-18 2024-05-14 中国空气动力研究与发展中心计算空气动力研究所 一种基于结霜相似律的平板结霜厚度预测方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4932610A (en) * 1986-03-11 1990-06-12 The United States Of America As Represented By The United States National Aeronautics And Space Administration Active control of boundary layer transition and turbulence
US10953979B2 (en) * 2015-11-11 2021-03-23 The Arizona Board Of Regents On Behalf Of The University Of Arizona Control of hypersonic boundary layer transition
CN108304600B (zh) * 2017-08-09 2020-03-31 北京空天技术研究所 一种高超声速飞行器转捩位置预测方法
CN108304601B (zh) * 2017-08-09 2019-12-31 北京空天技术研究所 一种高超声速飞行器边界层转捩的判断方法
CN112395694B (zh) * 2020-12-03 2023-05-02 中国人民解放军国防科技大学 一种超高速湍流边界层减阻控制方法
CN112861447B (zh) * 2021-02-09 2022-03-04 天津大学 一种基于稳定性理论的水下回转体首部线型优选设计方法
CN113657051B (zh) * 2021-08-19 2022-05-20 天津大学 一种水下航行体表面层流区壁面脉动压力频谱预测方法
CN114036873B (zh) * 2021-11-28 2022-05-27 天津大学 考虑前缘区修正的超疏水表面平板层流边界层的计算方法

Also Published As

Publication number Publication date
CN114528665A (zh) 2022-05-24

Similar Documents

Publication Publication Date Title
CN114528665B (zh) 一种壁面加热/冷却的水下平板边界层转捩位置的预测方法
CN113221350B (zh) 基于全局稳定性分析的高超声速飞行器转捩预测方法
CN108304600B (zh) 一种高超声速飞行器转捩位置预测方法
Jiao et al. Marangoni abnormal convection heat transfer of power-law fluid driven by temperature gradient in porous medium with heat generation
CN107832260B (zh) 一种平板冲击射流传热问题的数值模拟方法
Siddiqa et al. Natural convection flow of viscous fluid over triangular wavy horizontal surface
CN113657051B (zh) 一种水下航行体表面层流区壁面脉动压力频谱预测方法
Siddiqa et al. Double diffusive natural convection flow over a wavy surface situated in a non-absorbing medium
Li et al. Coupled simulation of fluid flow and propellant burning surface regression in a solid rocket motor
Zou et al. Control of the ventilated supercavity on the maneuvering trajectory
Fan et al. Prediction of hypersonic boundary layer transition with variable specific heat on plane flow
Mahmoodi et al. Cooling process of liquid propellant rocket by means of kerosene-alumina nanofluid
CN113792508B (zh) 考虑表面质量引射效应的气动热计算方法
Malinowski et al. Implementation of the axially symmetrical and three dimensional finite element models to the determination of the heat transfer coefficient distribution on the hot plate surface cooled by the water spray nozzle
Hossain et al. A numerical study on unsteady natural convection flow with temperature dependent viscosity past an isothermal vertical cylinder
Ahmed et al. UNSTEADY MHD CHEMICALLY REACTING FLUID THROUGH A POROUS MEDIUM BOUNDED BY A NON-ISOTHERMAL IMPULSIVELY-STARTED VERTICAL PLATE: A NUMERICAL TECHNIQUE.
Komurasaki A hydrothermal convective flow at extremely high temperature
Han et al. Buoyancy-driven convection heat transfer of copper–water nanofluid in a square enclosure under the different periodic oscillating boundary temperature waves
Peetala et al. Conjugate heat transfer analysis for a finite thickness cylinder in a hypersonic flow
CN107818197B (zh) 一种基于piv技术的超声速翼型的测力方法和装置
Khalid et al. Heat alleviation studies on hypersonic re-entry vehicles
Li et al. The Simulation of Wraparound Fins Aerodynamic Characteristics
Isaiah et al. MHD stokes heat and mass flow past an infinite permeable plate subjected to convective surface boundary conditions [J]
Illés et al. Thermal Characterization of Solid Structures during forced convection heating
Tangthieng et al. Thermosolutal transport and macrosegregation during freeze coating of a binary substance on a continuous moving object

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant