前两章讨论了对固体材料与静止水的加热，加热对象都是单一均匀物质。本章将讨
论含气多孔材料，即加热对象是固气共存的物体，或称固气多孔材料，如图4-1所示。
TaTP
TP Ta
空气
(a) 被上、下热板加热 (b)被上、下热空气流加热 图4-1 固气多孔平板加热
固气多孔材料是多孔材料大家庭中重要一员。例如，用于隔热与防震的泡沫塑料，
泡沫铝等粉末冶金材料, 加气混凝土与泡沫砖，绝干木材等都可视为固气多孔材料。这里的气体大都是单一物质，常指空气。发生在固气多孔材料内的传热传质过程视孔隙结构不同有较大差异。
本章先介绍反映多孔材料内渗流基本规律的达西定律及渗透率，推导描述固气多孔
材料内温度分布，气相速度及密度分布的质量与能量守恆方程。随后，由易及繁，分别讨论封闭型孔隙结构固气多孔材料传热，连通型孔隙结构固气多孔材料内一维与二维等温渗流以及在热板加热或热空气流加热条件下固气多孔材料内传热传质过程。针对每种情形，在认识其物理过程基础上建立数学模型，也即写出控制方程及定解条件，而后对其数值求解，最终详举实例，显示求解的结果。
顺便提及，本章讨论的内容与多孔介质传热传质中的干饱和多孔介质是十分近似
的，只是那里的气体是过热的水蒸气，这里是干空气。有关它们的联系，我们在第1章1.1.1小节中已有说明。
4.1 达西定律与渗透率
4.1.1 达西定律
达西定律是实验定律，最先由达西通过水在泥砂中渗流得到。按示意图4-2，其数学形式为
V=
•
K
μA
Δp
(4.1.1) H
式中，A与H是多孔材料矩形板试样的面积与高(厚)度，也是流体流经的截面积与长度，单位分别为m和m; V是牛顿流体以很低速度渗流通过的体积流率，单位为m/sec.;μ为流体的粘性系数，单位为Pa⋅sec.;Δp为流动压力降，即图4-2中p1−p2，单位为Pa;比例系数K为渗透系数，也称渗透率，单位为m。
1
2
2
•
3
p1 u 砂 •V p2A H 图4-2 达西定律示意图
x 与傅里叶定律式(2.1.10)相仿，将式(4.1.1) 表示成微分的形式，即 v=−或
→
K
μgradp (4.1.2)
−∇p−
μK
→
v=0 (4.1.3)
在直角坐标系内，常表示为
K∂pK∂pK∂p
， v=− ， w=− (4.1.4) μ∂xμ∂yμ∂z
式中，u,v与w分别是流体在多孔材料内沿坐标x,y与z方向的速度，常用单位为mm/sec.。需强调指出，这里的速度u,v与w并不是存在于多孔材料内孔隙间流体的真实
u=−
速度，而是一种折算速度，如图4-2所示，为 u=
V
(4.1.5) A
•
式(4.1.2)中的负号表示速度方向与压力梯度相反，也就是说，流体总是朝着压力下降的方向流动。关于这一折算速度，各类文献[1]中有不同的命名，如渗流(seepage)速度，过滤(filtration)速度，表面(superficial)速度，体积通量及达西速度。本书将使用达西速度。 根据体积孔隙率φ等于面积孔隙率的假定，达西速度与真实速度的关系可表示为： v=φV (4.1.6)
在引入达西速度后，多孔介质被视为连续介质，达西速度可对介质内任何一点定义，只是这个点要理解为足够的大，它是一个表征元容积，即表征元容积本身是体积微小的多孔体。正因为引入达西速度后可将多孔介质视为连续介质，人们在描述诸如多孔介质内速度、密度及温度分布时，直接用速度一词。而其真实含义是达西速度。读者在书写有关多孔材料质量、动量及能量守恆方程时，应特别注意达西速度的引入对它们的影响。这里顺便提及，分清达西速度与真实速度的差别如同在描述湿空气状态时要分清水蒸气容积密度与水蒸气密度的差别一样。
达西定律以其物理含义清晰及数学表示简明，与热传导中傅里叶定律，传质的费克定律，甚至与粘性流体中牛顿定律，弹性体中虎克定律及电学中欧姆定律一样，被后人广泛应用，并正在发挥着重要作用。尤其是，在讨论多孔材料传热传质问题时，常常用达西定律替代动量方程，使整个问题的描述与求解都大大简化。
→
→
2
为了更好理解与应用达西定律，现分析与讨论达西定律与粘性流体运动方程的关系，归结为如下二个方面：
(1)达西定律是粘性流体运动方程在特定条件下的结果
按粘性流体力学教科书[2,3]介绍，不可压缩的粘性流体其矢量形式的纳维-斯托克斯(N-S)方程为
→→→∂V
+(V⋅∇)V]=−∇p+μΔV (4.1.7) ρf[∂τ→
式(4.1.7)的基本物理意义为，在粘性流体流场内，单位体积物质在单位时间内动量的变化量等于作用在该单位体积上的力。这里，式(4.1.7)等号左边表示动量的变化，等号右边表示作用在单位体积上的力。式中的两项是压力项与粘性力项，都是作用在单位体积表面上的力，称作表面力。此处，体积力被忽略。
在如图4-3所示，间距为δ两平行平板间的定常流动，也称泊肃叶(Poiseuille)
流动，在这样特定条件下，式(4.1.7)被简化。在定常流动条件下，动量变化为0，即式(4.1.7)等号左边为0，并且，速度u不随x变化，也即，与坐标x无关。对于平行平板间的流体流动，意味着流动速度只能是坐标z的函数。这样，式(4.1.7)简化为
∂2u∂p
−+μ2=0 (4.1.8)
∂z∂x
在z=0与z=δ处，即流体在上、下平行板处，u=0。
Δp z
A,δ umux
H 图4-3平行槽道内低速流
根据己提供的边界条件求解方程(4.1.8)可得 u=
z(z−δ)∂p
(4.1.9)
2μ∂x
根据平均速度um的定义并将式(4.1.9)代入其中可得：
δ2∂pV1δ (4.1.10) um==∫udz=−
12μ∂xAδ0
将
•
∂pΔp
表示为，并代入式(4.1.10)，再整理成如式(4.1.1)所示达西定律的形式，则
H∂x
3
δ2Δp
A V= (4.1.11) 12μH
•
对比式(4.1.11)与式(4.1.1)，可得渗透率K为δ/12。 综上所述，式(4.1.11)是特定条件下粘性流体运动方程(4.1.8)的解，也即间接地说明达西定律是粘性流体运动方程在特定条件下的结果。这一结果反映多孔介质内粘性流体的压力差与粘性力的平衡，忽略流体动量变化。 (2)扩展的达西定律
当多孔介质内流体的动量变化不能被忽略时，用式(4.1.2)来计算流动速度就不行了。很多研究者提出很多模型对达西定律进行修正，这里仅扼要介绍Wooding[4]提出的方程。他从不可压缩粘性流体N-S方程(4.1.7)出发，对达西定律进行两项重大修正：
其一，将式(4.1.7)等号左边关于动量变化中的速度V通过式(4.1.6)用达西速度v代替，从而引入了孔隙率φ。
其二，将式(4.1.7)等号右边关于压力梯度项及粘性力项用式(4.1.3)的达西定律替代。
由此，使式(4.1.3)成为
→
1∂v1→μ→
ρf[+2(v⋅∇)v]=−∇p−v (4.1.12)
φ∂τφK
→
2
→→
4.1.2 渗透率
出现在式(4.1.1)中的渗透率K是反映多孔材料性质的重要参数。为了更好地理解
它的含义，将式(4.1.1) 改写成：
K=
Vμ (4.1.13)
AΔp/H
•
渗透率K基本计量单位为达西(darcy)，1Darcy是指粘性系数为1厘泊(cP)，即μ=1cP，的流体在1厘米具有1物理大气压的压降，即ΔP/H=1atm/cm，条件下，以
1秒渗流1立方厘米，即V=1cm/sec.时的渗透率。将上述单位代入式(4.1.13)，可得
•
3
1[cm3/sec]1[cP]
1[Darcy]= 2
1[cm]1[atm/cm]
再将粘性系数1P=0.1Pa⋅sec.，物理大气压1 atm=101325 Pa代入上式，可得：
1[cm3/sec]0.01⋅0.1[Pa⋅sec]
1[Darcy]=
1(cm2)101325(Pa/cm)
=
1[cm]0.01⋅0.1
=0.987⋅108cm2=0.987μm2=0.987⋅10−12m2
101325
2
2
在实际计算中所采用的渗透率单位常有Darcy、μm、cm与m多种，它们间相
差几个数量级，读者在使用时务必谨慎以免出错。
与诸多物理量同时具有相应的无量纲一样，渗透率也有相应的无量纲渗透率，称为达西数，记作Da，定义为
22
Da=
K 2L
4
这里，L 是定性尺寸。请注意，选什么几何尺度作为定性尺寸，会使Da的值有很大的差
别。例如，用多孔材料的整体(如平板)几何尺度，可能是cm或m的数量级。若用多孔材料孔隙的几何尺度，则可能是μm或mm的数量级。本书不采用无量纲渗透率Da来讨论问题。 大多数多孔材料的渗透率是由实验测量提供。表4-1(摘自[4])列出若干典型多孔材料的孔隙率与渗透率，表4-2(摘自[5])给出的是若干典型木材纵向(沿木材纹理方向)的渗透率，相应的横向渗透率，其值大约小3个数量级，在10-15~10-18m2这个范围内。
提醒读者注意，在不少文献[5,6]中，把K/μ定义为渗透率而把K定义为比渗透率。对于单相流体而言，为避免可能引起混淆，本书与很多著作[1,7]相一致，直接将K定义为渗透率。
表4-1常见多孔固体材料的孔隙率与渗透率
名称 2
孔隙率φ 渗透率K[m]
砖 烟 水泥(沥青) 铜粉 软木板 头发 皮革 石灰石 砂 砂岩(油砂) 硅土粉末 土壤
表4-2几种木材纵向(顺纹理方问)的渗透率
名称
渗透率
Dancy
m2
0.12~0.34
0.09~0.34
0.56~0.59 0.04~0.1 0.37~0.5 0.08~0.38 0.37~0.49 0.43~0.
10−15~2.2⋅10−13 1.1⋅10−9
10−13~2.3⋅10−11
3⋅10−10~1.5⋅10−9 2.4⋅10−11~5.1⋅10−11 8.3⋅10−10~1.2⋅10−9 9.5⋅10−14~1.2⋅10−13 2.0⋅10−15~4.5⋅10−14 2.0⋅10−11~1.8⋅10−10 5.0⋅10−16~3.0⋅10−12
1.3⋅10−14~5.1⋅10−14 2.9⋅10−13~1.4⋅10−11
红栎 1000~100 10-9~10-10 椴木 100~10 10-10~10-11 槭木、松木边材、花旗松边材(太平洋沿岸) 10~1 10-11~10-12
云松边材、松木边材 1~0.1 10-12~10-13 花旗松心材(太平洋海岸) 0.1~0.01 10-13~10-14 白栎心材、山毛榉心材、松木心材、花旗松心材(山地)
0.01~0.001 10-14~10-15
既然多孔材料的渗透率作为材料的一种特性，它由材料几何结构所决定。这样，对于结构形式比较简单的材料，应可以根据其几何特征用理论方法计算得到渗透率。如图4-3所示平行平板的结构，按式(4.1.11)与式(4.1.1)的对比中求得渗透率K为δ/12。若
2
δ=0.1mm，则K≈10−9m2。
确实，已有很多研究者对各类多孔材料建立数学模型以预测材料的渗透率[8]。其中比较成熟的应推Carman-Kozeny模型，它适用于由颗粒或纤维构成的固定床。模型中引入颗粒的有效平均直径dp，如果认为构成多孔材料的固体颗粒直径是均匀的，孔隙率是已知的
φ，则 Carman-Kozeny模型给出的理式为：
5
K=
23dpφ180(1−φ)2
(4.1.14)
式中，常数180是为使理论计算的数据与实验测量数据相一致而得到的。
对于结构复杂的众多多孔材料，Carman-Kozeny模型难以普遍适用。但该模型也给人们某种启示，在合理确定孔隙率、有效直径及相应常数有望使修正的Carman-Kozeny模型对材料的渗透率进行成功的预测。文献[9]在这方面做了有意义的尝试，用修正的Carman-Kozeny模型预测了木材人造板(定向刨花板)的渗透率。
从本节上述介绍的达西定律与渗透率中我们看到，像木材这类渗透率低的多孔材料，在其被加工等工作过程中材料内流体流动速度很低，用达西定律作为动量方程是可行的。虽然它是一种近似，但计算过程可大大简化。因此，在本书以后各章分析材料内流体的运动大都采用达西定律。
4.2 守恒方程及其封闭性
这一节推导固气多孔材料传热传质过程的质量守恆与能量守恆方程。 4.2.1质量守恒方程 (1)固相质量守恆
若孔材料内固相骨架是密实的，且不计固相骨架与气相之间化学变化等因素可能引
起质量交换，则材料被加热时，固相骨架的密度ρs不变(忽略体积膨胀)，表示为：
∂ρs
=0 (4.2.1) ∂τ 多孔材料内的气相物质，若处在封闭的孔隙内，其质量与密度也不发生变化。若孔隙连通，材料被加热时，相邻孔隙间气体就会流动，气体的密度ρa就会变化，其质量守恆的形式就不像式(4.2.1)那样简明。下面我们采用第2章中在特定坐标系内推导导热方程的方法来推导气体质量守恆方程。
(1)气相质量守恆 (ρaw)z+dzdxdy(ρau)xdydz
(ρav)ydxdz (ρav)y+dydxdz z
(ρau)x+dxdydz
y (ρaw)zdxdy
x
图4-4 推导固气多孔材料内气体质量守恒方程的单元体 在直角坐标系内取元体dv(=dxdydz)作为分析对象，如图4-4所示。对于该单元而
言，单位时间内通过边界面进入体系的质量应等于该体系质量的增加。其数学形式为
6
∂φρa
dxdydz=[(ρau)x−(ρau)x+dx]dydz∂τ+[(ρav)y−(ρau)y+dy]dxdz+[(ρaw)z−(ρaw)z+dz]dxdy
其中
(4.2.2)
∂ρau
dx (4.2.2a) ∂x∂ρv
(ρav)y+dy=(ρav)y+ady (4.2.2b)
∂y∂ρw
(ρaw)z+dz=(ρaw)z+adz (4.2.2c)
∂z
(ρau)x+dx=(ρau)x+
将式(4.2.2a)、(4.2.2b)和(4.2.2c)代入式(4.2.2)，经整理后，可得
∂φρa∂ρu∂ρv∂ρw
=−（a+a+a） (4.2.3) ∂τ∂x∂y∂z
式(4.2.3)即为固气多孔材料气相质量守恒方程。上式的物理含义是，单位体积单位时间内
气体质量的增加等于单位时间内流体通过体积边界流入体积的质量。特别提醒读者注意，式中的速度u、v与w为达西速度，密度ρa则是空气真实密度，体孔隙率与面孔隙率相等。 固气多孔材料的密度ρ可定义为
ρ=(1−φ)ρs+φρa (4.2.4)
现将式(4.2.1)改写成：
∂(1−φ)ρs
=0 (4.2.5)
∂τ再将式(4.2.5)与式(4.2.3)相加，并代入式(4.2.4)，可得：
∂ρ∂ρu∂ρv∂ρw=−（a+a+a） (4.2.6) ∂τ∂x∂y∂z
这是固气多孔材料质量守恆的另一种表示形式。 4.2.2能量守恒方程
与推导导热方程式(2.2.20)与上述气相质量守恆方程式(4.2.6)相仿，推导如下固相、
气相及固气多孔材料的能量守恆方程。
(1) 固相能量守恆方程
与4.2.1小节一样，针对以图4-4所示的元体dv(=dxdydz)为一热力学体系，可写出元体内固相骨架的能量平衡关系如下：
∂Ts∂2Ts∂2Ts∂2Ts
(1−φ)ρscs=(1−φ)ks(2+2+2)+hsgSsg(Ta−Ts) (4.2.7)
∂τ∂x∂y∂z
式(4.2.7)的物理含义可表示为
单位体积内固相骨架单位时间内以导热方 单位时间单位体积内在单位时间内能量的= 式通过体积边界导入 + 气相向固相的传热的 增加量 单位体积的能量 热量
式(4.2.7)中(1−φ)ρs是单位体积固气多孔材料中固相骨架的质量，此处的φ为体积孔隙
7
率。
式中(1−φ)ks中的φ为面孔隙率，认为两者的φ相同。式中hsg与Ssg分别为固气多孔材料内固气相间的热交换系数与单位体积内固气相之间接触表面积。这个表面积被称为比面，与孔隙率φ被视为多孔材料特性参数一样，比面Ssg也是多孔材料的重要特征量。在化学工业中作为催化剂的材料就追求高的比面。
(2) 气相能量守恆方程
气相的能量平衡关系为
∂ρacaTa∂2Ta∂2Ta∂2Ta
=φka(2+2+2)φ∂τ∂x∂y∂z
−c（a
∂ρauTa∂ρavTa∂ρawTa
++）+hsgSsg(Ts−Ta)∂x∂y∂z
(4.2.8)
式(4.2.8)的物理含义可解释为
以导热的方式在单位 固气多孔材料单位空 间内气体的能量在单= 时间内从单位空间边+ 界面传入的热量 位时间内的增加量 (3) 固气多孔材料能量守恆方程
以对流的方式在单位时间内从单位空间边界面传入的热量 单位时间内单位空间+ 内部固相传给气相的热量 为推导固气多孔材料能量守恆方程，再引入如下三项假定条件：
其一，在固气多孔材料内各处，固相骨架的温度Ts与气相的温度Ta始终相同，即
Ts=Ta (4.2.9)
应该说，上述假定是有条件的。条件是气相的流动速度足够的低，固气多孔材料的比面足够的大。式(4.2.9)的成立也被称为达到局部热平衡。 其二，固气多孔材料的导热系数k为：
k=(1−φ)ks+φka (4.2.10)
同样，式(4.2.10)这一假定也是有条件的。条件是固相骨架中的导热与气相中导热相互平行。更详尽的说明在4.3.1小节。
其三，固气多孔材料的比热容ρc为
ρc=(1−φ)ρscs+φρaca (4.2.11)
上式是按比热容的基本定义得到的，并不包含特定的条件。
与推导式(4.2.6)相仿，根据式(4.2.9)、(4.2.10)与(4.2.11)的假定，将式(4.2.7)
与式(4.2.8)相加，可得固气多孔材料的能量平衡关系为
∂ρcT∂2T∂2T∂2T∂ρcuT∂ρcvT∂ρcwT
=k(2+2+2)−（aa+aa+aa） (4.2.12) ∂τ∂x∂y∂z∂x∂y∂z
式(4.2.12)的物理含义为：单位体积的固气多孔材料在单位时间内能量的增加(等号左边)
等于单位时间内以导热方式通过体积边界导入体积的能量与以对流方式通过体积边界传入体积的能量之和(等号右边)。提醒读者注意，上式中，孔隙率不再出现，反映固相与气相间的热交换也不出现，这是由于上述式(4.2.9) 、(4.2.10)与(4.2.11)引入所致。也就是说，式(4.2.12)的成立是以式(4.2.9) 、(4.2.10)与(4.2.11)成立为前提。 4.2.3 动量方程与气体状态方程
8
在固气多孔材料有效孔隙率与渗透率很低，气体流动速度很小条件下，材料内气
体流动的动量方程可取达西定律的形式，式(4.1.4)即为固气多孔材料内气体的动量方程。也即，式(4.1.4)是关于固气多孔材料内寻求气体速度的控制方程。
由式(4.1.4)可知，气体流动的速度受压力梯度支配，而气体的压力取决于当地的温度与气体密度。若将气体视为理想气体，则气体压力pa可表示为：
Ma
式中，R是通用气体常数(8.314J/(mol⋅K))，Ma是空气的摩尔质量(28.97g/mol)，温度T以绝对温度(K)表示。
4.2.4 方程组的封闭性
回顾式(4.2.3)、(4.1.4)、(4.2.13)、(4.2.12)及(4.2.4)共计7个方程中有7个因变量。它们分别是由式(4.2.3)控制的气体密度ρa、由式(4.1.4)控制的气体速度u，v及w、由式(4.2.13)控制的气体压力pa、由式(4.2.12)控制的固气多孔材料的温度T以及由式(4.2.4)所示的多孔材料密度ρ与固相密度ρs及气相密度ρa的联系。这7个方程所包含的7个因变量的内在联系表示在如下图4-5：
ρ paρa
Tu,v,w
图4-5 7个因变量的内在联系框图
下面结合此图对固气多孔材料被加热后发生在材料内部传热传质现象，也即7个因
变量的相互影响作进一步解说。当材料被加热后，温度T就会上升，在材料内部会形成温度场，相应地，气体压力pa会变化，形成气体压力场，材料内各处气体压力不同。气体压力的差别就会推动气体的流动，从而在材料内部形成气体速度场u，v及w。气体的流动，一方面，引起了热对流，导致温度的变化，另一方面，会导致气体密度ρa的变化。气体密度的变化，一方面，导致材料密度ρ的变化，另一方面，则又引起气体压力的变化。
图4-5所示7个因变量的内在联系表明这7个物理量是相互影响的，或者说，往往
互为因果。从数学上讲，包含7个因变量的7个方程是相互耦合的，其中任何一个方程都不可能单独求解。
这种方程式个数与因变量个数相等，且能保证有唯一解的方程组为一封闭方程组。作者在此郑重提及，在构造任何用方程组形式表示的数学模型时，方程组的封闭性是正确数学模型的必要条件。
在求数学模型唯一解时，除方程组的封闭性是必要条件之外，还必须有保证唯一解的定解条件。关于固气多孔材料加热问题的定解条件将在以下几节按不同情形再作介绍。
pa=
ρaRT
(4.2.13)
4.3闭合孔隙多孔材料加热
孔隙闭合的固气多孔材料是一种常见的多孔材料，它们是固气多孔材料的一种极端
9
情形，例如，用于隔热泡沫塑料。
由于材料内的孔隙是闭合的，孔隙内的气体就不会在孔隙间流动，则上述式(4.1.4)
不再存在，u，v及w从因变量中消失。由于气体不流动，孔隙内气体的密度ρa不再变化，相应地，材料密度ρ也不再变化。这样，上述图4-5所示的7个因变量中只保留温度T与气体压力pa二个因变量。也就是，7个方程只保留式(4.2.12)与式(4.2.13)两个方程，且两个方程不再耦合，也即，气体压力pa唯一取决于温度T，而温度并不依赖于气体压力。这样，式(4.2.12)可不依赖于式(4.2.13)而单独求解，且式(4.2.12)退化为导热方程，其形式与密实固体导热方程相同，求解方法也相同。然而，两种导热方程的内容并不相同，为此引入等效导热系数。下面将就闭合孔隙固气多孔材料导热方程及等效导热系数展开讨论。
4.3.1 等效导热系数
闭合孔隙多孔材料导热方程是上述式(4.2.12)的简化，省略其中的对流项，表示为：
∂T∂2T∂2T∂2T
ρc=k(2+2+2) (4.3.1) ∂τ∂x∂y∂z
式中的比热容ρc需满足式(4.2.11)，导热系数k需满足式(4.2.10)。
关于比热容ρc式(4.2.11)的等效性与固气多孔材料结构形式无关，只取决于固相、气相各自的比热容及孔隙率。关于导热系数k式(4.2.10)的等效性不仅取决于固相、气相各自的导热系数及孔隙率，还取决于固气多孔材料结构形式。为此，把导热系数转化为等效导热系数的同时，隐含着以何种理想化结构形式来等效实际的多孔材料。人们为了等效不同形式的实际材料而提出了不同的结构模型[10-14]，对木材结构也不乏其人[15-18]。下面讨论固气多孔材料结构形式与等效导热系数的关系。 (1)简单串联、并联模型
若将固气多孔材料的结构视为简单的串联或并联结构，则与计算电阻网络的串并联相同，计算热阻网络。
QT QT 固相 气相 1−φ QLQL QL固相 QL 气相 φ ‘ QT QT
(a)真实材料 (b)等效热阻材料
图4-6简单串(并)联结构多孔材料 图4-6所示为孔隙率为φ的一单位立方体(例如，1mm)多孔材料，固相与气相的导热系数分别为ks与ka。从分析材料传热而言，图4-6(a)的结构可等效成图4-6(b)的结构。当热流QL纵向通过该材料时，该材料视为并联结构，其热阻网络如图4-7(a)所示的并联网络。图中Ra与Rs分别为气相与固相的热阻。热阻的计算公式与电阻相似，一段导线的电阻，它与导线长度成正比，与导线截面积成反正，比例系数为导电材料的电阻率。同样，
10
3
热阻的表达式为，热阻与热流流经的材料长度成正比，与热流流经的面积Sa成反比，比例等于孔隙率φ，系数为热阻率(导热系数的倒数)。而面积Sa在总截面积为1(1mm)条件下，则
2
Ra=
同理，固相热阻Rs为
1111
= (4.3.2) kaSakaφRs=
1111
(4.3..3) =
ksSsks1−φ
Rs
Rs Ra (a) 并联网络 (b)串联网络
图4-7 简单串(并)联结构热阻网络 根据并联热阻的运算规则，可得
Ra 11RaRs1
== (4.3.4) kL12Ra+Rsks(1−φ)+kaφ式(4.3.4)中第1个等号是图4-6所示材料纵向热阻定义式，其中kL为材料纵向等效导热系数，第2个等号是并联热阻表达式，第3个等号则是将式(4.3.2)的Ra与将式(4.3.3)的Rs代入式(4.3.4)所示的结果。由式(4.3.4)可得到材料纵向等效导热系数kL为
kL=(1−φ)ks+φka (4.3.5)
而当热流QT横向通过该材料时，该材料视为串联结构，其热阻网络如图4-7b所示
的串联网络。其中气相热阻Ra为
1φ⋅1Ra= (4.3.6)
ka12
固相热阻Rs
1(1−φ)⋅1
(4.3.7) Rs=2
ks1
材料横向热阻RT为
φ1−φ11
=+=+ (4.3.8) RT=RRas2
kT1kaks
由式(4.3.8)可得到材料横向等效导热系数kT为
kska
kT= (4.3.9)
ks(1−φ)+kaφ RL=
串联与并联模型的区别在于，并联模型中的固相是连通的，而串联模型中，固相被气相隔开，固相的贡献仅仅在于减少热流通过气相的厚度。虽然，简单串、并联模型过于简化，与实际多孔材料还有较大差距，但运用这种简单模型仍能有效地帮助我们分析实际材料的导热性能，例如，象木材这样的材料。
众所周知，木材在导热性能上有很强取向性，顺木材纹理方向的导热系数kL大于垂直纹理方向的kT。有诸多文献报导了实验结果，kL/kT≈2.5。这里，我们不妨用串、并联模型对它进行分析。视木材顺纹理方向为并联结构，而垂直纹理方向为串联结构，则由式(4.3.5)与式(4.3.9)之比为：
11
kL[(1−φ)ks+φka]2
= (4.3.10) kTkska
据木材学介绍，木材细胞壁固相物质的导热系数ks=0.3~0.7W/(m⋅K)，细胞腔内空气导热系数ks=0.026W/(m⋅K)，两者相差一个数量级。作为一种估算，不妨假定ka=0.1ks，φ=0.5。将它们代入式(4.3.10)，可得kL/kT=2.5。
(2)木材细胞胞管的等效导热系数
这一节，继续用串、并联模型来讨论木材纤维细胞管的等效导热系数。为便于讨论，突出物理概念，把木材纤维细胞壁视为正方形截面细长管。图4-8所示为木材纤维细胞胞管的一个截段，体积为单位立方体，细胞壁厚为δ，孔隙率为φ。
细胞壁 QTz 细胞腔
y
x δQ L
图4-8木材细胞胞管传热示意图 也即
Sa=φ⋅1 (4.3.11) 则纵向热流QL通过细胞壁(固相)的截面积Ss即为孔隙率1−φ，也即 Ss=(1−φ)⋅1 (4.3.12)
图4-8中细胞壁厚δ与孔隙率φ之间有着一一对应关系。因为截面积Ss还可表示为4个梯形面积之和，即
22
QLQT
由图4-8所示结构可知，纵向热流QL通过细胞腔(气相)的截面积Sa即为孔隙率φ，
Ss=4×
1+(1−2δ)
δ=(1−φ)⋅12 (4.3.13) 2
式(4.3.13)经整理后为
4δ2−4δ+1−φ=0 (4.3.14)
式(4.3.14)为关于δ的二次代数方程，其解为
δ=
1−φ (4.3.15) 2
1
⋅1，表示整个截面积被固相2
式(4.3.15)就是δ与φ的对应关系。当φ=0，细胞壁厚δ=占有。当φ=1，δ=0，表示整个截面积被气相占有。
当纵向热流QL通过图4-8结构，则此立方体被视为一并联结构，其热阻网络如图4-9(a)所示。
12
Rs Ra Rs1Rs2RaRs1Rs2(a) 并联网络 (b)串并联网络 图4-9 细胞胞管结构热阻网络
对照图4-9(a)与图4-8(a)可知，两者完全相同。相应地，关于Ra与Rs的表达式(4.3.11)与式(4.3.12)，它们与式(4.3.2)和式(4.3.3)完全相同，这样，纵向等效导热系数也完全相同，由式(4.3.5)表示，这里不再重复。
当横向热流QT通过图4-8结构，则此立方体被视为一串并联结构，其热阻网络如图4-9(b)所示。横向热流QT管胞的顶部壁、底部壁及左、右侧壁与气体夹层分三个分支通过表示结构顶(底)部壁沿y方向的热阻。图4-8结构。其中Rs1是第一分支与第三分支的热阻，其长度为1，截面积为δ⋅1，即
Rs1=
11
(4.3.16)
ksδ⋅1
图4-9(b)中的Rs2是第二分支中管胞壁部分的热阻，表示结构左(右)侧壁厚(沿y
方向)的热阻。热流通过的长度为δ，截面积为(1−2δ)⋅1，即
Rs2=
1δ (4.3.17)
ks(1−2δ)⋅1
图4-9(b)中的Ra是第二分支中管胞腔内气体部分的热阻，表示热流横向通过细胞腔气体部分呈现的热阻。热流通过的长度为(1−2δ)，截面积为(1−2δ)⋅1，即
Ra=
11−2δ1
= (4.3.18)
ka(1−2δ)⋅1ka⋅1
横向热流QT第二分支的热阻，即两层细胞管壁与气体串联的热阻Rsa为
Rsa=2Rs2+Ra=
2δ12kδ+ks(1−2δ)
+=a (4.3.19)
ks(1−2δ)⋅1ka⋅12kaks(1−2δ)⋅1
最终，横向热流QT通过图4-8所示结构，该结构呈现的热阻是三支热阻的并联，表示为：
2kaδ+ks(1−2δ)1
⋅
11RsaRs1/22kaks(1−2δ)⋅12ksδ⋅1
==RT=
kT12Rsa+Rs1/22kaδ+ks(1−2δ)+1 2kaks(1−2δ)⋅12ksδ⋅1 (4.3.20)
=
2kaδ+ks(1−2δ)
4kaksδ2+2ks2(1−2δ)δ⋅1+2kaks(1−2δ)⋅12
从式(4.3.20)可得图4-8所示结构横向等效导热系数kT为
4kaksδ2+(2ks2δ+2kaks⋅1)(1−2δ)
(4.3.21) kT=
[2kaδ+ks(1−2δ)]⋅1
13
由于细胞腔内气体的导热系数很低，作为一种估算，不妨假定ka≈0。这样，式(4.3.5) 与式(4.3.21)分别简化为
kL=(1−φ)ks (4.3.22)
与
kT=
将式(4.3.15)代入上式，可得：
2δks 1
kT=(1−φ)ks (4.3.23)
对比式(4.3.22)与式(4.3.23)可得：
kL1−φ=1+φ (4.3.24) =
kT1−φ对木材而言，孔隙率范围大致为φ=0.4~0.8，则kL/kT=1.6~1.9。 (3)其它模型
这里介绍最早由马克斯威尔(Maxwell,J.C.)于1904针对含有球型孔隙连续固体的导电性能研究[19]，而后由欧克金(Eucken,A)于1940年应用到传热研究的被称为Maxwell-Eucken等效导热系数ke模型，它的数学形式为
3φ] (4.3.25)
ka/ks+2)−φ(
ka/ks−1
不妨将φ=0.1~0.5，与ka/ks=0.1代入上式与式(4.3.5)，即可得到，Maxwell-Eucken模型的ke/ks=0.88~0.47，而并联模型的kL/ks=0.91~0.55。虽说两者有明显差别，
ke=ks[1+
但在一定条件下，作为一种估算，仍可认为有相当程度的接近。 至今，在传热，材料，化工等诸多领域学术刊物上，甚至从网上，可直接查阅到大量有关多孔材料等效导热系数数学模型的文献。它们往往针对材料固相骨架结构形式不同，孔隙形式及其中填充物不同，传热条件不同，运用不同基础理论，提出不同的预测等效导热系数的数学模型。
尽管各种数学模型的形式不同，但共同点是，在引入等效导热系数后把多孔材料简化为均质材料，把存在于多孔材料内复杂的传热现象简化为单纯的热传导现象。
4.3.2 一维闭合孔隙固气多孔材料平板加热
这一节仅以举例形式介绍用等效导热系数方法对固气多孔材料平板进行传热计算。 (1) 问题与已知条件
现有一固气多孔材料平板，例如，绝干木材，初始温度T0为20C，接受如图4-1(b)所示的空气对流加热，环境空气温度Ta为200C，速度为10m/sec.空气流对板的对流换热系数h约为12W/(m⋅K)(参见表3-5)。板厚H为12mm，板长与板宽远大于板厚，视平板为一维(厚度方向)传热。绝干木材细胞壁密度，也即固相密度ρs取1530kg/m，导热系数ks为0.5W/(m⋅K)，孔隙率φ为0.6，绝干木材比热容cs为1.3kJ/(kg⋅K)。孔隙内
14
3
2
0
0
充满干空气，密度ρa为1.16kg/m3，导热系数ka为0.03W/(m⋅K)，比热容为
1.0kJ/(kg⋅K)。 在上述已知条件下，求平板温度变化
(2) 求解过程 运用等效导热系数法就是将不均匀的多孔材料视为均匀的材料，相应地得到等效后均匀材料的热物性参数如下： 密度ρ： ρ=(1−φ)ρs+φρa≈(1−φ)ρs=612kg/m 热容ρc： ρc=(1−φ)ρscs+φρaca≈(1−φ)ρscs=826kJ/(m⋅K)
导热系数：木材是各相异性材料，若假定平板上、下表面平行于木材纹理方向，也即木材纤维细胞胞管与平板上、下表面平行，则板内一维热流垂直于细胞胞管。为此，取式(4.3.23)作为等效导热系数，也即 k=(1−33φ)ks=0.225ks=0.112W/(m⋅K) 在确定热容ρc及等效导热系数k后，导热方程(4.3.1)在图4-1(b)第三类边界条件下的解即为式(2.3.20-3)。现将式(2.3.2-3)代入式(2.3.20-3)，经整理后可得： 4sin(μm)ρcH2
cos(μ2z/H)e(4.3.26) ∑m2μsin(2μ)+m=1mm式(4.3.26)中，时间τ与地点z是自变量，温度T(z,τ)是因变量，其它各量均为参变量。这些参变量中，板厚H是几何参数，T0是初始参数，Ta是环境参数，等效导热系数k，热
T(z,τ)=Ta−(Ta−T0)∞−4
2μmkτ容ρc是材料物性参数，在本算例中作为已知条件被给出且为常数。μm则是超越方程(2.3.12-3)的根。 μtgμ=Bi (2.3.12-3) 式中Bi数已由式(2.3.21)给出定义，本算例中为：
hH12W/(m2⋅K)12⋅10−3m=Bi==0.3 2k2⋅0.112W/(m⋅K)关于超越方程(2.3.12-3)的求解在历史上采用作图的办法，如图4-10所示，在确定的Bi数
条件下，求得一组μm。针对不同的Bi值将求得的μm制作成数据表格供使用。 1008060402000-20-40-60-80-100123y y2=ctgμ y1=μBiμ1 μ2456μ378μ
图4-10 求超越方程(4.3.27)的根
15
在计算技术发达的今天，求解式(2.3.12-3)已变得很容易了，现介绍如下： 将式(2.3.12-3)改写成： 令：y1=
μBi
=ctgμ (4.3.27)
。这里y1是一条直线而y2是以π为周期的余切曲线，
μBi
，y2=ctgμ在一个周期内，y1与y2只有一个交点。μ从0出发时，y1 和图4-11b。将它们与图2-11和图2-12相比，可见不同边界条件下温度响应的一致性与差异。 2001801601401201008060402000100200300400500time(sec.)6007008009001000z/H=0.0z/H=0.1z/H=0.2z/H=0.3z/H=0.4z/H=0.5temperature(C) 图4-11a 第三类边界条件下平板内若干地点温度随时间变化 20018016014012010080604020-0.51 4 17 68 136 272 4 1088 T(c)-0.4-0.3-0.2-0.10z/H0.10.20.30.40.5 图4-11b 第三类边界条件下平板内温度分布 4.4 气体等温渗流 4.3节所讨论的问题是固气多孔材料在加热时，把存在于孔隙内的气体视为静止，把材料内传热过程简化为单纯热传导。而本节则着重讨论固气多孔材料内气体的流动，不计传热过程。也就是，材料与所处环境温度始终相同，只是由于固气多孔材料内气体压力与环境的压力不同所引起的气体流动及相应的气体压力的变化。先讨论一维稳态等温渗流，再讨伦一维与二维非稳态等温渗流。 4.4.1一维稳态等温渗流 针对一维稳态等温渗流这一特定情形，4.2节图4-5所示的关于固气多孔材料传热传质 16 问题的7个因变量中只保留关于空气速度u与空气压力pa两个因变量。它们相应的控制方程变成：(1)由于是一维稳态过程，空气速度u成为常数;(2)在空气速度u为常数条件下，由式(4.1.4)求空气压力pa的分布。 一块固气多孔材料平板置于两侧气体压力分别为p1与p2之间，求板内气体压力pa分布及气流速度u的大小，是典型一维稳态等温渗流问题，如图4-2所示。 为求板内压力分布，将式(4.1.4)中关于u的方程改写为 dpa=− 对式(4.4.1)等号两边进行积分， uμa dx (4.4.1) KuμaK ∫ 式(4.4.2)经运算后， pa p1 dpa=− ∫ x 0 dx (4.4.2) uμa x (4.4.3) K 式(4.4.3)显示了气体压力pa在板内的变化。它从板的左侧(x=0)的压力p1开始，随x线性地下降，直至右侧(x=H)的压力p2。 让式(4.4.3)中的x=H，pa=p2，经整理后可得 Kp1−p2 (4.4.4) u= μaH pa=p1− 若己知板厚H=1cm，p1与p2之差为一个大气压，即，p1−p2=10Pa，多孔材料的渗透率为K=10 5 m2，空气粘性系数为μa=18.6Pa⋅sec.。将这些已知条件代入式 (4.4.4)即可求得气流速度u=5.376mm/sec. −13 4.4.2一维非稳态等温渗流 一维非稳态等温渗流所要讨论的问题示意如图4-12。一块固气多孔材料内气体压力pa0大于环境压力pe的平板，其上、下表面盖有覆盖板，使气体无法从上、下表面与环境连通，压力高于pe的气体只能通过板的侧面流入环境。 zoxW 图4-12 一维固气多孔材料非稳态等温渗流示意图 为便于分析，不计沿板长度L方向(垂直纸面y方向)的流动，只计沿板宽度W方向 流动，则问题简化为一维非稳态等温渗流。对此，先写出控制方程。 (1)问题的控制方程 对于一维非稳态等温渗流这一情形，4.2节图4-5所示的关于固气多孔材料传热传质 17 问题的7个因变量中只保留关于空气速度u、空气压力pa与空气密度ρa三个因变量。它们相应的控制方程为： 1.关于空气速度u的式(4.1.4)，即 u=− K∂pa (4.1.4) μa∂x 2.关于空气压力pa的式(4.2.13) ，即 Pa= ρaRT (4.2.13) Ma 3.关于空气密度ρa的质量守恆方程，即 ∂φρa∂ρu =−a (4.4.5) ∂τ∂x 式(4.4.5)是式(4.2.3)在一维条件下的简化形式。 式(4.1.4)、式(4.2.13)及(4.4.5)构成本问题的控制方程组。相应的定解条件为： 1.几何条件：考虑到本问题所展示的现象对坐标轴z对称，故求解的几何域限定在从x=0到x=W/2平板的右半部分内。 2.初始条件：空气速度u的初始值为0，空气压力pa的初始值pa=pa0，空气密度ρa的初始值可由pa0按式(4.2.13)算得，ρa=对空气密度ρa无需给出相应的边界条件。 4.物理条件：本问题中的物理参数，即，渗透率K与空气粘性系数μa均为常数。 综上所述，本问题涉及3个因变量，对应有3个控制方程，由它们构成的方程组是封闭的。根据相应的定解条件即可求得一组关于u、ρa与pa的唯一解，也即多孔平板内确定的空气速度场、密度场与压力场。 应当指出，关于u、ρa与pa的三个方程，它们之间是相互耦合的，由图示意如下： 在这个耦合的方程组中，任何一个方程都不能单独求解。正是这种耦合关系使它们难以用解析的方法求解。差分方程的方法将比较容易解决这一难点。 要得到对应于上述控制方程及定解条件的差分方程，首先需将求解区域离散化。 (2) 区域离散化 与讨论一维非稳态导热问题相似，将求解域(0~W/2)划分成等步长Δx的有限个单元IX，即，W/2=IX⋅Δx。对于每个单元可定义压力pa与密度ρa分别为pa,i与ρa,i， ρa0。 3.边界条件：空气速度u在x=0处为0，空气压力pa对应的环境空气压力为pe， ρa paui=1,2,...IX。每个单元对应一个节点，节点置于单元的中心，用符号圆圈标出，如图4-13 所示。 对空气速度u离散化ui的处理有别于压力pa,i与密度ρa,i，它被定义在单元的边界 18 上，其节点用小方框标出，如图4-13所示。从式(4.1.4)看到，速度u与压力pa对x的一阶导数正成比，当一阶导数写成差分形式之后(见图2-15)，速度u应该在两个压力值pa,i−1与pa,i之间。因此，将ui定义在两个压力单元的交界面上也是非常合理的。如果将压力与密度对应的单元所建立的网格系统称为主网格的话，则速度所对应的网格与主网格相差半个空间步长，称这种两套网格为交错网格[20,21,22]。 uiui+1uIX+1 u 1 x pppa,1ppa,ia,i−1Δx a,IXa,i+1 ρa,i−1ρa,iρa,i+1ρa,IX ρa,1x =0x=L/2 图4-13一维非稳态等温渗流的单元网格 引入交错网格是保证计算多孔材料内气体速度取得成功的重要条件。 对时间区域离散化的处理遵循如下原则。在待解方程组的每个方程中，只有一个待求因变量有上标n+1，表示待求其下一时刻的值，其它因变量的上标均为n，即当前时刻的值，它们均已知。这样，使整个方程组中每个方程都呈显式。就方程组整体而言，相互是耦合的，每一个因变量依赖其它因变量，但求解的每一步，方程间的耦合被分开了。这就是使用数值求解的优越性。 (3)差分方程 速度u的式(4.1.4)相应的差分方程为 nnp−pKa,ia,i−1n+1 内单元 ui=− μaΔx i=2,3,...,IX (4.4.6) 左边界 ui右边界 u n+1 =0 i=1 (4.4.6a) nKpa0−pa,i−1 i=IX+1 (4.4.6b) =− μaΔx n+1 i 关于密度ρa的式(4.4.5)相应的差分方程可通过用差商代替微商的方法获得，即 n+1nnnnn −u−uρaρρρa,ia,ii+1,i φ (4.4.7) =a,i−1i ΔτΔx 经整理后即可为 n+1n ρa,i=ρa,i+ Δτnnnn (ρa,i−1ui−ρa,iui+1) i=1,2,...,IX (4.4.7a) φΔx 关于压力pa的式(4.2.13)相应的差分方程为： P n+1a,i = n+1ρa,iRT Ma (4.4.8) 上述式(4.4.6)、(4.4.7)和(4.4.8)中各物理量右上角的n与n+1表示各物理量在当 前时刻(n)与经过一个时间步长Δτ后的下一时刻(n+1)的值。 综上所述，式(4.4.6,a,b)、式(4.4.7a)和式(4.4.8)即为本小节一维非稳态等温渗流 19 问题的差分方程组。在初始条件给定后，即可由已知的当前时刻各物理量的值求得在区域(0~W/2)内下一时刻各物理量的值。也就是，通过差分方程可求得固气多孔材料内空气速度场、密度场与压力场。 这里需指出，在上述差分方程都表示成显式的格式，即求任何一个物理量在n+1时 刻的值时，将其它变量都用n时刻的值。这样，把原本相互耦合的问题容易地被解耦了，使计算顺利进行。但是，这种显式格式，对时间步长Δτ与空间步长Δx的选取都很严格。 这里还需指出，上述式(4.4.7)差分方程是通过用差商代替微商由微分方程(4.4.5)得到。其实，该差分方程还可通过控制容积法来建立。在很多情形下，控制容积法有显著的优点。见图4-13中单元i，该单元在Δτ时间内空气质量的增加(下式等号左边)等于Δτ时间内通过单元边界流入的空气质量(下式等号右边)，表示如下： n+1nnnnnφ(ρa,i−ρa,i)Δx=(ρa,i−1ui−ρa,iui+1)Δτ (4.4.9) 将上式稍加整理即可获得式(4.4.7)。这里，ρa,i−1uiΔτ和ρa,iui+1Δτ分别表示在Δτ时间内通过单元i左边界流入与通过单元i右边界流出单元的质量。鉴于本节讨论的问题，是多孔 材料平板内含有较高压力的空气通过平板右侧面向环境等温渗流，在等温渗流过程中，空气速度u始终沿x方向流动。在流动方向清楚的前提下，比较容易获得如式(4.4.7)的差分方程。然而，很多实际问题中，流动的方向不能事先确定，尽管仍可用差商代替微商的方法从数学上获得差分方程，但其完全可能违背实际物理过程。用控制容积法则比较容易建立适应流动方向不确定的差分方程。对此，在本章4.5节中有详尽介绍。 (4)控制方程的另一种形式 上述式(4.1.4)、式(4.2.13)及式(4.4.5)所构成的方程组也可合并成一个方程。为 此，将式(4.4.5)改写成 φnnnn ∂ρa∂u∂ρ=−ρa−ua (4.4.10) ∂τ∂x∂x 将式(4.1.4)与式(4.2.13)代入式(4.4.10)，稍加整理后可得 ∂paK∂2p∂pa2 [pa2+()] (4.4.11) φ= ∂τφμa∂x∂x 这个关于压力的方程与上述3个方程组成的方程组是等价的。对它进行差分处理所得到的差 分方程与上述差分方程(4.4.6)、(4.4.7)与(4.4.8)也是等价的。 (5)算例及计算结果讨论 本算例的已知条件有：板宽W=1m，材料内初始气体压力pa0=2⋅10Pa，多孔材 5 m2，空气粘性系数为μa=18.6Pa⋅sec.，并取空间 步长Δx=W/(2IX)=5cm，即右半平板内离散的单元数IX=10。取时间步长Δτ=0.01sec.。 按式(4.4.6,a,b)、(4.4.7a)和(4.4.8)计算得到空气速度ui、空气压力pa,i及空气 料的孔隙率φ=0.5，渗透率K=10 密度ρa,i。现将计算结果现示在如下图4-14, 图4-15与图4-16。 图4-14显示平板内不同地点空气速度ui随时间变化规律。由图可见，当初始气体 压力大于环境压力的平板突然放入该环境后，除侧表面处(x=0.5m)的速度从最大值迅速下降外，处于内部各处的气体速度都经历由0逐渐加大待达到最大值后再缓慢下降。各处所能达到的速度最大值从侧表面到中心由大变小。 若把图4-14中速度变化规律与第2章中图2-13所示一维非稳态导热中热流密度的 变化规律进行对比，可以看到，它们两者曲线的走向非常相似。 20 −13 111098u(mm/sec.)763210020406080100time(sec.)120140160180200x=0.05mx=0.2mx=0.35mx=0.5m 图4-14 平板内不同地点空气速度随时间变化规律 图4-15显示平板内不同地点空气压力pa,i随时间变化规律。由图可见，除表面处 (x=0.5m)的压力迅速下降外，处于内部各处的气体压力都是先缓慢而后逐渐加快最后又趋于缓慢下降的过程。 21.91.81.7p(100kPa).)1.61.51.41.31.21.11020406080100time(sec.)120140160180200x=0.05mx=0.2mx=0.35mx=0.5m 图4-15平板内不同地点空气压力随时间变化规律 图4-16是将地点作为横坐标，时间作为参变量绘制的气体压力分布随时间的变化。由图可见，气体压力的变化在开始阶段主要在表面附近，而逐渐地表现在板的中心部位。 2.121.91.81.7p(100kPa).)1.61.51.41.31.21.11-0.5-0.4-0.3-0.2-0.10x(m)0.10.20.30.40.50 5 10 20 40 80 120 200 图4-16平板内空气压力分布随时间变化 4.4.3二维非稳态等温渗流 本节所要讨论的二维非稳态等温渗流问题与4.4.2小节相似，只是将图4-12中盖在 平板上、下表面的板拿掉，如图4-17。一块固气多孔材料内气体压力pa0大于环境压力pe的平板，置于环境后，压力高于pe的气体通过板的上、下表面与左、右侧面流入环境，仍不计沿板长度L方向(垂直纸面y方向)的流动，只计沿板宽度W与厚度H方向流动，则问题归结为二维非稳态等温渗流。对此，先写出控制方程。 21 zoxWH 图4-17 二维非稳态等温渗流示意图 (1)问题的控制方程 对于二维非稳态等温渗流这一情形，上述图4-5中所示的关于固气多孔材料传热传质问题的7个因变量中需保留关于空气速度u、w、空气压力pa与空气密度ρa4个因变量。它们相应的控制方程为，关于空气速度u与w，即将式(4.1.4)改写为 u=− 和 Kx∂pa (4.4.9) μa∂x Kz∂pa (4.4.10) μa∂z w=− 关于空气压力pa的式(4.2.13)，即， Pa= ρaRT Ma (4.2.13) 以及关于空气密度ρa的质量守恆方程： ∂φρa∂ρu∂ρw =−a−a (4.4.11) ∂τ∂x∂z 式(4.4.11)是式(4.2.3)在二维条件下的简化形式。 式(4.4.9)、式(4.4.10)、式(4.2.13)及(4.4.11)构成本问题的封闭的控制方程组。相应的定解条件为： 1.几何条件：考虑到本问题所展示的现象对坐标轴z与坐标轴x对称，故求解的几何域限定在从x=0到x=W/2，从z=0到x=H/2，四分之一的平板内。 2.初始条件：空气速度u与w的初始值为0，空气压力pa的初始值pa=pa0，大于环境压力pe，空气密度ρa的初始值可由p0按式(4.2.13)算得，ρa=境空气压力为pe，对空气密度ρa无需给出相应的边界条件。 4.物理条件：本问题中的物理参数，即，沿x方向渗透率Kx，沿z方向渗透率Kz与空气粘性系数μa均为常数。 综上所述，本问题涉及4个因变量，对应有4个控制方程，由它们构成的方程组是封闭的。根据相应的定解条件即可求得一组关于u、w、ρa与pa的唯一解，也即多孔平板内确定的空气速度场、密度场与压力场。 关于u、w、ρa与pa的四个方程，它们之间是相互耦合的，由图示意如下： ρa0。 3.边界条件：空气速度u在x=0处为0，w在z=0处为0，空气压力pa对应的环 22 (2)差分方程 对应于上述控制方程及定解条件，可得到相应的差分方程。先讨论求解区域的离散化。 与讨论一维非稳态等温渗流问题相似，将区域(0~W/2)划分成等步长Δx的有限个 单元IX，即，W/2=IX⋅Δx基础上，再将区域(0~H/2)划分成等步长Δz的有限个单元JZ，即，H/2=JZ⋅Δz。这样，在1/4平板的求解域内总计有单元数为IX⋅JZ个，对于每个单元可定义压力pa与密度ρa分别为pa,i,j与ρa,i,j，i=1,2,...IX，j=1,2,...JZ。 与用交错网络法来定义沿x方向速度ui相同，可定义沿z方向速度wi如图4-18。 ρa,i,j+1,pa,i,j+1 wi,j+1 ρa,i,j,pa,i,jρa,i−1,j ui+1,jui,j pa,i−1,j wi,j ρa,i,j−1,pa,i,j−1 图4-18 单个单元中密度、压力与速度的节点 式(4.4.9) 相应的差分方程为 ρa pau,wρa,i+1,j pa,i+1,j u n+1i,j nn Kxpa,i,j−pa,i−1,j=− i=2,3,...,IX; j=1,2,...,JZ (4.4.12) μaΔx 1 uin,+j=0 i=1; j=1,2,...,JZ (4.4.12a) np−pKa0a,i−1,j1x uin,+ i=IX+1; j=1,2,...,JZ (4.4.12b) j=− μaΔx 式(4.4.10) 相应的差分方程为 w n+1 i,j nn Kzpa,i,j−pa,i,j−1 i=1,2,...,IX; j=2,3,...,JZ (4.4.13) =− μaΔz 1 win,+j=0 i=1,2,...,IX; j=1 (4.4.13a) np−pK0,i,j−1aa1z i=1,2,...,IX; j=JZ+1 (4.4.13b) win,+=−j μaΔz 式(4.4.11)相应的差分方程为 23 n+1nnnnnnnnnρaρaρa,i,j−ρa,i,j,i−1,jui,j−ρa,i,jui+1,j,i,j−1,wi−ρa,i,jwi,j+1 φ (4.4.14) += ΔzΔxΔτ式(4.4.14)经整理后也可为 ρn+1 a,i,j =ρna,i,j nnnnnnnn Δτρa,i−1,jui,j−ρa,i,jui+1,jρa,i,j−1wi,j−ρa,i,jwi,j+1 () ++φΔzΔx i=1,2,...,IX j=1,2,...,JZ (4.4.14a) 关于空气压力pa的式(4.2.13)相应的差分方程(4.4.8)仍然有效，保持不变。 综上所述，式(4.4.12,a,b)、式(4.4.13,a,b)、(4.4.14b)和(4.4.8)即为本小节二 维非稳态等温渗流问题的差分方程组。在初始条件给定后，即可由已知的当前时刻各物理量的值求得在区域0 本算例的已知条件有：板宽W=1m，板厚H=2cm，材料内初始气体压力 pa0=2⋅105Pa，多孔材料的孔隙率φ=0.5，渗透率K=10−13m2，空气粘性系数为μa=18.6Pa⋅sec.，并取空间步长Δx=W/(2IX)=5cm与Δz=H/(2JZ)=2mm，取时间步长Δτ=0.01sec.。根据本节由4个控制方程(4.4.9)、(4.4.10)、(4.2.13)及(4.4.11) 构成的定解问题而得到的由式(4.4.13)、(4.4.13a)、(4.4.13b)、(4.4.14a)与(4.4.8)组成的差分方程组可算得因变量空气速度u、w、空气压力pa与空气密度ρa随时间的变化。解的离散形式是，ui,j、wi,j、pa,i,j与ρa,i,j。其中，pa,i,j与ρa,i,j在等温条件下一一对应，它们的变化规律相似。这里通过图4-19及图4-20显示ui,j、vi,j与pa,i,j在平板中心处(x=0，z=0，也即，i=1，j=1)的变化规律，扼要说明二维非稳态等温渗流的特点。 图4-19是平板中心处u和w随时间变化规律。当多孔材料平板为各向同性时，z方向的渗透率等于x方向的渗透率，Kz=Kx，则z方向的速度w远大于x方向的速度u。而当多孔材料平板为各向异性，z方向的渗透率远小于x方向的渗透率，Kz=0.01Kx，则 n n n n n n n n n z方向的速度w远小于x方向的速度u。由此可见，对于二维问题，除了几何特征(宽度远 大于厚度)之外，更重要的是材料的性能的方向性。对于各向异性的木材，顺纹理方向与垂直纹理方向的渗透率的大小相差达3个数量级。所以，流体还是更容易沿纹理方向流动。 0.30.25u,w(mm/sec.)0.20.150.1u(kz=kx)w(kz=kx)u(kz=0.01kx)w(kz=0.01kx)0.050050100time(sec.)150200 图4-19平板中心处u和w随时间变化 图4-20是平板中心处压力变化曲线。由图可见，当Kz=Kx，中心压力迅速下降。 24 这是因为，气体大量沿z方向流向环境。当Kz 图4-20 平板中心处空气压力pa随时间变化 4.5固气多孔材料热板加热 本节是在4.3节讨论多孔材料热传导及在4.4节讨论等温渗流基础上，把两者结合起来，讨论孔隙连通的固气多孔材料的加热，如图4-20所示。 z Tp H o x Tp W 图4-20 二维固气多孔材料热板加热示意图 当上、下两块热板同时与固气多孔材料平板紧紧相接触时，热板对材料加热。与此同时，引起材料内各处孔隙中的空气压力有不同程度的升高，进而引起空气由高压处向低压处流动。空气的流动使材料内的传热过程在热传导的基础上又增添了热对流。由此，材料内部既有温度场又有空气速度场。本节介绍这一问题的控制方程与定解条件，相应差分方程及其求解，最终举例演示温度场与空气压力、密度及速度场。 4.5.1控制方程及定解条件 本问题涉及的因变量，除了4.4.3小节中关于二维等温非稳态渗流中的空气速度u、w、空气压力pa与空气密度ρa4个因变量之外，再加上因变量温度T。它们相应的控制方程罗列如下： 关于空气速度u与w，即式(4.4.9)与式(4.4.10)，即 u=− 和 Kx∂pa (4.4.9) μa∂x Kz∂pa (4.4.10) μa∂z 25 w=− 关于空气压力pa的式(4.2.13) ，即 Pa= ρaRT Ma (4.2.13) 关于空气密度ρa的质量守恆方程为 ∂ρu∂ρw∂φρa =−a−a (4.4.11) ∂z∂x∂τ以及关于温度T的能量守恆方程： ∂ρcuT∂ρcwT∂2T∂2T∂ρcT （aa+aa） (4.5.1) =k(2+2)− ∂z∂x∂z∂x∂τ式(4.5.1)是式(4.2.12)在二维条件下的简化，其中 k=(1−φ)ks+φka (4.2.10) ρc=(1−φ)ρscs+φρaca (4.2.11) 式(4.4.9)、式(4.4.10)、式(4.2.13)、(4.4.11)及(4.5.1)构成本问题的封闭的控 制方程组。相应的定解条件为： 1.几何条件：本问题坐标原点设在材料平板的左下角，故求解的几何域限定在0≤x≤W，0≤z≤H整个平板内。这里是一块二维平板，不计各物理量沿平板长度y方向的变化。 2.初始条件：空气速度u与w的初始值为0，平板材料温度T的初始值T0与平板周围环境的温度Te相同，即T=T0=Te。空气压力pa的初始值pa=pa0，它与平板周围环境的压力pe相同，即pa=pa0=pe。空气密度ρa的初始值可由pa0与T0按气体状态方程(4.2.13)求得，ρa= ρa0。 3.边界条件：对于空气速度u，无需提供边界条件，只提供空气环境压力pe即可。对于空气速度w，在z=0与z=H处w=0。对于空气压力pa也无需单独提供边界条件。对于空气密度ρa相应的边界条件为环境空气密度ρa0。对于温度T，其边界条件为 ∂T (4.5.1a) ∂x∂T x=W h(Te−T)=k (4.5.1b) ∂x z=0 T=TP (4.5.1c) z=H T=TP (4.5.1d) 这里，TP是上、下热板温度，它为常数。Te是环境温度，h是材料平板侧面与相邻环境空 x=0 h(Te−T)=−k 气之间的自然对流传热系数。 4.物理条件：沿x方向渗透率Kx与导热系数kx，沿z方向渗透率Kz与导热系数kz，它们均为常数。空气比热容ca与粘性系数μa均为温度的函数，详见4.5.4小节。 综上所述，本问题涉及5个因变量，对应有5个控制方程，由它们构成的方程组是封闭的。根据相应的定解条件即可求得一组关于u、w、ρa、pa与T的唯一解，也即多孔平板内确定的空气速度场、密度场、压力场以及温度场。 因变量u、w、ρa、pa与T对应的5个方程，它们之间相互耦合，再加上平板材料的密度ρ随空气密度ρa的变化有稍微的变化，则总计6个因变量。这6个因变量相互间的 26 耦合关系已由图4-5所示。 4.5.2差分方程 对应于上述控制方程及定解条件，可得到相应的差分方程。在推导差分方程前，先讨论求解区域的离散化，然后逐个讨论速度u与w，密度ρa，温度T以及压力pa的差分方程。 (1)区域离散化 将求解区域0≤x≤W，0≤z≤H划分成几何尺寸完全相同的有限个单元，如图4-21所示。单元的宽与高均为Δx与Δz，且平板的宽度W=IX⋅Δx，高度H=JZ⋅Δz。这样，平板总计有单元数为IX×JZ个，单元对应的节点都位于单元的，在图中用黑色小圆圈标出。对于每个单元可定义压力pa、密度ρa与温度T，分别为pa,i,j、ρa,i,j与Ti,j， i=1,2,...IX，j=1,2,...JZ。 i,j j i 图4-21 二维固气多孔材料平板离散化单元与节点 如图4-21a所示。这样，ui,j对 用交错网络法来定义沿x方向速度ui,j与沿z方向速度wi,j， 应的节点数为IX+1，wi,j对应的节点数为JZ+1。节点的位置在单元的边界上，用小方框标出。 ρa,i,j+1,pa,i,j+1,Ti,j+1 wi,j+1 ρa,i−1,jρa,i+1,jρa,i,j,pa,i,j,Ti,j pa,i−1,j pa,i+1,jΔz uui,j i+1,jTTi+1,ji−1,j wi,j ρa,i,j−1,pa,i,j−1,Ti,j−1 Δx 图4-21a 单个离散单元中温度、压力、密度与速度的节点 (2)关于速度u与w的差分方程 关于速度u的方程(4.4.9) 所相应的差分方程为 内部单元： 27 nnp−pKa,i,ja,i−1,j1x i=2,3,...,IX; j=1,2,...,JZ (4.5.2) uin,+j=− μa,i,jΔx 左边界： u n+1 i,j n Kxpa,i,j−pe i=1; j=1,2,...,JZ (4.5.2a) =− μa,i,jΔx 和右边界： u n+1 i,j n Kxpe−pa,i−1,j i=IX+1; j=1,2,...,JZ (4.5.2b) =− μa,i,jΔx 应该指出，由于传热的缘故，空气粘性系数μa不再是常数，它是温度的函数。为 此，有μa=μa,i,j，而μa,i,j可为 1 [μa(Ti,j)+μa(Ti−1,j)] 2 表示流体从单元i−1流到单元i时所具有的粘性系数是温度为Ti−1,j与Ti−1,j的粘性系 μa,i,j= 数的算术平均值。 关于速度w的方程(4.4.10)所相应的差分方程为： 内部单元： nnp−pKaija,,,i,j−11z i=1,2,...,IX; j=2,3,...,JZ (4.5.3) win,+=−j μa,i,jΔz 其中μa,i,j可为 μa,i,j=下边界： 1 win,+j=0 i=1,2,...,IX; j=1 (4.5.3a) 1 [μa(Ti,j)+μa(Ti,j−1)] 2 和 上边界： wi,j=0 i=1,2,...,IX; j=JZ+1 (4.5.3b) (2)关于空气密度ρa的差分方程 根据微分方程式(4.4.11)，可得 n+1 ΔρuΔρwΔφρa =−a−a (4.5.4) ΔzΔxΔτ式(4.5.4)经整理后为 ΔτΔρauΔρaw (+) (4.5.5) φΔxΔz i=1,2,...,IX j=1,2,...,JZ Δρau 式(4.5.5)中的 −表示单位时间内沿x方向流入单位空间的空气质量。其差分形式应 Δx 视u的方向不同而不同，可按如图4-22所示(a)、(b)、(c)、(d)四种情形进行分析。 n+1nρa,i,j=ρa,i,j− 28 ρa,i+1,j ρρρa,i−1,jρa,i,ja,i−1,jρa,i+1,ja,i,juui+1,jui,j u i,ji+1,j (a)左、右两侧都向右 (b)左侧向左，右侧向右 ρ ρρa,i−1,jρ ρa,i−1,jui,jρa,i,jui+1,ji+1,jui,j a,i,jui+1,j a,i+1,j (c) 左侧向右，右侧向左 (d) 左、右两侧都向左 图4-22 水平方向对流项的差分处理 情形(a) −Δρau=情形(b) −Δρau=情形(c) −Δρau=情形(d) −Δρau=式(4.5.5)中的 − nnnnρau−ρua,i,ji+1,j (4.5.6a) ,i−1,ji,jnnnnρa,i,jui,j−ρa,i,jui+1,j (4.5.6b) nnnnρau−ρua,i+1,ji+1,j (4.5.6c) ,i−1,ji,j nnnnρa,i,jui,j−ρa,i+1,jui+1,j (4.5.6d) Δρaw Δx ，表示单位时间内沿z方向流入单位空间的空气质量。其差 分形式应视w的方向不同而不同，可按如图4-23所示四种情形进行分析。 ρa,i,j+1 ρa,i,j+1 wi,j+1wi,j+1 ρa,i,j ρa,i,j wi,j wi,j ρa,i,j−1 ρa,i,j−1 (a) 上、下两侧都向上 (b) 下侧向下，上侧向上 ρa,i,j+1 ρa,i,j+1 wi,j+1 wi,j+1 ρa,i,jρa,i,j wi,jwi,j ρa,i,j−1 ρa,i,j−1 (c) 下侧向上，上侧向下 (d)上、下两侧都向下 图4-23垂直方向对流项的差分处理 情形(a)： −Δρaw=情形(b)： −Δρaw= nnnnρaw−ρwa,i,ji,j+1 (4.5.7a) ,i,j−1i,jnnnnρa,i,jwi,j−ρa,i,jwi,j+1 (4.5.7b) 29 情形(c)： −Δρaw=情形(d)：−Δρaw= nnnnρa,i,j−1wi,j−ρa,i,j+1wi,j+1 (4.5.7c) nnnnρa,i,jwi,j−ρa,i,j+1wi,j+1 (4.5.7d) 综上所述，空气密度ρa的差分方程(4.5.5)及(4.5.6)与(4.5.7)是通过对微分方程(4.4.11)进行差分处理获得的。其实，这些差分方程也可用控制容积法对图4-21a所示的单 元直接建立关于气体质量平衡的代数方程来获得，所得最终结果相同。 (3)关于温度T的差分方程 关于温度的差分方程可以通过对微分方程(4.5.1)进行差分处理来获得，也可用控制容积的方法。这里，我们用控制容积法来推导对应于微分方程(4.5.1)的差分方程。如图4-21所示，全部IX×JZ个单元可分为两类。一类是内单元，这些单元相邻的都是平板内其它单元。另一类是边界单元，这些单元的边界至少有一边不与其它单元相邻。边界单元中又可分两类，一类是上、下边界单元，它们与热板直接接触。另一类是左、右两侧边界单元，它们与环境相邻。下面分别就内单元，上、下边界单元及左、右边界单元介绍它们的差分方程。 内单元： 视图4-21a为一控制容积，或称控制单元。对它可建立能量平衡关系如下： nnnn Ti−Ti+Ti,nj−1−Ti,njTi,nj+1−Ti,nj1,j−Ti,j1,j−Ti,j Δz+kxkxΔz+kzΔx+kzΔxΔxΔxΔzΔz (4.5.8) n+1n (ρcTi,j)−(ρcTi,j) +Δ(ρacauT)Δz+Δ(ρacawT)Δx=ΔxΔz Δτi=2,3,⋅⋅⋅,IX−1 j=2,3,⋅⋅⋅,JZ−1 上式等号左边前四项，表示该单元在单位时间内通过单元的四个边界以热传导方式到达该单元的热量。其中前两项与后两项分别为沿x与沿z方向进入单元的。上式等号左边第5、6项，表示该单元在单位时间内通过单元的四个边界以热对流方式到达该单元的热量。第5与第6项分别为沿x与沿z方向进入单元的。由于热量以对流方式进入，它们的差分表示与气体速度u与w有关，下面将给予详细说明。现将式(4.5.8)重新整理，可得： Ti,j n+1 ={[kx ni,j+1 nnTi−1,j−Ti,j Δx ni,j Δz+kx nnTi+1,j−Ti,j Δx Δz+kz Ti,nj−1−Ti,nj Δz Δx (4.5.9) +kz T−T Δz Δx−Δ(ρacauT)Δz−Δ(ρacawT)Δx] Δτn +(ρc)nTi,j}/(ρc)n+1 ΔxΔz 式中，表示对流项的−ΔρacauT与−ΔρacawT其具体表示与上述式(4.5.5a,b,c,d) 与(4.5.6a,b,c,d)相仿，应按图4-22与图4-23所示的各四种情形予以准确表示。 1.关于−ΔρacauT 情形(a) −ΔρacauT=情形(b) −ΔρacauT=情形(c) −ΔρacauT=情形(d) −ΔρacauT=情形(a) −ΔρacawT=情形(b) −ΔρacawT=情形(c) −ΔρacawT= nnnnnnnnρa,i−1,jca,i−1,jTa,i−1,jui,j−ρa,i,jca,i,jTa,i,jui+1,j (4.5.9a) nnnnnnnnρa,i,jca,i,jTa,i,jui,j−ρa,i,jca,i,jTa,i,jui+1,j (4.5.9b) nnnnnnnnρa,i−1,jca,i−1,jTa,i−1,jui,j−ρa,i+1,jca,i+1,jTa,i+1,jui+1,j(4.5.9c) nnnnnnnnρa,i,jca,i,jTa,i,jui,j−ρa,i+1,jca,i+1,jTa,i+1,jui+1,j (4.5.9d) 2.关于−ΔρacawT nnnnnnnnρa,i,j−1ca,i,j−1Ta,i,j−1wi,j−ρa,i,jca,i,jTa,i,jwi,j+1 (4.5.9e) nnnnnnnnρacTw−ρcTwa,i,ja,i,ja,i,ji,j+1 (4.5.9f) ,i,ja,i,ja,i,ji,jnnnnnnnnρa,i,j−1ca,i,j−1Ta,i,j−1wi,j−ρa,i,j+1ca,i,j+1Ta,i,j+1wi,j+1 (4.5.9g) 30 情形(d) −ΔρacawT=nnnnnnnnρa,i,jca,i,jTa,i,jwi,j−ρa,i,j+1ca,i,j+1Ta,i,j+1wi,j+1 (4.5.9h) 上、下边界单元： 上边界单元是材料平板上表面与热板紧密接触的单元，它们与内单元差分方程的差别 只是用Tp替代了Ti,j+1，且相应地用Δz/2替代了Δz，因为热板与上表面单元所对应节点的间距为Δz/2。这样，上边界单元差分方程为 nn TinTinTi,nj−1−Ti,nj−1,j−Ti,j+1,j−Ti,jTi,j Δz+kxΔz+kzΔxTi,j={[kx ΔxΔxΔz (4.5.10) n Tp−Ti,jΔτn+kzΔx−Δ(ρacauT)Δz−Δ(ρacawT)Δx]+(ρc)nTi,j}/(ρc)n+1 Δz/2ΔxΔz 下边界单元与上边界单元相仿，其与内单元差分方程的差别只是用Tp替代Ti,j−1，且相应地用Δz/2替代了Δz，所得差分方程为 nnnnTi−Ti+Tp−Ti,njn+11,j−Ti,j1,j−Ti,j Δz+kxΔz+kzΔxTi,j={[kx ΔxΔxΔz/2 (4.5.11) nn −TTΔτn +kzi,j+1i,jΔx−Δ(ρacauT)Δz−Δ(ρacawT)Δx]+(ρc)nTi,j}/(ρc)n+1 ΔzΔxΔzn+1 边界单元方程中所包含的对流项与内单元中相同，仍可用上述式()表示。 左、右边界单元 左边界单元是材料平板左侧面与环境相邻的单元，其与内单元差分方程的差别只是来自左边的热传导kx Ti−1,j−Ti,j Δz用对流传热h(Te−Ti,j)Δz替代，则差分方程为 Δx nnTi+Ti,nj−1−Ti,njn+11,j−Ti,jn Ti,j={[h(Te−Ti,j)Δz+kxΔz+kzΔx ΔxΔz (4.5.12) nnT−TΔτn +(ρc)nTi,j}/(ρc)n+1+kzi,j+1i,jΔx−Δ(ρacauT)Δz−Δ(ρacawT)Δx] ΔxΔzΔz 右边界单元是材料平板右侧面与环境相邻的单元，其与内单元差分方程的差别只是来 自右边的热传导kx n+1 Ti−1,j−Ti,j Δx Δz用对流传热h(Te−Ti,j)Δz替代，则差分方程为 nnnn−TTTi−1,j−Ti,ji,j−1i,j Δz+h(Te−Ti,nj)Δz+kzΔxTi,j={[kx ΔxΔz (4.5.13) nnTi,j+1−Ti,jΔτn+kzΔx−Δ(ρacauT)Δz−Δ(ρacawT)Δx]+(ρc)nTi,j}/(ρc)n+1 ΔzΔxΔz 上述由式(4.5.9)、(4.5.10)、(4.5.11)、(4.5.12)与(4.5.13)联同式(4.5.9a-h)构 成温度T的差分方程。它们可统一写成如下形式： Ti,nj+1=(AW⋅Ti−1,j+AE⋅Ti+1,j+AD⋅Ti,j−1+AU⋅Ti,j+1+CA⋅Ti,j+DA)/EA(4.5.14) 式(4.5.14)表明，材料平板每个单元在某一时刻的温度都可表示为它周围四边相邻单元(包括热板与环境)温度及该单元本身前一时刻温度的线性组合。 对照式(4.5.14)与前述差分方程，可得到不同条件下式(4.5.14)中各项系数AW,AE,AD,AU,CA与EA的表达式，进而得到它们的数值。 (4) 关于空气压力pa的差分方程 空气压力pa的式(4.2.13)相应的差分方程(4.4.8)仍然有效，保持不变。 31 综上所述，式(4.5.2,a,b)、式(4.5.3,a,b)、(4.5.7)、(4.5.14)和(4.4.8)即为本小节固气多孔材料热板加热问题的差分方程组。在初始条件给定后，即可由已知的当前时刻各物理量的值求得在区域0≤x≤W，0≤z≤H内下一时刻各物理量的值。也就是，通过差分方程可求得固气多孔材料内空气速度场、密度场、压力场及温度埸。 4.5.3算例及计算结果讨论 本算例的已知条件有：板宽W=1m，板厚H=2cm，多孔材料的孔隙率φ=0.5，渗透率K=10 −13 m2，空气粘性系数μa与比热容ca为温度的函数， μa=17.326+0.0488T，单位为Pa⋅sec.T为0C，ca=1005.3−0.0026T−0.0006T2，单位为J/(kg⋅K)。多孔材料上、下表面在时间τ≥0与热板紧密接触，热板温度Tp为常数， Tp=2000C。环境温度Te=200C，材料初始温度与环境温度相同T0=200C，环境空气 的压力pe为大气压，pe=100kPa，环境空气密度ρe可根据T0与pe按理想气体状态方程算得。材料内初始气体压力pa0与密度ρa0为平板周围的环境压力pe与密度ρe。环境与多孔材料平板侧表面为自然对流冷却，自然对流传热系数h=5W/(m⋅K)。取空间步长 2 Δx=W/IX=1000/12=83.333mm与Δz=H/JZ=20/10=2mm，取时间步长Δτ=0.01sec.。 根据本节由5个控制方程(4.4.9)、(4.4.10)、(4.2.13) 、(4.4.11) 及(4.5.1)构成 的定解问题而得到的由式(4.5.2,a,b)、式(4.5.3,a,b)、(4.5.7)、(4.5.14)和(4.4.8)组成的差分方程组可算得因变量空气速度u、w、空气压力pa、空气密度ρa与温度T随时间的变化。解的离散形式是，ui,j、wi,j、pa,i,j、ρa,i,j与Ti,j(n=0,1,2,⋅⋅⋅)。计算结果以2D、1D与0D形式给出，它们的含义由图4-24说明。 b a 图4-24 二维平板若干特定离散单元与节点 n n n n n n n n n n d c 2D形式：2D表示ui,j、wi,j、pa,i,j、ρa,i,j与Ti,j这5个因变量在整个二维(2D)平板每个单元随时间的变化。它们将通过界面以图示形式给出。 1D形式：1D表示ui,j、wi,j、pa,i,j、ρa,i,j与Ti,j这5个因变量在整个二维平板内以特定的一维(1D)形式给出。这一维形式对pa,i,j、ρa,i,j与Ti,j这3个因变量而言为平板的水平中心轴与垂直中心轴，图4-24中用粗虚实线表示。对于速度ui,j与wi,j则兼用粗实线与点划线表示。 0D形式：0D表示ui,j、wi,j、pa,i,j、ρa,i,j与Ti,j这5个因变量在整个二维平板内仅以4个代表性单元的值随时间的变化。这4个代表性单元在平板中的位置由图4-24中4个深色方框a,b,c,d表示。 1D与0D的输出结果都以曲线形式表示。 下面我们通过图4-25至图4-29以2D形式以及图4-30、图4-31以0D形式给出ui,j、wi,j、 nnnpa,i,j、ρa,i,j与Ti,j这5个因变量的计算结果并进行分析讨论。 n n n n n n n n n n n n n n n n n 32 如前所述，固气多孔材料平板加热问题是4.3节平板热传导与4.4节等温渗流的结合。所以，在分析本算例计算结果时将对比前二节的结果。为便于讨论，将ui,j、wi,j、pa,i,j、 nnρa,i,j与Ti,j这5个因变量的计算结果逐个进行分析。 n n n (1)关于温度Ti,j 图4-25系列表示固气多孔材料受热板加热条件下材料的温度响应。图4-25a与图4-25b是垂直中心轴上单元的温度响应。图4-25a为平板从被加热后若干时刻(τ=0、5、15、50、100、150、200、250及300sec.)垂直中心轴温度分布。 200180160temperature(C)1401201008060402000510distance(mm)152005 15 50 100 150 200 250 300 n 图4-25a 垂直中心轴温度分布随时间变化 如图所示，由于垂直中心轴单元的温度分布对水平中心轴(z=10mm)呈现对称特征，则图4-25b给出的是垂直中心轴上半部分(z>10mm)节点(z=11,13,15,17,19mm)处温度随时间的变化。 250200temperature(C)150100500050100150200time(sec.)2503003501113151719 图4-25b垂直中心轴上各单元温度随时间变化 图4-25c是平板水平中心轴的温度影响。由于水平中心轴单元的温度分布对垂直中 心轴(x=500mm)呈现对称特征，则图4-25c给出的是水平中心轴右半部分(x>500mm)节点(x=1,624,707,791,874,957mm)温度随时间的变化。如图所示，从x=500mm到 x=900mm五个单元始终几乎相同，只有从x=900mm到x=1000mm，即在平板边缘 附近，温度呈现稍许的下降。这说明，在平板加热这个非稳态时间阶段内平板水平方向传热是微弱的。 180160140120100806040200050100150200time(sec.)2503003501624707791874957temperature(C) 图4-25c水平中心轴上各单元温度随时间变化 33 根据图4-25a、4-25b与4-25c提供的曲线，我们可得到这样的结论：在固气多孔材料平板被上、下热板加热后，传热主要发生在平板的厚度方向(垂直方向)，沿宽度方向(水平方向)是极微弱的。正是这样的缘故，在分析平板瞬态加热问题时，也可用沿厚度方向的一维瞬态传热问题来近似。 若在本算例的计算机程序中，设ui,j与wi,j始终为0，则运行程序所得温度场的结果即为固气多孔材料经历单纯热传导过程的温度响应。拿这些结果与图4-25系列相比，可说相差无几。若将图4-25a与b与第2章图2-11及图2-12相比，可见它们间的同一性。总之，就温度场而言，本算例提供的固气多孔材料加热过程十分近似单纯热传导过程。就其原因，一是固气多孔材料中空气相的热容远小于固相热容;二是气相的流动速度很小。 (2)关于空气压力pa,i,j 图4-26系列表示固气多孔材料受热板加热条件下材料内空气压力变化的情形。图4-26a给出的是垂直中心轴上半部分(z>10mm)单元空气压力随时间的变化。由图可见，在本算例条件下，垂直中心轴上各单元的压力响应似乎相同。然而，当仔细阅读pa,i,j的数据，仍可发现，在垂直中心轴上，平板表面与平板中心还有微小的压力差别。在加热过程开始的短期内，表面处压力稍高于中心处，而后又转为中心处压力稍高于表面处。这一特征与图4-25b所示的温度响应表面始终高于中心很不一样。 1.61.51.41.31.21.11050100150200time(sec.)2503003501113151719nn n n air pressure(100kPa) 图4-26a垂直中心轴上各单元空气压力随时间变化 图4-26b与图4-26c是水平中心轴上单元的压力响应。图4-26b为平板从被加热后若干时刻(τ=0、5、15、50、100、150、200、250及300sec.)水平中心轴压力分布。如图所示，由于水平中心轴单元的压力分布对垂直中心轴(x=500mm)呈现对称特征，则图4-26c给出的是右半水平中心轴上(x>500mm)各节点(x=1,624,707,791,874,957mm)压力随时间的变化。由图可见，在本算例条件下，水平中心轴上各单元的压力不尽相同，从中心到边缘由高而低。这一特征也与温度响应中沿水平方向几乎相同不一样。 1.61.5air pressure(100kpa)1.41.31.21.110.90200400600distance(mm)800100005 15 50 100 150 200 250 300 图4-26b 水平中心轴压力分布随时间变化 34 1.61.51.41.31.21.11050100150200time(sec.)2503003501624707791874957air pressure(100kPa) 图4-26c水平中心轴上各单元温度随时间变化 纵观本算例中的温度响应与压力响应，前者在垂直中心轴上呈现分布而在水平中心轴上几乎相同，只在边缘处稍有变化。而后者在垂直中心轴上几乎相同，在水平中心轴上呈现分布。两者的特征很不相同。 为了更好地认识本算例压力响应的特点，不妨将它与固气多孔材料单纯热传导条件下压力响应进行对比。图4-26d与图4-26e分别给出了单纯热传导条件下平板垂直中心轴与水平中心轴的压力随时间的变化。 1.71.6air pressure(100kPa)1.51.41.31.21.11050100150200time(sec.)2503003501113151719 图4-26d垂直中心轴上各单元空气压力随时间变化 1.61.51.41.31.21.110501001624707791874957air pressure(100kPa) 图4-26e水平中心轴上各单元空气压力随时间变化 一般而言，气体的压力与气体温度与密度有关，在单纯热传导条件下，气相速度为0，气体密度不发生变化，使气体的压力唯一取决于温度。这样，图4-26d与图4-26e中曲线走向与图4-25b与图4-25c完全相似。在气体流动的条件下，尽管是微弱的渗流，压力响应与温度响应就很不一样，气体的流动造成密度的变化。因此，要认识压力响应的行为还必须联系到气体密度的行为。 (3)关于气体密度ρa,i,j 图4-27a与图4-27b表示固气多孔材料受热板加热条件下材料内空气密度变化的情形。图4-27a是平板上半部分垂直中心轴上各节点(x=11,13,15,17,19mm)单元密度的变化过程。由图可见，在加热过程开始前，平板各单元气体密度相同，而当平板上、下表面被热板加热后，气体密度发生变化。表面处的密度先降而后转为升，而中心处的密度先升而后降。 35 n 150200time(sec.)250300350气体密度这种变化规律还需与气体流动的方向以及当地温度变化联系起来考虑。由于加热开始后表面处温度较中心处高，使表面处的压力较中心处高，引起气体由表面向中心流动，从而造成表面处气体密度下降而中心处上升，同时把表面与中心处的压力差别减小。这样，在加热过程前期，表现为表面处温度高与气体密度低而中心处温度低与气体密度高，两者的压力相差不大。 随着加热过程的延续，表面处的气体密度在下降一段时间后又出现回升。这是因为，中心处的温度逐渐升高再加之中心处气体密度升高，造成中心处气体压力升高并超过表面处的压力，造成气体由中心向表面流动，使中心处的气体密度由升高而下降，表面处密度由下降而上升。 1.41.3air density(kg/m^3)1.21.110.90.8050100150200time(sec.)2503003501113151719 图4-27a垂直中心轴上各单元空气密度随时间变化 再看图4-27b水平中心轴上各单元空气密度随时间变化。这里有两个明显的特征。其一是，各单元的气体密度都是先升而后转为下降。密度上升的原因是当地的压力上升。对照图4-26c压力响应曲线，可见两者有相当的类同。但压力响应在这个加热期间几乎一直在上升，只有靠*板侧面的单元出现缓慢的下降趋势，而密度则在较早时就出现由上升而转会下降。这是由于气体向平板两侧边缘流动的缘故。第二个特征是，平板边缘处的气体密度总是低于中心处，这同样可以从图4-26b压力分布的特点得到解释。 1.41.3air density(kg/m^3)1.21.110.90.8050100150200time(sec.)2503003501624707791874957 图4-27b水平中心轴上各单元空气密度随时间变化 (4)关于水平速度ui,j 图4-28a表示水平中心轴上各节点(x=500,583,667,750,833,917,1000mm)处 空气水平速度随时间变化。由图可见，从平板中心处(x=500mm)到平板边缘(x=1000mm)，水平速度不断加快，也就是，气体在水平方向有一个加速过程。 n 36 0.12air velocity(mm/sec.)0.10.080.060.040.020050100150200time(sec.)250300350500 583 667 750 833 917 1000 图4-28a右半水平中心轴上各单元空气水平速度随时间变化 图4-28b是平板右边界上各单元空气水平速度随时间变化。这里有两个明显的特 征。一是，速度都是由0迅速升高到最大值后而缓解减小。二是，中心处的速度总是大于表面处，最终趋于相同。 0.12air velocity(mm/sec.)0.10.080.060.040.020050100150200time(sec.)2503003501113151719 图4-28b平板右边界上各单元空气水平速度随时间变化 (5)关于速度wi,j 图4-29a给出的是下半块平板垂直中心轴上各节点(x=2,4,6,8mm)处空气垂直速度 n 随时间变化。这里也有两个明显的特点。其一，在加热过程短暂的前期(τ≤10sec)，速度大于0，为正，方向向上。而经过这个前期后，速度由正转为负，方向向下，且逐渐趋于0。另一个特点是，就绝对值而言，表面处的速度总是大于中心处。 0.030.025air velocity(mm/sec.)0.020.0150.010.0050-0.0050501002468 图4-29a下半垂直中心轴上各单元空气垂直速度随时间变化 图4-29b中的曲线正好与图4-29a曲线对水平中心轴呈对称分布。 150200time(sec.)250300350 37 0.005air velocity(mm/sec.)0-0.005-0.01-0.015-0.02-0.025-0.0305010012141618 图4-29b上半垂直中心轴上各单元空气垂直速度随时间变化 150200time(sec.)250300350关于图4-29a与图4-29b中曲线走向的物理解释，我们在关于气体压力行为的解释中已有一定的说明，即，在加热前期，表面空气压力高于中心，而后中心又高于表面，造成速度在加热前期为正而后又转为负。 图4-29c是图4-24中对应于j=3以点划线表示的单元上速度随时间变化的特征。由图可见，在水平方向不同位置上，速度的行为几乎相同。它们随时间的变化都是先正而后转为负，而后逐渐趋于0 0.030.025air velocity(mm/sec.)0.020.0150.010.0050-0.0050501001624707791874957 图4-29c右半水平轴上各单元空气垂直速度随时间变化 (6)关于温度Ti,j与压力pa,i,j的稳态趋向 前面(1)至(5)讨论的是5个因变量在加热过程的瞬态阶段的变化。这种瞬态过程在经历一段长时间后就趋于稳态。图4-30与图4-31显示温度Ti,j与压力pa,i,j由瞬态趋于稳态的过程。每个图中四条曲线对应于图4-24中四个标有a,b,c与d的代表性单元。由图4-30可见，温度响应在瞬态过程表现为表面(b与d)总是高于中心(a与c)，而趋于稳态后表现为a、b与c都相同，仅仅边缘中心处(c)的温度稍底。 250200temperature(C)1501005000200040006000time(sec.)800010000abcd150200time(sec.)250300350nn nn 图4-30温度响应的瞬态与稳态 38 由图4-31可见，压力响应无论是瞬态或是稳态过程，总是表现为平板中部(a与b)总是高于平板边缘(c与d)，而a与b总是几乎相同，c与d也总是几乎相同，。 1.61.51.41.31.21.11020004000time(sec)6000800010000air pressure(100kPa)abcd 图4-31气体压力响应的瞬态与稳态 关于密度与速度的变化这里没给出图示，但可以根据温度与压力响应曲线，联想它们的行为。对速度w，从瞬态到稳态过程中保持为0。气体密度的变化与压力变化趋势相似。速度的变化也随压力趋于稳态而趋于0。 4.6 固气多孔材料空气流加热 本节讨论固气多孔材料平板在空气流加热条件下其内部温度、空气压力、密度与速度的响应。与4.5节相比，只是改变了边界条件，由热板加热改为热空气流加热。如前讨论所知，对于薄板而言，其传热传质过程主要发生在沿板的厚度方向。因此，这里假定为一维传热传质问题，如图4-32所示，平板的两侧与底部都绝热。这样，热流从平板上表面的空气垂直向下传递，而固气多孔材料内的空气由平板内部垂直向上传递。 图4-32一维固气多孔材料空气流加热示意图 4.6.1控制方程及定解条件 本问题涉及的因变量是，空气速度w、空气压力pa、空气密度ρa及温度T。它们相应的控制方程如下： 关于空气速度w，即式(4.4.10)，也即 Taw=− K∂pa (4.4.10) μa∂z 关于空气压力pa的式(4.2.13)，也即 Pa= ρaRT Ma (4.2.13) 关于空气密度ρa的质量守恆方程为 39 ∂ρw∂φρa =−a (4.6.1) ∂z∂τ以及关于温度T的能量守恆方程为 ∂2T∂ρacawT∂ρcT (4.6.2) =k2− ∂z∂z∂τ其中导热系数k与热容ρc服从式(4.2.10)与式(4.2.11)。 式(4.4.10)、式(4.2.13)、(4.6.1)及(4.6.2)构成本问题的封闭的控制方程组。相应的定解条件为： 1.几何条件：本问题坐标原点设在材料平板的底部，故求解的几何域限定在从z=0到z=H。这里是一块一维平板，不计各物理量沿平板宽度x与长度y方向的变化。 2.初始条件：空气速度w的初始值为0，空气压力pa的初始值pa=pa0，它与平板周围环境的压力pe相同，即pa=pa0=pe。平板材料温度T的初始值T0，T=T0，与平板环境的温度Te不相同。空气密度ρa的初始值可由pa0与T0按气体状态方程(4.2.13)求得，ρa= ρa0。 3.边界条件：对于空气速度w，在z=0处，w=0。对于空气压力pa无需单独提供边界条件，环境空气压力pe供计算平板上表面处空气速度w用。对于空气密度ρa相应的边界条件为环境空气密度ρa0，它是对应于环境温度Te及压力pe的密度。对于温度T，其边界条件为 ∂T =0 (4.6.2a) ∂z ∂T z=H h(Te−T)=k (4.6.2b) ∂z 这里，环境温度Te，材料平板表面与环境空气之间的受迫对流传热系数h，都是常数。 4.物理条件：渗透率K与导热系数k为常数，空气比热容ca与粘性系数μa均为温度 z=0 的函数(见4.5.3小节)。 综上所述，本问题涉及4个因变量，对应有4个控制方程，由它们构成的方程组是封闭的。根据相应的定解条件即可求得一组关于w、ρa、pa与T的唯一解，也即多孔平板内确定的空气速度场、密度场、压力场以及温度场。 4.6.2差分方程 对应于上述控制方程及定解条件，可得到相应的差分方程。其离散的求解区域如图4-33所示，平板高度H=JZ⋅Δz。Δz是单元高度，常数。单元数为JZ个，单元对应的节点都位于。对于每个单元可定义压力pa,j、密度ρa,j与温度Tj，j=1,2,...JZ。 用交错网络法来定义沿z方向速度wj，由图4-33小方框所示。wj对应的节点数为 JZ+1。 w1 pa,1wj Δz pa,j−1wj+1 pa,jwJZ+1 pa,JZpa,j+1z z=H ρa,1 ρa,j−1ρa,jTj ρa,j+1Tj+1 ρa,JZTJZ z=0T1Tj−1 图4-33 求解区域离散化趋 40 应该指出，由于传热的缘故，空气粘性系数μa不再是常数，它是温度的函数。为 此，有μa=μa,j，且μa,j为 μa,j= 1 [μa(Tj)+μa(Tj−1)] 2 关于速度w的方程(4.4.10)所相应的差分方程为 w n+1i,j nn Kpa,j−pa,j−1 j=2,3,...,JZ (4.6.3) =− μa,jΔz 和 wi,j=0 j=1 (4.6.3a) 以及 w n+1 i,j n Kpe−pa,j−1 j=JZ+1 (4.6.3b) =− μa,jΔz n+1 关于空气密度ρa的差分方程为 n+1nρa,j=ρa,j+ ΔτΔρaw () j=1,2,...,JZ (4.6.4) φΔz 其中 Δρaw 的差分表示同式(4.5.5a,b,c,d) Δz 关于温度的差分方程也可用控制容积的方法获得。 内单元： +ΔρacawT=+kz ΔzΔz j=2,3,⋅⋅⋅,JZ−1 kz Tj−1−TjTj+1−Tj (ρcTj)n+1−(ρcTj)n ΔτΔz (4.6.5) 上式等号左边前二项，表示该单元在单位时间内通过其边界以热传导方式到达该单元的热量。等号左边第3项，表示该单元在单位时间内通过边界以热对流方式到达该单元的热量。现将式(4.6.5)重新整理，可得 Tj n+1 ={[kz Tj−TjT−TjΔτn +kzj−ΔρacawT]+(ρc)nTi,j}/(ρc)n+1 (4.6.6) ΔzΔzΔz 式中，表示对流项的ΔρacawT其具体表示与上述式(4.5.9e,f,g,h)相同。 上界单元： 上边界单元是材料平板上表面与环境相邻的单元，其与内单元差分方程的差别只是来 自上边的热传导k 用对流传热h(Te−Tj)替代，则差分方程为 ΔzTj−1−TjΔτn+1n Tj={[k+h(Te−Tj)−ΔρacawT]+(ρc)nTi,j}/(ρc)n+1(4.6.6a) ΔzΔz Tj+1−Tj 下界单元： 下边界单元是处于材料平板底部的单元，其与内单元差分方程的差别只是来自下边的 热传导为0，则差分方程为 41 Tj n+1 ={[k Tj+1−TjΔτn −ρa,jca,jwj+1T]+(ρc)nTi,j}/(ρc)n+1 (4.6.6b) ΔzΔz 上述由式(4.6.6)、(4.6.6a)与(4.6.6b)构成温度T的差分方程可统一写成 Tjn+1=(AD⋅Tj−1+AU⋅Tj+1+CA⋅Tj+DA)/EA (4.5.14) (4.6.7) 式(4.6.7)表明，材料平板每个单元在某一时刻的温度都可表示为它二边相邻单元(包括环境) 温度及该单元本身前一时刻温度的线性组合。 对照式(4.6.7)与前述差分方程，可得到不同条件下式(4.6.7)中各项系数AD,AU,CA与EA的表达式，进而得到它们的数值。 空气压力pa的式(4.2.13)相应的差分方程(4.4.8)仍然有效，保持不变。 综上所述，式(4.6.3,a,b)、式(4.6.4)和(4.6.7)即为本小节空气流加热条件下一维固气多孔材料非稳态渗流问题的差分方程组。在初始条件给定后，即可由已知的当前时刻各物理量的值求得在区域(0~H)内下一时刻各物理量的值。也就是，通过差分方程可求得固气多孔材料内空气速度场、密度场、压力场及温度埸。 4.6.3算例及计算结果 本算例的已知条件有：板厚H=2cm，多孔材料的孔隙率φ=0.5，渗透率 K=10−13m2，空气粘性系数μa与比热容ca为温度的函数，μa=17.326+0.0488T，单 位为Pa⋅sec.T为C，ca=1005.3−0.0026T−0.0006T，单位为J/(kg⋅K)。在时间 0 2 τ≥0，温度为Ta=2000C的空气流过固气多孔材料上表面，形成受迫对流加热，受迫对流 传热系数h=15W/(m⋅K)。材料初始温度T0=20C，环境空气的压力pe为大气压， 2 0 pe=100kPa，材料内初始气体压力pa0为环境压力pe。环境与材料内空气密度ρe与ρa0分别可根据Ta与T0以及pe按理想气体状态方程算得。取空间步长Δz=2mm，取时间步长Δτ=0.01sec.。 图4-34至图4-37显示了计算结果。 图4-34显示材料的温度响应，图中标注的1、7、13、19为单元节点距平板底部(x=0)的距离，单位为mm。由图曲线可见，它们与单纯热传导的结果非常接近。这再次显示，在一定条件下，固气多孔材料传热传质过程中，渗流的存在对温度响应的影响甚微。 250200temperature(C)15010050001020304050time(min)60708090171319 图4-34 温度响应 图4-35是压力响应，图中曲线是材料内空气压力pa,j与环境pe之差。这些曲线的特征表现为：1.对时间而言，每条曲线(即板内沿厚度不同地点)所示的压力都是先急速上升而后缓慢下降。上升是当地温度上升所致，而下降为空气在材料内渗流所致。曲线最终趋于 42 0表示材料内空气压力与环境压力相一致。2.对板内的地点而言，平板底部(x=0)的压力总是高于表面处(x=19mm)的压力。它们之间的差值，在接受加热前期被加大而后逐渐缩小并趋于0。 180160140air pressure(Pa)12010080604020001020304050time(min)60708090171319 图4-35 空气压力响应 图4-36为固气多孔材料内的空气密度在材料被加热后的响应。图中曲线都是单调地下降，表示空气向环境渗出。表面处下降得快，底部处下降得慢，但最终都平衡到环境热空气的密度。 1.31.2air pressure(Pa)1.110.90.80.701020304050time(min)60708090171319 图4-36空气密度响应 图4-37为固气多孔材料内的空气在材料被加热后引起流动而形成的速度分布随时间的变化。对照图4-37与图4-35中曲线走向，两者非常相似。这一点很容易理解，因为速度是压力差决定的。需注意的是，对压力而言，从平板底部到表面由大而小，是一个压力缓解的过程，但对速度而言，从平板底部到表面由小而大，是加速的过程。 6air velocity(1e-6*m/sec)321001020304050time(min)60708090281420 图4-37空气速度响应 43 参考文献 [1] Nied,D.A.& Bejan,A. Convection in Porous Media 2nd ed. Springer(2006) [2]章梓雄等，粘性流体力学.北京：清华大学出版社,2005 [3]南京工学院等，粘性流体力学.北京：高等教育出版社，1987 [4]Wooding,R.A. Steady state free thermal convection of liquid in a saturated permeable medium. J.Fluid Mech. 2:273-285,1957 [5]Siau,J.F., Transport Processes in wood. Springer-Verlag,1984 [6]高建民. 木材干燥学, 北京：科学出版社，2008 [7]王补宣，工程传热传质学：下册，北京：科学出版社，1998。 [8]Dullien,F.A. Porous Media:Fluid Transport and Pore Structure. 2nd ed. Academic Press, New York,1992 [9]Dai,C.,Yu,C.and Zhou,X. Heat and mass transfer in wood composite panels during hot pressing:Part2. Modeling void formation and mat permeability. Wood Fiber Sci. 37(2): 242-257 ,2005 [10]Abdallah Zerroug,Khaled Zehar and Larbi Refoufi. Thermal Conductivity Model of Porous Matwrials. Journal of Engineering and Applied Science, 2(4);722-727.2007. [11]Wang,J.,Carson,J.K.,North,M.F.,Cleland,D.J. A new structural model of effective thermal conductivity for heterogeneous materials with co-continuous phases. Int.J.Heat Mass Transfer 51(2008)23-2397. [12]Xue,Q. and Xu,W.M. A model of thermal conductivity of nanofluids with interfacial shells. Materials Chemistry and Physics 90(2005)298-301. [13]凤仪，朱震刚，陶宁，郑海务.闭孔泡沫铝的导热性能.金属学报，第39卷第8期， 817-820.2003. [14]Hadjov,K. Differential medium theory to predict the thermal conductivity of composites finite particle concentration. J.University of chemical technology and metallurgy, 43,4,438-442, 2008 [15]徐旭等.木材导热系数非线性拟合的神经网络模型.浙江大学学报(工学版)，第41(2007) 卷第7期，1201-1204. [16]俞自涛，胡亚才，田甜，范利武.木材横纹有效导热系数的分形模型. 浙江大学学报(工学 版)，第41(2007)卷第2期，351-355. [17]Harada,T.,Hata,T.and Ishihara,S. Thermal constants of wood during the heating process measured with the laser flash method. J.Wood Sci. 44, 425-431,1998 [18]Asako,Y et al. Effective conductivity of compressed woods. Int. J. Heat Mass Transf. 51, 23-2397,2008 [19]Maxwell,J.C.A treatise on electricity and magnetism,third ed., Dover Publications Inc., New York,19(chapter 9) [20]钱壬章，俞昌铭，林文贵，传热分析与计算，北京：高等教育出版社，1987。 [21]帕坦卡,S.V., 传热与流体流动的数值计算，张政译，北京：科学出版社，1984。 [22]陶文铨，数值传热学，西安：西安交通大学出版社，2001 44 第4章程序 4.3第3类边界条件一维导热解析解 C Last change: C 8 Mar 2009 10:04 pm c analytical solution for one dimension conduction with third BC REAL L,K DIMENSION RT(10),SETA2(10),TZ(10),Q(20) OPEN(1,FILE='101.TXT') PI=3.14159 T0=20 TF=200 L=12*0.001 ! thickness of plate m DEN=612 ! plate density kg/m^3 C=1350 ! specific heat J/kgK K=0.112 ! thermal conductivity W/mK H=12 ! ´Óvective heat transfer coefficient W/m^2K FO=0.001 C find roots of transcendental equation N=0 X=0.0001 BI=H*L/(2*K) ! Bi number WRITE(*,*)BI DO I=1,8 ! number of roots 10 Y1=X/BI Y2=1/TAN(X) IF (Y1.GT.Y2) THEN WRITE(*,*)N+1,X RT(I)=X N=N+1 X=3.14159*N+0.0001 ELSE X=X+0.0001 GOTO 10 END IF END DO WRITE(*,100)(RT(I),I=1,8) 100 FORMAT(20F10.4) PAUSE 20 XL=0 DO I=1,6 ! 6 location in half plate SETA=0 DO M=1,8 ! 8 iterm in a series AM=4*SIN(RT(M)) BM=2*RT(M)+SIN(2*RT(M)) XL2=COS(RT(M)*2*XL) F2=EXP(-RT(M)*RT(M)*4*FO) SETA1=AM/BM*XL2*F2 SETA=SETA+SETA1 END DO SETA2(I)=SETA TZ(I)=TF-(TF-T0)*SETA WRITE(*,100)FO,XL,SETA XL=XL+0.1 END DO TIME=FO*L*L*DEN*C/K WRITE(*,100)FO,(SETA2(I),I=1,6) WRITE(1,100)FO,(SETA2(I),I=1,6),TIME,(TZ(I),I=1,6) FO=FO*2 IF(FO.LT.1.1) GOTO 20 END 45 4.4(1)一维渗流 C Last change: C 15 Mar 2009 7:47 pm c onedf.for c it is used to cal. pressue, density and velosity of gas c in 1-D system IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=24,NY=12,NZ=20) DIMENSION PGAS(NX),LGASP(NX),LGASC(NX),UXC(NX) OPEN(1,FILE='OUTP.TXT') OPEN(2,FILE='OUTL.TXT') OPEN(3,FILE='OUTU.TXT') DDG=5.0 IG1=0 RGAS=8.315 ! J/(mole*K) T=273+20 WAIR=25.6 ! g/mole MU=18.6*1E-6 ! viscosity Pa*sec. K=1E-13/MU F=0.5 L=0.5 ! m IX=10 DX=L/REAL(IX) G=0 DG=0.01 DG0=DG DO I=1,IX PGAS(I)=2*1E5 LGASP(I)=PGAS(I)*WAIR/(RGAS*T)*0.001 !kg/m^3 END DO WRITE(*,100)(PGAS(I)*1E-5,I=1,IX) WRITE(*,100)(LGASP(I),I=1,IX) P0=1E5 LGASP0=P0*WAIR/(RGAS*T)*0.001 34 G=G+DG UXC(1)=0 c UXC(1)=K*(P0-PGAS(1))/DX DO I=2,IX UXC(I)=K*(PGAS(I-1)-PGAS(I))/DX END DO UXC(IX+1)=K*(PGAS(IX)-P0)/DX DO I=1,IX AUXPL=LGASP(I)*UXC(I+1) AUXLN=LGASP(I-1)*UXC(I) c IF(i.eq.1) AUXLN=LGASP0*UXC(I) LGASC(I)=LGASP(I)+DG/DX*(AUXLN-AUXPL)/F C WRITE(*,*)I,DG/DX*(AUXLN-AUXPL)/F PGAS(I)=LGASC(I)*RGAS*T/WAIR*1000 LGASP(I)=LGASC(I) END DO IG=INT(G/DDG) IF(IG.EQ.IG1) THEN WRITE(*,100)(UXC(I)*1000,I=1,IX+1) WRITE(1,100)G,(PGAS(I)*1E-5,I=1,IX) WRITE(2,100)G,(LGASC(I),I=1,IX) WRITE(3,100)G,(UXC(I)*1000,I=1,IX+1) IG1=IG1+1 END IF C STOP IF(G.LT.200) GOTO 34 100 FORMAT(200F8.3) 101 FORMAT(F8.3,100E10.2) END 46 4.4(2)二维渗流 C Last change: C 16 Mar 2009 3:50 pm c onedf.for c it is used to cal. pressue, density and velosity of gas c in 1-D system IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=24,NY=12,NZ=20) DIMENSION PG(NX,NZ),LGP(NX,NZ),LGC(NX,NZ),UXC(NX,NZ),VZC(NX,NZ) OPEN(1,FILE='OUTP.TXT') OPEN(2,FILE='OUTL.TXT') OPEN(3,FILE='OUTU.TXT') OPEN(4,FILE='OUTV.TXT') OPEN(7,FILE='OUT.TXT') DDG=5.0 IG1=0 RGAS=8.315 ! J/(mole*K) T=273+20 WAIR=25.6 ! g/mole MU=18.6*1E-6 ! viscosity Pa*sec. KX=1E-13/MU KZ=0.01*KX F=0.5 L=0.5 ! m IX=10 DX=L/REAL(IX) H=0.1 JZ=10 DZ=H/REAL(JZ) G=0 DG=0.01 DG0=DG DO I=1,IX DO J=1,JZ PG(I,J)=2*1E5 LGP(I,J)=PG(I,J)*WAIR/(RGAS*T)*0.001 !kg/m^3 END DO END DO DO J=1,JZ WRITE(*,100)(PG(I,J)*1E-5,I=1,IX) WRITE(*,100)(LGP(I,J),I=1,IX) END DO PAUSE P0=1E5 LGP0=P0*WAIR/(RGAS*T)*0.001 34 G=G+DG DO J=1,JZ UXC(1,J)=0 END DO DO I=1,IX VZC(I,1)=0 END DO c UXC(1)=K*(P0-PG(1))/DX DO I=2,IX UXC(I,1)=KX*(PG(I-1,1)-PG(I,1))/DX END DO DO J=2,JZ VZC(1,J)=KZ*(PG(1,J-1)-PG(1,J))/DZ END DO DO J=2,JZ DO I=2,IX UXC(I,J)=KX*(PG(I-1,J)-PG(I,J))/DX VZC(I,J)=KZ*(PG(I,J-1)-PG(I,J))/DZ END DO END DO 47 DO J=1,JZ UXC(IX+1,J)=KX*(PG(IX,J)-P0)/DX END DO DO I=1,IX VZC(I,JZ+1)=KZ*(PG(I,JZ)-P0)/DZ END DO DO J=1,JZ DO I=1,IX AUXPL=LGP(I,J)*UXC(I+1,J)/DX AUXLN=LGP(I-1,J)*UXC(I,J)/DX AVZPL=LGP(I,J)*VZC(I,J+1)/DZ AVZLN=LGP(I,J-1)*VZC(I,J)/DZ c IF(i.eq.1) AUXLN=LGP0*UXC(I,J) LGC(I,J)=LGP(I,J)+DG*(AUXLN-AUXPL+AVZLN-AVZPL)/F C WRITE(*,*)I,DG/DX*(AUXLN-AUXPL)/F PG(I,J)=LGC(I,J)*RGAS*T/WAIR*1000 LGP(I,J)=LGC(I,J) END DO END DO IG=INT(G/DDG) IF(IG.EQ.IG1) THEN WRITE(*,100)(UXC(I,1)*1000,I=2,IX+1),G WRITE(1,100)G WRITE(2,100)G WRITE(3,100)G WRITE(4,100)G WRITE(4,100)(VZC(I,JZ+1)*1000,I=1,IX) DO J=JZ,1,-1 WRITE(1,100)(PG(I,J)*1E-5,I=1,IX) WRITE(2,100)(LGC(I,J),I=1,IX) WRITE(3,100)(UXC(I,J)*1000,I=1,IX+1) WRITE(4,100)(VZC(I,J)*1000,I=1,IX) END DO WRITE(7,100)G,UXC(2,1)*1000,UXC(IX+1,1)*1000,UXC(2,JZ)*1000, # UXC(IX+1,JZ)*1000, VZC(1,2)*1000,VZC(IX,2)*1000, # VZC(1,JZ+1)*1000,VZC(IX,JZ+1)*1000, # LGC(1,1),LGC(IX,1),LGC(1,JZ),LGC(IX,JZ), # PG(1,1)*1E-5,PG(IX,1)*1E-5,PG(1,JZ)*1E-5,PG(IX,JZ)*1E-5 IG1=IG1+1 END IF C STOP IF(G.LT.200) GOTO 34 100 FORMAT(200F8.3) 101 FORMAT(F8.3,100E10.2) END 48 4.5.热板加热与渗流(二维) C Last change: C 18 Apr 2009 4:36 pm C HOT.2D from OSB082 11 JUN c================================================================== c c c main1 program c c c================================================================== IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=20,NZ=20) COMMON /BAS/G,DG,DX,DZ,IX,JZ,XW,ZH CALL INPUT(GMAX) GMAX=500.0 200 G=G+DG IF(G.GT.GMAX+DG+0.00001) GOTO 300 CALL VELOCITY CALL AIRDP CALL TEM CALL OUTPUT(GMAX) GOTO 200 300 WRITE(*,*)G END C----------------------------------------------------------------- C C s01 input C C inputing initial parameters C----------------------------------------------------------------- SUBROUTINE INPUT(GMAX) IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=20,NZ=20) COMMON /BAS/G,DG,DX,DZ,IX,JZ,XW,ZH COMMON /INP/TINI,TPTNT,TPTNB,TAM,HBCX,KPX,KPZ,PRST,KPTH,KPTV COMMON /TEMP/TP(NX,NZ),TC(NX,NZ) COMMON /PAIR/PAIRC(NX,NZ),P0AIR COMMON /LOAIR/LOAIRP(NX,NZ),LOAIRC(NX,NZ) COMMON /VELO/UX(NX,NZ),WZ(NX,NZ) COMMON/CONST/WWATER,WAIR,RGAS,DG00,FF1,AMV,BMV,CMV FCPA(X)= 0.0006*X*X - 0.026*X + 1005.3 ! specific heat of air J/kgK FVIS(X)= 0.0448*X + 17.326 ! viscosity Pa*sec XW=1.0 !0.4 !m THIP=0.02 TPTNT=200 TPTNB=200 TINI=20 IX=12 JZ=10 DG=0.03 ! 0.01 C mole mass of air (wair)28.97 g/mole C RGAS IS GAS CONSTANT 8.315 j/(mol*k) WAIR=28.97 RGAS=8.315 G=0.0 DX=XW/IX DZ=THIP/REAL(JZ) C POAIR P0AIR=100*1000 ! ambient pressure Pa (N/m^2) LOAIR= P0AIR*WAIR/(1000*RGAS*(TINI+273.15)) DO 41 I=1,IX DO 39 J=1,JZ TP(I,J)=TINI TC(I,J)=TINI LOAIRP(I,J)=LOAIR LOAIRC(I,J)=LOAIR 49 PAIRC(I,J)=P0AIR 39 CONTINUE 41 CONTINUE DO I=1,IX+1 DO J=1,JZ+1 UX(I,J)=0.0 WZ(I,J)=0.0 END DO END DO HBCX=10.0 ! convective heat transfer coefficient TAM=TINI ! ambient temperature PRST=0.6 KSX=0.5 ! W/m*k solid conductivity KAIR=0.026 ! air conductivity KPX=(1-PRST)*KSX+PRST*KAIR KPZ=KPX PERM=1E-13 ! permeability m^2 MUGAS=FVIS(TINI)*1E-6 !Pa*sec air viscosity KPTH=PERM/MUGAS KPTV=KPTH !*0.01 CALL OUTPUT(GMAX) RETURN END c----------------------------------------------------------------- c c s02 velociity c c calculating velocity (m/sec) of vapour flow acoording to c Darcy's law and permeability m^2/(Pa*sec) c----------------------------------------------------------------- SUBROUTINE VELOCITY IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=20,NZ=20) COMMON /BAS/G,DG,DX,DZ,IX,JZ,XW,ZH COMMON /INP/TINI,TPTNT,TPTNB,TAM,HBCX,KPX,KPZ,PRST,KPTH,KPTV COMMON /PAIR/PAIRC(NX,NZ),P0AIR COMMON /LOAIR/ LOAIRP(NX,NZ),LOAIRC(NX,NZ) COMMON /TEMP/TP(NX,NZ),TC(NX,NZ) COMMON /VELO/UX(NX,NZ),WZ(NX,NZ) FVIS(X)= 0.0448*X + 17.326 DO 52 I=1,IX+1 DO 50 J=1,JZ+1 IF(I.GT.1.AND.I.LT.IX+1) THEN TCL=TC(I-1,J) TCC=TC(I,J) VISCX=0.5*(FVIS(TCL)+FVIS(TCC)) UX(I,J)=KPTH/VISCX*(PAIRC(I-1,J)-PAIRC(I,J))/DX ELSE IF(I.EQ.1) THEN TCC=TC(I,J) VISCX=FVIS(TCC) UX(I,J)=2.0*KPTH/VISCX*(P0AIR-PAIRC(I,J))/DX END IF IF(I.EQ.IX+1) THEN TCL=TC(I-1,J) VISCX=FVIS(TCL) UX(I,J)=2.0*KPTH/VISCX*(PAIRC(IX,J)-P0AIR)/DX END IF END IF IF((J.EQ.JZ+1).OR.(J.EQ.1)) THEN WZ(I,J)=0.0 ELSE TCD=TC(I,J-1) TCC=TC(I,J) VISCY=0.5*(FVIS(TCD)+FVIS(TCC)) WZ(I,J)=KPTV/VISCY*(PAIRC(I,J-1)-PAIRC(I,J))/DZ IF(ABS(WZ(I,J)).LT.1E-9) WZ(I,J)=0.0 50 END IF C UX(I,J)=0 C WZ(I,J)=0 50 CONTINUE 52 CONTINUE 100 FORMAT(20F8.4) RETURN END c-----------------------------------------------------------------c c s03 AIRdp c c calculating pressure (Pa; N/m^2) and density (kg/m^3) of air based on c continuty equation and state equation c----------------------------------------------------------------- SUBROUTINE AIRDP IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=20,NZ=20) COMMON /BAS/G,DG,DX,DZ,IX,JZ,XW,ZH COMMON /INP/TINI,TPTNT,TPTNB,TAM,HBCX,KPX,KPZ,PRST,KPTH,KPTV COMMON /TEMP/TP(NX,NZ),TC(NX,NZ) COMMON /PAIR/PAIRC(NX,NZ),P0AIR COMMON /LOAIR/ LOAIRP(NX,NZ),LOAIRC(NX,NZ) COMMON /VELO/UX(NX,NZ),WZ(NX,NZ) COMMON/CONST/WWATER,WAIR,RGAS,DG00,FF1,AMV,BMV,CMV DO 39 I=1,IX DO 37 J=1,JZ UXLLR=LOAIRP(I-1,J)*UX(I,J)/DX !left side from L to R UXLRL=LOAIRP(I,J)*UX(I,J)/DX !left side from R to L UXRLR=LOAIRP(I,J)*UX(I+1,J)/DX UXRRL=LOAIRP(I+1,J)*UX(I+1,J)/DX WZDDU=LOAIRP(I,J-1)*WZ(I,J)/DZ !left side from L to R WZDUD=LOAIRP(I,J)*WZ(I,J)/DZ !left side from R to L WZUDU=LOAIRP(I,J)*WZ(I,J+1)/DZ WZUUD=LOAIRP(I,J+1)*WZ(I,J+1)/DZ IF(UX(I+1,J).GE.0) THEN IF(UX(I,J).GE.0) THEN LUX=UXLLR-UXRLR ELSE LUX=UXLRL-UXRLR END IF ELSE IF(UX(I,J).GE.0) THEN LUX=UXLLR-UXRRL ELSE LUX=UXLRL-UXRRL END IF END IF IF(WZ(I,J+1).GE.0) THEN IF(WZ(I,J).GE.0) THEN LWZ=WZDDU-WZUDU ELSE LWZ=WZDUD-WZUDU END IF ELSE IF(WZ(I,J).GE.0) THEN LWZ=WZDDU-WZUUD ELSE LWZ=WZDUD-WZUUD END IF END IF LOAIRC(I,J)=LOAIRP(I,J)+(LUX+LWZ)*DG/PRST IF(LOAIRC(I,J).LT.0.001) LOAIRC(I,J)=0.001 T71=TC(I,J) PAIRC(I,J)=LOAIRC(I,J)*1000*RGAS*(T71+273.15)/WAIR 37 CONTINUE 39 continue 51 RETURN END c----------------------------------------------------------------- c c s04 temperature c c calculating temperatures in mat according to energy conservation c----------------------------------------------------------------- SUBROUTINE TEM IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=20,NZ=20) COMMON /BAS/G,DG,DX,DZ,IX,JZ,XW,ZH COMMON /INP/TINI,TPTNT,TPTNB,TAM,HBCX,KPX,KPZ,PRST,KPTH,KPTV COMMON /TEMP/TP(NX,NZ),TC(NX,NZ) COMMON /LOAIR/ LOAIRP(NX,NZ),LOAIRC(NX,NZ) COMMON /VELO/UX(NX,NZ),WZ(NX,NZ) FCPA(X)= 0.0006*X*X - 0.026*X + 1005.3 ! specific heat of air J/kgK CPS=13 ! J/(kg*K) solid specific heat LOWS=1500 ! kg/m^3 solid density DO 51 I=1,IX DO 49 J=1,JZ CPAIR=FCPA(TP(I,J)) LOWCP=(1-PRST)*LOWS*CPS+PRST*LOAIRP(I,J)*CPAIR ULLR=CPAIR*LOAIRP(I-1,J)*UX(I,J) URLR=CPAIR*LOAIRP(I,J)*UX(I+1,J) URRL=CPAIR*LOAIRP(I+1,J)*UX(I+1,J) ULRL=CPAIR*LOAIRP(I,J)*UX(I,J) WDDU=CPAIR*LOAIRP(I,J-1)*WZ(I,J) WUDU=CPAIR*LOAIRP(I,J)*WZ(I,J+1) WUUD=CPAIR*LOAIRP(I,J+1)*WZ(I,J+1) WDUD=CPAIR*LOAIRP(I,J)*WZ(I,J) KDX=KPX/(DX*DX) KDZ=KPZ/(DZ*DZ) CA=-2*(KDX+KDZ)*DG+LOWCP DA=0 EA=LOWCP AW1=KDX AE1=KDX AD1=KDZ AU1=KDZ IF(I.EQ.1.OR.I.EQ.IX) THEN IF(I.EQ.1) AW1=0 IF(I.EQ.IX) AE1=0 CA=-(KDX+HBCX/DX+2*KDZ)*DG+LOWCP DA=HBCX*TAM/DX*DG END IF IF(J.EQ.1.OR.J.EQ.JZ) THEN IF(J.EQ.1) THEN AD1=0 DA=2*KDZ*TPTNT*DG END IF IF(J.EQ.JZ) THEN AU1=0 DA=2*KDZ*TPTNT*DG END IF CA=-(2*KDX+3*KPZ/(DZ*DZ))*DG+LOWCP END IF IF(UX(I+1,J).GE.0) THEN IF(UX(I,J).GE.0) THEN AW=(AW1+ULLR/DX)*DG AE=AE1*DG CA=CA-URLR/DX*DG ELSE AW=AW1*DG AE=AE1*DG CA=CA+(ULRL-URLR)/DX*DG END IF 52 ELSE IF(UX(I,J).GE.0) THEN AW=(AW1+ULLR/DX)*DG AE=(AE1-URRL/DX)*DG CA=CA ELSE AW=AW1*DG AE=(AE1-URRL/DX)*DG CA=CA+ULRL/DX*DG END IF END IF IF(WZ(I,J+1).GE.0) THEN IF(WZ(I,J).GE.0) THEN AD=(AD1+WDDU/DZ)*DG AU=AU1*DG CA=CA-WUDU/DZ*DG ELSE AD=AD1*DG AU=AU1*DG CA=CA+(WDUD-WUDU)/DZ*DG END IF ELSE IF(WZ(I,J).GE.0) THEN AD=(AD1+WDDU/DZ)*DG AU=(AU1-WUUD/DZ)*DG CA=CA ELSE AD=AD1*DG AU=(AU1-WUUD/DZ)*DG CA=CA+WDUD/DZ*DG END IF END IF TC(I,J)=(AW*TP(I-1,J)+AE*TP(I+1,J) # +AD*TP(I,J-1)+AU*TP(I,J+1)+CA*TP(I,J)+DA)/EA 49 CONTINUE 51 CONTINUE RETURN END c--------------------------------------------------------------- c c s05 output c c--------------------------------------------------------------- SUBROUTINE OUTPUT(GMAX) IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NX=20,NZ=20) COMMON /BAS/G,DG,DX,DZ,IX,JZ,XW,ZH COMMON /TEMP/TP(NX,NZ),TC(NX,NZ) COMMON /INP/TINI,TPTNT,TPTNB,TAM,HBCX,KPX,KPZ,PRST,KPTH,KPTV COMMON /PAIR/PAIRC(NX,NZ),P0AIR COMMON /LOAIR/ LOAIRP(NX,NZ),LOAIRC(NX,NZ) COMMON /VELO/UX(NX,NZ),WZ(NX,NZ) COMMON /DEM/NDIMENSION IF(G.LT.0.00001) THEN C time internal of output DDG=5.0 C DDG=DG IG1=0 END IF c o-dimensional output data OPEN(13,FILE='OUTPUT\\0D\\TLPUV.TXT') c one-dimensional output data OPEN(21,FILE='OUTPUT\\1D\\2TX.TXT') OPEN(22,FILE='OUTPUT\\1D\\2TZ.TXT') OPEN(23,FILE='OUTPUT\\1D\\2PAIRX.TXT') OPEN(24,FILE='OUTPUT\\1D\\2PARZ.TXT') 53 OPEN(25,FILE='OUTPUT\\1D\\2DENX.TXT') OPEN(26,FILE='OUTPUT\\1D\\2DENZ.TXT') OPEN(31,FILE='OUTPUT\\1D\\2UX.TXT') OPEN(32,FILE='OUTPUT\\1D\\2UZ.TXT') OPEN(33,FILE='OUTPUT\\1D\\2WX.TXT') OPEN(34,FILE='OUTPUT\\1D\\2WZ.TXT') c two-dimensional output data OPEN(51,FILE='OUTPUT\\2D\\3T.TXT') OPEN(52,FILE='OUTPUT\\2D\\3P.TXT') OPEN(53,FILE='OUTPUT\\2D\\3L.TXT') OPEN(,FILE='OUTPUT\\2D\\3U.TXT') OPEN(55,FILE='OUTPUT\\2D\\3V.TXT') IG=INT(G/DDG) IF(IG.EQ.IG1) THEN WRITE(*,1002)G,TC(IX/2,JZ/2),TC(IX,JZ),PAIRC(IX/2,JZ/2)*1E-5 WRITE(13,1002)G,TC(IX/2,JZ/2),TC(IX/2,JZ),TC(IX,JZ/2),TC(IX,JZ) # ,PAIRC(IX/2,JZ/2)*1E-5,PAIRC(IX/2,JZ)*1E-5 # ,PAIRC(IX,JZ/2)*1E-5,PAIRC(IX,JZ)*1E-5 WRITE(21,1002)G,(TC(I,JZ/2+1),I=1,IX) WRITE(22,1002)G,(TC(IX/2+1, J),J1,JZ) WRITE(24,1002)G,(PAIRC(IX/2+1,J)*1E-5,J=1,JZ) WRITE(23,1002)G,(PAIRC(I,JZ/2+1)*1E-5,I=1,IX) WRITE(25,1002)G,(LOAIRC(I,JZ/2+1),I=1,IX) WRITE(26,1002)G,(LOAIRC(IX/2+1,J),J=1,JZ) WRITE(31,1002)G,(UX(I,JZ/2+1)*1000,I=1,IX+1) WRITE(32,1002)G,(UX(IX+1,J)*1000,J=1,JZ) WRITE(33,1002)G,(WZ(I,3)*1000,I=1,IX) WRITE(34,1002)G,(WZ(IX/2+1,J)*1000,J==1,JZ+1) WRITE(51,1002)G WRITE(52,1002)G WRITE(53,1002)G WRITE(,1002)G WRITE(55,1002)G DO 71 J=JZ,1,-1 WRITE(51,1002)(TC(I,J),I=1,IX) WRITE(52,1002)(PAIRC(I,J)*1E-5,I=1,IX) WRITE(53,1002)(LOAIRC(I,J),I=1,IX) 71 CONTINUE DO J=JZ,1,-1 WRITE(,1002)(UX(I,J)*1000,I=1,IX+1) END DO DO J=JZ+1,1,-1 WRITE(55,1002)(WZ(I,J)*1000,I=1,IX) END DO IG1=IG1+1 END IF 1002 FORMAT(50F10.4) if(g.lt.0.00001) goto 300 DO 74 I=1,IX+1 DO 72 J=1,JZ+1 TP(I,J)=TC(I,J) LOAIRP(I,J)=LOAIRC(I,J) 72 CONTINUE 74 CONTINUE 300 RETURN END 4.6空气流一维加热 C Last change: C 18 Apr 2009 4:28 pm C HOT.2D from OSB082 11 JUN c================================================================== c c c main1 program c c c================================================================== IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NZ=20) COMMON /BAS/G,DG,DZ,JZ,ZH CALL INPUT(GMAX) GMAX=500.0 200 G=G+DG IF(G.GT.GMAX+DG+0.00001) GOTO 300 CALL VELOCITY CALL AIRDP CALL TEM CALL OUTPUT(GMAX) GOTO 200 300 WRITE(*,*)G END C----------------------------------------------------------------- C C s01 input C C inputing initial parameters C----------------------------------------------------------------- SUBROUTINE INPUT(GMAX) IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NZ=20) COMMON /BAS/G,DG,DZ,JZ,ZH COMMON /INP/TINI,TAM,HBC,KPZ,PRST,KPTV COMMON /TEMP/TP(NZ),TC(NZ) COMMON /PAIR/PAIRC(NZ),P0AIR COMMON /LOAIR/LOAIRP(NZ),LOAIRC(NZ) COMMON /VELO/WZ(NZ) COMMON/CONST/WWATER,WAIR,RGAS,DG00,FF1,AMV,BMV,CMV FCPA(X)= 0.0006*X*X - 0.026*X + 1005.3 ! specific heat of air J/kgK FVIS(X)= 0.0448*X + 17.326 ! viscosity Pa*sec THIP=0.02 !m TINI=20 JZ=10 DG=0.02 ! 0.01 C mole mass of air (wair)28.97 g/mole C RGAS IS GAS CONSTANT 8.315 j/(mol*k) WAIR=28.97 RGAS=8.315 G=0.0 DZ=THIP/REAL(JZ) C POAIR P0AIR=100*1000 ! ambient pressure Pa (N/m^2) LOAIR= P0AIR*WAIR/(1000*RGAS*(TINI+273.15)) DO 39 J=1,JZ TP(J)=TINI TC(J)=TINI LOAIRP(J)=LOAIR LOAIRC(J)=LOAIR PAIRC(J)=P0AIR 39 CONTINUE DO J=1,JZ+1 WZ(J)=0.0 END DO HBC=15.0 ! convective heat transfer coefficient 55 TAM=200 ! ambient temperature PRST=0.6 KSX=0.5 ! W/m*k solid conductivity KAIR=0.026 ! air conductivity KPX=(1-PRST)*KSX+PRST*KAIR KPZ=KPX PERM=1E-13 ! permeability m^2 MUGAS=FVIS(TINI)*1E-6 !Pa*sec air viscosity KPTH=PERM/MUGAS KPTV=KPTH !*0.01 CALL OUTPUT(GMAX) RETURN END c----------------------------------------------------------------- c c s02 velociity c c calculating velocity (m/sec) of vapour flow acoording to c Darcy's law and permeability m^2/(Pa*sec) c----------------------------------------------------------------- SUBROUTINE VELOCITY IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NZ=20) COMMON /BAS/G,DG,DZ,JZ,ZH COMMON /INP/TINI,TAM,HBC,KPZ,PRST,KPTV COMMON /PAIR/PAIRC(NZ),P0AIR COMMON /LOAIR/ LOAIRP(NZ),LOAIRC(NZ) COMMON /TEMP/TP(NZ),TC(NZ) COMMON /VELO/WZ(NZ) FVIS(X)= 0.0448*X + 17.326 DO 50 J=1,JZ+1 IF((J.EQ.JZ+1).OR.(J.EQ.1)) THEN IF(J.EQ.1) WZ(J)=0.0 IF(J.EQ.JZ+1) THEN TCD=TC(J-1) VISCY=FVIS(TCD) WZ(J)=KPTV/VISCY*(PAIRC(J-1)-P0AIR)/DZ END IF ELSE TCD=TC(J-1) TCC=TC(J) VISCY=0.5*(FVIS(TCD)+FVIS(TCC)) WZ(J)=KPTV/VISCY*(PAIRC(J-1)-PAIRC(J))/DZ IF(ABS(WZ(J)).LT.1E-9) WZ(J)=0.0 END IF c WZ(J)=0 50 CONTINUE 100 FORMAT(20F8.4) RETURN END c-----------------------------------------------------------------c c s03 AIRdp c c calculating pressure (Pa; N/m^2) and density (kg/m^3) of air based on c continuty equation and state equation c----------------------------------------------------------------- SUBROUTINE AIRDP IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NZ=20) COMMON /BAS/G,DG,DZ,JZ,ZH COMMON /INP/TINI,TAM,HBC,KPZ,PRST,KPTV COMMON /TEMP/TP(NZ),TC(NZ) COMMON /PAIR/PAIRC(NZ),P0AIR COMMON /LOAIR/ LOAIRP(NZ),LOAIRC(NZ) COMMON /VELO/WZ(NZ) COMMON/CONST/WWATER,WAIR,RGAS,DG00,FF1,AMV,BMV,CMV 56 DO 37 J=1,JZ WZDDU=LOAIRP(J-1)*WZ(J)/DZ ! down side from D to U WZDUD=LOAIRP(J)*WZ(J)/DZ !up side from D to U WZUDU=LOAIRP(J)*WZ(J+1)/DZ WZUUD=LOAIRP(J+1)*WZ(J+1)/DZ IF(WZ(J+1).GE.0) THEN IF(WZ(J).GE.0) THEN LWZ=WZDDU-WZUDU ELSE LWZ=WZDUD-WZUDU END IF ELSE IF(WZ(J).GE.0) THEN LWZ=WZDDU-WZUUD ELSE LWZ=WZDUD-WZUUD END IF END IF LOAIRC(J)=LOAIRP(J)+LWZ*DG/PRST IF(LOAIRC(J).LT.0.001) LOAIRC(J)=0.001 T71=TC(J) PAIRC(J)=LOAIRC(J)*1000*RGAS*(T71+273.15)/WAIR 37 CONTINUE RETURN END c----------------------------------------------------------------- c c s04 temperature c c calculating temperatures in mat according to energy conservation c----------------------------------------------------------------- SUBROUTINE TEM IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NZ=20) COMMON /BAS/G,DG,DZ,JZ,ZH COMMON /INP/TINI,TAM,HBC,KPZ,PRST,KPTV COMMON /TEMP/TP(NZ),TC(NZ) COMMON /LOAIR/ LOAIRP(NZ),LOAIRC(NZ) COMMON /VELO/WZ(NZ) FCPA(X)= 0.0006*X*X - 0.026*X + 1005.3 ! specific heat of air J/kgK CPS=13 ! J/(kg*K) solid specific heat LOWS=1500 ! kg/m^3 solid density DO 49 J=1,JZ CPAIR=FCPA(TP(J)) LOWCP=(1-PRST)*LOWS*CPS+PRST*LOAIRP(J)*CPAIR WDDU=CPAIR*LOAIRP(J-1)*WZ(J) WUDU=CPAIR*LOAIRP(J)*WZ(J+1) WUUD=CPAIR*LOAIRP(J+1)*WZ(J+1) WDUD=CPAIR*LOAIRP(J)*WZ(J) KDZ=KPZ/(DZ*DZ) CA=-2*(KDZ)*DG+LOWCP DA=0 EA=LOWCP AD1=KDZ AU1=KDZ IF(J.EQ.1.OR.J.EQ.JZ) THEN IF(J.EQ.1) THEN AD1=0 CA=-(KDZ)*DG+LOWCP END IF IF(J.EQ.JZ) THEN AU1=0 CA=-(KDZ+HBC/DZ)*DG+LOWCP DA=(HBC*TAM/DZ)*DG 57 END IF END IF IF(WZ(J+1).GE.0) THEN IF(WZ(J).GE.0) THEN AD=(AD1+WDDU/DZ)*DG AU=AU1*DG CA=CA-WUDU/DZ*DG ELSE AD=AD1*DG AU=AU1*DG CA=CA+(WDUD-WUDU)/DZ*DG END IF ELSE IF(WZ(J).GE.0) THEN AD=(AD1+WDDU/DZ)*DG AU=(AU1-WUUD/DZ)*DG CA=CA ELSE AD=AD1*DG AU=(AU1-WUUD/DZ)*DG CA=CA+WDUD/DZ*DG END IF END IF TC(J)=(AD*TP(J-1)+AU*TP(J+1)+CA*TP(J)+DA)/EA 49 CONTINUE RETURN END c--------------------------------------------------------------- c c s05 output c c--------------------------------------------------------------- SUBROUTINE OUTPUT(GMAX) IMPLICIT REAL*8(A-H,K,M,L,O-Z) PARAMETER(NZ=20) COMMON /BAS/G,DG,DZ,JZ,ZH COMMON /TEMP/TP(NZ),TC(NZ) COMMON /INP/TINI,TAM,HBC,KPZ,PRST,KPTV COMMON /PAIR/PAIRC(NZ),P0AIR COMMON /LOAIR/ LOAIRP(NZ),LOAIRC(NZ) COMMON /VELO/WZ(NZ) IF(G.LT.0.00001) THEN C time internal of output DDG=5.0 C DDG=DG IG1=0 END IF OPEN(22,FILE='TZ.TXT') OPEN(24,FILE='PARZ.TXT') OPEN(26,FILE='DENZ.TXT') OPEN(34,FILE='WZ.TXT') IG=INT(G/DDG) IF(IG.EQ.IG1) THEN G1=G/60 WRITE(*,1002)G1,TC(1),TC(JZ),PAIRC(1)*1E-5 ,PAIRC(JZ)*1E-5 WRITE(22,1002)G1,(TC(J),J=1,JZ) WRITE(24,1002)G1,((PAIRC(J)-P0AIR),J=1,JZ) WRITE(26,1002)G1,(LOAIRC(J),J=1,JZ) WRITE(34,1002)G1,(WZ(J)*1E+6,J=1,JZ+1) IG1=IG1+1 58 END IF 1002 FORMAT(50F10.4) if(g.lt.0.00001) goto 300 DO 72 J=1,JZ+1 TP(J)=TC(J) LOAIRP(J)=LOAIRC(J) 72 CONTINUE 300 RETURN END 59 因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- yrrf.cn 版权所有 赣ICP备2024042794号-2
违法及侵权请联系：TEL:199 1889 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务