《武汉工程大学学报》 2012年9期
54-58
出版日期：2012-10-10
ISSN:1674-2869
CN:42-1779/TQ
针束型回热器熵产率的分析
0引言热声发动机和制冷机没有或者很少有运动部件，并且使用环境友好气体作为工质，以上突出优点使其在近年来备受关注\[1\].回热器是热声热机的关键部件\[2-3\]，实际热声装置中使用了多种结构形式的回热器包括平行板结构、圆孔结构和钢丝型结构等\[4-5\].Swift等\[6-7\]引入了一种被称为针束结构的新型回热器：针束型结构如图1所示，在回热器轴向位置平行布置多根细针，单根针半径为rd，每根针处于一个正六边形流道中心.理论与实验均已证明针束型回热器结构优于其它规则几何结构的回热器\[8\].图1针束结构回热器横截面局部示意图
Fig.1A schematic diagram of the pinarray stack工质在回热器内振荡时的换热和黏性损耗会导致不可逆损失，这种不可逆损失可以通过熵产率\[9-10\]来分析.文献\[11-16\]对热声装置不可逆性进行了分析，但这些研究仅仅集中于单板、平行板或钢丝型回热器以及换热器.本文对针束回热器中由换热和黏性耗散导致的熵产率进行了分析，研究了间距、尺寸、角频率、温度梯度和阻抗比等参数对熵产率的影响.1理论分析根据针束型结构特点利用柱面坐标系(x,r,θ)进行分析，其中坐标x表示声波传播方向（针束沿此方向排列）.为便于计算，如图1所示：Swift等\[6-7\]将每根针对应的正六边形流道等效为一个面积相等的半径为r0的圆形流道，满足：r0=23/πy0≈1.05y0，其中y0为相邻针束1/2中心距.线性化动量守恒方程为\[4\] iωρ0u1=-dp1dx+μ2u1(1)
中i=-1，p1为工质的复压力振幅，ω为角频率，ρ0为工质平均密度，u1为工质的复速度振幅，μ为工质动力黏度.在无滑边界条件下，有：u1(rd)=0(2)在针的表面，对称边界条件下：⊥u1(hexagon)=u1(r0)r=0(3)
其中⊥表示在图1所示单元上与该六边形单元边界垂直方向求导.连续性方程为：ρ1t+ρ0u1x=0(4)
其中cp为工质的复密度振幅.线性化能量方程为：
ρ0cp(iωT1+u1dTdx)-iωT0βp1=k\[1rr(rT1r)\](5)
其中cp为比定压热容，T1为复温度振幅，dT/dx为温度梯度，T0为平均温度，β为等温膨胀系数，k为工质导热系数.线性化状态方程为：p1=ρ1a2+ρ0βT0(1-γ)T1(6)
其中a为工质绝热声速，γ为比热比.由公式（1）-（6）可得：
T1=T0βρ0cpp1(x)(1-fk)+iu1(x)ω×
dT0dx(1+σ1-σfμ-11-σfk)(7)
u1=u1(x)\[1-Y1(z0μ)J0(zμ)-J1(z0μ)Y0(zμ)Y1(z0μ)J0(zdμ)-J1(z0μ)Y0(zdμ)\]=
-1iωρ0dp1dx\[1-Y1(z0μ)J0(zμ)-J1(z0μ)Y0(zμ)Y1(z0μ)J0(zdμ)-J1(z0μ)Y0(zdμ)\](8)
其中
fk=Y1(z0k)J0(zk)-J1(z0k)Y0(zk)Y1(z0k)J0(zdk)-J1(z0k)Y0(zdk)(9)
fμ=Y1(z0μ)J0(zμ)-J1(z0μ)Y0(zμ)Y1(z0μ)J0(zdμ)-J1(z0μ)Y0(zdμ)(10)
其中J0和J1分别表示零阶和一阶Bessel函数，Y0和Y1分别表示零阶和一阶Neumann函数，zμ=r·-iωρ0/μ=(i-1)r/δμ,z0μ=(i-1)r0/δμ,zdμ=(i-1)rd/δμ,zk=(i-1)r/δk,z0k=(i-1)rk/δk,zdk=(i-1)rd/δk,δμ为工质的黏性渗透深度,δk为工质的热渗透深度.单位质量工质能量方程为
dEdt+pdVdt=V-V·(11)
其中E为比内能，V为比容，为黏性耗散函数，为热流向量.由傅里叶导热定律
=-k⊥T(12)
可得：
·=T··(-k⊥TT)=
-\[T·（k⊥TT）+kT(⊥T)2\](13)第9期杨志春，等：针束型回热器熵产率的分析
武汉工程大学学报第34卷
将公式（11）代入（5）得：
dEdt+pdVdt-VT·（k⊥TT）=V+VkT(⊥T)2(14)将Gibbs 方程Tds=dE+pdV代入公式（14）并忽略一阶小量得：
ρ0dsdt-·(k⊥T1T0)=T0+k(⊥T1T0)2(15)
其中s为单位体积熵产.公式（15）左边第一项为单位体积熵产率，第二项为熵流.公式（15）右边第一项为黏性不可逆性引起的熵产率，第二项为传热引起的熵产率.由此可得单位体积熵产率：Sit=T0+k(⊥TT0)2(16)在层流和小扰动条件下，黏性耗散函数为：
=μ(⊥u1)2(17)将公式（17）代入公式（16）并取时间平均得：
Sit=<Sit>t=k2T20(⊥T·⊥T*)+
μ2T0(⊥u·⊥u*)(18)
其中*表示复共轭.将公式（7）和公式（8）代入公式（18）并取实部可得针束型回热器时均熵产率：Sit=12T20\[T20β2a2u21Z2ωρ0cpRe(M11)+
(dT0dx)2u21ρ0cpω(1-σ)2Re(M21)+
u21T0βZρ0aσ1-σdT0dxIm(M31)\]+
ωρ02T0·u21·Fμ·F*μ(19)
其中
Re(M11)=AjkJ1(zk)A*jkJ*1(zk)+
AykY1(zk)A*ykY*1(zk)-
2Re\[AjkJ1(zk)A*ykY*1(zk)\]
Re(M21)=AjkJ1(zk)A*jkJ*1(zk)+
AykY1(zk)A*ykY*1(zk)-
2Re\[AjkJ1(zk)A*ykY*1(zk)\]+
σ{AjμJ1(zμ)A*jμJ*1(zμ)+
AyμY1(zμ)A*yμY*1(zμ)-
2Re\[AyμY1(zμ)A*jμJ*1(zμ)\]}-
2σRe\[AjμJ1(zμ)A*jkJ*1(zk)+
AyμY1(zμ)A*ykY*1(zk)-
AyμY1(zμ)A*jkJ*1(zk)-
AjμJ1(zμ)A*ykY*1(zk)\]
Im(M31)=2Im\[AykY1(zk)A*jμJ*1(zμ)+
AjkJ1(zk)A*yμY*1(zμ)+
AjμJ1(zμ)A*jkJ*1(zk)+
AyμY1(zμ)A*ykY*1(zk)\]
Fμ=Y1(z0μ)J1(zμ)-J1(z0μ)Y1(zμ)Y1(z0μ)J0(zdμ)-J1(z0μ)Y0(zdμ)
Ajk=Y1(z0k)Y1(z0k)J0(zdk)-J1(z0k)Y0(zdk)
Ayk=J1(z0k)Y1(z0k)J0(zdk)-J1(z0k)Y0(zdk)
Ajμ=Y1(z0μ)Y1(z0μ)J0(zdμ)-J1(z0μ)Y0(zdμ)
Ayμ=J1(z0μ)Y1(z0μ)J0(zdμ)-J1(z0μ)Y0(zdμ)
其中σ为工质Prandtl数，阻抗比Z=p1(x)ρ0au1(x).由公式（19）可知Sir与r有关，即Sir沿流道横向截面有一分布，取其截面平均得：
Si=(Sir)r=1π(r20-r2d)∫r0rd2πrSirdr=
12T20\[T20β2C20μ21Z2ωρ0cpF1+(dT0dx)2u21ρ0cpω(1-σ)2F2-u21T0βZρ0C0σ1-σdT0dxF3\]+ωρ02T0u21F4(20)
其中
F1=1π(r20-r2d)∫r0rd2πrRe(M11)dr=
2r20-r20\[AjkA*jkM+J1KJ*1K+AykA*ykMY1kY*1k-
Re(AjkA*ykMJ1kY*1k)\]
F2=1π(r20-r2d)∫r0rd2πrRe(M21)dr=
2r20-r2d{AjkA*jkMJ1kJ*1k+AykA*ykMY1kY*1k-
2Re(AjkA*ykMJ1kY*1k)+σ\[AjμA*jμMJ1μJ*1μ+
AyμA*yμMY1μY*1μ-2Re(AyμA*jμMY1μJ*1μ)\]
-2σRe(AjμA*jkMJ1μJ*1k+AyμA*ykMY1μY*1k-
AyμA*jkMY1μJ*1k-AjμA*ykMJ1μY*1k)}
F3=1π(r20-r2d)∫r0rd2πrRe(M31)dr=
4r20-r2dIm(AykA*jμMY1kJ*1μ+AjkA*yμMJ1kY*1μ+
AjμA*jkMJ1μJ*1k+AyμA*jkMY1μY*1k)
F4=1π(r20-r2d)∫r0rd2πrFμ·F*μdr=
2r20-r2d(AjμA*jμMJ1μJ*1μ-AjμA*yμMJ1μY*1μ-
AyμA*jμMY1μJ*1μ+AyμA*yμMY1μY*1μ)
MN1ηR1ζ=∫r0rdrN1ηR*1ζdr=
BηrdN0(Bη·rd)R1(B*ζ·rd)B2η-B*2ζ -
B*ζ·rdN1(Bη·rd)R0(B*ζ·rd)B2η-B*2ζ
(N,R=J,Y; η,ζ=μ,k)
Bk=i-1δk
Bμ=i-1δμ
2数值结果分析根据公式（19）和（20），可分析得出不同参数对熵产率的影响：2.1截面分布熵产率公式（19）表明截面分布不可逆熵产率Sir是针束距中心点距离r的函数.计算时所取参数T0=600 K,T0β=1,f=ω/2π=400 Hz,dT/dx=1 000 K/m,Z=1,rd=0.000 1 m,ρ0=1.9 kg/m3,σ=0.68和k=0.13 W/(m2·K).由图2可知：在不同的r0/rd下Sir随着r的增大而减小.这是因为离针束越远，换热与黏性耗散的程度越小，回热器的不可逆性越小.图2熵产率的截面分布
Fig.2Crosssection entropy generation rate distribution2.2温度梯度对截面平均分布熵产率的影响取参数T0=600 K,T0β=1,f=ω/2π=400 Hz,Z=1,rd=0.000 1 m，ρ0=1.9 kg/m3,σ=0.68和k=0.13 W/(m2·K)，如图3所示：在不同r0/rd时Sir随温度梯度dT/dx的增大而增大，这是由于换热导致的不可逆性受温度梯度影响.图3截面平均熵产率随温度变化规律
Fig.3Cross section average entropy generation
rate versus temperature gradient 2.3针束型回热器的优化频率给定特征尺寸r0/δk和rd/r0时，当角频率ω变化时，存在最小截面平均熵产率Si.由公式（20），当Si/ω=0时可得：
ω0=\[(dT0dx)2·c2pF2(T0cpF4+T20β2C20Z2F1)(1-σ)2\]12(21)ω0为截面平均不可逆熵产率最小时的优化角频率.计算时取参数T0=600 K,T0β=1,dT/dx=1 000 K/m,Z=1,rd=r0/5,ρ0=1.9 kg/m3,σ=0.68和k=0.13 W/(m2·K)，在不同r0/rd时Sir随ω的变化规律如图4所示.图4截面平均熵产率随角频率的变化规律
Fig.4Cross section average entropy generation
rate versus the circular frequency优化频率f0经计算可得：
f0=ω02π=434 Hz (r0δk=3)
425 Hz (r0δk=4)
423 Hz (r0δk=5)2.4针束型回热器的优化阻抗比取参数T0=600 K,T0β=1,f=ω/2π=400 Hz,dT/dx=1 000 K/m,rd=0.000 1 m,ρ0=1.9 kg/m3,σ=0.68和k=0.13 W/(m2·K)，在不同r0/rd时Sir随阻抗比Z的变化规律如图5所示.当Z变化时存在最小Si，对应的最优Z为：Z0=0.34 (r0δk=3)
0.35 (r0δk=4)
0.36 (r0δk=5)图5截面平均熵产率随阻抗比的变化规律
Fig.5Cross section average entropy generation
rate versus impedance ratio 3结语本文对针束型回热器中由于换热和黏性耗散的不可逆性导致的熵产率进行了研究，计算获得了截面熵产率的分布.在给定参数下，截面平均分布熵产率随温度梯度的增大而增大，随阻抗比变化时存在最小值.给定特征尺寸时，获得了最小熵产率对应的最优频率.本文结果有利于选择针束型回热器的优化尺寸.