外掠平板层流流动边界层微分方程积分解的实现
Implementation of Integration Solution to Boundary Layer Equation of Laminar Flow over Plane

作者: 林志敏 * , 王良璧 :兰州交通大学机电工程学院,甘肃 兰州; 张永恒 , 张旭耀 :兰州交通大学新能源与动力工程学院,甘肃 兰州;

关键词: 边界层微分方程流体力学传热学Boundary Layer Differential Equation Fluid Mechanics Heat Transfer

摘要: 边界层理论及其在典型流动问题中的应用是流体力学和传热学的重要内容,通过对外掠平板层流流动和换热等问题的分析,可以揭示流体流动和传热的基本规律,是认识更为复杂流体流动和传热现象的基础。在分析和介绍边界层微分方程建立方法、相似变量表达式以及边界层方程典型解法的基础上,重点探讨了外掠平板边界层相似变换微分方程的积分解解法,并给出了相应的MATLAB计算程序。计算结果表明,边界层相似变换微分方程的积分解法简单、易行,该方法对于教师讲解和学生理解边界层流动和传热特性有积极作用。

Abstract: Boundary layer theory and its application in typical flow problems are important contents of fluid mechanics and heat transfer courses. Through the analysis of laminar flow and heat transfer over a flat plate, the basic laws of fluid flow and heat transfer can be revealed, which is the basis of understanding more complex fluid flow and heat transfer phenomena. Based on the analysis and introduction of the establishment method of the boundary layer differential equation, the expression of similar variables and the typical solution of the boundary layer equation, the integral solution of the similar transformation differential equation of the boundary layer formed by fluid flowing over the flat plate was introduced, and the corresponding MATLAB program was given. The results showed that the integral method for solving boundary layer similarity transformation differential equation was simple and easy to use. This method has a positive effect on teachers’ explanation and students’ understanding of the characteristics of boundary layer flow and heat transfer.

1. 引言

流体流动涉及的动量和能量传递由连续性方程、动量方程和能量方程来描述,对于粘性流体的流动与传热,在靠近壁面附近存在着流体速度和温度沿壁面法线方向急剧变化的薄层,即流动边界层和热边界层。由于实际流体流动的复杂性,自1821年和1845年分别由纳维和斯托克斯建立起描述流体运动的N-S方程的半个多世纪里,除少数几种简单流动情况外,在解决实际流体的流动中N-S方程未能得到广泛的应用。自1904年普朗特提出边界层理论后,才使人们借助数学工具求解实际流体的流动问题成为可能,使得基于数学分析的理论流体动力学与主要基于实验的水力学有机结合,极大地促进了现代流体力学的发展 [1]。基于边界层理论和混合长度理论获得的典型内流、外流流动和传热的结果成为流体力学 [2] 和传热学 [3] 的重要内容。边界层理论的一个典型应用是求解外掠平壁的流动和传热,通过相似变换边界层方程简化为含远端渐进边界条件的非线性常微分方程,该方程最初由布拉休斯给出级数解,但该级数解中含有未定的起始端边界条件( s = f ( 0 ) ),随后由霍华斯(Howarth)通过数值解得出该未定起始端边界条件 [4]。尽管随着计算机和数值方法的发展,人们已能够通过差分法、打靶法、优化设计法等获得外掠平板和钝体的相似变换方程的数值解 [5] [6] [7] [8] [9],但人们对采用不同的方法获得相似变换方程的解及其性质的探索仍保持着浓厚的兴趣 [10]。

本文以外掠平板层流边界层为例,从边界层微分方程建立、相似变量引入、相似变换方程积分解等方面进行探讨,补充和完善流体力学和传热学的教学内容。

2. 外掠平板层流流动边界层动量方程和能量方程

描述流体外掠平板稳定流动的连续性方程和动量方程为:

u x + v y = 0 (1)

u u x + v u y = 1 ρ p x + μ ρ ( 2 u x 2 + 2 u y 2 ) (2)

