原理:由物体内能方程推导,原理是物体内能的变化=进出物体的内能+物体本身的能量。如下所示
其中,方程左侧代表的是物体内能的变化,右侧代表的是六个方向的进出能量,与物体本身的内能变化
由于是二维问题,所以Z方向上的B与F皆为0,所以可以化解为0,其中导热项由傅里叶定律可以变换,故最终推到式如下所示
然后进行合并同类项,可以将方程化解成以下格式 ,其中四个方向的1时刻的温度是我们所不知道,也就是我们所有计算的最终的结果,所以这里采用迭代的方法进行推算,当然当四个方向的温度采用原时刻也就是0时刻的时候,就变成显示求解了。
所以实际的工作分为两步,第一步进行网格划分,这里该问题仅针对矩形网格,调用的方法在如下
def fit_rectangle(self,n_length,m_width,length,width):
grid=Grid
grid.fit_rectangle(grid,5,5,3,2)
这里的n_length与m_width代表得是y方向以及x方向得网格长度,length与width代表的是实际区域的长和宽物理长度,其中length是Y方向,width是X方向。所以上面这行代码得意思就是一块长宽为3m*2m得矩形区域,划分成5*5得网格。
在进行完网格设定以后,接下来就是物理环境得设定,这里需要输入一个字典形式的变量,这里设定密度ρ为100,比热容c=1000,导热系数k=10,时间步长🔺t=100s,物体内能s=10
solver=Implicit_format_conduction
base={'rou':100,'c':1000,'k':10,'dett':100,'s':10}
solver.fit(solver,grid,base)
solver.boundary_condition(solver,'even','up',4)
solver.boundary_condition(solver,'even','down',2)
solver.boundary_condition(solver,'even','right',5)
solver.boundary_condition(solver,'even','left',3)
再然后就是边界条件的设定,这里采用第一类边界条件,也就是给定边界的温度定值,这里设定中‘even’表示的是均值,即整个边界条件的温度是不变的,然后‘up’等指的是方向,随后跟着的是具体的温度值,这里表示上方温度为4度,下方为2度,左边为5度,右边为3度。
同时隐式问题需要考虑到收敛问题,所以这里残差的计算方式是绝对值的均值变化,当小于10-6次方时停止运算,认为以及达到收敛。
然后是第二部,开始运算,这里输出三个变量,分别是最终的温度分布场result,以及temperature_variation是一个列表形式,中间每一个参数同样是列表形式代表的是每一步的运算过程,最后的residual是跟temperature_variation一样形式的残差变化过程,这里的100表示的是进行的一共的时间步长即100次🔺t,500表示的是单时间步长的残差收敛步数。
result,temperature_variation,residual=solver.Implicit_format_solver(solver,100,500)
下图所示的是第一个时间步长的残差变化收敛图,其实可以看出收敛得是非常快的。
![](https://i-blog.csdnimg.cn/blog_migrate/b9b182761f6f613a7bf1d826c6c7ef0a.png)
接下来得是最终得温度场结果,如下所示
[[5.0744004 5.70420753 6.16562093 6.6097002 7.1587784 ]
[4.90984364 5.63663855 6.25242643 6.83333912 7.42180063]
[5.00152196 5.81323045 6.46776248 7.03610846 7.55408465]
[5.38599562 6.37444432 6.96941806 7.38613703 7.71776144]
[6.48302421 7.76646494 7.86408545 7.87814045 7.82019123]]
整体类得代码如下所示
#Grid
import numpy as npclass
Grid:
def return_temp(self, temp):
return temp
def fit_rectangle(self,n_length,m_width,length,width):
grid=np.zeros((m_width,n_length))
self.grid=grid
self.n=n_length
self.m=m_width
self.length=length
self.width=width
self.detx=self.length/self.n
self.dety=self.width/self.m
return self
#Implicit_format_solver
import numpy as np
from Grid import Grid
def f_residual( grid1, grid2):
temp = abs(grid1 - grid2)
return np.mean(temp)
class Implicit_format_conduction:
def fit(self,grid,base):
self.grid = grid.grid
self.n = grid.n
self.m = grid.m
self.length = grid.length
self.width = grid.width
self.detx = grid.detx
self.dety = grid.dety
self.rou=base['rou']
self.c=base['c']
self.k=base['k']
self.dett=base['dett']
self.s=base['s']
def boundary_condition(self,model,seat,T):
if model=='even':
if seat=='up':
self.upgrid=np.zeros((self.m,1))
for i in range(self.m):
self.upgrid[i,0]=T
if seat=='down':
self.downgrid=np.zeros((self.m,1))
for i in range(self.m):
self.downgrid[i,0]=T
if seat=='left':
self.leftgrid=np.zeros((self.n,1))
for i in range(self.n):
self.leftgrid[i,0]=T
if seat=='right':
self.rightgrid=np.zeros((self.n,1))
for i in range(self.n):
self.rightgrid[i,0]=T
def Implicit_format_solver(self,Iteration,Iteration_steps):
#内部
self.ap1=self.rou*self.c*self.detx*self.dety/self.dett
self.ae0=self.k*self.dety/self.detx
self.aw0=self.k*self.dety/self.detx
self.an0=self.k*self.detx/self.dety
self.as0=self.k*self.detx/self.dety
self.ap0=self.rou*self.c*self.detx*self.dety/self.dett-self.ae0-self.aw0-self.an0-self.as0
self.list=[]
self.grid0=self.grid
self.residual=[]
for time in range(Iteration):
self.list_one = []
residual_list=[]
for Iteration_steps_0 in range(Iteration_steps):
if Iteration_steps_0>0:
if residual<0.000001:
break
self.list_one.append(self.grid)
newgrid=np.zeros((self.m,self.n))
for i in range(1,self.m-1):
for j in range(1,self.n-1):
newgrid[i,j]=self.ae0*self.grid[i+1,j]/self.ap1+self.aw0*self.grid[i-1,j]/self.ap1+\
self.an0*self.grid[i,j+1]/self.ap1+self.as0*self.grid[i,j-1]/self.ap1+self.ap0*self.grid0[i,j]/self.ap1+self.s*self.detx*self.dety/self.ap1
#边界
for i in range(1,self.m-1):
newgrid[i,self.n-1]=self.ae0*self.grid[i+1,self.n-1]/self.ap1+self.aw0*self.grid[i-1,self.n-1]/self.ap1+\
2*self.an0*self.upgrid[i,0]/self.ap1+self.as0*self.grid[i,self.n-2]/self.ap1+self.ap0*self.grid0[i,self.n-1]/self.ap1+self.s*self.detx*self.dety/self.ap1
#下边界
newgrid[i, 0] = self.ae0 * self.grid[i + 1, 0] / self.ap1 + self.aw0 * self.grid[i - 1, 0] / self.ap1 + \
self.an0 * self.grid[i, 1] / self.ap1 + 2*self.as0 * self.downgrid[i, 0] / self.ap1 +self.ap0*self.grid0[i,0]/self.ap1+ self.s * self.detx * self.dety / self.ap1
for j in range(1,self.n-1):
# 左边界
newgrid[0, j] = self.ae0 * self.grid[1, j] / self.ap1 + 2*self.aw0 * self.leftgrid[j,0] / self.ap1 + \
self.an0 * self.grid[0, j-1] / self.ap1 + self.as0 * self.grid[0, j+1] / self.ap1 +self.ap0*self.grid0[0,j]/self.ap1+ self.s * self.detx * self.dety / self.ap1
#右边界
newgrid[self.n-1, j] = 2*self.ae0 * self.rightgrid[j, 0] / self.ap1 + self.aw0 * self.grid[self.n-2, j] / self.ap1 + \
self.an0 * self.grid[self.n-1, j - 1] / self.ap1 + self.as0 * self.grid[self.n-1, j + 1] / self.ap1 +self.ap0*self.grid0[self.n-1,self.n-1]/self.ap1+ self.s * self.detx * self.dety / self.ap1
#左上
newgrid[0,self.n-1]=self.ae0 * self.grid[1, self.n-1] / self.ap1 + 2*self.aw0 * self.leftgrid[self.n-1,0] / self.ap1 + \
2*self.an0 * self.upgrid[0,0] / self.ap1 + self.as0 * self.grid[0, self.n-2] / self.ap1 +self.ap0*self.grid0[0,self.n-1]/self.ap1+ self.s * self.detx * self.dety / self.ap1
#左下
newgrid[0, 0] = self.ae0 * self.grid[1, 0] / self.ap1 + 2*self.aw0 * self.leftgrid[0, 0] / self.ap1 + \
self.an0 * self.grid[0, 1] / self.ap1 + 2*self.as0 * self.downgrid[0, 0] / self.ap1 +self.ap0*self.grid0[0, 0]/self.ap1+ self.s * self.detx * self.dety / self.ap1
#右上
newgrid[self.m-1,self.n-1]=2*self.ae0 * self.rightgrid[self.n-1,0] / self.ap1 + self.aw0 * self.grid[self.m-2,self.n-1] / self.ap1 + \
2*self.an0 * self.upgrid[self.m-1,0] / self.ap1 + self.as0 * self.grid[0, self.n-2] / self.ap1 +self.ap0*self.grid0[self.m-1,self.n-1]/self.ap1+ self.s * self.detx * self.dety / self.ap1
#右下
newgrid[self.m-1, 0] = 2*self.ae0 * self.rightgrid[0, 0] / self.ap1 + self.aw0 * self.grid[self.m-2, 0] / self.ap1 + \
self.an0 * self.grid[self.m-1, 1] / self.ap1 + 2*self.as0 * self.downgrid[self.m-1, 0] / self.ap1+self.ap0*self.grid0[self.m-1, 0]/self.ap1 + self.s * self.detx * self.dety / self.ap1
residual = f_residual(newgrid,self.grid)
residual_list.append(residual)
self.grid=newgrid
self.list.append(self.list_one)
self.residual.append(residual_list)
self.grid0=self.grid
return self.grid,self.list,self.residual