过程记录

于 2021-05-09 11:15:26 发布
阅读量212 收藏
点赞数
分类专栏： Python 文章标签： python
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Guilmots/article/details/116562335
版权

题目

在这里插入图片描述

显式格式

#边界条件是第一类和第二类最终
import pandas as pd
import numpy as np

#已知条件
length=3
height=2
TL=100;TD=50; #左边界和下边界
TR=5;TU=4; #右边界和上边界（测试）
den=100;cv=1000;k=10; #den密度，k导热系数
s=10;#内热源为10
dt=10;
qB=k*10#第二类边界条件热流密度
Tf=25
h=10

#迭代设置
n=10;m=10;#网格数
maxstep=200;
tr=0;#格式选择
miter=1;#迭代步数
maxiter=20;#最大迭代步
maxresi=1e-7;#最大残差

#划分网格
dx=length/n;
dy=height/m;


#准备变量
ae0=np.ones((n+2,m+2))#n代表y轴，m代表x轴，numpy默认转置
aw0=np.ones((n+2,m+2))
an0=np.ones((n+2,m+2))
as0=np.ones((n+2,m+2))
ap0=np.ones((n+2,m+2))
ap1=np.ones((n+2,m+2))
bp=np.ones((n+2,m+2))
T0=np.ones((n+2,m+2))
T=np.ones((n+2,m+2))
for i in range(1,n+1):
            for j in range(1,m+1): 
                T0[i][j]=3
#显示格式

Q=['time']
for i in range(1,n+1):
            for j in range(1,m+1): 
                Q.append("T"+str(i)+str(j))
                res = pd.DataFrame(columns=Q)
#迭代开始
for tstep in range(1,maxstep+1): 
#设定系数
    for i in range(2,n):
        for j in range(2,m):
            ae0[i][j]=k*dy/dx;
            aw0[i][j]=k*dy/dx;
            an0[i][j]=k*dx/dy;
            as0[i][j]=k*dx/dy;
            ap0[i][j]=den*cv*dx*dy/dt
            ap1[i][j]=ap0[i][j]+ae0[i][j]+aw0[i][j]+an0[i][j]+as0[i][j];
            bp[i][j]=s*dx*dy;
    
    for i in range(2,n):
#左面
        ae0[i][1]=k*dy/dx;
        aw0[i][1]=k*dy/(dx/2);
        an0[i][1]=k*dx/dy;
        as0[i][1]=k*dx/dy;
        ap0[i][1]=den*cv*dx*dy/dt
        ap1[i][1]=ap0[i][1]+ae0[i][1]+aw0[i][1]+an0[i][1]+as0[i][1];
        bp[i][1]=s*dx*dy;
                
#右面
        ae0[i][m]=0;
        aw0[i][m]=k*dy/dx;
        an0[i][m]=k*dx/dy;
        as0[i][m]=k*dx/dy;
        ap0[i][m]=den*cv*dx*dy/dt
        ap1[i][m]=ap0[i][m]+ae0[i][m]+aw0[i][m]+an0[i][m]+as0[i][m];
        bp[i][m]=s*dx*dy;
            
    for j in range(2,m):   
#上面
        ae0[1][j]=k*dy/dx;
        aw0[1][j]=k*dy/dx;
        an0[1][j]=0;
        as0[1][j]=k*dx/dy;
        ap0[1][j]=den*cv*dx*dy/dt
        ap1[1][j]=ap0[1][j]+ae0[1][j]+aw0[1][j]+an0[1][j]+as0[1][j]+dx/(1/h+dy/k);
        bp[1][j]=s*dx*dy;
                
#下面
        ae0[n][j]=k*dy/dx;
        aw0[n][j]=k*dy/dx;
        an0[n][j]=k*dx/dy;
        as0[n][j]=k*dx/(dy/2);
        ap0[n][j]=den*cv*dx*dy/dt
        ap1[n][j]=ap0[n][j]+ae0[n][j]+aw0[n][j]+an0[n][j]+as0[n][j];
        bp[n][j]=s*dx*dy;
                  
#左上角                    
    ae0[1][1]=k*dy/dx;
    aw0[1][1]=k*dy/(dx/2);
    an0[1][1]=0;
    as0[1][1]=k*dx/dy;
    ap0[1][1]=den*cv*dx*dy/dt
    ap1[1][1]=ap0[1][1]+ae0[1][1]+aw0[1][1]+an0[1][1]+as0[1][1]+dx/(1/h+dy/k);
    bp[1][1]=s*dx*dy;
                
