液体激光器热流场的基本理论

一、热传递及流体力学原理：
热量传递由三种方式：热传导、热对流和热辐射。

1. 热传导
热量有物体的一部分传到另一部分，或由一个物体传到与其接触的另一个物体的传热现象叫热传导。

其特点是，热量流动而物体各部分仍保持着宏观的静止。

a. 热传导定律——傅里叶假设
傅里叶假设：在单位时间内通过微元等温面dA 的热量d Q 与成正比，即
d dA Q n T
k
n
∂=-∂ (1-1) 式中：“-”号——热量是从温度高的地方流向温度低的地方，其方向与温度梯度相反；
k ——传热系数，它是表征材料导热能力的物理参数，单位是()
2W/m K ⋅；
T
n
∂∂——温度场梯度值； n ——等温面法线方向的单位矢量。

b. 热流强度
单位时间通过单位等温面积的热量叫热流强度。

按照傅里叶假设，热传导热流强度的计算公式为
d /dA=q Q n T
k
n
∂=-∂ (1-2) 式中q 为热传导热流强度，单位是2W/m 。

热传导热流强度q 是一个矢量。

2. 对流换热
通过流体的流动而换热的现象叫对流换热。

对流换热是一个复杂的过程，它既包括流体的对流作用，也包含流体的导热作用。

按牛顿公式对流换热的热流强度的计算公式为
()q n c e s h T T =- (1-3)
式中： q c ——对流热流强度，单位是2W/m ； h ——表面传热系数，单位是()
2W/m K ⋅； e T ——流体的平均温度； s T ——固体表面的温度；
n ——指向固体表面的单位法向矢量。

对流热流强度q
是一个矢量，其方向垂直于固体表面，以指向固体为
c
正。

对流系数h是考虑了影响对热换热的各种因素的一个综合系数，这些因素是流体的速度、流体的物理特征、换热表面的形状和尺寸等。

热辐射在本研究中不做考虑。

3. 管内流体对流换热
a. 对流速度边界层
当具有粘性的流体流过壁面时，就会在壁面上产生粘滞力。

粘滞力阻碍了流体的运动，使靠近壁面的速度降低，使直接贴附于壁面的流体近于停滞不动。

如图1所示，沿壁面发现y方向不同点的速度v，就是速度分布图。

速度分布图表明：从y=0处v=0开始，速度v随着y方向的增加而迅速增大；。

把所有y=δ的点连接成线当经过厚度为δ的薄层时，v就接近主流速度v
∞
即为速度边界层，δ为边界层厚度。

边界层
图1. 速度边界层
b. 层流、湍流
流体进入管口一段距离后，管壁两侧的边界层才能在管中心回合，这时管断面流速分布和流动状态才达到定型，这段距离通常称为入口段。

之后的流速分布不再改变，流态定型，流动达到充分发展，称为流动充分段。

如图
2、图3所示。

图2层流状态的管内流动流速分布
图3紊流状态的管内流动流速分布
层流和紊流的流态判断由雷诺准则(1-4)式判定
ud
Re v
= (1-4) 式中 Re ——雷诺数；
v ——运动粘度，单位2m ； u ——流体速度；
d ——圆管内径(
2)，单位m 当2320Re <时为层流；当42320<<10Re 时为过渡时期；当410Re <时为旺盛紊流。

层流时，速度分布曲线为抛物线形状；紊流时，速度分布曲线呈对数曲线形状。

c. 管子入口段的长度、弯曲管道对流场的影响
入口段管内流动的流速分布不稳定，并且随着流距而改变；当流体在弯管内流动时，因离心力而产生二次环流，增加了流动中的扰动，使流动场变为非充分发展流动，并且弯曲半径越小，影响越大。

介于上面的情况，如果在分析中希望达到充分发展以保持流场的稳定，必须增长管的长度(弯曲后管
的长度)。

在紊流状态下管长与管内径之比50L
d
>时(层流要小些)，可以认为
管内流体已经达到充分发展，即为稳定流场。

