¹兰州交通大学机电工程学院，甘肃 兰州

²兰州交通大学新能源与动力工程学院，甘肃 兰州

收稿日期：2021年2月3日；录用日期：2021年4月2日；发布日期：2021年4月9日

摘要

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

关键词

边界层，微分方程，流体力学，传热学

Implementation of Integration Solution to Boundary Layer Equation of Laminar Flow over Plane

Zhimin Lin^1*, Yongheng Zhang², Xuyao Zhang², Liangbi Wang¹

¹School of Mechatronic Engineering, Lanzhou Jiaotong University, Lanzhou Gansu

²School of New Energy and Power Engineering, Lanzhou Jiaotong University, Lanzhou Gansu

Received: Feb. 3^rd, 2021; accepted: Apr. 2^nd, 2021; published: Apr. 9^th, 2021

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.

Keywords:Boundary Layer, Differential Equation, Fluid Mechanics, Heat Transfer

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

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

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

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

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

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (1)

$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}u}{\partial x^2}\text{+}\frac{{\partial }^{2}u}{\partial y^2}\right)$ (2)

$u\frac{\partial v}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}v}{\partial x^2}\text{+}\frac{{\partial }^{2}v}{\partial y^2}\right)$ (3)

能量方程为：

$u\frac{\partial t}{\partial x}+v\frac{\partial t}{\partial y}=\frac{\lambda }{\rho c}\left(\frac{{\partial }^{2}t}{\partial x^2}\text{+}\frac{{\partial }^{2}t}{\partial y^2}\right)$ (4)

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

连续性方程：

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (5)

动量方程：

$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{\text{d}p}{\text{d}x}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial y^2}$ (6)

在边界层外缘：

$p_s+\frac{1}{2}\rho u_s^2=const$ (7)

能量方程为：

$u\frac{\partial t}{\partial x}+v\frac{\partial t}{\partial y}=a\frac{{\partial }^{2}t}{\partial y^2}$ (8)

边界条件为：