#右上角
    ae0[1][m]=0;
    aw0[1][m]=k*dy/dx;
    an0[1][m]=0;
    as0[1][m]=k*dx/dy;
    ap0[1][m]=den*cv*dx*dy/dt
    ap1[1][m]=ap0[1][m]+ae0[1][m]+aw0[1][m]+an0[1][m]+as0[1][m]+dx/(1/h+dy/k);
    bp[1][m]=s*dx*dy;
#左下角
    ae0[n][1]=k*dy/dx;
    aw0[n][1]=k*dy/(dx/2);
    an0[n][1]=k*dx/dy;
    as0[n][1]=k*dx/(dy/2);
    ap0[n][1]=den*cv*dx*dy/dt
    ap1[n][1]=ap0[n][1]+ae0[n][1]+aw0[n][1]+an0[n][1]+as0[n][1];
    bp[n][1]=s*dx*dy;
                
#右下角
    ae0[n][m]=0;
    aw0[n][m]=k*dy/dx;
    an0[n][m]=k*dx/dy;
    as0[n][m]=k*dx/(dy/2);
    ap0[n][m]=den*cv*dx*dy/dt
    ap1[n][m]=ap0[n][m]+ae0[n][m]+aw0[n][m]+an0[n][m]+as0[n][m];
    bp[n][m]=s*dx*dy;
                
#边界条件设定
    for j in range(1,m+1): 
            T0[0][j]=TU
            T[0][j]=TU;
            T0[n+1][j]=TD; 
            T[n+1][j]=TD;
    for i in range(1,n+1): 
            T0[i][0]=TL;
            T[i][0]=TL;
            T0[i][n+1]=TR; 
            T[i][n+1]=TR;
#显式格式当前时间步计算
    for i in range(1,n+1):
            for j in range(1,m+1):
                T[i][j]=(ae0[i][j]*T0[i][j+1]+aw0[i][j]*T0[i][j-1]+an0[i][j]*T0[i-1][j]+as0[i][j]*T0[i+1][j]+ap0[i][j]*T0[i][j]+bp[i][j])/ap1[i][j];
    #第二类边界条件
    for i in range(1,n+1):
            T[i][m]=(ae0[i][j]*T0[i][j+1]+aw0[i][j]*T0[i][j-1]+an0[i][j]*T0[i-1][j]+as0[i][j]*T0[i+1][j]+ap0[i][j]*T0[i][j]+qB*dy+bp[i][j])/ap1[i][j]; 
    #第三类边界条件
    for j in range(1,m+1):
            T[1][j]=(ae0[i][j]*T0[i][j+1]+aw0[i][j]*T0[i][j-1]+an0[i][j]*T0[i-1][j]+as0[i][j]*T0[i+1][j]+ap0[i][j]*T0[i][j]+(Tf/(1/h+dy/k))*dx+bp[i][j])/ap1[i][j]; 
    for i in range(1,n+1):
            for j in range(1,m+1):
                T0[i][j]=T[i][j]     
##把矩阵写成一个横行    
    dict={}
    for i in range(1,n+1):
            for j in range(1,m+1):
                X='T'+str(i)+str(j)
                dict['time']=tstep*dt
                dict[X] = T[i][j] 
    res = res.append([dict],ignore_index=True)
#输出矩阵到excel
res.to_csv('显式输出21.csv') 
print('成功')

隐式格式

import pandas as pd
import numpy as np

#已知条件
length=3
height=2
TL=3.5;TD=4; #左边界和下边界
TR=5;TU=4; #右边界和上边界（测试）
den=100;cv=1000;k=10; #den密度，k导热系数
s=10;#内热源为10
dt=10;
qB=k*10#第二类边界条件热流密度
Tf=25#第三类边界条件
h=10

#迭代设置
n=10;m=10;#网格数
maxstep=200;#时间步
tr=0;#格式选择
iter=1;#迭代步数
maxiter=1000;#最大迭代步
maxresi=1e-7;#最大残差

#划分网格
dx=length/n;
dy=height/m;


#准备变量
ae1=np.ones((n+2,m+2))#n代表y轴，m代表x轴，numpy默认转置
aw1=np.ones((n+2,m+2))
an1=np.ones((n+2,m+2))
as1=np.ones((n+2,m+2))
ap0=np.ones((n+2,m+2))
ap1=np.ones((n+2,m+2))
bp=np.ones((n+2,m+2))
T0=np.ones((n+2,m+2))
Tn0=np.ones((n+2,m+2))#上一个迭代步的值
resi=np.ones((n+2,m+2))
T=np.ones((n+2,m+2))
for i in range(1,n+1):
            for j in range(1,m+1): 
                T0[i][j]=3