v
二、有限元
1. 原理
将一个表示结构或连续体的求解域离散为若干个单元，并通过它们的边界上的节点相互联结成为组合体。

用每个单元内所假设的近似函数来分片地表示全求解域内待求的未知场变量。

而每个单元内的近似函数由未知场函数(或及其导数，为叙述方便，后面略去此加注)在单元各个结点上的数值和与其对应的插值函数来表达。

由于联结相邻单元的结点上，场函数应具有相同的数值，银耳将它们用作数值求解的基本未知量。

这样一来，求解原来待求场函数的无穷多自由度问题转换为求解场函数结点值的有限自由度问题。

通过和原问题数学模型(基本方程、边界条件)等效的辩分原理或加权余量法，建立求解基本未知量(场函数的结点值)的代数方程组或常微分方程组。

2. 有限元解决热传导问题 a. 引言
在一般的三维问题中，瞬态温度场的场变量(),,,x y z t φ在直角坐标中应满足的微分方程是
0x y z c
k k k Q t x x y y z z φφφφρρ⎛⎫∂∂∂∂∂∂∂⎛⎫⎛⎫----= ⎪ ⎪ ⎪∂∂∂∂∂∂∂⎝⎭⎝⎭
⎝⎭ (在Ω内) (2-1) 边界条件是
φφ
= (在1Γ边界上) (2-2)
x
x y y z z k n k n k n q x y z
φφφ∂∂∂++=∂∂∂ (在2Γ边界上) (2-3) ()x
x y y z z a k n k n k n h x y z
φφφφφ∂∂∂++=-∂∂∂ (在3Γ边界上) (2-4) 式中： ρ——材料的密度，单位是3kg/m ； c ——材料的比热容，单位是()J/kg K ⋅； (),t φφ=Γ——在1Γ边界上的给定温度；
(),q q t =Γ——在2Γ边界上的给定热流密度，单位是2W/m ；
(),a a t φφ=Γ——对于3Γ边界，在自然对流条件下，a φ是外界环境温度；在强迫对流条件下，a φ是边界层的绝热壁温度。

微分方程(2-1)是热量平衡方程。

式中第1项是微体升温需要的热量；第2，3，4项是由,x y 和z 方向传入微体的热量；最后一项是微体内热源产生的
热量。

微分方程表明：微体升温所需的热量应与传入微体的热量以及微体内热源产生的热量相平衡。

(2-2)式是在1Γ边界上给定温度(),t φφ=Γ，称为第一类边界条件，它是强制边界条件。

(2-3)式是在2Γ边界上给定热流量(),q q t =Γ，称为第二类边界条件，当0q =时就是绝热边界条件。

(2-4)式是在3Γ边界上给定对流换热的条件，称为第三类边界条件。

第二、三类边界条件是自然边界条件。

123Γ+Γ+Γ=Γ，Γ是域Ω的全部边界。

当在一个方向上，例如z 方向温度变化为零时，方程(2-1)就退化为二维问题的热传导微分方程：
0x y c k k Q t x x y y φφφρρ⎛⎫
∂∂∂∂∂⎛⎫---= ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭ (在Ω内) (2-5) 这时当场变量(),,x y t φ不再是z 的函数，场变量应满足的边界条件是：
(),t φφ=Γ (在1Γ边界上) (2-6)
(),x
x y y k n k n q t x y
φφ∂∂+=Γ∂∂ (在2Γ边界上) (2-7) ()x
x y y a k n k n h x y
φφφφ∂∂+=-∂∂ (在3Γ边界上) (2-8) 在(2-5)~(2-8)式中，各项符号意义与在(2-1)~(2-4)式中的相同。

求解瞬态温度场问题是求解在初始条件下，即在
()()00,,,,,x y z x y z φφ= (2-9)
条件下满足瞬态热传导方程及边界条件的场函数φ，φ应是坐标和时间的函数。

