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

一、热传递及流体力学原理：热量传递由三种方式：热传导、热对流和热辐射。1.热传导热量有物体的一部分传到另一部分，或由一个物体传到与其接触的另一个物体的传热现象叫热传导。其特点是，热量流动而物体各部分仍保持着宏观的静止。热传导定律——傅里叶假设傅里叶假设：在单位时间内通过微元等温面dA的热量dQ与成正比，即dQ二_dQ二_k—ndA

dn(1-1)式中：“-”号——热量是从温度高的地方流向温度低的地方，其方向与温度梯度相反；传热系数，它是表征材料导热能力的物理参数，单位是W/(m2•K)；冬 温度场梯度值；dnn——等温面法线方向的单位矢量。热流强度单位时间通过单位等温面积的热量叫热流强度。按照傅里叶假设，热传导热流强度的计算公式为dTq=dQ/dA=一k——n (1-2)dn式中q为热传导热流强度，单位是W/m2。热传导热流强度q是一个矢量。2.对流换热通过流体的流动而换热的现象叫对流换热。对流换热是一个复杂的过程，它既包括流体的对流作用，也包含流体的导热作用。按牛顿公式对流换热的热流强度的计算公式为q=h(T_T)n (1-3)ces式中：q——对流热流强度，单位是W/m2；ch——表面传热系数，单位是W/(m2•K)；T——流体的平均温度；eT——固体表面的温度；sn——指向固体表面的单位法向矢量。对流热流强度q是一个矢量，其方向垂直于固体表面，以指向固体为c正。对流系数h是考虑了影响对热换热的各种因素的一个综合系数，这些因素是流体的速度、流体的物理特征、换热表面的形状和尺寸等。热辐射在本研究中不做考虑。3.管内流体对流换热a.对流速度边界层当具有粘性的流体流过壁面时，就会在壁面上产生粘滞力。粘滞力阻碍了流体的运动，使靠近壁面的速度降低，使直接贴附于壁面的流体近于停滞不动。如图1所示，沿壁面发现y方向不同点的速度v，就是速度分布图。速度分布图表明：从y=0处v=0开始，速度v随着y方向的增加而迅速增大;当经过厚度为8的薄层时，v就接近主流速度v。把所有y=8的点连接成线即为速度边界层，8为边界层厚度。边界层图1.速度边界层b.层流、湍流流体进入管口一段距离后，管壁两侧的边界层才能在管中心回合，这时管断面流速分布和流动状态才达到定型，这段距离通常称为入口段。之后的流速分布不再改变，流态定型，流动达到充分发展，称为流动充分段。如图2、图3所示。

