之前我说过的准确的说90%都是有限差分法。虽然简单,但是他致命的缺点是守恒难以保证。所以想说一下有限体积法。目前,只讨论扩散部分!
能量守恒:
在一个时间段dt,控制体内能增加量=进入控制体的能量-离开控制体的能量。
对于一个微元控制体,假设是标准的正方体,它的能量守恒定律表示为:(U表示控制体的内能,Q表示通过表面的热通量,S表示热源量。)
【在0时刻到1时刻】
如果用温度来表示内能,并且考虑到各处的温度可能不相等的情况,我们引入内能密度平均值:
那么方程(1)的左边项,可以写成:
因为它的意义很明显就是表示不同时间的差值,所以称它为:非稳态项。
再来看(1)式的右边中的通量Q项,如果我们仍然考虑可能发生热通量在时间尺度上也不是均匀值得情况,利用导热的基本方程傅里叶定律:
左边面w:
这个方程很明显是与温度的梯度有关,符合扩散项的形式,可以看作是扩散项。
最后来看源项,并考虑源项的平均值:
现在,来分析一下1D情况,给出几点假设:1.常物性;2.单位体的边长为单位1;这样方程化简为:
现在把一系列的控制体1-n连起来,就是一维导热模型了。
于是,这个物理模型就转化成了n个类似于(5)的方程,每个方程代表一个控制体。
这里,我们做几个合理的假设,第一:如果说把微元体做的无穷小,小到一个近似的质点,那么这个质点就可以代表整个控制体,点的物理量,也可以合理的表示整个控制体的平均值。
于是(5)式的左侧就可以写成:
等式右侧,考虑一下的模型:
第二:利用中心差分的思想,如果我们的微元体足够的小,那么相邻微元体的物理量可以近似看成是线性分布,这样的话就可以利用平均值的思想来表示中间界面的物理量微分:
观察这个式子,对于空间项,我们利用了中心差分的思想,所以他又叫做中心差分格式。
第三:此时,我们在做一个关于时间项的合理假设,如果说时间步长无限小,那么就可以用其中一个时间点的值代替两个时间点之间的平均值。于是对于上式每个时间点可以写:
? 同样观察这个式子,对于时间项,我们取两个时间点的平均值,并将平均值赋给了前的面已知时间点,这种操作叫做时间项的显式;相反如果我们赋值给了后面的时间节点,那么就叫做隐式。
此时(5)式考虑上源项就变成了:
这样(6)就叫做中心差分显式格式的离散方程。
现在反过来想这么一个问题,这中间考虑了很多假设,一旦有假设就意味着它和现实多少式有些出入的,有出入就意味着有误差,这个误差就叫做离散方程的误差,它会随着控制体的体积和时间步长的减少而不断减少。
现在对(6)进行整理化简:
现在对(7)进行进一步化简,将系数简写:
通过观察,可以发现一个有意思的现象:
言外之意就是下一个时刻P处的温度是上一时层W、P、E三点温度加权平均值。这里的系数可以理解为影响因子,物理意义也很明确:ap表达式本身就有热容项,aw和aE是热阻的倒数。于是可以这么理解:如果两边控制体的热阻越小那么他们对在同一时层上的中间控制体的影响就越大;如果中间控制体的热容越大,那么中间控制体在同一时层上受两边控制体的影响就越小!
为方便后续处理,可将(8)式改写为:
为了说明FVM中处理边界控制体和中间部分控制体的区别,举例一个简单的三控制体模型。
中间的控制体,按照前面的结论有:
于是,得:
(符合物理意义,中间控制体的影响因子很大在0.98,两边控制体因为距离相等,所以它两的影响因子相等且均等于0.01.)
? 左侧边界的控制体,类似于引入虚拟点的思路,我们引入一个厚度为0的两个边界“控制体”。他们的温度由边界条件所决定。
于是,参考之前的思路,有:
于是有:
(物理意义也很明确,对于1来说,左边的T0和右边的T2,由于左边的距离是右边的0.5,所以T0影响值是T2的2倍。)
右侧边界的控制体:
得:
(物理意义和左侧得类似。)
于是,这个问题就转化成了如下方程组问题:
解这个方程组,可以手动一个一个算直到解的误差达到要求,也可以利用线性代数的方法,视情况。关于方程组的解法,后面可以单独聊一聊。
可以想到的是,这里我只分了3个网格,如果再细分,那么解也会更加的精确。
*************************************************************************
如果我把网格增加到50个,将速度值导出Excel表格,观察:
出现了一堆“无法显示”的数值,这就是因为我们之前在FDM中说过的:因为50个网格情况下,dx很小,不能满足CFL条件,此时,我们需要适当的减少dt的值;在FVM中“P”控制体的系数一定是小于0了,不符合物理意义了!解决的办法仍然是需要减少dt。减少dt后50个控制体的温度分布如下:
? 随着控制体数量的增加,同一个点的温度变化将趋于恒定。比如:从50个控制体增加到100个控制体相对误差为10的-5次方数量级。