b. 瞬态热传导问题
1) 瞬态热传导有限元的一般格式
首先建立三维瞬态热传导问题的微分方程(2-1)式和边界条件(2-2)，(2-3)和(2-4)式的等效积分形式，即
()()123123d d d d 0x y z x x y y z z x x y y z z a w c k k k Q t x x y y z z w w k n k n k n q x y z w k n k n k n h x y z φφφφρρφφφφφφφφ
φφΩΓΓΓ⎡⎤⎛⎫∂∂∂∂∂∂∂⎛⎫⎛⎫----Ω+⎢
⎥ ⎪ ⎪ ⎪∂∂∂∂∂∂∂⎝⎭⎝⎭⎝⎭⎣⎦
⎛⎫∂∂∂-Γ+++-Γ+ ⎪∂∂∂⎝⎭⎛⎫∂∂∂++--Γ= ⎪∂∂∂⎝⎭
⎰⎰⎰⎰ (2-10) 其中w ，1w ，2w ，3w 是任意函数。

按伽辽金方法选择任意函数，设1Γ上已满足φφ=，则10w =，并且不失一般地可令
23w w w δφ=== (2-11)
将上式代入(2-10)式，并对其中第1个域Ω内积分的第2，3和4项进行分部积分，则可得到
()2
3
d d d 0
x y z a c k k k Q t x x y y z z q h φδφφδφφδφφδφρδφρδφδφφφΩΓΓ⎡⎤⎛⎫∂∂∂∂∂∂∂⎛⎫⎛⎫⎛⎫+++-Ω⎢⎥ ⎪ ⎪ ⎪ ⎪∂∂∂∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦+Γ--Γ=⎰⎰⎰ (2-12) 利用(2-12)式可以建立瞬态温度场有限元的一般格式。

首先将空间域Ω离散为有限个单元体，在典型单元内温度φ仍可以近似地用结点温度i φ插值得到，但是注意此时的结点温度是时间的函数，即
()()1
,,e
n i i i N x y z t φφφ===∑ (2-13)
插值函数i N 只是空间域的函数，它应具有插值函数的基本性质。

将上式代入(2-12)式，并且考虑i δφ的任意性，就可以得到用来确定n 个结点温度i φ的有限元求解方程
φφ+=C K P (2-14)
这是一组以时间t 以独立变量的线性常微分方程组。

式中： C ——热容矩阵；
K ——热传导矩阵，C 和K (在引入给定温度以后)都是对称正定
矩阵；
P ——温度在和列阵；
φ——结点温度列阵；
φd d t φ⎛⎫
=
⎪⎝⎭
——结点温度对时间的导数列阵。

矩阵C ，K 和P 的元素由单元的相应的矩阵元素集成，即
e e
ij ij ij e
e
K K H =+∑∑
e ij ij e
C C =∑ (2-15)
i i i
e e e i Q q H e
e
e
P P P P =++∑∑∑ 式中： e
ij K ——单元对热传导矩阵的贡献；
e
ij H ——单元热交换边界对热传导矩阵的修改； e ij C ——单元对热容矩阵的贡献； i
e Q P ——单元热源产生的温度载荷；
i
e q P ——单元给定热流边界的温度载荷； i
e H
P ——单元的对流换热边界的温度载荷。

这些单元的矩阵元素由下列各式给出：
d e j j j e
i i i ij x y z N N N N N N K k k k x x
y y z z Ω∂∂∂⎛⎫∂∂∂=++Ω ⎪∂∂∂∂∂∂⎝⎭⎰ (2-16) 3
d e e ij i j H hN N Γ=Γ⎰ (2-17)
d e e ij i j C cN N ρΩ
=Ω⎰ (2-18)
d e i
e Q i P QN ρΩ
=Ω⎰ (2-19) 2
d e i e q i P qN Γ=Γ⎰ (2-20)
3
d e i
e H a i P h N φΓ=Γ⎰ (2-21) 至此，已经将时间域和空间域的偏微分方程问题在空间域内离散为N 个结点温度()i t φ的常微分方程的初值问题。

