粘性流体力学上机实践作业
姓名：xx 学号：1439110xx
一、边界层方程的数值解
1问题描述
自编程序完成边界层方程：2f '''+f ⋅f ''=0 的数值解。 此即布拉修斯方程。对于实壁，f w =0，边界条件成为
'=0 η=0:f =f w =0；f '=f w η=∞:f '=1
其中：η=
y
ψ=u e (x ) ⋅g (x ) ⋅f (η) . =g (x ) 引入相似变量η和流函数ψ，使偏微分方程式变成常微分方程式。 求解η在[0,10]的范围内f , f ', f ''的数值解。
2求解过程
由于上述方程是非线性方程，可采用MATLAB 7.0软件来求解，步骤如下：
>> syms x1 x2 x3; y1=x1;y2=x2;y3=-x3*x2/2; u21 ;
[x,y,z]= DELGKT4_kuta121017003(y1,y2,y3,100,0,10,0,0,u21) function [x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,z0) format long ; h = (b-a)/N;
x = zeros(N+1,1);y = zeros(N+1,1);z = zeros(N+1,1);%二维全零矩阵, 长为n+1,宽为1
x(1)=f0;y(1) = y0;z(1) = z0; for i=2:N+1
k1 = subs(sym(f),findsym(sym(f)),y(i-1)); l1 = subs(sym(g),findsym(sym(g)),z(i-1));
=
newton_DELGK4_kuta_12101901(y1,y2,y3,100,0,10,0,0,0.2,0.4,0.00000001)
m1=subs(sym(e),[findsym(sym(e))],[z(i-1),x(i-1)]);%输出变量的先后顺序与定义的先后顺序一致
k2 = subs(sym(f),findsym(sym(f)),y(i-1)+h*l1/2); l2 = subs(sym(g),findsym(sym(g)),z(i-1)+h*m1/2);
m2=subs(sym(e),[findsym(sym(e))],[z(i-1)+h*m1/2,x(i-1)+h*k1/2]); k3 = subs(sym(f),findsym(sym(f)),y(i-1)+h*l2/2); l3 = subs(sym(g),findsym(sym(g)),z(i-1)+h*m2/2);
m3=subs(sym(e),[findsym(sym(e))],[z(i-1)+h*m2/2,x(i-1)+h*k2/2]); k4 = subs(sym(f),findsym(sym(f)),y(i-1)+h*l3); l4 = subs(sym(g),findsym(sym(g)),z(i-1)+h*m3);
m4=subs(sym(e),[findsym(sym(e))],[z(i-1)+h*m3,x(i-1)+h*k3]); x(i) = x(i-1)+h*(k1+2*k2+2*k3+k4)/6; y(i) = y(i-1)+h*(l1+2*l2+2*l3+l4)/6; z(i) = z(i-1)+h*(m1+2*m2+2*m3+m4)/6; end end
function u21 =
newton_DELGK4_kuta_12101901(f,g,e,N,a,b,f0,y0,u10,u20,eps) [x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u10);v10=y(N+1); [x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u20);v20=y(N+1); u21=u10+(1-v10)*(u20-u10)/(v20-v10);[x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u21); while abs(newton1(f,g,e,N,a,b,f0,y0,u21)-1)>eps v10=v20;u10=u20;u20=u21;
[x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u21); v20=y(N+1);
u21=u10+(1-v10)*(u20-u10)/(v20-v10); end
function f = newton1(f,g,e,N,a,b,f0,y0,u21)
[x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u21); f=y(N+1);
3结果分析
边界层方程解散点图如下：
二、Fluent 软件数值计算
1．问题描述
12. 950℃的铝液流经二维通道, 尺寸如图所示。试应用数值计算方法求解铝液流经二维通道时的流速分布和静压分布，并获得流经300cm 时的阻力损失。
2.
2.1建模和网格划分
（1）采用gambit 前处理软件进行建模和网格划分，铝液通道建模如图1所示。
图1 铝液通道建模
（2）用gambit 进行网格划分，由于模型较简单，直接对面采用四边形结构化网格划分。网格如图2所示。
图2 通道网格划分
（3）设置边界条件。
左边进口设置为inlet 名称，velocity-inlet 类型；右边出口设置名称为outlet ，pressure-outlet 类型；上边界是自由表面，设置名称为top ，pressure-inlet 类型，其余边界设置为wall 。
2.2 fluent求解
2.2.1 网格导入、检查
在fluent 中导入gambit 生成的mesh 文件，用scale 修改网格单位，检查网格的属性，是否有负数体积出现。
图3 网格展示和检查
2.2.2fluent 求解定义
定义材料，流动环境，求解器，选取模型根据题目示意，流动属于湍流流动，打开能量方程。选取k-e 双方程模型。设置速度入口为0.15m/s，初始温度为950℃。
打开残差控制，采取用inlet 进行初始化。进行迭代计算，设置迭代次数为5000次。
2.2.3迭代计算
大约在迭代次数到2350的时候，求解收敛。 残差控制图如下：
2.3 数据后处理及分析
静压分布云图如下：
总压云图：
速度云图：
流函数云图：
流函数等值线图：
速度矢量图：
绘制出口的沿y 方向速度剖面图：
2.4 阻力损失的计算
在fluent report中查询进出口截面的平均速度和平均静压，得到V 0=0.15m /s 1=0.04m /s ，P 0=0.43Pa ，P 1=0Pa
2则由伯努利方程知： +=h ρg 2g 0+sun
3ρg +122g 其中ρ=2.3⨯10kg /m , g =9.8
代入方程，则可以求出阻力损失
-3=1.085⨯10m
h sun 3
粘性流体力学上机实践作业
姓名：xx 学号：1439110xx
一、边界层方程的数值解
1问题描述
自编程序完成边界层方程：2f '''+f ⋅f ''=0 的数值解。 此即布拉修斯方程。对于实壁，f w =0，边界条件成为
'=0 η=0:f =f w =0；f '=f w η=∞:f '=1
其中：η=
y
ψ=u e (x ) ⋅g (x ) ⋅f (η) . =g (x ) 引入相似变量η和流函数ψ，使偏微分方程式变成常微分方程式。 求解η在[0,10]的范围内f , f ', f ''的数值解。
2求解过程
由于上述方程是非线性方程，可采用MATLAB 7.0软件来求解，步骤如下：
>> syms x1 x2 x3; y1=x1;y2=x2;y3=-x3*x2/2; u21 ;
[x,y,z]= DELGKT4_kuta121017003(y1,y2,y3,100,0,10,0,0,u21) function [x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,z0) format long ; h = (b-a)/N;
x = zeros(N+1,1);y = zeros(N+1,1);z = zeros(N+1,1);%二维全零矩阵, 长为n+1,宽为1
x(1)=f0;y(1) = y0;z(1) = z0; for i=2:N+1
k1 = subs(sym(f),findsym(sym(f)),y(i-1)); l1 = subs(sym(g),findsym(sym(g)),z(i-1));
=
newton_DELGK4_kuta_12101901(y1,y2,y3,100,0,10,0,0,0.2,0.4,0.00000001)
m1=subs(sym(e),[findsym(sym(e))],[z(i-1),x(i-1)]);%输出变量的先后顺序与定义的先后顺序一致
k2 = subs(sym(f),findsym(sym(f)),y(i-1)+h*l1/2); l2 = subs(sym(g),findsym(sym(g)),z(i-1)+h*m1/2);
m2=subs(sym(e),[findsym(sym(e))],[z(i-1)+h*m1/2,x(i-1)+h*k1/2]); k3 = subs(sym(f),findsym(sym(f)),y(i-1)+h*l2/2); l3 = subs(sym(g),findsym(sym(g)),z(i-1)+h*m2/2);
m3=subs(sym(e),[findsym(sym(e))],[z(i-1)+h*m2/2,x(i-1)+h*k2/2]); k4 = subs(sym(f),findsym(sym(f)),y(i-1)+h*l3); l4 = subs(sym(g),findsym(sym(g)),z(i-1)+h*m3);
m4=subs(sym(e),[findsym(sym(e))],[z(i-1)+h*m3,x(i-1)+h*k3]); x(i) = x(i-1)+h*(k1+2*k2+2*k3+k4)/6; y(i) = y(i-1)+h*(l1+2*l2+2*l3+l4)/6; z(i) = z(i-1)+h*(m1+2*m2+2*m3+m4)/6; end end
function u21 =
newton_DELGK4_kuta_12101901(f,g,e,N,a,b,f0,y0,u10,u20,eps) [x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u10);v10=y(N+1); [x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u20);v20=y(N+1); u21=u10+(1-v10)*(u20-u10)/(v20-v10);[x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u21); while abs(newton1(f,g,e,N,a,b,f0,y0,u21)-1)>eps v10=v20;u10=u20;u20=u21;
[x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u21); v20=y(N+1);
u21=u10+(1-v10)*(u20-u10)/(v20-v10); end
function f = newton1(f,g,e,N,a,b,f0,y0,u21)
[x,y,z]= DELGKT4_kuta121017003(f,g,e,N,a,b,f0,y0,u21); f=y(N+1);
3结果分析
边界层方程解散点图如下：
二、Fluent 软件数值计算
1．问题描述
12. 950℃的铝液流经二维通道, 尺寸如图所示。试应用数值计算方法求解铝液流经二维通道时的流速分布和静压分布，并获得流经300cm 时的阻力损失。
2.
2.1建模和网格划分
（1）采用gambit 前处理软件进行建模和网格划分，铝液通道建模如图1所示。
图1 铝液通道建模
（2）用gambit 进行网格划分，由于模型较简单，直接对面采用四边形结构化网格划分。网格如图2所示。
图2 通道网格划分
（3）设置边界条件。
左边进口设置为inlet 名称，velocity-inlet 类型；右边出口设置名称为outlet ，pressure-outlet 类型；上边界是自由表面，设置名称为top ，pressure-inlet 类型，其余边界设置为wall 。
2.2 fluent求解
2.2.1 网格导入、检查
在fluent 中导入gambit 生成的mesh 文件，用scale 修改网格单位，检查网格的属性，是否有负数体积出现。
图3 网格展示和检查
2.2.2fluent 求解定义
定义材料，流动环境，求解器，选取模型根据题目示意，流动属于湍流流动，打开能量方程。选取k-e 双方程模型。设置速度入口为0.15m/s，初始温度为950℃。
打开残差控制，采取用inlet 进行初始化。进行迭代计算，设置迭代次数为5000次。
2.2.3迭代计算
大约在迭代次数到2350的时候，求解收敛。 残差控制图如下：
2.3 数据后处理及分析
静压分布云图如下：
总压云图：
速度云图：
流函数云图：
流函数等值线图：
速度矢量图：
绘制出口的沿y 方向速度剖面图：
2.4 阻力损失的计算
在fluent report中查询进出口截面的平均速度和平均静压，得到V 0=0.15m /s 1=0.04m /s ，P 0=0.43Pa ，P 1=0Pa
2则由伯努利方程知： +=h ρg 2g 0+sun
3ρg +122g 其中ρ=2.3⨯10kg /m , g =9.8
代入方程，则可以求出阻力损失
-3=1.085⨯10m
h sun 3