基于python的热传导_差分


Euler法原理

Euler法是一般数值分析课程中最基础也最重要的一种数值方法,用于求解微分方程的初值问题。一阶常微分方程的初值问题一般形式是



数值方法就是要求问题(1)的解


在若干离散点的近似值


,(n=1,2,...)的方法,建立数值解法,首先要将微分方程离散化,转化为差分方程,用数值积分公式或泰勒展开法




由此得到迭代公式



用差分方程初值问题的解来近似微分方程初值问题(1)的解, 即由(3)依次算出


的近似值


,(n=1,2,...)。这组计算式称为

向前 Euler 法。

向前 Euler 法是显式一阶方法,本质上是利用第N点的函数值和一阶导数(斜率)去近似求解第n+1点的函数值,因此它的精度不高,好处是不需要迭代求解。

如果能够利用第N+1点的一阶导数来修正误差,就能够得到更高的精度。直观上用梯形公式计算数值积分要比矩形公式好,梯形公式利用第N+1点的斜率,与第N点的斜率取平均得到等效的斜率,具有二阶精度,但是它是隐式格式,一般需用迭代法解方程



如果实际计算时精度要求不太高,用公式(3)求解时,每步可以只迭代一次,则得到改进的Euler 法,相当于将Euler公式与梯形公式结合使用:先用 Euler 公式求


的一个初步近似值


(预测值),然后用梯形公式校正求得近似值


。为了编程方便,把斜率单独提取出来计算等效斜率,迭代式如下



改进Euler法具有二阶精度,但收敛率略小于梯形法。


二维热传导方程

二维热传导方程(常物性)具有如下形式:



是热扩散系数,


为热源项,初值条件


。采用Euler法对时间项进行离散化:



Euler法只处理了时间偏导项(向前差分),进行接下来对空间二阶偏导进行离散化。用


表示第i行(y方向)第j列(x方向)的温度值,t时刻的温度场以(Ny+1,Nx+1)大小的矩阵存储,


代表x,y方向的分段数量,每一段的等分长度为





基于python的热传导_差分_02

矩阵的每个值代表一个网格节点温度,注意y方向代表行标,x方向代表列标


每一点对


,


方向均采用中心差分格式:



把二阶偏导组装成矩阵形式,维度与


相同




写成矩阵乘积形式,记为




称为差分系数矩阵,均为三对角对称矩阵,但是大小不同,


为(Nx+1,Nx+1),


为(Ny+1,Ny+1)。这样迭代公式变为



这里需要注意式(9-1)和(9-2)中求得的二阶偏导矩阵的边界点是无效的,因为差分系数矩阵第一行和最后一行是不完整的[-2 1]和[1 -2],相当于边界点的温度恒为0,因此需要单独给定边界条件才能确定矩阵边界的温度值。


Python计算代码

因为给出了矩阵化的迭代式,所以编程就很简单了,直接采用numpy的向量化运算。计算10秒,前5秒加热,后5秒不加热,边界条件采用第一类边界,即给定温度值。


import


基于python的热传导_二维稳态热传导 代码实现_03


基于python的热传导_微分方程的数值解法与程序实现 pdf_04


由于向前Euler法只有一阶精度,当


太大或


太小均会造成函数不收敛,得到错误数值结果。当Nt<2000,即dt>0.005时出现不收敛现象。


基于python的热传导_差分_05


由于Euler法的向量化运算十分简单快速,实际操作中可以用较大的时间分段Nt,去增加总的时间精度而不会明显降低求解速度,当需要更高精度和速度时可以使用改进Euler法或Runge-Kutta法,工程上常用四阶精度的Runge-Kutta法求解。


总结

Euler法是显式一阶方法,可以快速求解ODE初值问题。

热传导方程本质上是PDE方程,在采用差分离散化空间变量后转换为时间变量的ODE方程,方可采用Euler法求解,其最终格式与热传导显式差分格式完全相同。

在处理温度U节点及其空间二阶导时采用向量化和矩阵化形式,避免编程中使用多重for-loop,而是采用向量化运算提高效率。

一般小规模问题,可用较大的时间分段Nt去提高收敛性和精度,而不必采用复杂的迭代式。

下一篇文章将讲解如何处理第二类边界条件和非矩形的边界形状。