对于给定温度值的边界1Γ上1n 个结点，方程组(2-14)式中相应项应引入的条件是
()112i
i i n φφ==⋅⋅⋅,,, (2-22)
式中，i 是1Γ上1n 个结点的编号。

2) 一阶常微分方程的求解
瞬态热传导问题的有限元求解方程(2-14)式是一阶常微分方程组。

以下
讨论有限元的直接积分法求解该方程。

直接积分法是指求解常微分方程组前不对它的形式进行变换，而是直接对它进行数值积分。

它通常基于两个基本概念，一是求解的时间域0t T ≤≤划分为若干个时间步长t ∆，在一定数目的t ∆时间区域内，假设φ和φ的函数形式来近似方程的精确解；二是以仅在相隔t ∆的离散时间点上满足微分方程来代替时间域内任何时刻t 都满足微分方程。

在以下的讨论中首先将时间域0～T 等分(不限于等分，这里是为了讨论方便)为M 个时间步长，/t T M ∆=，并认为0t =时的初始温度列阵0φ是已知的。

进一步假设0t =，()1t t =∆，()22t t =∆，···，()n t nt =∆时刻的解n φ已经求得。

下一步要计算的是()1n n t t t +=+∆时刻的1n φ+。

从该求解过程建立起求解所有离散时间点场函数的一般算法步骤。

2.1） 用加权余量法建立两点循环公式
对于只有一阶导数的常微分方程组(2-14)式，可以再两个时间点n t 和1
n t +之间的t ∆时间区域内，假设φ和φ采取如图2所示的线性插值形式，即
()11n n n n n t t N N φξφφ+++∆=+ (2-23)
其中
11n n t t N N ξξξ+⎫
=
⎪
∆⎪
=-⎬⎪
=⎪⎭
()01ξ≤≤
(2-23)式对时间t 求导，可得
11n n n n N N φφφ++=+ (2-24)
其中
111n n N t N t +⎫=-
⎪⎪∆⎬
⎪=
⎪∆⎭
图1 建立两点循环公式的插值函数及权函数
由于采用(2-23)式的近似插值，在时间区域t ∆内，方程(2-14)式将产生余量。

对于这一时间区域，典型的加权余量格式可以表示为如下形式：
1
11110d 0C K P n n n n n n n n w N N N N φφφφξ++++⎡⎤⎛⎫⎛⎫+++-= ⎪ ⎪⎢
⎥⎝⎭⎝⎭⎣⎦
⎰ (2-25) 当求解初值问题时，如果已知一组参数n φ，则可以利用(2-25)式近似确定另一组参数1n φ+。

将(2-23)和(2-24)式代入(2-25)式就可以得到时间区域
1
0,w θ==
12θ=
1
θ=12
θ=
13
θ=
23θ=
(a)
(b)
(d)
(c)
(f)
(e)
t ∆前后结点上两组参量的关系式，即
()1111100001
d d d 1-d d 0
K C K C P n n w w w w t t w ξξξξφξξφξ+⎛⎫⎛
⎫++- ⎪ ⎪∆∆⎝
⎭⎝
⎭-=⎰⎰⎰⎰⎰ (2-26) 式中可以代入不同的权函数。

在以上讨论中假定热传导矩阵K 和热容矩阵C
不随时间t 而变化。

(2-26)式可以表示为任何权函数都使用的一般形式，即
()()11C/K C/K P n n t t θφθφ+⎡⎤∆++-∆+-=⎣⎦ (2-27)
式中
11
011
00d d d d P P w w w w θξξ
ξξξ⎫=⎪⎬⎪
=⎭
⎰⎰⎰⎰ (2-28) 一种方便的做法就是假定P 采用与未知场函数φ相同的插值表达式，这时将得到
()11P P P n n θθ+=+- (2-29)
当n φ和P n 都已知时，就可以由(2-27)式求得下一时刻的1n φ+。

这就是连点循环公式。