$\{\begin{array}{l}y=0,u=0,v=0,t=t_w\\ y=\infty ,u=u_s,v=0,t=t_s\end{array}$ (9)

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

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

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

$\eta =\frac{y}{\sqrt{\nu x/u_s}}$ (10)

$f\left(\eta \right)=\frac{\psi }{\sqrt{u_s\nu x}}$ (11)

取外流速度 $u_s=c{x}^m$ [13]，动量方程(6)变为：

$f^{‴}+\frac{m+1}{2}f{f}^{″}+m\left(1-{f^{\prime }}^2\right)=0$ (12)

能量方程(8)变为：

${{\theta }^{″}}_1+\frac{1}{2}\left(m+1\right)Prf{{\theta }^{\prime }}_1=0$ (13)

其中， ${\theta }_1=\frac{t-t_w}{t_s-t_w}$。

边界条件为：

$\{\begin{array}{l}\eta =0,f\left(0\right)=f^{\prime }\left(0\right)={\theta }_1\left(0\right)=0\\ \eta \to \infty ,f^{\prime }\left(\infty \right)={\theta }_1\left(\infty \right)=1\end{array}$ (14)

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

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

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

$\eta =B{x}^{p}y,\Psi =A{x}^{q}f\left(\eta \right)$ (15)

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

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

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

$f^{‴}+\frac{1}{2}f{f}^{″}=0$ (16)

${{\theta }^{″}}_1+\frac{1}{2}Prf{{\theta }^{\prime }}_1=0$ (17)

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

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

(1) 给定初始试探值 $f^{″}\left(0\right)=a_1$，$f^{″}\left(0\right)=a_2$ 和 ${{\theta }^{\prime }}_1\left(0\right)=b_1$，${{\theta }^{\prime }}_1\left(0\right)=b_2$

(2) 通过求解初值问题得到 $f_1^{*}{}^{\prime }\left(\infty \right)$ 、 $f_2^{*}{}^{\prime }\left(\infty \right)$ 和 ${\theta }_{11}^{*}\left(\infty \right)$ 、 ${\theta }_{12}^{*}\left(\infty \right)$ ；

(3) 根据以下两式： $f^{″}\left(0\right)=a_2+\frac{a_2-a_1}{f_2^{*}{}^{\prime }\left(\infty \right)-f_1^{*}{}^{\prime }\left(\infty \right)}\left(f^{\prime }\left(\infty \right)-f_2^{*}{}^{\prime }\left(\infty \right)\right)$ 和 ${{\theta }^{\prime }}_1\left(0\right)=b_2+\frac{b_2-b_1}{{\theta }_{12}^{*}\left(\infty \right)-{\theta }_{11}^{*}\left(\infty \right)}\left({\theta }_1\left(\infty \right)-{\theta }_{12}^{*}\left(\infty \right)\right)$ 修正初始试探值；

(4) 判断误差 $\left|f^{\prime }\left(\infty \right)-f_2^{*}{}^{\prime }\left(\infty \right)\right|$ 、 $\left|{\theta }_1\left(\infty \right)-{\theta }_{12}^{*}\left(\infty \right)\right|$ 是否小与给定误差限，若小与给定的误差限，则停止迭代，否则转第(2)步。

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

求 $\stackrel{\to }{x}=\left[f^{″}\left(0\right),{{\theta }^{\prime }}_1\left(0\right)\right]$ 使

$\min F\left(f^{″}\left(0\right),{{\theta }^{\prime }}_1\left(0\right)\right)$ (18)

目标函数可表示为： $F\left(f^{″}\left(0\right),{{\theta }^{\prime }}_1\left(0\right)\right)={\left(f^{\prime }\left(\infty \right)-1\right)}^2+{\left({\theta }_1\left(\infty \right)-1\right)}^2$。

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

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

$f^{″}\left(\eta \right)=\frac{\exp \left(-\frac{1}{2}{\displaystyle {\int }_0^{\eta }f\text{d}\eta }\right)}{{\displaystyle {\int }_0^{\eta }\left[\exp \left(-\frac{1}{2}{\displaystyle {\int }_0^{\eta }f\text{d}\eta }\right)\right]\text{d}\eta }}$ (19)

$f^{\prime }\left(\eta \right)={\displaystyle {\int }_0^{\eta }{f}^{″}\left(\eta \right)\text{d}\eta }$ (20)

$f\left(\eta \right)={\displaystyle {\int }_0^{\eta }{f}^{\prime }\left(\eta \right)\text{d}\eta }$ (21)

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

(1) 预先选定一个足够大的无量纲坐标 $\eta$ 的值，例如取 ${\eta }_{\max }=10$，并将其分为n等分，得到

${\eta }_i=\frac{{\eta }_{\max }}{n-1}i,\left(i=0,1,2,\cdots ,n-1\right)$ ；

(2) 初始假设 $f_i={\eta }_i$ ；

(3) 计算式(19)分母的积分，并记为S；

(4) 依次计算式(19)~(21)分子的积分，结合步骤(3)的结果得到近似解 $f^{″}\left({\eta }_i\right),f^{\prime }\left({\eta }_i\right),f\left({\eta }_i\right)$ ；

(5) 根据误差 ${\displaystyle \underset{i=1}{\overset{n}{\sum }}\left|f_{i,old}-f_{i,new}\right|}$ 是否小于给定的误差限，决定是否中止计算；或以新算出的 $f\left({\eta }_i\right)$ 重复

步骤(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)，解法简单，学生易于理解。

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

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

5. 结束语

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

基金项目

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

文章引用

林志敏,张永恒,张旭耀,王良璧. 外掠平板层流流动边界层微分方程积分解的实现
Implementation of Integration Solution to Boundary Layer Equation of Laminar Flow over Plane[J]. 创新教育研究, 2021, 09(02): 285-292. https://doi.org/10.12677/CES.2021.92045

参考文献