python环境下的CFD运算（二维隐式导热问题，已写成类，可轻松调动）

最新推荐文章于 2024-03-22 09:48:59 发布

wwy2014wwy

最新推荐文章于 2024-03-22 09:48:59 发布

阅读量1.1k

点赞数 1

分类专栏： python环境下得CFD 文章标签： python 经验分享课程设计 Powered by 金山文档

本文链接：https://blog.csdn.net/wwy2014wwy/article/details/129368771

版权

python环境下得CFD 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

原理：由物体内能方程推导，原理是物体内能的变化=进出物体的内能+物体本身的能量。如下所示

$\text{[math]}$

其中，方程左侧代表的是物体内能的变化，右侧代表的是六个方向的进出能量，与物体本身的内能变化

由于是二维问题，所以Z方向上的B与F皆为0，所以可以化解为0，其中导热项由傅里叶定律可以变换，故最终推到式如下所示

$\text{[math]}$

然后进行合并同类项，可以将方程化解成以下格式，其中四个方向的1时刻的温度是我们所不知道，也就是我们所有计算的最终的结果，所以这里采用迭代的方法进行推算，当然当四个方向的温度采用原时刻也就是0时刻的时候，就变成显示求解了。

$\text{[math]}$

所以实际的工作分为两步，第一步进行网格划分，这里该问题仅针对矩形网格，调用的方法在如下

 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)

下图所示的是第一个时间步长的残差变化收敛图，其实可以看出收敛得是非常快的。

接下来得是最终得温度场结果，如下所示

[[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