可以更明显地表示成
11K Q n n φ++= ()012,,,,n M =⋅⋅⋅ (2-30)
其中
()()1111K C K
Q C K P P n n n n t t θθφθθ++⎫=∆+⎪
⎬⎡⎤=∆--+-+⎪⎣⎦⎭
(2-31)
利用上式，从0t =出发，提一次一次递推求得结点温度列阵()t φ的各个瞬时值12,,,M φφφ⋅⋅⋅。

2.2) 算法步骤
利用直接积分的两点循环公式递推求解瞬态热传导问题的有限元微分方程的算法步骤可以归结如下：
初始计算
○
1形成系统系数矩阵C 和K ； ○2给定0
φ； ○
3选择参数θ和时间步长t ∆； ○
4形成系统的有效系数矩阵K C K t θ=∆+；
○5三角分解K ，K LDL T
=.
对于每一时间步长：
○1形成向量1
P n +； ○2形成有效向量()()1111Q C K P P n n n n t θφθθ++⎡⎤=∆--+-+⎣⎦；
○3回代求解1n φ+，11LDL Q T n n φ++=。

2.3) 参数θ的选择
以上通过加权余量法得到的求解瞬态热传导方程的两点循环公式(2-27)本质上是一组加权的差分公式，因为在每个时间区域t ∆内，对于()n t t θ+∆点
()01θ≤≤建立的差分公式为
()()11n n n t t φθθφθφ++∆=-+ (2-32)
()()1n n n t t t φθφφ++∆=- (2-33)
该两式实际上就是利用时间区域t ∆内线性插值公式得到的(2-23)和(2-24)式，只是坐标参数ξ换成θ而已。

将(2-32)和(2-33)式代入(2-14)，并将P 表示成和φ相同的差分形式，则得到建立于()n t t θ+∆时刻的差分方程。

它和前面利用加权余量法的两点循环公式(2-27)式相同。

从此可以清楚地认识θ的物理意义，θ的取值决定了在t ∆时间区域内建立差分方程的具体地点。

从图2给出的一组权函数及相应的θ值
的关系也可以看出θ的物理意义。

图2中的(a)、(b)、(c)集中在点n ，1
2
n +，
1n +上加权，得到的是著名的前差分公式(欧拉差分公式)，中心差分公式(Crank-Nicholson 差分公式)和后差分公式。

(d)是在时间区域内权函数等于常数，其结果和中心差分的结果相同。

(e)和(f)为伽辽金型的权函数，其结果分
别和权函数集中在13n +和2
3
n +点上加权时相同。

加权余量法在何处集中加
权，也就是要求在该处微分方程得到满足，因而θ的取值将直接影响到解的精度和稳定性。

三、模型建立
1. 有限元模型建立
单面均匀泵浦横流抽运无机液体激光器的有限元模型如图3所示，其中z ，y ，x 方向分别为液体介质的长，宽，高。

其中泵浦方向是沿着y 轴向右，流场方向是沿着x 轴向上，出光方向则是沿着z 轴。

其中参数如表1所示。

图3 单面均匀泵浦横流抽运无机液体激光器模型
表1 无机液体介质参数
液体介质高液体介质宽
方向 x
y
z
光源中的一块液体介质来进行分析，如图4所示(由于泵浦光源的均匀性，所以取z α=其中()010mm α<<)。

图4 有限元分析中建立的液体介质模型
2. 问题描述 如图4所示，建立的有限元模型中矩形液体介质长宽分别60mm 、8mm ；单面均匀泵浦光强度为5610⨯ 2W m ,假设在达到液体介质面之前没有衰减；入口流体温度为20度，假设在进入泵浦光作用范围已经达到充分发展即流场为稳定流场；上下两边考虑空气冷却作用；每个周期泵浦时间为0.3ms ；
分别对○
1泵浦光源稳定(持续泵浦或者频率
为10Hz)，但流速的不同(0,0.5,1,1.5,2v = ()m s )；○2流速稳定控制在1v = ()m 的流场，但泵浦频率变化(1 Hz ， 10 Hz ，100 Hz)
泵浦方向
流场方向。