COMSOL中有限元解题过程示例——均匀介质中的热传递

我们大家在大学数学中学习过一些基本的偏微分方程,相信对它有一定的了解,它的一般形式如下: http://s11/middle/6cc1d6b1tc748e5115aea&690,
其中,前两项是与时间有关的项,分别为函数u随着时间变化的一阶导和二阶导,http://s11/middle/6cc1d6b1tc7490131fdda&690。
对于偏微分方程的求解方法有多种,随着计算机的发展,越来越多的人开始接触并运用数值的方法来求解偏微分方程的解,而数值解法有多种,我们的重点是用有限元法来求解。现在举个例子结合数理方程运用数值解法中的有限元方法来求解一个物理中常见的传热问题,并用仿真软件——COMSOL Multiphysics来实现它的求解过程。
例:有一个均匀的圆柱体结构,例如铜棒、铁棒等,在整个棒上分布一个恒定的热源Q=10J/(sm),棒的长度为L=1m,横截面积为A=1m2,穿过单位横截面积的热流是q=1.25J/(m2s),棒的左端给定了初始温度为T0=1K。在棒的热传导满足以下方程:http://s13/middle/6cc1d6b1tc749045fc3ec&690。
看到这个方程我们先想到什么呢,首先这是一个热传导方程,在求解域中也就是整个棒中满足以上控制条件,Q代表源项,由于在棒的横截面上的热流是均匀的,可以将其简化为只沿着x轴的一维传热模型,所以上式中的自变量为x,而不是空间坐标。
下面用有限元的方法进行分析和求解,有限元方法求解的步骤如下:
第一步,首先是进行变分,将方程化解为积分形式,这样就降低了原函数的连续性,即变为方程的“弱”形式,选择一个任意的函数w,作为试函数,在上述微分方程的两边同乘以w后进行积分,分布积分后可得到下式: http://s8/middle/6cc1d6b1tc74908a166d7&690这个式子为(1)式。
第二步,离散求解域,这里进行几何离散的依据是形函数(此步骤对应COMSOL中划分网格),将一维模型划分为三个单元,四个节点,由上所述,此一维模型中的相邻节点处的温度随着x的变化可以认为满足下式:http://s12/middle/6cc1d6b1tc7490d1ca1eb&690(T为形函数),
将其联立求解后可以得到a和b的值,将得到的a和b分别带入到T=a+bx中,可得到基元内任意点的温度变化:http://s10/middle/6cc1d6b1tc7491096b709&690,
令Ni=(1/l)(xj-x),Nj=(1/l)(x-xi)(Ni为基函数),所以整个结构上的温度为http://s13/middle/6cc1d6b1tc7491259384c&690,
将此式带入到(1)式中,可以得到下式
http://s3/middle/6cc1d6b1tc749142d2c62&690此式为(2)式。
第三步,进行单元分析,得到单元刚度矩阵,首先分析第一个单元,它与第一个和第二个节点有关,上式中的Ik可以用下式来表达:http://s4/middle/6cc1d6b1tc74916d794f3&690,
可以将其写成矩阵的形式也就是单元刚度矩阵了,同理,第二个单元与第二个和第三个节点有关,也可以写成类似上述单元刚度矩阵的形式。
第四步,将以上得到的三个单元刚度矩阵组合成整体刚度矩阵。
由于在一个单元和第二个单元处共用一个节点2,所以在节点2处的值是有加和关系的,也就是将矩阵中对角线上的元素在单元中有重复的部分做一个加和后就得到总体刚度矩阵。
而(2)式中的等号右边Ib和Is分别代表边界项和源项,它们的值可以通过数值计算得到,此时,将(2)式左侧的T提出来,可以写成两个矩阵乘积的形式,如下:http://s9/middle/6cc1d6b1tc749289d1ed8&690,
由于在左侧的温度值T1是已知的,所以没有必要求解未知的q值了,可以将上式中的第一行和第一列去掉,这样就变成了一个3×3矩阵了,而N的值也是可以求出来的,从而可以求解除每个节点处的T值,结果为T1=1,T2=1.7,T3=2.01,T4=2.11。
下面,我们进一步用COMSOL软件来进行以上求解过程。
首先,分析问题,选择模型,由上述过程我们知道,这个物理问题我们可以简化为一维模型,它的传热过程是在均匀介质中进行的,即选择模型为固体中的热传递,偏微分方程可知它与时间无关,是一个稳态问题。我们在软件中首先选择好一维稳态固体传热模型,之后进行创建几何体,一维模型中建立一个0到1的线段,在固体传热模型中加入一些参数值和条件,有热传导系数设定为3.3J/(Kms),左侧端点的温度值为T0=1,整个区域中的热源Q=10J/(ms),右侧端点的热流为q=1.25J/(m2s),这些步骤都是在例题中直接给出以及进行初步分析得到的,进行以上操作后,要进行网格剖分,对应有限元求解中的离散几何域,然后进行计算,得到结果如图所示,从结果中我们可以得到T1=1,T2=1.96,T3=2.59,T4=2.89。此结果比用有限元计算的更加精确,由于在有限元中我们只分为了三个单元,但是在软件中我们可以细化网格,也就是将划分的单元数增多,所以得到的结果更加精确。
http://s4/middle/6cc1d6b1tc7492ecd9d93&690