图2层流状态的管内流动流速分布v入口段图2层流状态的管内流动流速分布v入口段充分发展段(1-4)图3紊流状态的管内流动流速分布层流和紊流的流态判断由雷诺准则（1-4）式判定ud(1-4)Re=v式中 Re——雷诺数；v 运动粘度，单位m2Js；u——流体速度；d 圆管内径（方管则约为截面面积丘x2）,单位m当Re<2320时为层流；当2320vRevl04时为过渡时期；当104<Re时为旺盛紊流。层流时，速度分布曲线为抛物线形状；紊流时，速度分布曲线呈对数曲线形状。管子入口段的长度、弯曲管道对流场的影响入口段管内流动的流速分布不稳定，并且随着流距而改变；当流体在弯管内流动时，因离心力而产生二次环流，增加了流动中的扰动，使流动场变为非充分发展流动，并且弯曲半径越小，影响越大。介于上面的情况，如果在分析中希望达到充分发展以保持流场的稳定，必须增长管的长度（弯曲后管的长度）。在紊流状态下管长与管内径之比->50时（层流要小些），可以认为d管内流体已经达到充分发展，即为稳定流场。二、有限元1.原理将一个表示结构或连续体的求解域离散为若干个单元，并通过它们的边界上的节点相互联结成为组合体。用每个单元内所假设的近似函数来分片地表示全求解域内待求的未知场变量。而每个单元内的近似函数由未知场函数(或及其导数，为叙述方便，后面略去此加注)在单元各个结点上的数值和与其对应的插值函数来表达。由于联结相邻单元的结点上，场函数应具有相同的数值，银耳将它们用作数值求解的基本未知量。这样一来，求解原来待求场函数的无穷多自由度问题转换为求解场函数结点值的有限自由度问题。通过和原问题数学模型(基本方程、边界条件)等效的辩分原理或加权余量法，建立求解基本未知量(场函数的结点值)的代数方程组或常微分方程组。2.有限元解决热传导问题a.引言在一般的三维问题中，瞬态温度场的场变量0(x,y,z,)在直角坐标中应满足的微分方程是Pcdt dxPcdt dxvxdx丿dydzvzQz-PQ二0(在o内)(2-1)边界条件是Q0 边界条件是Q0 Q0 Q0kn+kn+kn=qxQxx yQyy zQzz(在r边界上) (2-2)1(在r边界上) (2-3)2Q0 Q0 Q0 Q0 Q0kn+kn+k nxQxx yQyyzQzz二h(0-0)a(在r边界上)3(2-4)式中：P 材料的密度，单位是kg/m3；c 材料的比热容，单位是J/(kg-K)；0=0(r,t)——在r边界上的给定温度；1q=q(r,t)——在r边界上的给定热流密度，单位是w/m2；20=0(r,t)——对于r边界，在自然对流条件下，0是外界环境温度;aa3a在强迫对流条件下，0是边界层的绝热壁温度。a微分方程(2-1)是热量平衡方程。式中第1项是微体升温需要的热量；第2,3,4项是由x,y和z方向传入微体的热量；最后一项是微体内热源产生的

热量。微分方程表明：微体升温所需的热量应与传入微体的热量以及微体内热源产生的热量相平衡。(2-2)式是在r边界上给定温度O=e(r,t),称为第一类边界条件，它是1强制边界条件。(2-3)式是在r边界上给定热流量q=q(r,t),称为第二类边2界条件，当q=0时就是绝热边界条件。(2-4)式是在r边界上给定对流换热3的条件，称为第三类边界条件。第二、三类边界条件是自然边界条件。r+r+r二r，r是域o的全部边界。123当在一个方向上，例如z方向温度变化为零时，方程(2-1)就退化为二维问题的热传导微分方程：pc也

apc也

ata匚a©)——k——axvxax丿k汨-PQ=0Vyay丿(在o内) (2-5)这时当场变量©(x,y,t)不再是z的函数，场变量应满足的边界条件是:帖®(r,t)+k竺n帖®(r,t)+k竺n=q(r,t)yayyk越nxax+k越n=h(©-©)yayya(在r边界上)1(在r边界上)2(在r边界上)3在(2-5)~(2-8)式中，各项符号意义与在(2-1)~(2-4)式中的相同。求解瞬态温度场问题是求解在初始条件下，即在©(x,y,z,0)=©0(x,y,z(2-6)(2-7)(2-8)(2-9)条件下满足瞬态热传导方程及边界条件的场函数©，©应是坐标和时间的函数。b.瞬态热传导问题1)瞬态热传导有限元的一般格式首先建立三维瞬态热传导问题的微分方程(2-1)式和边界条件(2-2)，(2-3)和(2-4)式的等效积分形式，即

