二维稳态对流扩散问题，无源项，QUICK格式，四周边界定壁温的Python程序

二维稳态对流扩散，QUICK格式

原理

模型

四面为定壁温

例子

初始情况是物体温度为0℃，西面表面温度突然变成60℃，北面表面温度突然变成600℃，温度保持不变，东面和南面表面温度恒定为0℃。水平速度为2.8m/s，方向由西向东，竖直方向速度为2.6m/s，方向由北向南。扩散系数为0.1，密度为1。求物体内部温度分布。

代码

#二维稳态对流扩散问题
#无源项，QUICK格式
from numpy import *
from numpy.linalg import *
import pandas as pd
X = 2.8  #物体x方向总长度
Y = 3.4  #物体y方向总长度
u = 2.8  #x方向速度
v = 2.6  #y方向速度
n = 40  #x方向划分网格数
m = 40  #y方向划分网格数
I = 0.1  #扩散系数
p = 1  #密度

#边界条件
tA = 60  #左边界A点，定壁温
tB = 0  #右边界B点，定壁温
tC = 600  #上边界C点，定壁温
tD = 0  #上边界D点，定壁温

dx = X/n  #x方向控制体宽度
dy = Y/m  #y方向控制体宽度
De, Dw, Ds, Dn = I/dx, I/dx, I/dy, I/dy
DA, DB, DC, DD = 2*I/dx, 2*I/dx, 2*I/dy, 2*I/dy
Fe, Fw, FA, FB = p*u, p*u, p*u, p*u
Fs, Fn, FC, FD = p*v, p*v, p*v, p*v

Pe = max(Fe/De, Fn/Dn, FA/DA, FC/DC)
if Pe > 8/3:
    print(f'Pe={Pe}>8/3,解不稳定')
else:
    print('解是稳定的')
    A = zeros((m*n,m*n)) #创建零数组
    B = zeros((m*n,1))
    #节点i=1, j=1
    A[0][0] = (7/8*Fe+De+9/6*DA)*dy+(7/8*Fs+Ds+9/6*DC)*dx
    A[0][1] = (3/8*Fe-1/6*DA-De)*dy
    A[0][1*n] = (3/8*Fs-1/6*DC-Ds)*dx
    B[0][0] = (2/8*Fe+FA+8/6*DA)*dy*tA + (2/8*Fs+FC+8/6*DC)*dx*tC
    #节点i=1, j=2
    A[1][0] = -(1/8*Fe+7/8*Fw+Dw)*dy
    A[1][1] = (6/8*Fe-3/8*Fw+De+Dw)*dy+(7/8*Fs+Ds+9/6*DC)*dx
    A[1][2] = (3/8*Fe-De)*dy
    A[1][(2-1)*n+1] = (3/8*Fs-Ds-1/6*DC)*dx
    B[1][0] = -2/8*Fw*dy*tA + (2/8*Fs+FC+8/6*DC)*dx*tC
    #节点i=1, 3≤j＜n
    for j in range(3,n):
        A[j-1][j-3] = 1/8*Fw*dy
        A[j-1][j-2] = -(1/8*Fe+6/8*Fw+Dw)*dy
        A[j-1][j-1] = (6/8*Fe-3/8*Fw+De+Dw)*dy+(7/8*Fs+Ds+9/6*DC)*dx
        A[j-1][j] = (3/8*Fe-De)*dy
        A[j-1][(2-1)*n+j-1] = (3/8*Fs-Ds-1/6*DC)*dx
        B[j-1][0] = (8/6*DC+2/8*Fs+FC)*dx*tC
    #节点i=1, j=n
    A[n-1][n-3] = 1/8*Fw*dy
    A[n-1][n-2] = -(6/8*Fw+1/6*DB+Dw)*dy
    A[n-1][n-1] = (-3/8*Fw+Dw+9/6*DB)*dy+(7/8*Fs+Ds+9/6*DC)*dx
    A[n-1][(2-1)*n+n-1] = (3/8*Fs-Ds-1/6*DC)*dx
    B[n-1][0] = (8/6*DB-FB)*dy*tB + (2/8*Fs+FC+8/6*DC)*dx*tC
    #节点i=2, j=1
    A[n][0] = -(1/8*Fs+7/8*Fn+Dn)*dx
    A[n][(2-1)*n] = (7/8*Fe+De+9/6*DA)*dy+(6/8*Fs-3/8*Fn+Ds+Dn)*dx
    A[n][(2-1)*n+1] = (3/8*Fe-De-1/6*DA)*dy
    A[n][(3-1)*n] = (3/8*Fs-Ds)*dx
    B[n][0] = (2/8*Fe+FA+8/6*DA)*dy*tA - 2/8*Fn*dx*tC
    #节点i=2, j=2
    A[n+1][1] = -(7/8*Fn+Dn+1/8*Fs)*dx
    A[n+1][(2-1)*n] = -(1/8*Fe+7/8*Fw+Dw)*dy
    A[n+1][(2-1)*n+1] = (6/8*Fe-3/8*Fw+De+Dw)*dy+(6/8*Fs-3/8*Fn+Ds+Dn)*dx
    A[n+1][(2-1)*n+2] = (3/8*Fe-De)*dy
    A[n+1][(3-1)*n+1] = (3/8*Fs-Ds)*dx
    B[n+1][0] = -2/8*Fw*dy*tA - 2/8*