#显示格式

Q=['time']
for i in range(1,n+1):
            for j in range(1,m+1): 
                Q.append("T"+str(i)+str(j))
                res = pd.DataFrame(columns=Q)
for i in range(2,n):
    for j in range(2,m):
        ae1[i][j]=k*dy/dx;
        aw1[i][j]=k*dy/dx;
        an1[i][j]=k*dx/dy;
        as1[i][j]=k*dx/dy;
        ap0[i][j]=den*cv*dx*dy/dt
        ap1[i][j]=ap0[i][j]+ae1[i][j]+aw1[i][j]+an1[i][j]+as1[i][j];
        bp[i][j]=s*dx*dy;

for i in range(2,n):
#左面
    ae1[i][1]=k*dy/dx;
    aw1[i][1]=k*dy/(dx/2);
    an1[i][1]=k*dx/dy;
    as1[i][1]=k*dx/dy;
    ap0[i][1]=den*cv*dx*dy/dt
    ap1[i][1]=ap0[i][1]+ae1[i][1]+aw1[i][1]+an1[i][1]+as1[i][1];
    bp[i][1]=s*dx*dy;
            
#右面
    ae1[i][m]=0;
    aw1[i][m]=k*dy/dx;
    an1[i][m]=k*dx/dy;
    as1[i][m]=k*dx/dy;
    ap0[i][m]=den*cv*dx*dy/dt
    ap1[i][m]=ap0[i][m]+ae1[i][m]+aw1[i][m]+an1[i][m]+as1[i][m];
    bp[i][m]=s*dx*dy;
        
for j in range(2,m):   
#上面
    ae1[1][j]=k*dy/dx;
    aw1[1][j]=k*dy/dx;
    an1[1][j]=0;
    as1[1][j]=k*dx/dy;
    ap0[1][j]=den*cv*dx*dy/dt
    ap1[1][j]=ap0[1][j]+ae1[1][j]+aw1[1][j]+an1[1][j]+as1[1][j]+dx/(1/h+dy/k);
    bp[1][j]=s*dx*dy;
            
#下面
    ae1[n][j]=k*dy/dx;
    aw1[n][j]=k*dy/dx;
    an1[n][j]=k*dx/dy;
    as1[n][j]=k*dx/(dy/2);
    ap0[n][j]=den*cv*dx*dy/dt
    ap1[n][j]=ap0[n][j]+ae1[n][j]+aw1[n][j]+an1[n][j]+as1[n][j];
    bp[n][j]=s*dx*dy;
              
#左上角                    
ae1[1][1]=k*dy/dx;
aw1[1][1]=k*dy/(dx/2);
an1[1][1]=0;
as1[1][1]=k*dx/dy;
ap0[1][1]=den*cv*dx*dy/dt
ap1[1][1]=ap0[1][1]+ae1[1][1]+aw1[1][1]+an1[1][1]+as1[1][1]+dx/(1/h+dy/k);
bp[1][1]=s*dx*dy;
            
#右上角
ae1[1][m]=0;
aw1[1][m]=k*dy/dx;
an1[1][m]=0;
as1[1][m]=k*dx/dy;
ap0[1][m]=den*cv*dx*dy/dt
ap1[1][m]=ap0[1][m]+ae1[1][m]+aw1[1][m]+an1[1][m]+as1[1][m]+dx/(1/h+dy/k);
bp[1][m]=s*dx*dy;
#左下角
ae1[n][1]=k*dy/dx;
aw1[n][1]=k*dy/(dx/2);
an1[n][1]=k*dx/dy;
as1[n][1]=k*dx/(dy/2);
ap0[n][1]=den*cv*dx*dy/dt
ap1[n][1]=ap0[n][1]+ae1[n][1]+aw1[n][1]+an1[n][1]+as1[n][1];
bp[n][1]=s*dx*dy;
            
#右下角
ae1[n][m]=0;
aw1[n][m]=k*dy/dx;
an1[n][m]=k*dx/dy;
as1[n][m]=k*dx/(dy/2);
ap0[n][m]=den*cv*dx*dy/dt
ap1[n][m]=ap0[n][m]+ae1[n][m]+aw1[n][m]+an1[n][m]+as1[n][m];
bp[n][m]=s*dx*dy;
            