3e3c3e] 3匚3ewpc——-——Ik——-——k——3e3c3e] 3匚3ewpc——-——Ik——-——k——9t 9x(x3x丿9y帥JwG-e)dT+JT11JwT3 3k理nx3xx3(帥)_k丄-pQdQ+3zIz3z丿 _7r73e73ek——n+k——n+k——n-qdr+x3xxy3yyz3zz丿wr22+k遊n+k遊n-h(e-0)]dr=0y3yyz3zza丿(2-10)其中w，w，12 3 1已满足©=&，则w=0，并且不失一般地可令1(2-11)W=W=W=8e(2-11)238efPc越〕+98e帥)98^(帥)38^k—+ k8efPc越〕+98e帥)98^(帥)38^k—+ k—+ k3e]-80pQdQ3z(z3z丿JQ+J8eqdr-J8eh(e-e)dr=o厂 厂 a丄2 丄3利用(2-12)式可以建立瞬态温度场有限元的一般格式。首先将空间域Q9t丿9x(x9x丿k——+3yIy9y丿(2-12)离散为有限个单元体，在典型单元内温度e仍可以近似地用结点温度e插值i得到，但是注意此时的结点温度是时间的函数，即e=e=*n(x,y,z)e(t) (2-13)iii=1插值函数N只是空间域的函数，它应具有插值函数的基本性质。将上式i代入(2-12)式，并且考虑8e的任意性，就可以得到用来确定n个结点温度eii的有限元求解方程ce+ke=p (2-14)这是一组以时间t以独立变量的线性常微分方程组。式中：C——热容矩阵；K——热传导矩阵，C和K(在引入给定温度以后)都是对称正定矩阵；P——温度在和列阵；e——结点温度列阵；ef=业］——结点温度对时间的导数列阵。Idt丿

矩阵C,K和P的元素由单元的相应的矩阵元素集成，即K丄Ke+ZHeijeC=工ije矩阵C,K和P的元素由单元的相应的矩阵元素集成，即K丄Ke+ZHeijeC=工ijeKe+He

ij ijeCeij(2-15)式中：P=工iPe+QieKeijPe+ Peqi Hiee单元对热传导矩阵的贡献Heij单元热交换边界对热传导矩阵的修改；Ceij单元对热容矩阵的贡献；PeQi单元热源产生的温度载荷；Pe

qi单元给定热流边界的温度载荷；厂6N6N 6N6N 6N6N)Ke=fij3(2-16)He=ihNNdrijreij3Ce=fpcNNdoij oe ijPe=fpQNdoQioe iPe=fqNdrqKe=fij3(2-16)He=ihNNdrijreij3Ce=fpcNNdoij oe ijPe=fpQNdoQioe iPe=fqNdrqi re2iPe=fhNdrH reaii re3(2-17)(2-18)(2-19)(2-20)(2-21)至此，已经将时间域和空间域的偏微分方程问题在空间域内离散为N个结点温度0(t)的常微分方程的初值问题。对于给定温度值的边界r上n个结i 1 1点，方程组(2-14)式中相应项应引入的条件是(2-22)0=0(i=1,2,…,n)(2-22)ii 1式中，i是r上n个结点的编号。112)一阶常微分方程的求解瞬态热传导问题的有限元求解方程(2-14)式是一阶常微分方程组。以下讨论有限元的直接积分法求解该方程。直接积分法是指求解常微分方程组前不对它的形式进行变换，而是直接对它进行数值积分。它通常基于两个基本概念，一是求解的时间域0<t<T划分为若干个时间步长At，在一定数目的At时间区域内，假设Q和$的函数形式来近似方程的精确解；二是以仅在相隔At的离散时间点上满足微分方程来代替时间域内任何时刻t都满足微分方程。在以下的讨论中首先将时间域0〜T等分(不限于等分，这里是为了讨论方便)为M个时间步长，At二TIM，并认为t=0时的初始温度列阵$是已知0的。进一步假设t=0，t(二At)，t(二2At),•••,t(=ntA)时刻的解$已经TOC\o"1-5"\h\z1 2 n n求得。下一步要计算的是t(=t+At)时刻的$ 。从该求解过程建立起求解n+1 n n+1所有离散时间点场函数的一般算法步骤。2.1)用加权余量法建立两点循环公式对于只有一阶导数的常微分方程组(2-14)式，可以再两个时间点t和tn n+1之间的At时间区域内，假设$和$采取如图2所示的线性插值形式，即$(t+S)At)=N$+N$ (2-23)n nn n+1n+1其中N=1—爲 (0<g<1)nN=£n+1(2-23)式对时间t求导，可得(2-24)$=N$+N$(2-24)nn n+1n+1其中1Atn+n+1AtJ(a)(b)(c)o=2(d)(e)图1建立两点循环公式的插值函数及权函数(a)(b)(c)o=2(d)(e)图1建立两点循环公式的插值函数及权函数由于采用(2-23)式的近似插值，在时间区域At内，方程(2-14)式将产生余量。对于这一时间区域，典型的加权余量格式可以表示为如下形式：n ，"‘丄 n ，"‘丄 n+1丿+kN©+N©n n+1)n ，"‘丄 n+1丿-Pdg二0(2-25)当求解初值问题时，如果已知一组参数©，则可以利用(2-25)式近似确n定另一组参数©。将(2-23)和(2-24)式代入(2-25)式就可以得到时间区域n+1At前后结点上两组参量的关系式，即fKJ1fKJ1诞dg+CJ1w朮'Vo oAt丿-J1wPd^=000+fKJ1w(1弋)dg-CJ1w朮'00At丿n+1(2-26)式中可以代入不同的权函数。在以上讨论中假定热传导矩阵K和热容矩阵C不随时间t而变化。(2-26)式可以表示为任何权函数都使用的一般形式，即(C/At+K0)0 +「-C/At+K(1-8)0]=P (2-27)n+1 n式中(2-28)J1诞d^/f1wdgJ1wPd^.J1wdg0■'0(2-28)一种方便的做法就是假定P采用与未知场函数0相同的插值表达式，这时将得到P=P0+P(1-0) (2-29)TOC\o"1-5"\h\zn+1 n当0和P都已知时，就可以由(2-27)式求得下一时刻的0 。这就是连点循环nn n+1公式。可以更明显地表示成(2-30)(2-31)K0=Q (n=0,1,2,M)(2-30)(2-31)\o"CurrentDocument"n+1 n+1其中K=C/At+0KQ=「C/At-(l-0)K]0+(1-0)P+0P>n+1「 」n n n+1丿利用上式，从t=0出发，提一次一次递推求得结点温度列阵0(t)的各个瞬时值0,0，…,0。1 2M2.2)算法步骤利用直接积分的两点循环公式递推求解瞬态热传导问题的有限元微分方程的算法步骤可以归结如下：初始计算Q形成系统系数矩阵C和K；Q给定0；0Q选择参数0和时间步长At；Q形成系统的有效系数矩阵K=C/At+0K；

