﻿ 外掠平板层流流动边界层微分方程积分解的实现

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

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. 引言

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)

$\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)

$\left\{\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)

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

$\eta =\frac{y}{\text{ }\sqrt{\nu x/{u}_{s}}}$ (10)

$f\left(\eta \right)=\frac{\psi }{\text{ }\sqrt{{u}_{s}\nu x}}$ (11)

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

${{\theta }^{″}}_{1}+\frac{1}{2}\left(m+1\right)Pr\text{ }f{{\theta }^{\prime }}_{1}=0$ (13)

$\left\{\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)

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

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

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

${{\theta }^{″}}_{1}+\frac{1}{2}Pr\text{ }f{{\theta }^{\prime }}_{1}=0$ (17)

(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) 判断误差 $|{f}^{\prime }\left(\infty \right)-{f}_{2}^{*}{}^{\prime }\left(\infty \right)|$$|{\theta }_{1}\left(\infty \right)-{\theta }_{12}^{*}\left(\infty \right)|$ 是否小与给定误差限，若小与给定的误差限，则停止迭代，否则转第(2)步。

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

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

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

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

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

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

${\eta }_{i}=\frac{{\eta }_{\mathrm{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) 根据误差 $\underset{i=1}{\overset{n}{\sum }}|{f}_{i,old}-{f}_{i,new}|$ 是否小于给定的误差限，决定是否中止计算；或以新算出的 $f\left({\eta }_{i}\right)$ 重复

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;

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

5. 结束语

NOTES

*通讯作者。

[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