流体力学——平板边界层编程

合集下载

相关主题

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

对于本次编程编程作业，小组运用matlab 和c++两种程序对平板边界层问题和绕过楔形体边界层流动问题进行分析研究。以下是运用matlab 解决问题的过程。一、平板边界层问题

该问题可以归结为在已知边界层条件下解一个高阶微分方程,即解0''5.0'''=+ff f 。Matlab 提供了解微分方程的方法，运用换元法将高阶微分方程降阶，然后运用“ode45”函数进行求解。函数其难点在于如何将边界条件中1',→∞→f η运用好，由四阶龙格-库塔方法知其核心是换元试算匹配，故在运用函数时通过二分法实现1',→∞→f η是可行的。程序如下：

第一问m 函数

function dy = rigid(t,y) dy = zeros(3,1); dy(1) = y(2); dy(2) = y(3);

dy(3) = -0.5*y(1)*y(3);

%第一问main 程序

[T,Y] = ode45('rigid',[0 5],[0 0 0]);

plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') %二分法试算f ’’的初始值以满足f ’趋向无穷时的边界条件，图像上可以清晰看出f ’无穷时的结果 >> [T,Y] = ode45('rigid',[0 5],[0 0 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> [T,Y] = ode45('rigid',[0 5],[0 0 0.5]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> [T,Y] = ode45('rigid',[0 5],[0 0 0.25]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')

>> [T,Y] = ode45('rigid',[0 5],[0 0 0.375]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> grid on

>> [T,Y] = ode45('rigid',[0 5],[0 0 0.3125]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') >> grid on

>> [T,Y] = ode45('rigid',[0 5],[0 0 0.34375]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on

>> [T,Y] = ode45('rigid',[0 5],[0 0 0.328125]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on

>> [T,Y] = ode45('rigid',[0 10],[0 0 0.328125]);%当f ’’为0.328125时，逼近结果已经很好，在0到5的变化范围内已经非常接近精确解 plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on

>> [T,Y] = ode45('rigid',[0 5],[0 0 0.335975]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') grid on

选取f’’=0.335975时的数据展示，T代表η的变化，Y的第一二三列代表f

f

,f

,'

''

为形象展示所得结果，如图-1所示：

二、楔形绕流边界层问题

解决绕过楔形的边界层流动问题与平板问题的不同之处在于微分方程变为0

α，其中α可取1，并不失其普遍性，由于β

+f

fβ

+

ff

'''2=

)

'

1(

''

-

小于-0.199时会发生分离，该微分方程不再适用，故β取值大于-0.199。其解法与平板绕流类似，在平板绕流基础上增加β的变化即可。

程序如下：

第二问m函数

function dz = rigid1(t,z)

dz = zeros(3,1);

t=-0.199; %其中t表示β的变化，可由-0.199至无穷大

dz(1) = z(2);

dz(2) = z(3);

dz(3) = -z(1)*z(3)-t*(1-z(2)^2);

第二问main函数

[T,Y] = ode45('rigid1',[0 5],[0 0 0]);

plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') %二分法试算f’’的初始值以满足f’趋向无穷时的边界条件，图像上可以清晰看出f’无穷时的结果

>> [T,Y] = ode45('rigid1',[0 5],[0 0 1]);