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

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

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

由于是二维问题,所以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)

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

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

[[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
  • 1
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
隐式方法是求解一维热传导方程的一种数值方法,它与显式方法相比更加稳定,但计算量较大。下面是使用隐式方法求解一维热传导方程的示例代码: ```python import numpy as np import matplotlib.pyplot as plt # 定义参数 L = 1 # 物体长度 T = 1 # 总时间 N = 100 # 空间步数 M = 1000 # 时间步数 alpha = 0.01 # 热扩散系数 # 计算步长 dx = L / N dt = T / M # 初始化温度矩阵 u = np.zeros((N+1, M+1)) # 设置初始条件 u[:, 0] = 0 # 初始温度为0 u[0, :] = 1 # 左边界温度为1 # 构建系数矩阵 A = np.zeros((N-1, N-1)) for i in range(N-1): A[i, i] = 1 + 2 * alpha * dt / dx**2 if i > 0: A[i, i-1] = -alpha * dt / dx**2 if i < N-2: A[i, i+1] = -alpha * dt / dx**2 # 迭代求解 for j in range(0, M): b = u[1:N, j] u[1:N, j+1] = np.linalg.solve(A, b) # 绘制温度分布曲线 x = np.linspace(0, L, N+1) t = np.linspace(0, T, M+1) X, T = np.meshgrid(x, t) fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_surface(X, T, u.T, cmap='hot') ax.set_xlabel('Position') ax.set_ylabel('Time') ax.set_zlabel('Temperature') plt.show() ``` 以上示例代码使用了隐式方法来近似求解一维热传导方程,通过构建系数矩阵A和迭代计算,得到了温度在空间和时间上的分布,并使用matplotlib库绘制了温度分布曲线。隐式方法中的求解步骤使用了线性方程组的求解函数`np.linalg.solve`,通过解线性方程组来计算下一个时间步的温度值。请注意,隐式方法计算量较大,但稳定性更好,适合处理较大的时间步长和热扩散系数。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值