利用有限体积算法三阶迎风型QUICK 离散格式求解
二维不可压缩黏性流体方腔流动问题
1.二维不可压缩黏性流体方腔流动问题
二维不可压缩黏性流体方腔流动（cavity flow ）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为=1.0=1.0p ρ、它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。
2. 基本方程组、初始条件和边界条件
设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动
N - S
()()u v x y x x φρφρφμ??????
+= ?????????
） 其中u 为变量φ在水平x 方向的流速，v 为φ在垂直y 方向的流速，μ为黏度，S φ为源项。源项中不仅包含压力梯度项，也包含时间导数项。
初始条件：方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动。 边界条件：
流动速度u v 、均可采用无滑移边界条件，压力p 采用自由输出边界条件。
3.计算网格划分和控制体单元与节点定义
采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。
图F.1二维不可压缩黏性
方腔流问题示意图
节点P 所在主控制单元如图F.2中有阴影部分所示。在x 方向与节点P 相邻的节点为W 和E ，在y 方向与节点P 相邻的节点为S 和N ，主控制单元界面分
别为w s e n 、、、
。压力p 和速度u v 、分别在三套不同网格中如图F.3中有阴影部分所示。
4.有限体积算法三阶迎风型QUICK 离散格式
对方程（F.1）在图F.2所示节点P 所在控制体单元内积分，有：
()()V V V V V u dV v dV dV dV S dV x y x x y y φφφρφρφμμ???????
????????+=++ ? ????????????????（F.2） 由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积V ?仅是面积，而它的边界是长度。设 ,w e s n A A y A A x ==?==?，利用Gauss 定理，可将方程（F.2）改写成如下有限体积算法离散格式：
()()()()e w n s e w n s u A u A v A v A A A A A S x y
x x y y φρφρφρφρφφφφφμμμμ????-+-????
?????????
???=-+-+?? ? ? ? ??????
??????? （F.3） 对上式中e w n s
x x y y φφφφ????
???????? ? ? ? ?????????????、、
、采用一阶向前差分近似，则有： 图F.2方腔流动计算网格、 控制体单元和节点示意图
图F.3计算采用的交错网格示意图
,,P W E P e w
N P P S
n s
x x x x y y y y φφφφφφφφφφφφ--??????
== ? ?
????????????--??=
= ? ????????? （F.4）
同时记：
()()()(),,e e w w
e w n n s s
n s F u A F u A F v A F v A ρρρρ==== （F.5）
,,,e e w w e w PE PW
n n s s
n s PN PS
A A
D D x x A A D D y y μμδδμμδδ=
===
（F.6）
则可由式（F.2）写成：
()()()()e e w w n n s s e E P w P W n N P s P S F F F F D D D D S x y
φφφφφφφφφφφφφ-+-=---+
---+?? （F.7）
式中P E W N S e w n s D D D D φφφφφ、、、、、、、、都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量e e w w n n s s F F F F φφφφ、、、，就可以求出节点上未知量P E W N S φφφφφ、、、、。
为了便于讨论，现对一维对流扩散方程的三阶迎风型QUICK 离散格式进行分析：在三阶迎风型QUICK 离散格式中，计算主控制单元界面上流动量φ需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。
图F.4三阶迎风型QUICK 离散格式示意图
当0,0e w u u >>时，通过WW 、W 和P 三个节点值拟合曲线来计算主控制单元左侧界面参数w φ。通过节点W 、P 和E 三个节点值拟合曲线来计算主控制单元右侧界面参数e φ。当0,0e w u u <<，则分别通过节点W 、P 、E 和P 、E 、EE 三个节点值计算主控制单元左、右两侧界面参数w φ和e φ。根据上述计算原则，可以得到界面参数w φ计算公式如下：
当0w u >时，界面参数w φ计算公式为：
636888
ww wW wP wWW φφφφ=+- （F.8a ）
当0e u >时，界面参数w φ计算公式为：
636
888
e P E W φφφφ=+- （F.8b ）
对于一维无源项一维对流扩散方程三阶迎风型QUICK 离散格式： 当0,0e w u u >>时，三阶迎风型QUICK 离散格式为：
P P W W E E WW WW a a a a φφφφ=++ （F.9） 其中
()
31
883
81
8
W w w e
E e e
WW w
P W E WW e w a D F F a D F a F a a a a F F =++=-=-=+++- （F.9a ） 同理，若0,0e w u u <<，三阶迎风型QUICK 离散格式为：
P P W W E E EE EE a a a a φφφφ=++ （F.10） 其中
()
3
,
861
,
8818
W w w E e e w WW e
P W E EE e w a D F a D F F a F a a a a F F =+=--==+++- （F.10a ） 将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型QUICK 离散格式：
P P W W E E EE EE EE EE a a a a a φφφφφ=+++
（F.11） 其中
()()()()()
613
18881
8361
118881
18
W w w w e e w w
WW w w
E e e e e e w w EE e e
P W E WW EE e w a D F F F a F a D F F F a F a a a a a F F αααααααα=+++-=-=-----=-=++++- （F.11a ）
式中
01000100w w w w e e e e F F F F αααα>=<=>=<=当时，；
当时，；
当时，；当时，。
（F.11b ）
同理，可以得到带有源项的二维对流扩散方程三阶迎风型QUICK 离散格式为：
P P W W E E EE EE EE EE S S N N SS SS NN NN a a a a a a a a a S V
φφφφφφφφφφ=++++
++++? （F.12）
其中S φ为有限体积算法中源项平均值。式中各个系数为：
()()()613
18881
8361
11888W w w w e e w w
WW w w
E e e e e e w w
a D F F F a F a D F F F ααααααα=+++-=-=----- ()()()()()1
18
613
188818361
118881
18EE e e S s s s n n s s
SS s s N n n n n n s s
NN n n
a F a D F F F a F a D F F F a F ααααααααα=
-=+++-=-=-----=-
()P W E WW EE S N SS NN e w n s a a a a a a a a a F F F F =++++++++-+- （F.12a ）
式中
0100,,,k k k k F F k w e s n
αα>=<==当时，；
当时，。 （F.12b ）
源项S φ为：
()u p
S t x
φρ??=-
-?? （F.13） 若把()n
u ρ表示n t 时刻动量，()1
n u ρ+表示1n t +时刻动量，则可以得到源项S φ离散
格式为：
()()()1
n n
P
P
e w V
u u S dV x y p p y t
φρρ+?-?=-
??--??? （F.14）
最后，得到有限体积算法二维对流扩散方程三阶迎风型QUICK 离散格式：
()()()1
n n n n n
P P W W E E S S N N n n
n
n P
P
e
w
a u a u a u a u a u u u x y p p
y
t
ρρ+=+++-
-??--?? （F.15）
式中系数k a 为一阶迎风格式中各对应系数。
5.计算结果分析
利用三阶迎风型QUICK 离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。图F.5是不同雷诺数Re 条件下采用三阶迎风型QUICK 离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。
计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。这表明有限体积算法三阶迎风型QUICK 离散格式具有相当高的计算精度。
图F.5不同雷诺数Re 条件下采用三阶迎风型QUICK
离散格式计算二维不可压缩黏性方腔流动的计算结果
Re =100 Re =1 000
Re =5 000 Re =10 000
从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并不对称。这是因为，方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动时，在方腔上壁面两侧的两个顶角处不再满足边界条件，这是一个带奇性的方腔流动。
计算结果表明，方腔流动和雷诺数有关，随着雷诺数Re增加，计算精度在降低。当雷诺数Re较低时，方腔下壁面的两个底角附近的二次小涡十分清晰，随着雷诺数Re的增加，二次小涡变得越来越模糊。由于三阶迎风型QUICK离散格式计算精度较高，因此三阶迎风型QUICK离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。
using System;
using System.Collections.Generic;
using https://www.doczj.com/doc/4711934516.html,ponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.IO;
using System.Windows.Forms;
using https://www.doczj.com/doc/4711934516.html,;
using System.Collections;
namespace FVM_计算方腔流场
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
public int MX= 30;
public int MY =30;
public double Re =100.00;
public double Ha = 10.00;
public double Ri = 1.0;
public double Pr = 0.733;
public double dt =0.0005;
public double c2 =2.25;
public double[,] u, v, p, t,uo, vo, po,to, un, vn, pn,tn,streamline;
//public double dx, dy;
FileStream fileResult, file1,f4;
StreamWriter fileResultWriter, fileW1,fw4;
public double dx, dy, err, value,x,y;
private void Form1_Load(object sender, EventArgs e)
{
u = new double[MX + 1, MY + 2];
v = new double[MX + 2, MY + 1];
p = new double[MX + 2, MY + 2];
t = new double[MX + 2, MY + 1];
uo = new double[MX + 1, MY + 2];
vo = new double[MX + 2, MY + 1];
po = new double[MX + 2, MY + 2];
to = new double[MX + 2, MY + 1];
un = new double[MX + 1, MY + 2];
vn = new double[MX + 2, MY + 1];
pn = new double[MX + 2, MY + 2];
tn = new double[MX + 2, MY + 1];
streamline = new double[MX + 1, MY + 1];
dx = 1.0 / MX;
dy = 1.0 / MY;
fileResult = new FileStream(" X_Matrix.txt", FileMode.OpenOrCreate); fileResultWriter = new StreamWriter(fileResult);
file1 = new FileStream(" file1.txt", FileMode.OpenOrCreate);
fileW1 = new StreamWriter(file1);
f4 = new FileStream(" f4.txt", FileMode.OpenOrCreate);
fw4 = new StreamWriter(f4);
//double dx1, dy1,
//init(u, v, p, dx, dy);
//MessageBox.Show(u[5, MY + 1].ToString());
}
public double max(double a, double b)
{
if (a < b)
return b;
else
return a;
}
public double alfa(double x)
{
if (x >= 0)
return 1.0;
else
return 0.0;
}
void init(double[,] u, double[,] v, double[,] p, double[,] t, double dx, double dy) {
//u = new double[MX + 1, MY + 2];
//v = new double[MX + 2, MY + 1];
//p = new double[MX + 2, MY + 2];
int i,j;
//dx=1.0/MX;
//dy=1.0/MY;
for(i=0;i<=MX;i++)
{
for(j=0;j<=MY+1;j++)
{
u[i,j]=0.0;
if(j==MY+1) u[i,j]=4.0/3.0;
if(j==MY) u[i,j]=2.0/3.0;
}
}
for(i=0;i<=MX+1;i++)
for(j=0;j<=MY;j++)
v[i,j]=0.0;
for(i=0;i<=MX+1;i++)
for(j=0;j<=MY+1;j++)
p[i,j]=1.0;
for (i = 0; i <= MX+1; i++)
{
for (j = 0; j <= MY; j++)
{
t[i, j] = 0.0;
if (i == MX + 1) t[i, j] = 4.0 / 3.0;
if (i == MX) t[i, j] = 2.0 / 3.0;
}
}
}
//
// 一阶迎风型离散格式
// 二维的三阶迎风型离散格式为9点格式，因此有两层边界网格需要
//处理，本程序采用一阶迎风型离散格式处理内层，用物理边界条件处理外层。
// 入口:u、v、p、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号；
// 出口：un，新的x方向速度。
//
void upwind_u(double[,] u, double[,] v, double[,] p, double[,] t, double[,] un, double dx, double dy, int i, int j)
{
//u = new double[MX + 1, MY + 2];
//v = new double[MX + 2, MY + 1];
//p = new double[MX + 2, MY + 2];
//un = new double[MX + 1, MY + 2];
double aw, ae, a_s, an, df, ap, miu;
miu = 1.0 / Re;
aw = miu + max(0.5 * (u[i - 1, j] + u[i, j]) * dy, 0.0);
ae = miu + max(0.0, -0.5 * (u[i, j] + u[i + 1, j]) * dy);
a_s = miu + max(0.5 * (v[i, j - 1] + v[i + 1, j - 1]) * dx, 0.0);
an = miu + max(0.0, -0.5 * (v[i, j] + v[i + 1, j]) * dx);
df = 0.5 * (u[i + 1, j] - u[i - 1, j]) * dy + 0.5 * (v[i, j] + v[i + 1, j] - v[i, j - 1] - v[i + 1, j - 1]) * dx;
ap = aw + ae + a_s + an + df;
un[i, j] = u[i, j] + dt / dx / dy * (-ap * u[i, j] + aw * u[i - 1, j] + ae * u[i + 1, j] + a_s * u[i, j - 1] + an * u[i, j + 1]) - dt * (p[i + 1, j] - p[i, j]) / dx;
}
//入口: u、v、p、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号；
//出口：vn，新的y方向速度。
void upwind_v(double[,] u, double[,] v, double[,] p, double[,] t, double[,] vn, double dx, double dy, int i, int j)
{
//u = new double[MX + 1, MY + 2];
//v = new double[MX + 2, MY + 1];
//p = new double[MX + 2, MY + 2];
//vn = new double[MX + 2, MY + 1];
double aw, ae, a_s, an, df, ap, miu;
miu = 1.0 / Re;
aw = miu + max(0.5 * (u[i - 1, j] + u[i - 1, j + 1]) * dy, 0.0);
ae = miu + max(0.0, -0.5 * (u[i, j] + u[i, j + 1]) * dy);
a_s = miu + max(0.5 * (v[i, j - 1] + v[i, j]) * dx, 0.0);
an = miu + max(0.0, -0.5 * (v[i, j] + v[i, j + 1]) * dx);
df = 0.5 * (u[i, j] + u[i, j + 1] - u[i - 1, j] - u[i - 1, j + 1]) * dy + 0.5 * (v[i, j + 1] - v[i, j - 1]) * dx;
ap = aw + ae + a_s + an + df;
vn[i, j] = v[i, j] + dt / dx / dy * (-ap * v[i, j] + aw * v[i - 1, j] + ae * v[i + 1, j] + a_s * v[i, j - 1] + an * v[i, j + 1])
- dt * (p[i, j + 1] - p[i, j]) / dy + dt * Ri * t[i, j] - dt * Ha * Ha * v[i, j] / Re;
}
//入口: u、v、p、t、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号；
//出口：tn，新的温度。
void upwind_t(double[,] u, double[,] v, double[,] p, double[,] t, double[,] tn, double dx, double dy, int i, int j)
{
//u = new double[MX + 1, MY + 2];
//v = new double[MX + 2, MY + 1];
//p = new double[MX + 2, MY + 2];
//vn = new double[MX + 2, MY + 1];
double aw, ae, a_s, an, df, ap, miu;
miu = 1.0 / Re/Pr;
aw = miu + max(0.5 * (u[i - 1, j] + u[i - 1, j + 1]) * dy, 0.0);
ae = miu + max(0.0, -0.5 * (u[i, j] + u[i, j + 1]) * dy);
a_s = miu + max(0.5 * (v[i, j - 1] + v[i, j]) * dx, 0.0);
an = miu + max(0.0, -0.5 * (v[i, j] + v[i, j + 1]) * dx);
df = 0.5 * (u[i, j] + u[i, j + 1] - u[i - 1, j] - u[i - 1, j + 1]) * dy + 0.5 * (v[i, j + 1] - v[i, j - 1]) * dx;
ap = aw + ae + a_s + an + df;
tn[i, j] = t[i, j] + dt / dx / dy * (-ap * t[i, j] + aw * t[i - 1, j] + ae * t[i + 1, j] + a_s * t[i, j - 1] + an * t[i, j + 1]);
}
//三阶迎风型离散格式
//入口:u、v、p、dx、dy，当前速度、压强，空间步长；
//出口：un、vn，新的速度。
void quick(double[,] u, double[,] v, double[,] p, double[,] t, double[,] un, double[,] vn, double[,] tn, double dx, double dy)
{
//u = new double[MX + 1, MY + 2];
//v = new double[MX + 2, MY + 1];
//p = new double[MX + 2, MY + 2];
//un = new double[MX + 1, MY + 2];
//vn = new double[MX + 2, MY + 1];
double miu,miu1, fw, fe, fs, fn, df, aw, ae, a_s, an, aww, aee, ass, ann, ap;
int i, j;
miu = 1.0 / Re;
miu1 = 1.0 / Re/Pr ;
for(i=2;i<=MX-2;i++)
{
for(j=2;j<=MY-1;j++)
{
fw = 0.5 * (u[i - 1, j] + u[i, j]) * dy;
fe = 0.5 * (u[i, j] + u[i + 1, j]) * dy;
fs = 0.5 * (v[i, j - 1] + v[i + 1, j - 1]) * dx;
fn = 0.5 * (v[i, j] + v[i + 1, j]) * dx;
df = fe - fw + fn - fs;
aw = miu + 0.750 * alfa(fw) * fw + 0.125 * alfa(fe) * fe + 0.375 * (1.0 - alfa(fw)) * fw;
aww = -0.125 * alfa(fw) * fw;
ae = miu - 0.375 * alfa(fe) * fe - 0.750 * (1.0 - alfa(fe)) * fe - 0.125 * (1.0 - alfa(fw)) * fw;
aee = 0.125 * (1.0 - alfa(fe)) * fe;
a_s = miu + 0.750 * alfa(fs) * fs + 0.125 * alfa(fn) * fn + 0.375 * (1.0 - alfa(fs)) * fs;
ass = -0.125 * alfa(fs) * fs;
an = miu - 0.375 * alfa(fn) * fn - 0.750 * (1.0 - alfa(fn)) * fn - 0.125 * (1.0 - alfa(fs)) * fs;
ann = 0.125 * (1.0 - alfa(fn)) * fn;
ap = aw + ae + a_s + an + aww + aee + ass + ann + df;
//aw、ae、as、an...均为有限体积算法中各项系数，详见前文三阶迎风型QUICK离散格式。
un[i, j] = u[i, j] + dt / dx / dy * (-ap * u[i, j] + aw * u[i - 1, j] + ae * u[i + 1, j] + a_s * u[i, j - 1]+ an * u[i, j + 1] + aww * u[i - 2, j]
+ aee * u[i + 2, j] + ass * u[i, j - 2] + ann * u[i, j + 2])- dt * (p[i + 1, j] - p[i, j]) / dx;
}
}
;
//j=1;
for (i = 2; i <= MX - 2; i++)
{
upwind_u(u, v, p,t, un, dx, dy, i, 1);
}
// upwind_u(u, v, p, un, dx, dy, i, j);
//j=MY;
for (i = 2; i <= MX - 2; i++)
{
upwind_u(u, v, p,t, un, dx, dy, i, MY);
}
// upwind_u(u, v, p, un, dx, dy, i, j);
//i=1;
for (j = 1; j <= MY; j++)
{
upwind_u(u, v, p,t, un, dx, dy, 1, j);
}
// upwind_u(u, v, p, un, dx, dy, i, j);
//i=MX-1;
for (j = 1; j <= MY; j++)
{
upwind_u(u, v, p, t, un, dx, dy, (MX - 1), j);
}
//upwind_u(u, v, p, un, dx, dy, i, j);
//内层边界由一阶迎风型离散格式得到
for(i=1;i<=MX-1;i++)
{
un[i, 0] = -un[i, 1];
un[i, MY + 1] = 2.0 - un[i, MY];
}
for(j=0;j<=MY+1;j++)
{
un[0, j] = 0.0;
un[MX, j] = 0.0;
}
//外层边界条件按物理边界条件给出
for(i=2;i<=MX-1;i++)
{
for(j=2;j<=MY-2;j++)
{
fw = 0.5 * (u[i - 1, j] + u[i - 1, j + 1]) * dy; fe = 0.5 * (u[i, j] + u[i, j + 1]) * dy;
fs = 0.5 * (v[i, j - 1] + v[i, j]) * dx;
fn = 0.5 * (v[i, j] + v[i, j + 1]) * dx;
df = fe - fw + fn - fs;
aw = miu + 0.750 * alfa(fw) * fw + 0.125 * alfa(fe) * fe + 0.375 * (1.0 - alfa(fw)) * fw;
aww = -0.125 * alfa(fw) * fw;
ae = miu - 0.375 * alfa(fe) * fe - 0.750 * (1.0 - alfa(fe)) * fe - 0.125 * (1.0 - alfa(fw)) * fw;
aee = 0.125 * (1.0 - alfa(fe)) * fe;
a_s = miu + 0.750 * alfa(fs) * fs + 0.125 * alfa(fn) * fn + 0.375 * (1.0 - alfa(fs)) * fs;
ass = -0.125 * alfa(fs) * fs;
an = miu - 0.375 * alfa(fn) * fn - 0.750 * (1.0 - alfa(fn)) * fn - 0.125 * (1.0 - alfa(fs)) * fs;
ann = 0.125 * (1.0 - alfa(fn)) * fn;
ap = aw + ae + a_s + an + aww + aee + ass + ann + df;
vn[i, j] = v[i, j] + dt / dx / dy * (-ap * v[i, j] + aw * v[i - 1, j] + ae * v[i + 1, j] + a_s * v[i, j - 1]
+ an * v[i, j + 1] + aww * v[i - 2, j] + aee * v[i + 2, j] + ass * v[i, j - 2] + ann * v[i, j + 2])
- dt * (p[i, j + 1] - p[i, j]) / dy + dt * Ri * t[i, j] - dt * Ha * Ha * v[i, j] / Re;
}
}
//
//j=1;
for (i = 2; i <= MX - 1; i++)
{
upwind_v(u, v, p,t, vn, dx, dy, i, 1);
}
// upwind_v(u, v, p, vn, dx, dy, i, j);
//j=MY-1;
for (i = 2; i <= MX - 1; i++)
{
upwind_v(u, v, p, t, vn, dx, dy, i, (MY - 1));
}
// upwind_v(u, v, p, vn, dx, dy, i, j);
//i=1;
for (j = 1; j <= MY - 1; j++)
{
upwind_v(u, v, p, t, vn, dx, dy, 1, j);
}
// upwind_v(u, v, p, vn, dx, dy, i, j);
//i=MX;
for (j = 1; j <= MY - 1; j++)
{
upwind_v(u, v, p, t, vn, dx, dy, MX, j);
}
//upwind_v(u, v, p, vn, dx, dy, i, j);
//
for(i=1;i<=MX;i++)
{
vn[i, 0] = 0.0;
vn[i, MY] = 0.0;
}
for(j=0;j<=MY;j++)
{
vn[0, j] = -vn[1, j];
vn[MX + 1, j] = -vn[MX, j];
}
for (i = 2; i <= MX - 1; i++)
{
for (j = 2; j <= MY - 2; j++)
{
fw = 0.5 * (u[i - 1, j] + u[i - 1, j + 1]) * dy;
fe = 0.5 * (u[i, j] + u[i, j + 1]) * dy;
fs = 0.5 * (v[i, j - 1] + v[i, j]) * dx;
fn = 0.5 * (v[i, j] + v[i, j + 1]) * dx;
df = fe - fw + fn - fs;
aw = miu1 + 0.750 * alfa(fw) * fw + 0.125 * alfa(fe) * fe + 0.375 * (1.0 - alfa(fw)) * fw;
aww = -0.125 * alfa(fw) * fw;
ae = miu1 - 0.375 * alfa(fe) * fe - 0.750 * (1.0 - alfa(fe)) * fe - 0.125 * (1.0 - alfa(fw)) * fw;
aee = 0.125 * (1.0 - alfa(fe)) * fe;
a_s = miu1 + 0.750 * alfa(fs) * fs + 0.125 * alfa(fn) * fn + 0.375 * (1.0 - alfa(fs)) * fs;
ass = -0.125 * alfa(fs) * fs;
an = miu1 - 0.375 * alfa(fn) * fn - 0.750 * (1.0 - alfa(fn)) * fn - 0.125 * (1.0 - alfa(fs)) * fs;
ann = 0.125 * (1.0 - alfa(fn)) * fn;
ap = aw + ae + a_s + an + aww + aee + ass + ann + df;
tn[i, j] = t[i, j] + dt / dx / dy * (-ap * t[i, j] + aw * t[i - 1, j] + ae * t[i + 1, j] + a_s * t[i, j - 1]
+ an * t[i, j + 1] + aww * t[i - 2, j] + aee * t[i + 2, j] + ass * t[i, j - 2] + ann * t[i, j + 2]);
}
}
//
//j=1;
for (i = 2; i <= MX - 1; i++)
{
upwind_t(u, v, p, t, tn, dx, dy, i, 1);
}
// upwind_v(u, v, p, vn, dx, dy, i, j);
//j=MY-1;
for (i = 2; i <= MX - 1; i++)
{
upwind_t(u, v, p, t, tn, dx, dy, i, (MY - 1)); }
// upwind_v(u, v, p, vn, dx, dy, i, j);
//i=1;
for (j = 1; j <= MY - 1; j++)
{
upwind_t(u, v, p, t, tn, dx, dy, 1, j);
}
// upwind_v(u, v, p, vn, dx, dy, i, j);
//i=MX;
for (j = 1; j <= MY - 1; j++)
{
upwind_t(u, v, p, t, tn, dx, dy, MX, j);
}
//upwind_v(u, v, p, vn, dx, dy, i, j);
//
for (i = 1; i <= MX; i++)
{
tn[i, 0] = tn[i, 1];
tn[i, MY] = tn[i, MY-1];
}
for (j = 0; j <= MY; j++)
{
tn[0, j] = -tn[1, j];
tn[MX + 1, j] = 2-tn[MX, j];
}
}
//更新压强
//入口: un、vn、p、dx、dy，新的速度，当前压强，空间步长；
//出口: pn，新的压强。
void getp(double[, ] un , double[,] vn, double[,] p, double[,] pn, double dx, double dy) {
int i, j;
for (i = 1; i <= MX; i++)
for (j = 1; j <= MY; j++)
pn[i, j] = p[i, j] - dt * c2 / dx * (un[i, j] - un[i - 1, j] + vn[i, j] - vn[i, j - 1]);
for (i = 1; i <= MX; i++)
{
pn[i, 0] = pn[i, 1];
pn[i, MY + 1] = pn[i, MY];
}
for (j = 0; j <= MY + 1; j++)
{
pn[0, j] = pn[1, j];
pn[MX + 1, j] = pn[MX, j];
}
}
private void button1_Click(object sender, EventArgs e)
{
int i, j, step;
init(u, v, p,t, dx, dy);
MessageBox.Show(dx.ToString());
err = 100.0;
step = 0;
while (err > 0.0001 ) //err<1e-5，定常解判据；step，限制迭代次数
{
step++;
err = 0.0;
quick(u, v, p,t, un, vn,tn, dx, dy);
//for (int i0 = 0; i0 < MX + 1; i0++)
//{
// for (int j0 = 0; j0 < MY + 2; j0++)
// {
// fileResultWriter.Write(Math.Round(un[i0, j0], 3).ToString() + " ");
// }
// fileResultWriter.Write("\r\n");
//}
//fileResultWriter.Write(" " + "\r\n");
getp(un, vn, p, pn, dx, dy);
for (i = 0; i <= MX; i++)
{
for (j = 0; j <= MY + 1; j++)
{
value = Math.Abs(un[i,j] - u[i,j]) / dt;
if (value > err) err = value;
u[i,j] = un[i,j];
}
}
for (i = 0; i <= MX + 1; i++)
{
for (j = 0; j <= MY; j++)
{
value = Math.Abs(vn[i,j] - v[i,j]) / dt;
if (value > err) err = value;
v[i,j] = vn[i,j];
}
}
for (i = 0; i <= MX + 1; i++)
{
for (j = 0; j <= MY + 1; j++)
{
value = Math.Abs(pn[i,j] - p[i,j]) / c2 / dt;
if (value > err) err = value;
p[i,j] = pn[i,j];
}
}
for (i = 0; i <= MX + 1; i++)
{
for (j = 0; j <= MY; j++)
{
value = Math.Abs(tn[i, j] - t[i, j]) / dt;
if (value > err) err = value;
t[i, j] = tn[i, j];
}
}
fileW1.Write(step.ToString() + " "+ Math.Round(err, 6).ToString() + "\r\n");
//for (int i0 = 0; i0 < MX + 1; i0++)
//{
// for (int j0 = 0; j0 < MY + 1; j0++)
// {
// fileResultWriter.Write(Math.Round(pn[i0, j0], 3).ToString() + " ");
// }
// fileResultWriter.Write("\r\n");
//}
//fileResultWriter.Write(" " + "\r\n");
}
for (i = 0; i <= MX; i++)
{
for (j = 0; j <= MY; j++)
{
uo[i,j] = 0.5 * (u[i,j] + u[i,j + 1]);
vo[i,j] = 0.5 * (v[i,j] + v[i + 1,j]);
to[i, j] = 0.5 * (t[i, j] + t[i + 1, j]);
po[i,j] = 0.25 * (p[i,j] + p[i + 1,j] + p[i,j + 1] + p[i + 1,j + 1]);
}
}
for (int i0 = 0; i0 < MX + 1; i0++)
{
for (int j0 = 0; j0 < MY + 1; j0++)
{
fileResultWriter.Write(Math.Round(po[i0, j0], 3).ToString() + " ");
}
fileResultWriter.Write("\r\n");
}
fileResultWriter.Write(" " + "\r\n");
for (int i0 = 0; i0 < MX + 1; i0++)
{
for (int j0 = 0; j0 < MY + 1; j0++)
{
fileResultWriter.Write(Math.Round(uo[i0, j0], 3).ToString() + " ");
}
fileResultWriter.Write("\r\n");
}
fileResultWriter.Write(" " + "\r\n");
化工原理习题及答案 第五章传热 姓名 ____________ 班级 ____________ 学号 _____________ 成绩 ______________ 一、填空题： 1.（ 6 分）某大型化工容器的外层包上隔热层, 以减少热损失 , 若容器外表温度为500℃ ,而 环境温度为20℃ ,采用某隔热材料, 其厚度为240mm,λ =此时单位面积的热损失为_______。 ( 注 : 大型容器可视为平壁) *** 答案 ***1140w 2.（ 6 分）某大型化工容器的外层包上隔热层, 以减少热损失 , 若容器外表温度为500℃ ,而 环境温度为20℃ ,采用某隔热材料, 其厚度为120mm, λ =此时单位面积的热损失为 _______。 ( 注 : 大型容器可视为平壁) *** 答案 *** 1000w 3.（ 6 分）某大型化工容器的外层包上隔热层, 以减少热损失 , 若容器外表温度为150℃ ,而 环境温度为20℃ , 要求每平方米热损失不大于500w, 采用某隔热材料, 其导热系数λ =则其 厚度不低于 _______。 ( 注 : 大型容器可视为平壁) *** 答案 *** 91mm 4.（ 6 分）某间壁换热器中 , 流体被加热时 , 圆形直管内湍流的传热系数表达式为 ___________________. 当管内水的流速为0.5m.s时,计算得到管壁对水的传热系数α= .K). 若水的其它物性不变, 仅改变水在管内的流速, 当流速为 0.8m.s时,此时传热系数α =_____________. *** 答案 ***α =(λ /d)Re Pr α = .K) 5.（ 6 分）某间壁换热器中 , 流体被加热时 , 圆形管内湍流的传热系数表达式为 _____________________. 当管内水的流速为0.5m.s时,计算得到管壁对水的传热系数α= .K). 若水的其它物性不变, 仅改变水在管内的流速, 当流速为 1.2m.s时,此时传热系数α =________________. *** 答案 ***α =(λ /d)Re Pr
抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a t u +??=??，T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.1）的定解问题分两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： （1.2） ()()x x u ?=0,, ∞<<∞-x 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： ()13.1 ()()x x u ?=0,, l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。
现在考虑边值问题（1.1），（1.3）的差分逼近 取 N l h = 为空间步长，M T = τ为时间步长，其中N ，M 是 自然数， jh x x j ==, ()N j ,,1,0Λ=； τ k y y k ==, ()M k ,,1,0Λ= 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 ()j i y x ,表 示网格节点； h G 表示网格内点（位于开矩形G 中的网格节点）的集合； h G 表示位于闭矩形G 中的网格节点的集合； h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点()k i t x ,处的待求近似解，N j ≤≤0，M k ≤≤0。 注意到在节点()k i t x ,处的微商和差商之间的下列关系 （(,)k j k j u u x t t t ????≡ ? ????）： ()() ()ττ O t u t x u t x u k j k j k j +??? ????=-+,,1 ()() ()2112,,ττ O t u t x u t x u k j k j k j +??? ????=--+ ()()()h O x u h t x u t x u k j k j k j +??? ????=-+,,1 ()() ()h O x u h t x u t x u k j k j k j +??? ????=--,,1 ()() ()2112,,h O x u h t x u t x u k j k j k j +??? ????=--+ ()()() ()2 222 11,,2,h O x u h t x u t x u t x u k j k j k j k j +???? ????=+--+ 可得到以下几种最简差分格式
传热习题及答案 一、选择题： 1、关于传热系数K 下述说法中错误的是（ ）C A 、传热过程中总传热系数K 实际是个平均值； B 、总传热系数K 随着所取的传热面不同而异； C 、总传热系数K 可用来表示传热过程的强弱，与冷、热流体 的物性无关； D 、要提高K 值，应从降低最大热阻着手； 2、在确定换热介质的流程时，通常走管程的有（ ），走壳程 的有（ ）。Ａ、Ｃ、Ｄ；Ｂ、Ｅ、Ｆ Ａ、高压流体； Ｂ、蒸汽； Ｃ、易结垢的流 体； Ｄ、腐蚀性流体； Ｅ、粘度大的流体； Ｆ、被冷却的流 体； 3、影响对流传热系数的因素有( )。A 、B 、C 、D 、E A 、产生对流的原因； B 、流体的流动状况； C 、流体的物性； D 、 流体有无相变；E 、壁面的几何因素； 4、某套管换热器，管间用饱和水蒸气将湍流流动的空气加热 至指定温度，若需进一步提高空气出口温度，拟将加热管管径 增加一倍（管长、流动状态及其他条件均不变），你认为此措 施是：A A 、不可行的； B 、可行的； C 、可能行，也可能不行； D 、视具 体情况而定； 解：原因是：流量不变 2d u =常数 当管径增大时，a. 2/u l d ∝，0.80.2 1.8/1/u d d α∝= b. d 增大时，α增大，d α∝ 综合以上结果， 1.81/A d α∝，管径增加，A α下降 根据()21p mc t t KA -=m Δt 对于该系统K α≈∴ 21 12ln m t t KA t A T t T t α-?≈-- 即 12 1 ln p mc A T t T t α=-- ∵A α↓ 则12ln T t T t -↓-∴2t ↓
第1章前言 1.1问题背景 在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。诸如粒子扩散或神经细胞的动作电位。也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。热方程及其非线性的推广形式也被应用与影响分析。 在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。 科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。 解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。 1.2问题现状 近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。 目前，在欧美各国MATLAB的使用十分普及。在大学的数学、工程和科学系科，MATLAB
《计算机数学基础(2)》辅导六 第14章常微分方程的数值解法 一、重点内容 1.欧拉公式： (k＝0，1，2，…，n－1) 局部截断误差是O(h2)。 2. 改进欧拉公式： 或表示成： 平均形式： 局部截断误差是O(h3)。 3. 四阶龙格――库塔法公式： 其中κ1＝f(x k，y k)；κ2＝f(x k+ 0.5h，y k+ 0.5 hκ1)；κ3＝f(x k+ 0.5 h，y k+ 0.5 hκ2)； κ4＝f(x k+h，y k+hκ3) 局部截断误差是O(h5)。
二、实例 例1用欧拉法解初值问题 取步长h＝0.2。计算过程保留4位小数。 解h＝0.2，f(x，y)＝－y－xy2。首先建立欧拉迭代格式 ＝0.2y k(4－x k y k) (k＝0，1，2) 当k＝0，x1＝0.2时，已知x0＝0，y0＝1，有 y(0.2)≈y1＝0.2×1(4－0×1)＝0.8 当k＝1，x2＝0.4时，已知x1＝0.2，y1＝0.8，有 y(0.4)≈y2＝0.2×0.8×(4－0.2×0.8)＝0.6144 当k＝2，x3＝0.6时，已知x2＝0.4，y2＝0.6144，有 y(0.6)≈y3＝0.2×0.6144×(4－0.4×0.6144)＝0.4613 例2 用欧拉预报－校正公式求解初值问题 取步长h＝0.2，计算y(1.2)，y(1.4)的近似值，小数点后至少保留5位。 解步长h＝0.2，此时f(x，y)＝－y－y2sin x 欧拉预报－校正公式为： 有迭代格式：
当k＝0，x0＝1，y0＝1时，x1＝1.2，有 ＝y0(0.8－0.2y0sin x0)＝1×(0.8－0.2×1sin1)＝0.63171 y(1.2)≈y1 ＝1×(0.9－0.1×1×sin1)－0.1(0.63171＋0.631712sin1.2)＝0.71549 当k＝1，x1＝1.2，y1＝0.71549时，x2＝1.4，有 ＝y1(0.8－0.2y1sin x1)＝0.71549×(0.8－0.2×0.71549sin1.2) ＝0.47697 y(1.4)≈y2 ＝0.71549×(0.9－0.1×0.71549×sin1.2) －0.1(0.47697＋0.476972sin1.4) ＝0.52611 例3写出用四阶龙格――库塔法求解初值问题 的计算公式，取步长h＝0.2计算y(0.4)的近似值。至少保留四位小数。 解此处f(x，y)＝8－3y，四阶龙格――库塔法公式为 其中κ1＝f(x k，y k)；κ2＝f(x k+ 0.5h，y k+ 0.5 hκ1)；κ3＝f(x k+ 0.5 h，y k+ 0.5 hκ2)；
三 计算题 1 （15分）在如图所示的输水系统中，已知 管路总长度（包括所有当量长度，下同）为 100m ，其中压力表之后的管路长度为80m ， 管路摩擦系数为0.03，管路内径为0.05m ， 水的密度为1000Kg/m 3，泵的效率为0.85， 输水量为15m 3/h 。求： （1）整个管路的阻力损失，J/Kg ； （2）泵轴功率，Kw ； （3）压力表的读数，Pa 。 解：（1）整个管路的阻力损失，J/kg ； 由题意知， s m A V u s /12.2) 4 05.03600(15 2 =??==π 则kg J u d l h f /1.1352 12.205.010003.022 2=??=??=∑λ （2）泵轴功率，kw ； 在贮槽液面0-0′与高位槽液面1-1′间列柏努利方程，以贮槽液面为基准水平面，有： ∑-+++=+++10,1 21020022f e h p u gH W p u gH ρ ρ 其中， ∑=kg J h f /1.135, u 0= u 1=0， p 1= p 0=0（表压）, H 0=0， H=20m 代入方程得： kg J h gH W f e /3.3311.1352081.9=+?=+=∑ 又 s kg V W s s /17.410003600 15 =?= =ρ 故 w W W N e s e 5.1381=?=， η=80%， kw w N N e 727.11727===η 2 （15分）如图所示，用泵将水从贮槽送至敞口高位槽，两槽液面均恒定 不变，输送管路尺寸为φ83×3.5mm ，泵的进出口管道上分别安装有真空表和压力表，真空表安装位置离贮槽的水面高度H 1为4.8m ，压力表安装位置离贮槽的水面高度H 2为5m 。当输水量为36m 3/h 时，进水管道全部阻力损失为1.96J/kg ，出水管道全部阻力损失为4.9J/kg ，压力表读数为2.452×
文献综述 信息与计算科学 热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展. 计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”. 在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程. 有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计
传热部分习题答案 1-7 热电偶常用来测量气流温度。如附图所示，用热电偶来测量管道中高温气流的温度T f ,壁管温度f w T T <。试分析热电偶结点的换热方式。 解：具有管道内流体对节点的对流换热，沿偶丝到节点的导热和管道内壁到节点的热辐射 1-21 有一台气体冷却器，气侧表面传热系数1h ＝95W/(m 2 .K)，壁面厚δ＝2.5mm ， )./(5.46K m W =λ水侧表面传热系数58002=h W/(m 2.K)。设传热壁可以看成平壁，试计 算各个环节单位面积的热阻及从气到水的总传热系数。你能否指出，为了强化这一传热过程，应首先从哪一环节着手？ 解： ;010526.0111== h R ;10376.55.460025.052-?===λδR ; 10724.1580011423-?===h R 则λδ++= 21111 h h K ＝94.7 )./(2 K m W ，应强化气体侧表面传热。 1-22 在上题中，如果气侧结了一层厚为2mm 的灰，)./(116.0K m W =λ;水侧结了一层厚为1mm 的水垢)./(15.1K m W =λ。其他条件不变。试问此时的总传热系数为多少？ 解：由题意得 5800115.1001.05.460025.0116.0002.09511 111 2 3322111++++= ++++= h h K λδλδλδ ＝34.6)./(2 K m W 1-32 一玻璃窗，尺寸为60cm cm 30?，厚为4mm 。冬天，室内及室外温度分别为20℃及 -20℃，内表面的自然对流换热表面系数为W ，外表面强制对流换热表面系数为50)./(K m W 。玻璃的导热系数)./(78.0K m W =λ。试确定通过玻璃的热损失。 解： λδA Ah A h T + +?= Φ2111 ＝57.5W -2 一冷藏室的墙由钢皮矿渣棉及石棉板三层叠合构成，各层的厚度依次为0.794mm.,152mm 及9.5mm ，导热系数分别为45)./(K m W ,0. 07)./(K m W 及0.1)./(K m W 。冷藏室的有效换热面积为37.22 m ，室内外气温分别为-2℃及30℃，室内外壁面的表面传热系数可分别按 1.5)./(2K m W 及 2.5 )./(2 K m W 计算。为维持冷藏室温度恒定，试确定冷藏室内的冷却排管每小时需带走的热量。 解：由题意得 332211212 111λδλδλδ++++-? =Φh h t t A ＝2.371.00095.007.0152.045000794.05.215.11) 2(30?+ +++-- ＝357.14W 357.14×3600＝1285.6KJ
化工原理课后习题答案第4章传热习题解答
习 题 1. 如附图所示。某工业炉的炉壁由耐火砖λ1=1.3W/（m·K ）、绝热层λ2=0.18W/（m·K ）及普通砖λ3=0.93W/（m·K ）三层组成。炉膛壁内壁温度1100o C ，普通砖层厚12cm ，其外表面温度为50 o C 。通过炉壁的热损失为1200W/m 2，绝热材料的耐热温度为900 o C 。求耐火砖层的最小厚度及此时绝热层厚度。 设各层间接触良好，接触热阻可以忽略。 已知：λ1=1.3W/m·K ，λ2=0.18W/m·K ， λ3=0.93W/m·K ，T 1＝1100 o C ，T 2＝900 o C ，T 4＝50o C ，3 δ＝12cm ，q ＝ 1200W/m 2，Rc ＝0 求： 1 δ＝？2 δ＝？ 解： ∵δλT q ?= ∴1 δ＝m q T T 22.01200 900 11003.12 1 1 =-? =- λ 又∵3 3 224 23 4 33 2 3 22 λδλδδλδλ+-= -=-=T T T T T T q ∴W K m q T T /579.093 .012 .0120050900233422 2?=--=--= λδλ δ 得：∴m 10.018.0579.0579.022 =?==λδ
习 题1附图 习题2附图 2. 如附图所示。为测量炉壁内壁的温度，在炉外壁及距外壁1/3厚度处设置热电偶，测得t 2=300 o C ，t 3=50 o C 。求内壁温度t 1。设炉壁由单层均质材料组成。 已知：T 2=300o C ，T 3=50o C 求： T 1=？ 解： ∵δ λ δλ3 13 2 3 T T T T q -=-= ∴T 1－T 3＝3（T 2－T 3） T 1＝2(T 2－T 3)+T 3=3×(300-50)+50＝800 o C
第五章练习 1、棉花保温性能好，主要是因为______。 （A）棉纤维素导热系数小； （B）棉花中含有相当数量的油脂； （C）棉花中含有大量空气，而空气的运动又受到极为严重的阻碍； （D）棉花白色，因而黑度小。 2、在房间中利用火炉进行取暖时，其传热方式为； （A）传导和对流 （B）传导和辐射 （C）传导、对流和辐射，但对流和辐射是主要的。 3、一定流量的液体在一φ25×2.5mm的直管内作湍流流动，其对流传热系数α＝1000w/m2.℃。如流量与物性都不变，改用一 φ19×2mm的直管，则其α值将变为_______w/m2.℃。 （A）1059 （B）1496 （C）1585 （D）1678 4、有一列管换热器，由38根φ25×2.5无缝钢管组成，某有机蒸汽在管外冷凝，管内通冷却水，水流速0.5m/s，相应的Re>104及对流传热膜系数为α，当水流量减为一半时，对流传热膜系数改变为α'＝_____。 （A）α'=（1/2）α（B）α'>（1/2）α （C）α'<（1/2）α（D）不一定 5、蒸汽冷凝时的热阻_______。 （A）决定于汽膜厚度（B）决定于液膜厚度 （C）决定于汽膜和液膜厚度（D）主要决定于液膜厚度，但汽膜厚度也有影响 6、沸腾传热的的过热度增大，其传热系数_________。 （A）增大（B）减小 （C）只在某范围变大（D）沸腾传热系数与过热度无关 7、在一列管式换热器中，用冷却水冷凝酒精蒸汽，冷却水应走_______。 （A）管程（B）壳程 （C）易于清洗侧（D）抗压性大的一侧 8、红砖的黑度为0.93，若温度为300℃，则红砖的辐射能力为______W/m2。 （A）6112.3 （B）5684 （C）916.8 （D）1000 9、已知当温度为T时，耐火砖的辐射能力大于铝板的辐射能力，则铝的黑度______耐火砖的黑度。 （A）大于（B）等于 （C）不能确定是否大于（D）小于
一、选择题 1、关于传热系数K，下述说法中错误的是（） A、传热过程中总传热系数K实际是个平均值； B、总传热系数K随着所取的传热面不同而异； C、总传热系数K可用来表示传热过程的强弱，与冷、热流体的物性无关； D、要提高K值，应从降低最大热阻着手； C 2、在确定换热介质的流程时，通常走管程的有（），走壳程的有（）。 A、高压流体； B、蒸汽； C、易结垢的流体； D、腐蚀性流体； E、粘度大的流体； F、被冷却的流体； A、C、D； B、E、F 3、影响对流传热系数的因素有( )。 A、产生对流的原因； B、流体的流动状况； C、流体的物性； D、流体有无相变； E、壁面的几何因素； A、B、C、D、E 4、对下述几组换热介质，通常在列管式换热器中K值从大到小正确的排列顺序应是（）。 A、②>④>③>①； B、③>④>②>①； C、③>②>①>④； D、②>③>④>①； 冷流体热流体 ①水气体 ②水沸腾水蒸气冷凝 ③水水 ④水轻油 D 5、下述各种情况下对流传热系数由大到小的正确顺序应该是（）。 A、③>④>①>②； B、④>③>②>①； C、③>④>②>①； D、③>②>④>①； ①空气流速为30m/S时的a；②水的流速为1.5m/s时的a； ③蒸汽滴状冷凝时的a；④水沸腾时的a； C 6、传热过程中当两侧流体的对流传热系数都较大时，影响传热过程的将是（）。 A、管避热阻； B、污垢热阻； C、管内对流传热热阻； D、管外对流传热热阻； B 7、关于辐射传热，下述几种说法中错误的是（）。 A、除真空和大多数固体外，热射线可完全透过； B、热辐射和光辐射的本质完全相同，不同的仅仅是波长的范围； C、热射线和可见光一样，都服从折射定律； D、物体的温度不变，其发射的辐射能也不变； A 8、冷热水通过间壁换热器换热，热水进口温度为90?C，出口温度为50?C，冷
第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ， 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。） function rechuandaopde %以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;
h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差
化工原理习题第二部分热量传递 一、填空题： 1.某大型化工容器的外层包上隔热层，以减少热损失,若容器外表温度为500℃, 而环境温度为20℃, 采用某隔热材料，其厚度为240mm，λ=0.57w/m.K，此时单位面积的热损失为____ 1140w ___。(注:大型容器可视为平壁) 2.牛顿冷却定律的表达式为____ q=αA△t _____，给热系数（或对流传热系数）α的单位是__ w/m2.K _____。 3.某并流操作的间壁式换热器中,热流体的进出口温度为90℃和50℃,冷流体的进出口温度为30℃和40℃,此时传热平均温度差△t=____27.9K _____。 3. 某并流操作的间壁式换热器中,热流体的进出口温度为90℃和50℃，冷流体的进出口温度为15℃和30℃,此时传热平均温度差△t=____ 41.6K _____。 4.热量传递的方式主要有三种:__ 热传导___、___热对流 ____、热辐射。 5.对流传热中的努塞特准数式是＿＿Nu＝αｌ/λ＿＿＿＿, 它反映了对流传热过程几何尺寸对α的影响。 6.稳定热传导是指传热系统中各点的温度仅随位置变不随时间而改变。 7.两流体的间壁换热过程中，计算式Q＝α.A.△ｔ，A表示为α一侧的换热壁面面积_______。 8.在两流体通过圆筒间壁换热过程中，计算式Q=K.A.△t中，A表示为____________ A 泛指传热面, 与K 相对应________。 9.两流体进行传热,冷流体从10℃升到30℃,热流体从80℃降到60℃，当它们逆流流动时, 平均传热温差△tm=_____ 50℃_______，当并流时,△tm=___ 47.2℃______。 10.冷、热气体在间壁换热器中换热，热气体进口温度T＝400℃，出口温度T 为200℃，冷气体进口温度ｔ＝50℃，两股气体的质量流量相同,物性数据可视为相同,若不计热损失时，冷气体出口温度为＿250_＿℃；若热损失为5％时，冷气体出口温度为＿＿240℃＿。 11.一列管换热器，列管规格为φ38×3, 管长4m,管数127根，则外表面积F=__F1=127×4π×0.038=60.6m2，而以内表面积计的传热面积F____ F2=127×4π×0.032=51.1m2__________。
第一章 1-1 对于附图所示的两种水平夹层，试分析冷、热表面间热量交换的方式有何不同？如果要通过实验来测定夹层中流体的导热系数，应采用哪一种布置？ 解：（ a ）中热量交换的方式主要有热传导和热辐射。 （ b ）热量交换的方式主要有热传导，自然对流和热辐射。 所以如果要通过实验来测定夹层中流体的导热系数，应采用（ a ）布置。 1-7 一炉子的炉墙厚 13cm ，总面积为 20m 2 ，平均导热系数为 1.04w/m · k ，内外壁温分别是 520 ℃及 50 ℃。试计算通过炉墙的热损失。如果所燃用的煤的发热量是 2.09 × 10 4 kJ/kg ，问每天因热损失要用掉多少千克煤？ 解：根据傅利叶公式 每天用煤 1-9 在一次测定空气横向流过单根圆管的对流换热实验中，得到下列数据：管壁平均温度 t w = 69 ℃，空气温度 t f = 20 ℃，管子外径 d= 14mm ，加热段长 80mm ，输入加热段的功率 8.5w ，如果全部热量通过对流换热传给空气，试问此时的对流换热表面传热系数多大？ 解：根据牛顿冷却公式
1-14 宇宙空间可近似的看作 0K 的真空空间。一航天器在太空中飞行，其外表面平均温度为250K ，表面发射率为 0.7 ，试计算航天器单位表面上的换热量？ 解：航天器单位表面上的换热量 1-27 附图所示的空腔由两个平行黑体表面组成，孔腔内抽成真空，且空腔的厚度远小于其高度与宽度。其余已知条件如图。表面 2 是厚δ = 0.1m 的平板的一侧面，其另一侧表面 3 被高温流体加热，平板的平均导热系数λ =17.5w/m ? K ，试问在稳态工况下表面 3 的 t w3 温度为多少？ 解： 表面 1 到表面 2 的辐射换热量 = 表面 2 到表面 3 的导热量 第二章
第8章 常微分方程数值解法 本章主要内容： 1．欧拉法、改进欧拉法. 2．龙格-库塔法。 3．单步法的收敛性与稳定性。 重点、难点 一、微分方程的数值解法 在工程技术或自然科学中，我们会遇到的许多微分方程的问题，而我们只能对其中具有较简单形式的微分方程才能够求出它们的精确解。对于大量的微分方程问题我们需要考虑求它们的满足一定精度要求的近似解的方法，称为微分方程的数值解法。本章我们主要 讨论常微分方程初值问题?????==00 )() ,(y x y y x f dx dy 的数值解法。 数值解法的基本思想是：在常微分方程初值问题解的存在区间[a,b]内，取n+1个节点a=x 0＜x 1＜…＜x N =b （其中差h n = x n –x n-1称为步长，一般取h 为常数，即等步长），在这些节点上把常微分方程的初值问题离散化为差分方程的相应问题，再求出这些点的上的差分方程值作为相应的微分方程的近似值（满足精度要求）。 二、欧拉法与改进欧拉法 欧拉法与改进欧拉法是用数值积分方法对微分方程进行离散化的一种方法。 将常微分方程),(y x f y ='变为() *+=?++1 1))(,()()(n x n x n n dt t y t f x y x y 1．欧拉法（欧拉折线法） 欧拉法是求解常微分方程初值问题的一种最简单的数值解法。 欧拉法的基本思想：用左矩阵公式计算（＊）式右端积分，则得欧拉法的计算公式为：N a b h N n y x hf y y n n n n -= -=+=+)1,...,1,0(),(1 欧拉法局部截断误差 11121 )(2 ++++≤≤''=n n n n n x x y h R ξξ或简记为O （h 2）。
第五章 传热过程基础 1．用平板法测定固体的导热系数，在平板一侧用电热器加热，另一侧用冷却器冷却，同时在板两侧用热电偶测量其表面温度，若所测固体的表面积为0.02 m 2，厚度为0.02 m ，实验测得电流表读数为0.5 A ，伏特表读数为100 V ，两侧表面温度分别为200 ℃和50 ℃，试求该材料的导热系数。 解：传热达稳态后电热器的加热速率应与固体的散热（导热）速率相等，即 L t t S Q 2 1-=λ 式中 W 50W 1005.0=?==IV Q m 02.0C 50C 200m 02.0212=?=?==L t t S ，，， 将上述数据代入，可得 ()() )()C m W 333.0C m W 5020002.002 .05021??=??-??=-= t t S QL λ 2．某平壁燃烧炉由一层400 mm 厚的耐火砖和一层200 mm 厚的绝缘砖砌成，操作稳定后，测得炉的内表面温度为1500 ℃，外表面温度为100 ℃，试求导热的热通量及两砖间的界面温度。设两砖接触良好，已知耐火砖的导热系数为10.80.0006t λ=+，绝缘砖的导热系数为 20.30.0003t λ=+，W /(m C)??。两式中的t 可分别取为各层材料的平均温度。 解：此为两层平壁的热传导问题，稳态导热时，通过各层平壁截面的传热速率相等，即 Q Q Q ==21 （5-32） 或 2 32212 11b t t S b t t S Q -=-=λλ （5-32a ） 式中 115000.80.00060.80.0006 1.250.00032t t t λ+=+=+?=+ 21000.30.00030.30.00030.3150.000152 t t t λ+=+=+?=+ 代入λ1、λ2得 2.0100)00015.0315.0(4.01500)000 3.025.1(-+=-+t t t t 解之得 C 9772?==t t ()()()C m W 543.1C m W 9770003.025.10003.025.11??=???+=+=t λ 则 () 221 11 m W 2017m W 4 .0977 1500543.1=-? =-=b t t S Q λ 3．外径为159 mm 的钢管，其外依次包扎A 、B 两层保温材料，A 层保温材料的厚度为50 mm ，导热系数为0.1 W /(m·℃)，B 层保温材料的厚度为100 mm ，导热系数为1.0 W /(m·℃)，
常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate
化工原理—传热复习题及答案 一、选择题 1、关于传热系数K，下述说法中错误的是（） A、传热过程中总传热系数K实际是个平均值； B、总传热系数K随着所取的传热面不同而异； C、总传热系数K可用来表示传热过程的强弱，与冷、热流体的物性无关； D、要提高K值，应从降低最大热阻着手； C 2、在确定换热介质的流程时，通常走管程的有（），走壳程的有（）。 A、高压流体； B、蒸汽； C、易结垢的流体； D、腐蚀性流体； E、粘度大的流体； F、被冷却的流体； A、C、D； B、E、F 3、影响对流传热系数的因素有()。 A、产生对流的原因； B、流体的流动状况；
C、流体的物性； D、流体有无相变； E、壁面的几何因素； A、B、C、D、E 4、对下述几组换热介质，通常在列管式换热器中K值从大到小正确的排列顺序应是（）。 A、②>④>③>①； B、③>④>②>①； C、③>②>①>④； D、②>③>④>①； 冷流体热流体 ①水气体 ②水沸腾水蒸气冷凝 ③水水 ④水轻油 D 5、下述各种情况下对流传热系数由大到小的正确顺序应该是（）。 A、③>④>①>②； B、④>③>②>①； C、③>④>②>①； D、③>②>④>①； ①空气流速为30m/S时的a；②水的流速为1.5m/s时的a；
③蒸汽滴状冷凝时的a；④水沸腾时的a； C 6、传热过程中当两侧流体的对流传热系数都较大时，影响传热过程的将是（）。 A、管避热阻； B、污垢热阻； C、管内对流传热热阻； D、管外对流传热热阻； B 7、关于辐射传热，下述几种说法中错误的是（）。 A、除真空和大多数固体外，热射线可完全透过； B、热辐射和光辐射的本质完全相同，不同的仅仅是波长的范围； C、热射线和可见光一样，都服从折射定律； D、物体的温度不变，其发射的辐射能也不变； A 8、冷热水通过间壁换热器换热，热水进口温度为90?C，出口温度为50?C，冷水进口温度为15?C，出口温度为53?C，冷热水的流量相同，且假定冷热水的物性为相同，则热损失占传热量的（）。 A、5%； B、6%； C、7%； D、8%；
传热复习题1 （1）保温瓶在设计和使用过程中采取了哪些防止热损失的措施？ 答：首先，保温瓶瓶胆设计成玻璃夹层结构。夹层因空气被抽出接近真空，可防止对流散热损失。其次，瓶胆夹层内两表面均镀有银、铝等低黑度涂层，增加了辐射传热热阻大幅度降低了辐射散热量。举例说，如夹层内壁温度为98ο C ，外壁温度为28ο C ，黑度为0.95的玻璃表面镀上黑度为0.02的银层后,其辐射散热量可由原来的5502m W 降至6.152m W 。第三，在使用保温瓶时，瓶盖选用导热系数很小的软木制作，而且在灌水时还要在瓶颈处留出一段空气的导热系数比水在小得多，从而有效地降低了瓶口的导热热损失。 （2）计算蒸气在水平管外冷凝的凯恩（Kern ）公式为：=a 0.725(t d r g ??????02 3 μρ λ)。试定性说明各种因素对冷凝给热系数a 的影响。 答：应当指出，冷凝给热的热阻是凝液造成的，因此式中各物性常数都是凝液的物性，而非蒸气的物性。 当λ大时，液膜导热性能良好，a 自然就大；ρ大，液膜容易从壁面上往下滑，同样使a 增大，潜热r 大，a 也大。相反，若蒸气温度和壁温之间的温差t ?大，则意味着单位时间内凝液量增多，凝液膜增厚，这反而不利于传热，因此a 会变小；当粘度μ增大时，因流动阻力增大，液膜增厚，a 必然减小；至于水平管的直径d 0大了，会使管子下部液膜加厚，同样不利于传热， a 也要变小。 将上面定性分析结果与凯恩公式对照后可以发现：二者是完全一致的． （3）换热器的热负荷与传热速率有何不同？ 答：冷、热流体在单位时间内在换器中所交换的热量，称为换热器的热负荷．它是针对换热任务的需要提出的，是生产上要求换热器应具有的换热能力，热负荷可根据生产中物流量、进出口温度及状态变化求化求出。 而传热速率则是指换热器本身在一定操作条件下所具有的传递热量的能力，是换热器本身的特性，二者是不相同的。 容易混淆的是，实际生产中设计或标定换热器时，常把传热速率与热负荷视为相等，一般都是通过热负荷的计算，求得换热器应具有的传热速率，再依据传热基本方程求出所需换热器的传热面积，尽管二者在数值上常视为相等，但就其本质讲，含义是完全不同的。 （4）何谓换热器的控制热阻? 答:换热器的总热阻1/K 主要取决于冷、热流体的对流传热热阻，当然也和管壁的热阻及污垢热阻有关，即， λααb K i ∑++=0111 若忽略管壁及污垢热阻，则有 1 11αα+≈i K 如果i α和0α相接近，也就是两种流体的传热阻力差不多时，在谋求强化传热过程中，一般要考虑把i α、0α都增大。 但往往有这种情况，两者的α 值相差很大，例如i α>>0α,则 1 1 αα<< i 。 这时 1 1α≈K K ≈0α 即总传热系数K 值接近对流传热系数小的一侧流体的α 值，在本例条件下总热阻完全被管外的对流传热热阻所控制。1/0α被称为控制热阻。 控制热阻的概念对强化传热研究是非常重要的。它可以使我们在寻求强化的过程中紧紧抓住主要矛盾，提醒我们为了增大换热器的总传热系数，减少热阻，关键在于增大控制热阻一侧的α。 （5）室内暖气片为什么只把外表面制成翅片状？ 答：室内暖气片内部的热流体（水或水蒸气）与片外的冷流体（空气）之间的热递属对流—传导—对流过程。因为