u v x + v u y = 1 ρ p y + μ ρ ( 2 v x 2 + 2 v y 2 ) (3)

能量方程为:

u t x + v t y = λ ρ c ( 2 t x 2 + 2 t y 2 ) (4)

在方程(1)~(3)中,认为流体物性为常数,不计作用于流体的质量力,忽略能量方程中的耗散项。根据普朗特边界层理论和量级分析,可得边界层微分方程组为:

连续性方程:

u x + v y = 0 (5)

动量方程:

u u x + v u y = 1 ρ d p d x + μ ρ 2 u y 2 (6)

在边界层外缘:

p s + 1 2 ρ u s 2 = c o n s t (7)

能量方程为:

u t x + v t y = a 2 t y 2 (8)

边界条件为:

{ y = 0 , u = 0 , v = 0 , t = t w y = , u = u s , v = 0 , t = t s (9)

方程(5)~(8)除可采用量级分析法得出外,文献 [11] 和 [12] 分别用摄动法和尺度化分析法导出边界层微分方程。

3. 边界层动量方程和能量方程的相似变换

基于对边界层内惯性力与粘滞力具有相同数量级的认识,并引入流函数,可定义如下的相似变量和无量纲流函数:

η = y ν x / u s (10)

f ( η ) = ψ u s ν x (11)

取外流速度 u s = c x m [13],动量方程(6)变为:

f + m + 1 2 f f + m ( 1 f 2 ) = 0 (12)

能量方程(8)变为:

θ 1 + 1 2 ( m + 1 ) P r f θ 1 = 0 (13)

其中, θ 1 = t t w t s t w

边界条件为:

{ η = 0 , f ( 0 ) = f ( 0 ) = θ 1 ( 0 ) = 0 η , f ( ) = θ 1 ( ) = 1 (14)

方程(12)、(13)为流体外掠楔形体的相似动量方程和能量方程,当m = 0时,则变为Blasius方程。与

平板长度相比较,边界层厚度δ是一个小的量级,且与 x R e x 成正比,形如式(10)的相似变量η的表达式

是基于对边界层的认识,其构造形式存在一定的技巧性 [14],该文通过定义当量时间 τ = x / u s ,应用量纲分析导出相似变量的表达式。文献 [15] 依据群论,从数学上更为严密地导出相似变量的表达式。确定相似变量表达式的另一种方法是,令

η = B x p y , Ψ = A x q f ( η ) (15)

用流函数表示速度u和v,并将其代入动量方程(6),通过比较方程两边同类项的指数,并引入边界条件,也可得到相似变量的表达式。

4. 边界层动量方程和能量方程的积分解及实现

对于流体外掠平板的边界层相似变换的动量方程和能量方程分别为:

f + 1 2 f f = 0 (16)

θ 1 + 1 2 P r f θ 1 = 0 (17)

边界条件仍然为式(14)。

由方程(16)、(17)和边界条件(14)描述的问题是含渐进边界条件的高阶非线性常微分方程两点边值问题。采用打靶法求解该问题的原理是将边值问题转化为初值问题进行求解,具体步骤为:

(1) 给定初始试探值 f ( 0 ) = a 1 f ( 0 ) = a 2 θ 1 ( 0 ) = b 1 θ 1 ( 0 ) = b 2

(2) 通过求解初值问题得到 f 1 * ( ) f 2 * ( ) θ 11 * ( ) θ 12 * ( )

(3) 根据以下两式: f ( 0 ) = a 2 + a 2 a 1 f 2 * ( ) f 1 * ( ) ( f ( ) f 2 * ( ) ) θ 1 ( 0 ) = b 2 + b 2 b 1 θ 12 * ( ) θ 11 * ( ) ( θ 1 ( ) θ 12 * ( ) ) 修正初始试探值;

(4) 判断误差 | f ( ) f 2 * ( ) | | θ 1 ( ) θ 12 * ( ) | 是否小与给定误差限,若小与给定的误差限,则停止迭代,否则转第(2)步。

打靶法的求解过程也可用优化算法实现,对应的优化问题的数学描述为:

x = [ f ( 0 ) , θ 1 ( 0 ) ] 使

min F ( f ( 0 ) , θ 1 ( 0 ) ) (18)

目标函数可表示为: F ( f ( 0 ) , θ 1 ( 0 ) ) = ( f ( ) 1 ) 2 + ( θ 1 ( ) 1 ) 2

常微分方程初值问题和优化问题的求解可借助MATLAB软件的相关函数实现。

观察方程(16)、(17)可以发现,采用直接积分的方法也可获得方程的解。以方程(16)为例,对该方程积分三次,并引入边界条件,得:

f ( η ) = exp ( 1 2 0 η f d η ) 0 η [ exp ( 1 2 0 η f d η ) ] d η (19)

f ( η ) = 0 η f ( η ) d η (20)

f ( η ) = 0 η f ( η ) d η (21)

上述积分方程的求解步骤如下:

(1) 预先选定一个足够大的无量纲坐标 η 的值,例如取 η max = 10 ,并将其分为n等分,得到

η i = η max n 1 i , ( i = 0 , 1 , 2 , , n 1 )

(2) 初始假设 f i = η i

(3) 计算式(19)分母的积分,并记为S;

(4) 依次计算式(19)~(21)分子的积分,结合步骤(3)的结果得到近似解 f ( η i ) , f ( η i ) , f ( η i )

(5) 根据误差 i = 1 n | f i , o l d f i , n e w | 是否小于给定的误差限,决定是否中止计算;或以新算出的 f ( η i ) 重复

步骤(3)和(4)若干次,如重复计算10次,终止计算,并输出结果。

用积分法求解相似变换动量方程(16)的MATLAB程序如下:

functionboundary_inte

clc;

clear;

close all;

eta_infi=10;

det=0.1;

eta=0:det:eta_infi;

N1=length(eta);

ff=eta;

f1p=zeros(1,N1);

fk=zeros(1,N1);

etak=zeros(1,N1);

f2p=zeros(1,N1);

etaj=zeros(1,N1);

fi=zeros(1,N1);

etai=zeros(1,N1);

fenmu1=zeros(1,N1);

NN=10;

fornn=1:NN

for i=1:N1

fi=ff(1:i);

etai=eta(1:i);

fenmu1(i)=inte(fi,etai);

fenmu1(i)=exp(-fenmu1(i)/2);

end

fenmu=inte(fenmu1,eta);

for i=1:N1

for j=1:i

for k=1:j

fk=ff(1:k);

etak=eta(1:k);

f2p(k)=inte(fk,etak);

f2p(k)=exp(-f2p(k)/2)/fenmu;

end

etaj=eta(1:j);

f1p(j)=inte(f2p,etaj);

end

etai=eta(1:i);

ff(i)=inte(f1p,etai);

end

end

[eta' ff' f1p' f2p']

plot(eta,f1p,'-k','linewidth',1.5)

grid on

axis([0 10 0 1.2])

xlabel('\eta')

ylabel('f^\prime')

function f=inte(yy,xx)

N=length(xx);

if N==1

dx=0;

else

dx=xx(2)-xx(1);

end

sum=0;

for i=2:N-1

sum=sum+yy(i);

end

f=(sum+(yy(1)+yy(N))/2)*dx;

计算结果表明,即使数值积分采用精度较低的梯形积分也可得到较好的结果,所得无量纲速度随相似变量变化的曲线如图1所示。用积分法动量方程(16),解法简单,学生易于理解。

Figure 1. Distribution of the dimensionless velocity in the laminar boundary layer over a flat plate

图1. 外掠平板边界层内无因次速度分布

通过对边界层方程的数值求解,学生能够更直观和深入地理解边界层理论,能够自主地对典型外掠物体壁面附近剪切应力的变化以及对流换热问题进行分析。此外,壁面有抽吸、喷注和圆形自由射流等层流流动问题也具有形如式(12)的相似变换微分方程,方程(16)的数值解对于认识和了解这类问题的流动特性也有借鉴作用。由边界层相似变换得出的非线性常微分方程的求解问题还是引出摄动法的来源之一 [1],促进了相关数学方法的发展 [16] [17]。

5. 结束语

对工程问题的数学描述和量化表示,是深刻认识和把握问题内在规律的重要途径。流体外掠典型物体的边界层流动和传热问题是N-S方程和能量方程的具体应用实例,对于理解和掌握流体力学和传热学基本理论,了解复杂工程问题的分析方法具有重要作用。通过对外掠平板层流边界层相似变换微分方程积分解法的分析和运用可以看出,积分求解方法简单、易行,计算过程稳定。边界层理论及其应用是流体力学和传热学的重点和难点内容,应用数值方法分析边界层流动问题对于学生定量理解边界层的发展,分析边界层内流动及流体与壁面的换热特性具有积极作用。

基金项目

兰州交通大学教学改革项目“专业基础课中渐进培养学生解决复杂工程问题能力的探索与实践——以《流体力学》为例”(No. JGY201947)资助。

NOTES

*通讯作者。

文章引用: 林志敏 , 张永恒 , 张旭耀 , 王良璧 (2021) 外掠平板层流流动边界层微分方程积分解的实现。 创新教育研究, 9, 285-292. doi: 10.12677/CES.2021.92045

参考文献

[1] Schlichting, H. and Gersten, K. (2017) Boundary-Layer Theory. 9th Edition, Springer-Verlag, Berlin, Heidelberg.

[2] 孔珑. 工程流体力学[M]. 第4版. 北京: 中国电力出版社, 2014.

[3] 陶文铨. 传热学[M]. 第5版. 北京: 高等教育出版社, 2019.

[4] Liao, S.-J. (1999) An Explicit, Totally Analytic Approximate Solution for Blasius’ Viscous Flow Problems. International Journal of Non-Linear Mechanics, 34, 759-778.
https://doi.org/10.1016/S0020-7462(98)00056-0

[5] 许维德, 孔祥谦. 层流边界层的有限差分计算[J]. 力学与实践, 1981, 3(4): 43-48.

[6] 施天莫. 计算传热学[M]. 陈越南, 范正翘, 陈善年, 等, 译. 北京: 科学出版社, 1987.

[7] 王诚, 周振, 石少卿, 陈暲. 试射法求解二维层流边界层中的Falkner-Skan方程[J]. 重庆理工大学学报, 2016, 30(5): 53-56.

[8] 郭宽良, 孔祥谦, 阵善年. 计算传热学[M]. 合肥: 中国科学技术大学出版社, 1988.

[9] 张永恒, 王良璧, 梁士轩. 基于优化方法的边界层微分方程求解[J]. 科技导报, 2008, 26(2): 46-50.

[10] Varian, V.P. (2013) A Solution of the Blasius Problem. https://library.keldysh.ru/preprint.asp?id=2013-40&lg=e

[11] 王培光, 王振东. 用摄动法建立边界层方程[J]. 力学与实践, 1995, 17(5): 59.

[12] 樊福梅, 梁平, 龙新峰. Prandtl边界层方程推导中的尺度化分析[J]. 华南理工大学学报(自然科学版), 2003, 31(10): 61-64.

[13] 侯义元. 边界层方程相似性解法简述[J]. 北京水利电力经济管理学院学报, 1986(1): 93-102.

[14] 岳曾元. Blasius问题相似性解的一个简单证明[J]. 力学与实践, 1984, 6(1): 52-53.

[15] 杨祖樱. 用群论方法求粘性流动的相似性解[J]. 福州大学学报(自然科学版), 1986(4): 43-52.

[16] 王飞. 基于同伦分析方法的几类边界层问题研究[D]: [硕士学位论文]. 福州: 福州大学, 2016.

[17] 周冉. 重整化群方法在边界层问题中的应用[D]: [博士学位论文]. 长春: 吉林大学, 2019.

分享
Top