#边界条件设定
for j in range(1,m+1): 
        T[n+1][j]=TD;
        Tn0[n+1][j]=TD
for i in range(1,n+1): 
        T[i][0]=TL; 
        Tn0[i][0]=TL; 
#迭代开始
for tstep in range(1,maxstep+1): 
#设定系数

#隐式格式当前时间步计算
    for i in range(1,n+1):
        for j in range(1,m):
            Tn0[i][j]=T0[i][j]
#隐式格式迭代步计算            
    for iter in range(1,maxiter+1):
        #第一类边界条件
        for i in range(2,n+1):
            for j in range(1,m):
                T[i][j]=(ae1[i][j]*Tn0[i][j+1]+aw1[i][j]*Tn0[i][j-1]+an1[i][j]*Tn0[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+bp[i][j])/ap1[i][j];
        #第二类边界条件
            T[i][m]=(ae1[i][j]*Tn0[i][j+1]+aw1[i][j]*Tn0[i][j-1]+an1[i][j]*Tn0[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+qB*dy+bp[i][j])/ap1[i][j];     
       #第三类边界条件
        for j in range(1,m+1):
            T[1][j]=(ae1[i][j]*Tn0[i][j+1]+aw1[i][j]*Tn0[i][j-1]+an1[i][j]*Tn0[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+(Tf/(1/h+dy/k))*dx+bp[i][j])/ap1[i][j]; 
        #代入最新计算值继续迭代    
        for i in range(1,n+1):
                for j in range(1,m+1):
                    Tn0[i][j]=T[i][j]
#结束迭代，数值传给上一时间步    
    for i in range(1,n+1):
            for j in range(1,m+1):
                T0[i][j]=T[i][j]     
##把矩阵写成一个横行    
    dict={}
    for i in range(1,n+1):
            for j in range(1,m+1):
                X='T'+str(i)+str(j)
                dict['time']=tstep*dt
                dict[X] = T[i][j] 
    res = res.append([dict],ignore_index=True)
#输出矩阵到excel
res.to_csv('显式输出19.csv') 
print('成功')

G-S高斯-塞格尔迭代

import pandas as pd
import numpy as np

#已知条件
length=3
height=2
TL=3.5;TD=2; #左边界和下边界
TR=5;TU=4; #右边界和上边界（测试）
den=100;cv=1000;k=10; #den密度，k导热系数
s=10;#内热源为10
dt=20;
qB=k*10#第二类边界条件热流密度
Tf=25#第三类边界条件
h=10

#迭代设置
n=15;m=15;#网格数
maxstep=100;#时间步
tr=0;#格式选择
iter=1;#迭代步数
maxiter=200;#最大迭代步
maxresi=1e-7;#最大残差

#划分网格
dx=length/n;
dy=height/m;


#准备变量
ae1=np.ones((n+2,m+2))#n代表y轴，m代表x轴，numpy默认转置
aw1=np.ones((n+2,m+2))
an1=np.ones((n+2,m+2))
as1=np.ones((n+2,m+2))
ap0=np.ones((n+2,m+2))
ap1=np.ones((n+2,m+2))
bp=np.ones((n+2,m+2))
T0=np.ones((n+2,m+2))
Tn0=np.ones((n+2,m+2))#上一个迭代步的值
resi=np.ones((n+2,m+2))
T=np.ones((n+2,m+2))
for i in range(1,n+1):
            for j in range(1,m+1): 
                T0[i][j]=3
#显示格式

Q=['time']
for i in range(1,n+1):
            for j in range(1,m+1): 
                Q.append("T"+str(i)+str(j))
                res = pd.DataFrame(columns=Q)
#迭代开始
for tstep in range(1,maxstep+1): 
#设定系数
    for i in range(2,n):
        for j in range(2,m):
            ae1[i][j]=k*dy/dx;
            aw1[i][j]=k*dy/dx;
            an1[i][j]=k*dx/dy;
            as1[i][j]=k*dx/dy;
            ap0[i][j]=den*cv*dx*dy/dt
            ap1[i][j]=ap0[i][j]+ae1[i][j]+aw1[i][j]+an1[i][j]+as1[i][j];
            bp[i][j]=s*dx*dy;
    
    for i in range(2,n):
#左面
        ae1[i][1]=k*dy/dx;
        aw1[i][1]=k*dy/(dx/2);
        an1[i][1]=k*dx/dy;
        as1[i][1]=k*dx/dy;
        ap0[i][1]=den*cv*dx*dy/dt
        ap1[i][1]=ap0[i][1]+ae1[i][1]+aw1[i][1]+an1[i][1]+as1[i][1];
        bp[i][1]=s*dx*dy;
                
#右面
        ae1[i][m]=0;
        aw1[i][m]=k*dy/dx;
        an1[i][m]=k*dx/dy;
        as1[i][m]=k*dx/dy;
        ap0[i][m]=den*cv*dx*dy/dt
        ap1[i][m]=ap0[i][m]+ae1[i][m]+aw1[i][m]+an1[i][m]+as1[i][m];
        bp[i][m]=s*dx*dy;
            
    for j in range(2,m):   
#上面
        ae1[1][j]=k*dy/dx;
        aw1[1][j]=k*dy/dx;
        an1[1][j]=0;
        as1[1][j]=k*dx/dy;
        ap0[1][j]=den*cv*dx*dy/dt
        ap1[1][j]=ap0[1][j]+ae1[1][j]+aw1[1][j]+an1[1][j]+as1[1][j]+dx/(1/h+dy/k);
        bp[1][j]=s*dx*dy;
                
#下面
        ae1[n][j]=k*dy/dx;
        aw1[n][j]=k*dy/dx;
        an1[n][j]=k*dx/dy;
        as1[n][j]=k*dx/(dy/2);
        ap0[n][j]=den*cv*dx*dy/dt
        ap1[n][j]=ap0[n][j]+ae1[n][j]+aw1[n][j]+an1[n][j]+as1[n][j];
        bp[n][j]=s*dx*dy;
                  
#左上角                    
    ae1[1][1]=k*dy/dx;
    aw1[1][1]=k*dy/(dx/2);
    an1[1][1]=0;
    as1[1][1]=k*dx/dy;
    ap0[1][1]=den*cv*dx*dy/dt
    ap1[1][1]=ap0[1][1]+ae1[1][1]+aw1[1][1]+an1[1][1]+as1[1][1]+dx/(1/h+dy/k);
    bp[1][1]=s*dx*dy;
                
#右上角
    ae1[1][m]=0;
    aw1[1][m]=k*dy/dx;
    an1[1][m]=0;
    as1[1][m]=k*dx/dy;
    ap0[1][m]=den*cv*dx*dy/dt
    ap1[1][m]=ap0[1][m]+ae1[1][m]+aw1[1][m]+an1[1][m]+as1[1][m]+dx/(1/h+dy/k);
    bp[1][m]=s*dx*dy;
#左下角
    ae1[n][1]=k*dy/dx;
    aw1[n][1]=k*dy/(dx/2);
    an1[n][1]=k*dx/dy;
    as1[n][1]=k*dx/(dy/2);
    ap0[n][1]=den*cv*dx*dy/dt
    ap1[n][1]=ap0[n][1]+ae1[n][1]+aw1[n][1]+an1[n][1]+as1[n][1];
    bp[n][1]=s*dx*dy;
                
#右下角
    ae1[n][m]=0;
    aw1[n][m]=k*dy/dx;
    an1[n][m]=k*dx/dy;
    as1[n][m]=k*dx/(dy/2);
    ap0[n][m]=den*cv*dx*dy/dt
    ap1[n][m]=ap0[n][m]+ae1[n][m]+aw1[n][m]+an1[n][m]+as1[n][m];
    bp[n][m]=s*dx*dy;
                
#边界条件设定
    for j in range(1,m+1): 
            Tn0[0][j]=TU;
            Tn0[n+1][j]=TD;
            T[0][j]=TU;
            T[n+1][j]=TD;
    for i in range(1,n+1): 
            Tn0[i][0]=TL; 
            Tn0[i][n+1]=TR;
            T[i][0]=TL; 
            T[i][n+1]=TR;
#隐式格式当前时间步计算
    for i in range(1,n+1):
        for j in range(1,m):
            Tn0[i][j]=T0[i][j]
#隐式格式迭代步计算            
    for iter in range(1,maxiter+1):
        #第一类边界条件
        for i in range(2,n+1):
            for j in range(1,m):
                T[i][j]=(ae1[i][j]*Tn0[i][j+1]+aw1[i][j]*T[i][j-1]+an1[i][j]*T[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+bp[i][j])/ap1[i][j];
        #第二类边界条件
            T[i][m]=(ae1[i][j]*Tn0[i][j+1]+aw1[i][j]*T[i][j-1]+an1[i][j]*T[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+qB*dy+bp[i][j])/ap1[i][j];     
       #第三类边界条件
        for j in range(1,m+1):
            T[1][j]=(ae1[i][j]*Tn0[i][j+1]+aw1[i][j]*T[i][j-1]+an1[i][j]*T[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+(Tf/(1/h+dy/k))*dx+bp[i][j])/ap1[i][j]; 
        #代入最新计算值继续迭代    
        for i in range(1,n+1):
            for j in range(1,m+1):
                Tn0[i][j]=T[i][j]
#结束迭代，数值传给上一时间步    
    for i in range(1,n+1):
            for j in range(1,m+1):
                T0[i][j]=T[i][j]     
##把矩阵写成一个横行    
    dict={}
    for i in range(1,n+1):
            for j in range(1,m+1):
                X='T'+str(i)+str(j)
                dict['time']=tstep*dt
                dict[X] = T[i][j] 
#    res = res.append([dict],ignore_index=True)
#输出矩阵到excel
#res.to_csv('显式输出7.csv') 
print('成功')

ADI迭代（测试中）

import numpy as np

#已知条件
length=3
height=2
TL=3.5;TD=2; #左边界和下边界
TR=5;TU=4; #右边界和上边界（测试）
den=100;cv=1000;k=10; #den密度，k导热系数
s=10;#内热源为10
dt=200;
qB=k*10#第二类边界条件热流密度
Tf=25#第三类边界条件
h=10

#迭代设置
n=5;m=5;#网格数
maxstep=3;#时间步
tr=0;#格式选择
iter=1;#迭代步数
maxiter=2;#最大迭代步
maxresi=1e-7;#最大残差

#划分网格
dx=length/n;
dy=height/m;


#准备变量
ae1=np.ones((n+2,m+2))#n代表y轴，m代表x轴，numpy默认转置
aw1=np.ones((n+2,m+2))
an1=np.ones((n+2,m+2))
as1=np.ones((n+2,m+2))
ap0=np.ones((n+2,m+2))
ap1=np.ones((n+2,m+2))
bp=np.ones((n+2,m+2))
T0=np.ones((n+2,m+2))
Tn0=np.ones((n+2,m+2))#上一个迭代步的值
resi=np.ones((n+2,m+2))
T=np.ones((n+2,m+2))
for i in range(1,n+1):
            for j in range(1,m+1): 
                T0[i][j]=3
                Tn0[i][j]=T0[i][j]
a=np.ones((n+2))
b=np.ones((n+2))
c=np.ones((n+2))
d=np.ones((n+2))
x=np.ones((n+2))
p=np.ones((n+2))
q=np.ones((n+2))
#定义TDMA
def TDMA():    
    p[1]=-c[1]/a[1]
    q[1]=d[1]/a[1]
    for i in range(2,n+1):
        p[i]=-c[i]/(a[i]+b[i]*p[i-1]);
        q[i]=(d[i]-b[i]*q[i-1])/(a[i]+b[i]*p[i-1])
        x[n]=(d[n]-b[n]*q[n-1])/(a[n]+b[n]*p[n-1])
    for i in range(n-1,0,-1):
        x[i]=p[i]*x[i+1]+q[i]
 #设定系数
for i in range(2,n):
    for j in range(2,m):
        ae1[i][j]=k*dy/dx;
        aw1[i][j]=k*dy/dx;
        an1[i][j]=k*dx/dy;
        as1[i][j]=k*dx/dy;
        ap0[i][j]=den*cv*dx*dy/dt
        ap1[i][j]=ap0[i][j]+ae1[i][j]+aw1[i][j]+an1[i][j]+as1[i][j];
        bp[i][j]=s*dx*dy;

for i in range(2,n):
#左面
    ae1[i][1]=k*dy/dx;
    aw1[i][1]=k*dy/(dx/2);
    an1[i][1]=k*dx/dy;
    as1[i][1]=k*dx/dy;
    ap0[i][1]=den*cv*dx*dy/dt
    ap1[i][1]=ap0[i][1]+ae1[i][1]+aw1[i][1]+an1[i][1]+as1[i][1];
    bp[i][1]=s*dx*dy;

#右面
    ae1[i][m]=0;
    aw1[i][m]=k*dy/dx;
    an1[i][m]=k*dx/dy;
    as1[i][m]=k*dx/dy;
    ap0[i][m]=den*cv*dx*dy/dt
    ap1[i][m]=ap0[i][m]+ae1[i][m]+aw1[i][m]+an1[i][m]+as1[i][m];
    bp[i][m]=s*dx*dy;

for j in range(2,m):
#上面
    ae1[1][j]=k*dy/dx;
    aw1[1][j]=k*dy/dx;
    an1[1][j]=0;
    as1[1][j]=k*dx/dy;
    ap0[1][j]=den*cv*dx*dy/dt
    ap1[1][j]=ap0[1][j]+ae1[1][j]+aw1[1][j]+an1[1][j]+as1[1][j]+dx/(1/h+dy/k);
    bp[1][j]=s*dx*dy;

#下面
    ae1[n][j]=k*dy/dx;
    aw1[n][j]=k*dy/dx;
    an1[n][j]=k*dx/dy;
    as1[n][j]=k*dx/(dy/2);
    ap0[n][j]=den*cv*dx*dy/dt
    ap1[n][j]=ap0[n][j]+ae1[n][j]+aw1[n][j]+an1[n][j]+as1[n][j];
    bp[n][j]=s*dx*dy;

#左上角                    
    ae1[1][1]=k*dy/dx;
    aw1[1][1]=k*dy/(dx/2);
    an1[1][1]=0;
    as1[1][1]=k*dx/dy;
    ap0[1][1]=den*cv*dx*dy/dt
    ap1[1][1]=ap0[1][1]+ae1[1][1]+aw1[1][1]+an1[1][1]+as1[1][1]+dx/(1/h+dy/k);
    bp[1][1]=s*dx*dy;

#右上角
    ae1[1][m]=0;
    aw1[1][m]=k*dy/dx;
    an1[1][m]=0;
    as1[1][m]=k*dx/dy;
    ap0[1][m]=den*cv*dx*dy/dt
    ap1[1][m]=ap0[1][m]+ae1[1][m]+aw1[1][m]+an1[1][m]+as1[1][m]+dx/(1/h+dy/k);
    bp[1][m]=s*dx*dy;
#左下角
    ae1[n][1]=k*dy/dx;
    aw1[n][1]=k*dy/(dx/2);
    an1[n][1]=k*dx/dy;
    as1[n][1]=k*dx/(dy/2);
    ap0[n][1]=den*cv*dx*dy/dt
    ap1[n][1]=ap0[n][1]+ae1[n][1]+aw1[n][1]+an1[n][1]+as1[n][1];
    bp[n][1]=s*dx*dy;

#右下角
    ae1[n][m]=0;
    aw1[n][m]=k*dy/dx;
    an1[n][m]=k*dx/dy;
    as1[n][m]=k*dx/(dy/2);
    ap0[n][m]=den*cv*dx*dy/dt
    ap1[n][m]=ap0[n][m]+ae1[n][m]+aw1[n][m]+an1[n][m]+as1[n][m];
    bp[n][m]=s*dx*dy;

#边界条件设定
for j in range(1,m+1): 
        Tn0[0][j]=TU;
        Tn0[n+1][j]=TD;
        T[0][j]=TU;
        T[n+1][j]=TD;
for i in range(1,n+1): 
        Tn0[i][0]=TL; 
        Tn0[i][n+1]=TR;
        T[i][0]=TL; 
        T[i][n+1]=TR;
#迭代开始
for iter in range(1,maxiter+1):                  
#设定系数
    for i in range(2,n):
        for j in range(2,m):
            ae1[i][j]=k*dy/dx;
            aw1[i][j]=k*dy/dx;
            an1[i][j]=k*dx/dy;
            as1[i][j]=k*dx/dy;
            ap0[i][j]=den*cv*dx*dy/dt
            ap1[i][j]=ap0[i][j]+ae1[i][j]+aw1[i][j]+an1[i][j]+as1[i][j];
            bp[i][j]=s*dx*dy;
    
    for i in range(2,n):
    #左面
        ae1[i][1]=k*dy/dx;
        aw1[i][1]=k*dy/(dx/2);
        an1[i][1]=k*dx/dy;
        as1[i][1]=k*dx/dy;
        ap0[i][1]=den*cv*dx*dy/dt
        ap1[i][1]=ap0[i][1]+ae1[i][1]+aw1[i][1]+an1[i][1]+as1[i][1];
        bp[i][1]=s*dx*dy;
    
    #右面
        ae1[i][m]=0;
        aw1[i][m]=k*dy/dx;
        an1[i][m]=k*dx/dy;
        as1[i][m]=k*dx/dy;
        ap0[i][m]=den*cv*dx*dy/dt
        ap1[i][m]=ap0[i][m]+ae1[i][m]+aw1[i][m]+an1[i][m]+as1[i][m];
        bp[i][m]=s*dx*dy;
    
    for j in range(2,m):
    #上面
        ae1[1][j]=k*dy/dx;
        aw1[1][j]=k*dy/dx;
        an1[1][j]=0;
        as1[1][j]=k*dx/dy;
        ap0[1][j]=den*cv*dx*dy/dt
        ap1[1][j]=ap0[1][j]+ae1[1][j]+aw1[1][j]+an1[1][j]+as1[1][j]+dx/(1/h+dy/k);
        bp[1][j]=s*dx*dy;
    
    #下面
        ae1[n][j]=k*dy/dx;
        aw1[n][j]=k*dy/dx;
        an1[n][j]=k*dx/dy;
        as1[n][j]=k*dx/(dy/2);
        ap0[n][j]=den*cv*dx*dy/dt
        ap1[n][j]=ap0[n][j]+ae1[n][j]+aw1[n][j]+an1[n][j]+as1[n][j];
        bp[n][j]=s*dx*dy;
    
    #左上角                    
        ae1[1][1]=k*dy/dx;
        aw1[1][1]=k*dy/(dx/2);
        an1[1][1]=0;
        as1[1][1]=k*dx/dy;
        ap0[1][1]=den*cv*dx*dy/dt
        ap1[1][1]=ap0[1][1]+ae1[1][1]+aw1[1][1]+an1[1][1]+as1[1][1]+dx/(1/h+dy/k);
        bp[1][1]=s*dx*dy;
    
    #右上角
        ae1[1][m]=0;
        aw1[1][m]=k*dy/dx;
        an1[1][m]=0;
        as1[1][m]=k*dx/dy;
        ap0[1][m]=den*cv*dx*dy/dt
        ap1[1][m]=ap0[1][m]+ae1[1][m]+aw1[1][m]+an1[1][m]+as1[1][m]+dx/(1/h+dy/k);
        bp[1][m]=s*dx*dy;
    #左下角
        ae1[n][1]=k*dy/dx;
        aw1[n][1]=k*dy/(dx/2);
        an1[n][1]=k*dx/dy;
        as1[n][1]=k*dx/(dy/2);
        ap0[n][1]=den*cv*dx*dy/dt
        ap1[n][1]=ap0[n][1]+ae1[n][1]+aw1[n][1]+an1[n][1]+as1[n][1];
        bp[n][1]=s*dx*dy;
    
    #右下角
        ae1[n][m]=0;
        aw1[n][m]=k*dy/dx;
        an1[n][m]=k*dx/dy;
        as1[n][m]=k*dx/(dy/2);
        ap0[n][m]=den*cv*dx*dy/dt
        ap1[n][m]=ap0[n][m]+ae1[n][m]+aw1[n][m]+an1[n][m]+as1[n][m];
        bp[n][m]=s*dx*dy;
    
    #边界条件设定
    for j in range(1,m+1): 
            Tn0[0][j]=TU;
            Tn0[n+1][j]=TD;
            T[0][j]=TU;
            T[n+1][j]=TD;
    for i in range(1,n+1): 
            Tn0[i][0]=TL; 
            Tn0[i][n+1]=TR;
            T[i][0]=TL; 
            T[i][n+1]=TR;
    #x方向线迭代
    for i in range(1,n+1):
        if j==1:
            a[i]=ap1[i][j]
            c[i]=-ae1[i][j]
            d[i]=an1[i][j]*Tn0[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+bp[i][j]+aw1[i][j]*T[i][j-1]

        if j in range(2,m):
            a[i]=ap1[i][j]
            b[i]=-aw1[i][j]
            c[i]=-ae1[i][j]
            d[i]=an1[i][j]*Tn0[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+bp[i][j]

        if j==m:
            a[i]=ap1[i][j]
            b[i]=-aw1[i][j]
            d[i]=an1[i][j]*Tn0[i-1][j]+as1[i][j]*Tn0[i+1][j]+ap0[i][j]*T0[i][j]+bp[i][j]+ae1[i][j]*T[i][j+1]+qB*dy
        TDMA()
        for j in range(1,m+1):
            T[i][j]=x[i]
        print("第"+str(i)+"列线迭代：")
        print(ap1)
    print("第"+str(iter)+"次迭代步")
EULE
关注
专栏目录
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值