⑨三角分解K,K=LDLT.对于每一时间步长：Q形成向量P；n+19形成有效向量Qn+1Pn+1二「C/At-(1-O)k伸+(i-e)9形成有效向量Qn+1Pn+1TOC\o"1-5"\h\zn n\o"CurrentDocument"9回代求解0 ,LDL0 =Q。n+1 n+1 n+12.3)参数0的选择以上通过加权余量法得到的求解瞬态热传导方程的两点循环公式(2-27)本质上是一组加权的差分公式，因为在每个时间区域At内，对于(t+0At)点n(0<0<1)建立的差分公式为(2-32)0(t+0At)=(l-0)0+00(2-32)n n n+10(t+0At)=(0n n+1(2-33)该两式实际上就是利用时间区域At内线性插值公式得到的(2-23)和(2-24)式,只是坐标参数g换成0而已。将(2-32)和(2-33)式代入(2-14),并将P表示成和0相同的差分形式，则得到建立于(t+0At)时刻的差分方程。它和前面利用加权余量法的两点循环公n式(2-27)式相同。从此可以清楚地认识0的物理意义，0的取值决定了在At时间区域内建立差分方程的具体地点。从图2给出的一组权函数及相应的0值1的关系也可以看出0的物理意义。图2中的(a)、(b)、(c)集中在点n，n+丄，2n+1上加权，得到的是著名的前差分公式(欧拉差分公式)，中心差分公式(Crank-Nicholson差分公式)和后差分公式。(d)是在时间区域内权函数等于常数，其结果和中心差分的结果相同。(e)和⑴为伽辽金型的权函数，其结果分12别和权函