精华内容
下载资源
问答
  • 汽车轮胎二维稳态温度场的数值分析.pdf
  • 二维温度场matlab分析

    2018-10-23 14:06:11
    基于matlab的二维钢板冷却的温度场分析,运用差分方程
  • 教育精品资料
  • 二维稳态导热的数值计算主要采用了热平衡法。用差分法建立节点的热平衡方程, 将节 点所在的单元体的四个方向传递的热流密度,内热源在单元体产生的热流密度, 根据能量守 恒的原则建立方程,可以得到每一个节点...
  • 在MATLAB环境下的c语言编写的温度场模型,采用先进的的温度模型算法,全面真实 在MATLAB环境下的c语言编写的温度场模型,采用先进的的温度模型算法,全面真实
  • 二维稳态导热程序

    2013-04-01 23:59:31
    fortran90编写的二维常物性 钢坯加热通用程序
  • 数值传热的二维稳态导热实验。根据确定的边界条件,利用C语言编写程序对研究区域进行网格划分,获取各网格点的温度值,再利用MATLAB进行数值模拟,确定温度场分布云图
  • 根据能量守恒定律和傅立叶定律,采用有限体积法对二维稳态平面温度场进行求解,合理地确定控制体的大小,建立了平面稳态温度场问题的线性方程组及其离散模型。结果表现了平面温度场三角形单元变分计算公式的实际物理...
  • 用有限差分法和 Matlab 计算二维热加工温度场分析 Eg1 薄板焊接过程温度场分析 取焊件的一半为模型进行离散化起始点为 o 点以后以 v 速度沿 x 轴运动 根据题意为二维稳态导热二维稳态导热方程为 2 2 1 T T T Q 2...
  • 和一稳态扩散算例一样的初始和边界条件 ,采用乘方格式 时间步长为0.001s,初始温度场为200,速度为2m/s,长度为2cm,t=0s时刻东侧温度突然降至0C。时间差分采用全隐式格式
  • 二维稳态导热

    千次阅读 2020-03-17 18:27:41
    计算求解在有内热源和无内热源情况下,固体材料分别为铜体和保温材料(λ=0.54+0.00058{t}℃\lambda=0.54+0.00058 \left\{t\right\}℃λ=0.54+0....⋅K)h=5\sim 500W/(m^{2}\cdot K)h=5∼500W/(m2⋅K)情况下的温度分布...

    计算题目

    计算求解在有内热源和无内热源情况下,固体材料分别为铜体和保温材料( λ = 0.54 + 0.00058 { t } ℃ \lambda=0.54+0.00058 \left\{t\right\}℃ λ=0.54+0.00058{t})时的,对流换热系数从 h = 5 ∼ 500 W / ( m 2 ⋅ K ) h=5\sim 500W/(m^{2}\cdot K) h=5500W/(m2K)情况下的温度分布。

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Uv8DxZfk-1584440806238)(geo.png)]{#img width=".4\textwidth"}

    数学物理模型

    物理模型

    二维有内热源的稳态导热模型。

    数学模型

    如图1所示建立平面直角坐标系,上述物理模型可用下述数学模型表
    ∂ ∂ x ( λ ∂ T ∂ x ) + ∂ ∂ y ( λ ∂ T ∂ y ) + S = 0 \frac{\partial}{\partial x}(\lambda \frac{\partial T}{\partial x})+\frac{\partial}{\partial y}(\lambda \frac{\partial T}{\partial y})+S=0 x(λxT)+y(λyT)+S=0
    物理条件: λ = λ ( T ) \lambda = \lambda(T) λ=λ(T) 边界条件: KaTeX parse error: Undefined control sequence: \notag at position 148: …{x=1.5}-T_{f}) \̲n̲o̲t̲a̲g̲ ̲\\ & \left.…

    计算区域及方程离散

    区域离散

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-69lqkHFq-1584440806243)(node.png)]{#img width=".4\textwidth"}

    对题示二维导热问题物理区域采用内点法进行网格离散,如图2所示。

    方程离散

    使用控制体积积分法离散方程

    ∂ ∂ x ( λ ∂ T ∂ x ) + ∂ ∂ y ( λ ∂ T ∂ y ) + S = 0 \frac{\partial}{\partial x}(\lambda\frac{\partial T}{\partial x})+\frac{\partial}{\partial y}(\lambda\frac{\partial T}{\partial y}) + S=0 x(λxT)+y(λyT)+S=0
    将方程在控制体积内积分得
    ( λ e T E − T P ( δ x ) e − λ w T P − T W ( δ x ) w ) Δ y + ( λ n T N − T P ( δ y ) n − λ s T P − T S ( δ y ) s ) Δ x + S Δ x Δ y (\lambda_{e}\frac{T_{E}-T_{P}}{(\delta x)_{e}}-\lambda_{w}\frac{T_{P}-T_{W}}{(\delta x)_{w}})\Delta y + (\lambda_{n}\frac{T_{N}-T_{P}}{(\delta y)_{n}}-\lambda_{s}\frac{T_{P}-T_{S}}{(\delta y)_{s}})\Delta x + S\Delta x\Delta y (λe(δx)eTETPλw(δx)wTPTW)Δy+(λn(δy)nTNTPλs(δy)sTPTS)Δx+SΔxΔy
    通过源项处理并整理,最终可以得到 KaTeX parse error: Undefined control sequence: \notag at position 229: …\Delta y)T_{P} \̲n̲o̲t̲a̲g̲ ̲\\ =\frac{\…

    内部节点

    a P T P = a E T E + a W T W + a N T N + a S T S + b a_{P}T_{P}=a_{E}T_{E}+a_{W}T_{W}+a_{N}T_{N}+a_{S}T_{S}+b aPTP=aETE+aWTW+aNTN+aSTS+b 其中
    KaTeX parse error: Undefined control sequence: \notag at position 127: …w}/\lambda_{w}}\̲n̲o̲t̲a̲g̲ ̲\\ a_{N}=\…
    KaTeX parse error: Undefined control sequence: \notag at position 76: …lta x \Delta y \̲n̲o̲t̲a̲g̲ ̲\\ b=S_{c}…

    边界节点

    对于上边界节点 KaTeX parse error: Undefined control sequence: \notag at position 72: …{\partial y}=10\̲n̲o̲t̲a̲g̲ ̲\\ a_{P}^{…
    其中 q B = 10 ,      a N = 0    ( λ B = 0 ) q_{B}=10,~~~~a_{N}=0~~(\lambda_{B}=0) qB=10,    aN=0  (λB=0) 对于下边界节点
    KaTeX parse error: Undefined control sequence: \notag at position 66: …partial y} = 0 \̲n̲o̲t̲a̲g̲ ̲\\ a_{P}^{…
    其中 q B = 0 ,      a S = 0    ( λ B = 0 ) q_{B}=0,~~~~a_{S}=0~~(\lambda_{B}=0) qB=0,    aS=0  (λB=0) 对于左边界节点
    KaTeX parse error: Undefined control sequence: \notag at position 41: …0 ~~~~~&T=300K \̲n̲o̲t̲a̲g̲ ̲\\ a_{P}^{…
    对于右边界节点 KaTeX parse error: Undefined control sequence: \notag at position 85: …l x}=h(T-T_{f})\̲n̲o̲t̲a̲g̲ ̲\\ a_{P}^{…
    其中
    a E = 0    ( λ B = 0 ) ,      q B = T f − T P 1 h + ( δ x ) B λ B a_{E}=0~~(\lambda_{B}=0),~~~~q_{B}=\frac{T_{f}-T_{P}}{\frac{1}{h}+\frac{(\delta x)_{B}}{\lambda_{B}}} aE=0  (λB=0),    qB=h1+λB(δx)BTfTP
    整理得
    ( a P + S P , a d Δ V ) T P = a W T W + a N T N + a S T S + b + S C , a d Δ V (a_{P}+S_{P,ad}\Delta V)T_{P}=a_{W}T_{W}+a_{N}T_{N}+a_{S}T_{S}+b+S_{C,ad}\Delta V (aP+SP,adΔV)TP=aWTW+aNTN+aSTS+b+SC,adΔV
    其中 KaTeX parse error: Undefined control sequence: \notag at position 95: … \lambda_{B}]} \̲n̲o̲t̲a̲g̲ ̲\\ S_{C,ad…

    数值方法及程序流程

    交叉方向线迭代

    本文使用交叉方向线迭代法来计算二维情况下的稳态导热问题。该算法意旨通过 x x x方向与 y y y方向交叉计算三对角矩阵
    从图3可以看出,该算法诣在将二维导热问题的五对角矩阵转化为三对角矩阵求解。

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-S33DF4vn-1584440806245)(mesh_0.png)]{#img width=".4\textwidth"}

    由图可知扫描 x x x轴方向时,计算每一列的值,此时的 N N N S S S两点使用显示格式,故每一列都是一个三对角矩阵,以此重复计算所有列的值。然后再扫描 y y y方向

    扫描方向{#img width=".4\textwidth"}

    由图4可以看出在扫描 y y y轴方向时,需要计算每一行的值,此时的 W W W E E E两点使用显示格式,故每一行都是一个三对角矩阵,以此重复计算所有行的值,然后进行迭代。

    扫描方向{#img width=".4\textwidth"}

    有源项的二维热传导方程的代数形式

    若方程有源项,例如 S = 0.001 T 2 S=0.001T^{2} S=0.001T2,通过泰勒展开处理源项可得到
    S = 0.001 T 2 = 0.001 ( 2 T ∗ 2 − T 2 = 0.004 T ∗ 2 − 2 T ∗ T ) S = 0.001T^{2} = 0.001(2T^{*2}-T^{2} = 0.004T^{*2}-2T^{*}T) S=0.001T2=0.001(2T2T2=0.004T22TT)
    可以得到此时内部节点的代数方程 KaTeX parse error: Undefined control sequence: \notag at position 233: …\Delta y)T_{P} \̲n̲o̲t̲a̲g̲ ̲\\ =\frac{\…
    根据左边界条件,得到左边界节点的代数方程
    T P = ( a E T E + a W T W + a N T N + a S T S + 0.004 T ∗ 2 Δ x Δ y ) / a P T_{P} = (a_{E}T_{E}+a_{W}T_{W}+a_{N}T_{N}+a_{S}T_{S} + 0.004T^{*2}\Delta x\Delta y ) / a_{P} TP=(aETE+aWTW+aNTN+aSTS+0.004T2ΔxΔy)/aP
    根据右边界条件,得到右边界节点的代数方程
    KaTeX parse error: Undefined control sequence: \/ at position 133: … (\Delta x)_{B}\̲/̲lambda_{B}]}\De…
    根据上边界条件,得到上边界节点的代数方程
    T P = ( a W T W + a E T E + a N T N + a S T S + 10 Δ x Δ V Δ V + 0.004 T ∗ 2 Δ x Δ y ) / a P T_{P} = (a_{W}T_{W}+a_{E}T_{E}+a_{N}T_{N}+a_{S}T_{S}+\frac{10\Delta x}{\Delta V}\Delta V + 0.004T^{*2}\Delta x\Delta y)/a_{P} TP=(aWTW+aETE+aNTN+aSTS+ΔV10ΔxΔV+0.004T2ΔxΔy)/aP
    根据下边界条件,得到下边界节点的代数方程
    T P = ( a W T W + a E T E + a N T N + 0.004 T ∗ 2 Δ x Δ y ) / a P T_{P} = (a_{W}T_{W}+a_{E}T_{E}+a_{N}T_{N} + 0.004T^{*2}\Delta x\Delta y) / a_{P} TP=(aWTW+aETE+aNTN+0.004T2ΔxΔy)/aP

    无源项的二维热传导方程的代数形式

    可以得到此时内部节点的代数方程 KaTeX parse error: Undefined control sequence: \notag at position 205: …bda_{s}})T_{P} \̲n̲o̲t̲a̲g̲ ̲\\ =\frac{\…
    根据左边界条件,得到左边界节点的代数方程
    T P = ( a E T E + a W T W + a N T N + a S T S ) / a P T_{P} = (a_{E}T_{E}+a_{W}T_{W}+a_{N}T_{N}+a_{S}T_{S}) / a_{P} TP=(aETE+aWTW+aNTN+aSTS)/aP
    根据右边界条件,得到右边界节点的代数方程
    KaTeX parse error: Undefined control sequence: \/ at position 103: … (\Delta x)_{B}\̲/̲lambda_{B}]}\De…
    根据上边界条件,得到上边界节点的代数方程
    T P = ( a W T W + a E T E + a N T N + a S T S + 10 Δ x Δ V Δ V ) / a P T_{P} = (a_{W}T_{W}+a_{E}T_{E}+a_{N}T_{N}+a_{S}T_{S}+\frac{10\Delta x}{\Delta V}\Delta V)/a_{P} TP=(aWTW+aETE+aNTN+aSTS+ΔV10ΔxΔV)/aP
    根据下边界条件,得到下边界节点的代数方程
    T P = ( a W T W + a E T E + a N T N ) / a P T_{P} = (a_{W}T_{W}+a_{E}T_{E}+a_{N}T_{N}) / a_{P} TP=(aWTW+aETE+aNTN)/aP

    计算结果及网格独立性考核

    网格独立性考核

    网格数量对数值计算结果有着重要的影响。当网格足够细密以至于再进一步加密网格已对数值计算结果基本上没有影响时所得到的数值解称为网格独立解。
    本文通过保持几何形状的长宽比,即1.5:1的比例将网络进行缩放,分别将网格数理增大
    2, 4, 6, 8, 10, 12, 14, 16, 18, 20倍来分析网络独立解。

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5nQEIrpi-1584440806246)(indepn.png)]{#img width=".4\textwidth"}

    由图可得,在第四个点也就是将网格等比例放大8倍,即12x8的网格此后的解变化不明显,故认为此时的解为网格独立解。

    二维铜体导热温度分布

    材料为铜体时导热系数 λ = 400 W / m ⋅ K \lambda = 400 W/m\cdot K λ=400W/mK

    有内热源源二维铜体导热温度分布

    当铜体内存在内热源 S = 0.001 T 2 S=0.001T^{2} S=0.001T2时的,不同对流换热系数下温度分布

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rZmz4GFi-1584440806246)(cu.png)]{#img width=“15cm”}

    无内热源二维铜体导热温度分布

    无内热源时,不同对流换热系数下的温度分布

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0xuFsBsT-1584440806248)(cu_s.png)]{#img width=“15cm”}

    二维保温材料导热温度分布

    当材料为保温材料时,其导热系数 λ = 0.54 + 0.00058 t ℃ \lambda = 0.54 + 0.00058{t}℃ λ=0.54+0.00058t
    在这里使用调和平均法计算节点的导热系数 KaTeX parse error: Undefined control sequence: \notag at position 165: …}+\lambda_{P}} \̲n̲o̲t̲a̲g̲ ̲\\ \lambda_…

    有内热源源二维保温材料导热温度分布

    当保温材料内存在内热源 S = 0.001 T 2 S=0.001T^{2} S=0.001T2

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WJFQd9JL-1584440806248)(bw_s.png)]{#img
    width=“15cm”}

    无内热源源二维保温材料导热温度分布

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rbKknyO2-1584440806249)(bw.png)]{#img
    width=“15cm”}

    结果分析与讨论

    对于保温材料来说,由于导热系数很小所以在右端界面下由对流换热带走的热量很少,所以不同的换热系数对温度场的影响并不明显,而内热源的存在使得保温材料的温度场变化非常明显,由于保温材料导热系数低,
    所以热量难以散失所以在存在内热源的情况下,整体温度较高,升温非常明显。

    对于铜体来说,与保温材料不同,对比有热源与无内热源的图像可以看出,两者温度场的差别并不大,这是因为铜体材料的导热系数很大,所以热量散失的很多,
    ,所以内热源的存在对温度场的影响不明显;而在对流换热系数改变的情况下其温度场的变化非常明显,当对流换热系数很小时,即界面的换热能力较弱时,温度场的
    差异主要存在于上下两个边界条件,热流密度 q = 10 W / m 2 q=10W/m^{2} q=10W/m2处,由于与外界存在热交换所以温度较高,而对于绝热的边界条件的下边界来说,温度较低,在对流换热系数
    很大时,右端界面的换热能里逐步增强,较大的对流换热系数的情况下,在此界面带走了大量的热,所以相对而言整体温度在下降,温度变化打的地方也从上下方向转移到了
    左右方向。

    综上所述,对流换热系数的改变对铜体的温度场影响较大,而内热源的存在与否对保温材料的温度场具有更明显的贡献。

    参考文献

    [1]黄善波,刘忠良.计算传热学基础 [M].
    东营:中国石油大学出版社,2009.12

    [2] 同登科,周田生,张高民.计算方法(第二版)[M].
    东营:中国石油大学出版社,2009.7

    附录

    # 无内热源保温材料
    
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.axes3d import Axes3D
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    from matplotlib import cm
    
    M = 15   # the node number of x-axis 
    N = 10   # the node number of y-axis
    
    # h = 5000.
    Tf = 300.
    dx = 1.5 / M
    dy = 1. / N
    dv = dx * dy
    
    
    a = np.zeros(5)
    a[1] = M/(1.5 * N)    # the value of a_W 
    a[2] = M/(1.5 * N)    # the value of a_E
    a[3] = 1.5 * N / M    # the value of a_S
    a[4] = 1.5 * N / M    # the value of a_N
    a[0] = a[1] + a[2] + a[3] + a[4] # the value of a_P
    
    # T = np.zeros((M,N))
    Tw = 300.
    
    class TMDAwithADI(object):
        def __init__(self):
            super().__init__()
        def lam(self,T):
            return 0.54 + 0.00058 * (T-273.15)
        
        def lamb(self,T1,T2):
            return 2*self.lam(T1)*self.lam(T2) / (self.lam(T1)+self.lam(T2))
    
        def calc(self,h,T):
            for k in range(1000):
                # Scan x-axis
                temp = np.array(T)
                for i in range(M):
                    for j in range(N):
                        # 左边界
                        if i == 0:
                            # 上节点
                            if j == 0:
                                T[i,j] = (self.lamb(Tw,temp[i,j])*self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1] + 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0] - self.lamb(temp[i,j],temp[i,j])*a[4] )
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0] - self.lamb(temp[i,j],temp[i,j])*a[3] )
                            # 内节点
                            else:
                                T[i,j] = (self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] )
                    
                        # 右节点
                        elif i == M-1:
                            S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                            S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                            # 上节点
                            if j == 0:
                                T[i,j] = (self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(temp[i,j],temp[i,j+1])*a[3]*temp[i,j+1]+ 15 / M+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-a[2]-self.lamb(temp[i,j],temp[i,j])*a[4] )
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1]+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2]-self.lamb(temp[i,j],temp[i,j])*a[3] )
                            # 内节点
                            else:
                                T[i,j] = (self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1]+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2] )
                        
                        # 内节点
                        else:
                            # 上节点
                            if j == 0:
                                T[i,j] = (self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j],temp[i,j+1])*a[3]*temp[i,j+1]+ 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[4] )
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[3] )
                            # 内节点
                            else:
                                T[i,j] = (self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j],temp[i,j-1])*a[4]*temp[i,j-1] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] )
    
                # Scan y-axis
                for j in range(N):
                    for i in range(M):
                        # 上边界
                        if j == 0:
                            # 左节点
                            if i == 0:
                                T[i,j] = (self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1]+ 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[4] )
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1]+15./M +S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2]-self.lamb(temp[i,j],temp[i,j])*a[4] )
                            # 内节点
                            else:
                                T[i,j] = (self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1]+ 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[4] )
                        
                        # 下边界
                        elif j == N-1:
                            # 左节点
                            if i == 0:
                                T[i,j] = (self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*T[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[3] )
                            # 右节点
                            elif i == M-1:
                                S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i,j],T[i,j-1])*a[4]*T[i,j-1]+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2]-self.lamb(temp[i,j],temp[i,j])*a[3] )
                            # 内节点
                            else:
                                T[i,j] = (self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*T[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[3] )
                        # 内节点
                        else:
                            # 左节点
                            if i == 0:
                                T[i,j] = (self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j-1],temp[i,j])*a[4]*T[i,j-1] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] )
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*T[i,j-1] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1] + S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2] )
                            # 内节点
                            else:
                                T[i,j] = (self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j-1],temp[i,j])*a[4]*T[i,j-1] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] )
    
    
    # gs = TMDAwithADI()
    # T = np.zeros((M,N))
    # gs.calc(5000,T)
    # fig = plt.figure()
    # x = np.arange(N)
    # y = np.arange(M)
    # X,Y = np.meshgrid(x,y)
    # plt.contourf(x,y,T,cmap=cm.jet)
    # plt.colorbar()
    # plt.show()
    
    # 有内热源保温材料
    
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.axes3d import Axes3D
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    from matplotlib import cm
    
    M = 15   # the node number of x-axis 
    N = 10   # the node number of y-axis
    
    # h = 5000.
    Tf = 300.
    dx = 1.5 / M
    dy = 1. / N
    dv = dx * dy
    
    
    a = np.zeros(5)
    a[1] = M/(1.5 * N)    # the value of a_W 
    a[2] = M/(1.5 * N)    # the value of a_E
    a[3] = 1.5 * N / M    # the value of a_S
    a[4] = 1.5 * N / M    # the value of a_N
    a[0] = a[1] + a[2] + a[3] + a[4] # the value of a_P
    
    # T = np.zeros((M,N))
    Tw = 300.
    
    class TMDAwithADI(object):
        def __init__(self):
            super().__init__()
        def lam(self,T):
            return 0.54 + 0.00058 * (T-273.15)
        
        def lamb(self,T1,T2):
            return 2*self.lam(T1)*self.lam(T2) / (self.lam(T1)+self.lam(T2))
    
        def calc(self,h,T):
            for k in range(10000):
                # Scan x-axis
                temp = np.array(T)
                for i in range(M):
                    for j in range(N):
                        # 左边界
                        if i == 0:
                            # 上节点
                            if j == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(Tw,temp[i,j])*self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1] + 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0] - self.lamb(temp[i,j],temp[i,j])*a[4] + 1e-3*2*temp[i,j]*dv)
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0] - self.lamb(temp[i,j],temp[i,j])*a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] + 1e-3*2*temp[i,j]*dv)
                    
                        # 右节点
                        elif i == M-1:
                            S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                            S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                            # 上节点
                            if j == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(temp[i,j],temp[i,j+1])*a[3]*temp[i,j+1]+ 15 / M+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-a[2]-self.lamb(temp[i,j],temp[i,j])*a[4] + 1e-3*2*temp[i,j]*dv)
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1]+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2]-self.lamb(temp[i,j],temp[i,j])*a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1]+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2] + 1e-3*2*temp[i,j]*dv)
                        
                        # 内节点
                        else:
                            # 上节点
                            if j == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j],temp[i,j+1])*a[3]*temp[i,j+1]+ 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[4] + 1e-3*2*temp[i,j]*dv)
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*temp[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(T[i-1,j],temp[i,j])*a[1]*T[i-1,j] + self.lamb(T[i+1,j],temp[i,j])*a[2]*T[i+1,j] + self.lamb(temp[i,j],temp[i,j-1])*a[4]*temp[i,j-1] + self.lamb(temp[i,j+1],temp[i,j])*a[3]*temp[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] + 1e-3*2*temp[i,j]*dv)
    
                # Scan y-axis
                for j in range(N):
                    for i in range(M):
                        # 上边界
                        if j == 0:
                            # 左节点
                            if i == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1]+ 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[4] + 1e-3*2*temp[i,j]*dv)
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1]+15./M +S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2]-self.lamb(temp[i,j],temp[i,j])*a[4] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1]+ 15. / M)/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[4] + 1e-3*2*temp[i,j]*dv)
                        
                        # 下边界
                        elif j == N-1:
                            # 左节点
                            if i == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*T[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[3] + 1e-3*2*temp[i,j]*dv)
                            # 右节点
                            elif i == M-1:
                                S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i,j],T[i,j-1])*a[4]*T[i,j-1]+ S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2]-self.lamb(temp[i,j],temp[i,j])*a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*T[i,j-1])/(self.lamb(temp[i,j],temp[i,j])*a[0]-self.lamb(temp[i,j],temp[i,j])*a[3] + 1e-3*2*temp[i,j]*dv)
                        # 内节点
                        else:
                            # 左节点
                            if i == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(Tw,temp[i,j])*a[1]*Tw + self.lamb(T[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j-1],temp[i,j])*a[4]*T[i,j-1] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] + 1e-3*2*temp[i,j]*dv)
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i,j-1],temp[i,j])*a[4]*T[i,j-1] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1] + S_c/(N)) / (self.lamb(temp[i,j],temp[i,j])*a[0]+S_p/N-self.lamb(temp[i,j],temp[i,j])*a[2] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+self.lamb(temp[i-1,j],temp[i,j])*a[1]*temp[i-1,j] + self.lamb(temp[i+1,j],temp[i,j])*a[2]*temp[i+1,j] + self.lamb(T[i,j-1],temp[i,j])*a[4]*T[i,j-1] + self.lamb(T[i,j+1],temp[i,j])*a[3]*T[i,j+1])/(self.lamb(temp[i,j],temp[i,j])*a[0] + 1e-3*2*temp[i,j]*dv)
    
    
    # gs = TMDAwithADI()
    # T = np.zeros((M,N))
    # gs.calc(5000,T)
    # fig = plt.figure()
    # x = np.arange(N)
    # y = np.arange(M)
    # X,Y = np.meshgrid(x,y)
    # plt.contourf(x,y,T,cmap=cm.jet)
    # plt.colorbar()
    # plt.show()
    
    # 无内热源铜体
    
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.axes3d import Axes3D
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    from matplotlib import cm
    
    M = 15   # the node number of x-axis 
    N = 10   # the node number of y-axis
    
    # h = 5000.
    Tf = 300.
    
    a = np.zeros(5)
    a[1] = 400*M/(1.5 * N)    # the value of a_W 
    a[2] = 400*M/(1.5 * N)    # the value of a_E
    a[3] = 400*1.5 * N / M    # the value of a_S
    a[4] = 400*1.5 * N / M    # the value of a_N
    a[0] = a[1] + a[2] + a[3] + a[4] # the value of a_P
    
    # T = np.zeros((M,N))
    
    Tw = 300.
    
    class TMDAwithADI(object):
        def __init__(self):
            super().__init__()
    
        def lam(self,T):
            # return 0.54 + 0.00058 * (T-273.15)
            return 400 
    
        def calc(self,h,T):
            for k in range(100000):
                # Scan x-axis
                temp = np.array(T)
                for i in range(M):
                    for j in range(N):
                        # 左边界
                        if i == 0:
                            # 上节点
                            if j == 0:
                                T[i,j] = (a[1]*Tw + a[2]*T[i+1,j] + a[3]*temp[i,j+1] + 15. / M)/(a[0] - a[4])
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (a[1]*Tw + a[2]*T[i+1,j] + a[4]*temp[i,j-1])/(a[0] - a[3])
                            # 内节点
                            else:
                                T[i,j] = (a[1]*Tw + a[2]*T[i+1,j] + a[4]*temp[i,j-1] + a[3]*temp[i,j+1])/(a[0])
                    
                        # 右节点
                        elif i == M-1:
                            S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                            S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                            # 上节点
                            if j == 0:
                                T[i,j] = (a[1]*T[i-1,j] + a[3]*temp[i,j+1]+ 15 / M+ S_c/(N)) / (a[0]+S_p/N-a[2]-a[4])
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (a[1]*T[i-1,j] + a[4]*temp[i,j-1]+ S_c/(N)) / (a[0]+S_p/N-a[2]-a[3])
                            # 内节点
                            else:
                                T[i,j] = (a[1]*T[i-1,j] + a[4]*temp[i,j-1] + a[3]*temp[i,j+1]+ S_c/(N)) / (a[0]+S_p/N-a[2])
                        
                        # 内节点
                        else:
                            # 上节点
                            if j == 0:
                                T[i,j] = (a[1]*T[i-1,j] + a[2]*T[i+1,j] + a[3]*temp[i,j+1]+ 15. / M)/(a[0]-a[4])
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (a[1]*T[i-1,j] + a[2]*T[i+1,j] + a[4]*temp[i,j-1])/(a[0]-a[3])
                            # 内节点
                            else:
                                T[i,j] = (a[1]*T[i-1,j] + a[2]*T[i+1,j] + a[4]*temp[i,j-1] + a[3]*temp[i,j+1])/(a[0])
    
                # Scan y-axis
                for j in range(N):
                    for i in range(M):
                        # 上边界
                        if j == 0:
                            # 左节点
                            if i == 0:
                                T[i,j] = (a[1]*Tw + a[2]*temp[i+1,j] + a[3]*T[i,j+1]+ 15. / M)/(a[0]-a[4])
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (a[1]*temp[i-1,j] + a[3]*T[i,j+1]+15./M +S_c/(N)) / (a[0]+S_p/N-a[2]-a[4])
                            # 内节点
                            else:
                                T[i,j] = (a[1]*temp[i-1,j] + a[2]*temp[i+1,j] + a[3]*T[i,j+1]+ 15. / M)/(a[0]-a[4])
                        
                        # 下边界
                        elif j == N-1:
                            # 左节点
                            if i == 0:
                                T[i,j] = (a[1]*Tw + a[2]*temp[i+1,j] + a[4]*T[i,j-1])/(a[0]-a[3])
                            # 右节点
                            elif i == M-1:
                                S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (a[1]*temp[i-1,j] + a[4]*T[i,j-1]+ S_c/(N)) / (a[0]+S_p/N-a[2]-a[3])
                            # 内节点
                            else:
                                T[i,j] = (a[1]*temp[i-1,j] + a[2]*temp[i+1,j] + a[4]*T[i,j-1])/(a[0]-a[3])
                        # 内节点
                        else:
                            # 左节点
                            if i == 0:
                                T[i,j] = (a[1]*Tw + a[2]*temp[i+1,j] + a[4]*T[i,j-1] + a[3]*T[i,j+1])/(a[0])
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (a[1]*temp[i-1,j] + a[4]*T[i,j-1] + a[3]*T[i,j+1] + S_c/(N)) / (a[0]+S_p/N-a[2])
                            # 内节点
                            else:
                                T[i,j] = (a[1]*temp[i-1,j] + a[2]*temp[i+1,j] + a[4]*T[i,j-1] + a[3]*T[i,j+1])/(a[0])
    
    # 有内热源铜体
    
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.axes3d import Axes3D
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    from matplotlib import cm
    
    M = 15   # the node number of x-axis 
    N = 10   # the node number of y-axis
    
    # h = 5000.
    Tf = 300.
    dx = 1.5 / M
    dy = 1. / N
    dv = dx * dy
    
    
    a = np.zeros(5)
    a[1] = 0.54*M/(1.5 * N)    # the value of a_W 
    a[2] = 0.54*M/(1.5 * N)    # the value of a_E
    a[3] = 0.54*1.5 * N / M    # the value of a_S
    a[4] = 0.54*1.5 * N / M    # the value of a_N
    a[0] = a[1] + a[2] + a[3] + a[4] # the value of a_P
    
    # T = np.zeros((M,N))
    Tw = 300.
    
    class TMDAwithADI(object):
        def __init__(self):
            super().__init__()
        def lam(self,T):
            # return 0.54 + 0.00058 * (T-273.15)
            return 0.54
    
        def calc(self,h,T):
            for k in range(10000):
                # Scan x-axis
                temp = np.array(T)
                for i in range(M):
                    for j in range(N):
                        # 左边界
                        if i == 0:
                            # 上节点
                            if j == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*Tw + a[2]*T[i+1,j] + a[3]*temp[i,j+1] + 15. / M)/(a[0] - a[4] + 1e-3*2*temp[i,j]*dv)
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*Tw + a[2]*T[i+1,j] + a[4]*temp[i,j-1])/(a[0] - a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*Tw + a[2]*T[i+1,j] + a[4]*temp[i,j-1] + a[3]*temp[i,j+1])/(a[0] + 1e-3*2*temp[i,j]*dv)
                    
                        # 右节点
                        elif i == M-1:
                            S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                            S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                            # 上节点
                            if j == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*T[i-1,j] + a[3]*temp[i,j+1]+ 15 / M+ S_c/(N)) / (a[0]+S_p/N-a[2]-a[4] + 1e-3*2*temp[i,j]*dv)
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*T[i-1,j] + a[4]*temp[i,j-1]+ S_c/(N)) / (a[0]+S_p/N-a[2]-a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*T[i-1,j] + a[4]*temp[i,j-1] + a[3]*temp[i,j+1]+ S_c/(N)) / (a[0]+S_p/N-a[2] + 1e-3*2*temp[i,j]*dv)
                        
                        # 内节点
                        else:
                            # 上节点
                            if j == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*T[i-1,j] + a[2]*T[i+1,j] + a[3]*temp[i,j+1]+ 15. / M)/(a[0]-a[4] + 1e-3*2*temp[i,j]*dv)
                            # 下节点
                            elif j == N-1:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*T[i-1,j] + a[2]*T[i+1,j] + a[4]*temp[i,j-1])/(a[0]-a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*T[i-1,j] + a[2]*T[i+1,j] + a[4]*temp[i,j-1] + a[3]*temp[i,j+1])/(a[0] + 1e-3*2*temp[i,j]*dv)
    
                # Scan y-axis
                for j in range(N):
                    for i in range(M):
                        # 上边界
                        if j == 0:
                            # 左节点
                            if i == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*Tw + a[2]*temp[i+1,j] + a[3]*T[i,j+1]+ 15. / M)/(a[0]-a[4] + 1e-3*2*temp[i,j]*dv)
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*temp[i-1,j] + a[3]*T[i,j+1]+15./M +S_c/(N)) / (a[0]+S_p/N-a[2]-a[4] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*temp[i-1,j] + a[2]*temp[i+1,j] + a[3]*T[i,j+1]+ 15. / M)/(a[0]-a[4] + 1e-3*2*temp[i,j]*dv)
                        
                        # 下边界
                        elif j == N-1:
                            # 左节点
                            if i == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*Tw + a[2]*temp[i+1,j] + a[4]*T[i,j-1])/(a[0]-a[3] + 1e-3*2*temp[i,j]*dv)
                            # 右节点
                            elif i == M-1:
                                S_p =  1/((1/h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1/h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*temp[i-1,j] + a[4]*T[i,j-1]+ S_c/(N)) / (a[0]+S_p/N-a[2]-a[3] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*temp[i-1,j] + a[2]*temp[i+1,j] + a[4]*T[i,j-1])/(a[0]-a[3] + 1e-3*2*temp[i,j]*dv)
                        # 内节点
                        else:
                            # 左节点
                            if i == 0:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*Tw + a[2]*temp[i+1,j] + a[4]*T[i,j-1] + a[3]*T[i,j+1])/(a[0] + 1e-3*2*temp[i,j]*dv)
                            # 右节点
                            elif i == M-1:
                                S_p =  1./((1./h + (1.5/M)/self.lam(T[i,j])))
                                S_c = Tf / ((1./h + (1.5/M)/self.lam(T[i,j])))
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*temp[i-1,j] + a[4]*T[i,j-1] + a[3]*T[i,j+1] + S_c/(N)) / (a[0]+S_p/N-a[2] + 1e-3*2*temp[i,j]*dv)
                            # 内节点
                            else:
                                T[i,j] = (1e-3*3*temp[i,j]**2*dv+a[1]*temp[i-1,j] + a[2]*temp[i+1,j] + a[4]*T[i,j-1] + a[3]*T[i,j+1])/(a[0] + 1e-3*2*temp[i,j]*dv)
    
    # 绘制不同对流换热系数的温度场
    
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.axes3d import Axes3D
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    from matplotlib import cm
    import TDMA_ADI_bw
    import TDMA_ADI_S_bw
    import TDMA_ADI_Cu
    import TDMA_ADI_S_Cu
    
    M = 15
    N = 10
    
    hs = [5,50,500,5000]
    
    tdma_adi_bw = TDMA_ADI_bw.TMDAwithADI()
    i = 0
    for h in hs:
        i += 1
        plt.subplot(2,2,i)
        T = np.zeros((M,N))
        x = np.arange(N)
        y = np.arange(M)
        tdma_adi_bw.calc(h,T)
        plt.contourf(x,y,T)
        plt.title("h = "+str(h))
        plt.colorbar()
    plt.show()
    
    
    tdma_adi_bw_s = TDMA_ADI_S_bw.TMDAwithADI()
    i = 0
    for h in hs:
        i += 1
        plt.subplot(2,2,i)
        T = np.zeros((M,N))
        x = np.arange(N)
        y = np.arange(M)
        tdma_adi_bw_s.calc(h,T)
        plt.contourf(x,y,T)
        plt.title("h = "+str(h))
        plt.colorbar()
    plt.show()
    
    tdma_adi_cu = TDMA_ADI_Cu.TMDAwithADI()
    i = 0
    for h in hs:
        i += 1
        plt.subplot(2,2,i)
        T = np.zeros((M,N))
        x = np.arange(N)
        y = np.arange(M)
        tdma_adi_cu.calc(h,T)
        plt.contourf(x,y,T)
        plt.title("h = "+str(h))
        plt.colorbar()
    plt.show()
    
    tdma_adi_s_cu = TDMA_ADI_S_Cu.TMDAwithADI()
    i = 0
    for h in hs:
        i += 1
        plt.subplot(2,2,i)
        T = np.zeros((M,N))
        x = np.arange(N)
        y = np.arange(M)
        tdma_adi_s_cu.calc(h,T)
        plt.contourf(x,y,T)
        plt.title("h = "+str(h))
        plt.colorbar()
    plt.show()
    
    # 网格独立性考核
    
        import TDMA_ADI
        import numpy as np
        import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d.axes3d import Axes3D
        from matplotlib.ticker import LinearLocator, FormatStrFormatter
        from matplotlib import cm
        
        tdma_adi = TDMA_ADI.TMDAwithADI()
        
        ns = [2, 4, 6, 8, 10, 12,14, 16,18,20]
        
        ans = []
        
        for n in ns:
            M = int(1.5 * n)
            N = int(1 * n)
            T = np.zeros((M,N))
            tdma_adi.calc(M,N,T)
            res = T[int(0.4*n*1.5),int(0.4*n*1)]
            ans.append(res)
        
        print(ans)
        plt.plot(np.linspace(0,1,len(ans)), ans)
        plt.scatter(np.linspace(0,1,len(ans)),ans)
        for xy in zip(ns,ans):
            plt.annotate("(%d, %.2f)" % xy, xy=xy)
        plt.show()
    
    展开全文
  • 用有限元法计算了加热治癌时在人体肝癌截面内由同心线圈感应的二维电磁场、稳态温度场和瞬态温度场。为改善被加热人体内电磁功率分布,提出用充有电解液的水袋为屏蔽物的局部电磁屏蔽法。此法在某种程度上能消除在...
  • Fluent简单流场和温度场计算,简单易学的教程内容,对于初学者有帮助
  • 本文采用二维完全隐式有限差分法对铜铝方形件摩擦焊接的非稳态温度场进行了有限差分计算,计算结果表明,异种方形件与圆形件的温度场有明显不同。沿焊接界面径向对角线方向随着半径增大温度先增大后减小,在轴向上及...
  • 二维稳态导热数值计算C程序代码 本代码可操作性强 应用于C++程序
  • 为研究矿井风流温度周期性变化对巷道围岩温度场的影响,依据能量守恒定律及围岩散热的特点,建立了周期性变化边界条件下的一稳态导热控制方程。引入了无因次准数,对围岩导热控制方程进行了无因次化,并利用变量分离...
  • +∞), 内部无热源,边界柱面温度保持为F(x,y)F(x,y)F(x,y),求柱内稳态温度分布。 注意到圆柱的对称性及柱面温度分布与z无关,可设柱内温度为u(x,y)u(x,y)u(x,y)与z无关。u满足二维Laplace方程圆内第一边值问题,也...

    今有无限长圆柱体 ( x 2 + y 2 < a 2 , − ∞ < z < + ∞ ) (x^2+y^2<a^2,-\infty<z<+\infty) (x2+y2<a2,<z<+), 内部无热源,边界柱面温度保持为 F ( x , y ) F(x,y) F(x,y),求柱内稳态温度分布。

    注意到圆柱的对称性及柱面温度分布与z无关,可设柱内温度为 u ( x , y ) u(x,y) u(x,y)与z无关。u满足二维Laplace方程圆内第一边值问题,也称Dirichlet问题
    { ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = 0 , x 2 + y 2 < a 2 u ∣ x 2 + y 2 = a 2 = F ( x , y ) (5) \begin{cases} \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=0,\quad x^2+y^2<a^2 \\ u|_{x^2+y^2=a^2}=F(x,y) \tag{5} \end{cases} {x22u+y22u=0,x2+y2<a2ux2+y2=a2=F(x,y)(5)
    在第1章中曾就 F ( x , y ) F(x,y) F(x,y)为多项式利用调和多项式求解,对任意给定的 F ( x , y ) F(x,y) F(x,y)也可利用“复变函数”中给出的Poisson积分公式求解,这里用分离变量法求解。

    由于求解区域为平面圆域,应采用极坐标 ( r , θ ) : x = r ⋅ c o s θ , y = r ⋅ s i n θ (r,\theta):x=r·cos\theta,y=r·sin\theta (r,θ):x=rcosθ,y=rsinθ。记柱内温度分布为 u ( r , θ ) u(r,\theta) u(r,θ),极坐标下边值(5)式为
    { 1 r ∂ ∂ r ( r ∂ u ∂ r ) + 1 r 2 ∂ 2 u ∂ θ 2 = 0 u ∣ r = a = F ( a c o s θ , a s i n θ ) = f ( θ ) \begin{cases} \frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial u}{\partial r})+\frac{1}{r^2}\frac{\partial^2u}{\partial \theta^2}=0 \\ u|_{r=a}=F(acos\theta,asin\theta)=f(\theta) \end{cases} {r1r(rru)+r21θ22u=0ur=a=F(acosθ,asinθ)=f(θ)
    这是 ( r , θ ) (r,\theta) (r,θ)平面的矩形域 ( 0 ≤ r < a , − ∞ < θ < + ∞ ) (0\leq r<a,-\infty<\theta<+\infty) (0r<a,<θ<+)上的问题,可进行分离变量。

    u ( r , θ ) = R ( r ) Θ ( θ ) u(r,\theta)=R(r)\Theta(\theta) u(r,θ)=R(r)Θ(θ),代入方程,同除以 u u u,有
    1 r ( r R ′ ( r ) ) ’ R ( r ) + Θ ′ ′ ( θ ) r 2 Θ ( θ ) = 0 \frac{\frac{1}{r}(rR'(r))’}{R(r)}+\frac{\Theta''(\theta)}{r^2\Theta(\theta)}=0 R(r)r1(rR(r))+r2Θ(θ)Θ(θ)=0
    分离出两个常微分方程
    Θ ′ ′ ( θ ) + λ Θ ( θ ) = 0 ( 6 a ) 1 r ( r R ′ ( r ) ) ′ − λ r 2 R ( r ) = 0 ( 7 a ) \Theta''(\theta)+\lambda \Theta(\theta)=0 \quad (6a) \\ \frac{1}{r}(rR'(r))'-\frac{\lambda}{r^2}R(r)=0 \quad (7a) Θ(θ)+λΘ(θ)=0(6a)r1(rR(r))r2λR(r)=0(7a)
    Θ ( θ ) \Theta(\theta) Θ(θ)蕴含周期条件
    Θ ( θ + 2 π ) = Θ ( θ ) \Theta(\theta+2\pi)=\Theta(\theta) Θ(θ+2π)=Θ(θ)
    Θ ( θ ) \Theta(\theta) Θ(θ)的方程构成固有值问题
    { Θ ′ ′ ( θ ) + λ Θ ( θ ) = 0 Θ ( θ + 2 π ) = Θ ( θ ) (6b) \begin{cases} \Theta''(\theta)+\lambda\Theta(\theta)=0 \\ \Theta(\theta+2\pi)=\Theta(\theta) \tag{6b} \end{cases} {Θ(θ)+λΘ(θ)=0Θ(θ+2π)=Θ(θ)(6b)
    类似上例讨论可知

    λ = − ω 2 < 0 \lambda=-\omega^2<0 λ=ω2<0时,(6a)式的通解为 Θ ( θ ) = A e ω θ + B e − ω θ \Theta(\theta)=Ae^{\omega\theta}+Be^{-\omega\theta} Θ(θ)=Aeωθ+Beωθ,无非零周期解;

    λ = 0 \lambda=0 λ=0时,通解为 Θ ( θ ) = A + B θ \Theta(\theta)=A+B\theta Θ(θ)=A+Bθ,有周期解 Θ 0 ( θ ) ≡ 1 \Theta_0(\theta)\equiv 1 Θ0(θ)1;

    λ = ω 2 > 0 \lambda=\omega^2>0 λ=ω2>0时,(6a)式的通解为
    Θ ( θ ) = A c o s ω θ + B s i n w θ \Theta(\theta)=Acos\omega\theta+Bsinw\theta Θ(θ)=Acosωθ+Bsinwθ
    当且仅当 w = n ( n = 1 , 2 , ⋅ ⋅ ⋅ ) w=n(n=1,2,···) w=n(n=1,2,)时,该通解以 2 π 2\pi 2π为周期。

    综上所述,以 2 π 2\pi 2π为周期的周期条件下的固有值问题(6b)式的固有值是
    λ n = n 2 , n = 0 , 1 , 2 , ⋅ ⋅ ⋅ , \lambda_n=n^2,\quad n=0,1,2,···, λn=n2,n=0,1,2,,
    相应的固有函数为
    Θ n ( θ ) = A n c o s   n θ + B n s i n   n θ \Theta_n(\theta)=A_ncos\space n\theta+B_nsin\space n\theta Θn(θ)=Ancos nθ+Bnsin nθ
    或记为
    Θ n ( θ ) = { c o s   n θ s i n   n θ } \Theta_n(\theta)= \begin{Bmatrix} cos\space n\theta \\ sin \space n\theta \end{Bmatrix} Θn(θ)={cos nθsin nθ}
    注意,这里对 n ≠ 0 n\neq 0 n=0,一个固有值对应于两个线性无关的固有函数,物理上称为简并现象。

    当固有值 λ n = n 2 \lambda_n=n^2 λn=n2,相应的 R ( r ) R(r) R(r)的方程(7a)是欧拉(Euler)方程
    r 2 R ′ ′ ( r ) + r R ′ ( r ) − n 2 R ( r ) = 0 r^2R''(r)+rR'(r)-n^2R(r)=0 r2R(r)+rR(r)n2R(r)=0
    在自变量代换 r = e t r=e^t r=et下,变为常系数方程
    d 2 R d t 2 − n 2 R = 0 \frac{d^2R}{dt^2}-n^2R=0 dt2d2Rn2R=0
    其通解为
    R 0 = C 0 + D 0 t = C 0 + D 0 l n r R n = C n e n t + D n e − n t = C n r n + D n r − n , n = 1 , 2 , ⋅ ⋅ ⋅ , R_0=C_0+D_0t=C_0+D_0lnr \\ R_n=C_ne^{nt}+D_ne^{-nt}=Cnr^n+D_nr^{-n},\quad n=1,2,···, R0=C0+D0t=C0+D0lnrRn=Cnent+Dnent=Cnrn+Dnrn,n=1,2,
    由问题的物理意义,显然在中心 r = 0 r=0 r=0需有自然边界条件
    ∣ R ( 0 ) ∣ < + ∞ |R(0)|<+\infty R(0)<+

    D 0 = D n = 0 R 0 ( r ) = 1 , R n ( r ) = r n D_0=D_n=0 \\ R_0(r)=1,\quad R_n(r)=r^n D0=Dn=0R0(r)=1,Rn(r)=rn
    得园内二维Laplace方程变量分离形状的周期解
    u 0 ( r , θ ) = 1 , u n ( r , θ ) = { c o s   n θ s i n   n θ } r n , n = 1 , 2 , ⋅ ⋅ ⋅ . u_0(r,\theta)=1, \quad \\ u_n(r,\theta)= \begin{Bmatrix} cos\space n\theta \\ sin\space n\theta \end{Bmatrix} r^n,\quad n=1,2,···. u0(r,θ)=1,un(r,θ)={cos nθsin nθ}rn,n=1,2,.
    由叠加原理,设
    u ( r , θ ) = A 0 2 + ∑ n = 1 + ∞ ( r a ) n ( A n c o s   n θ + B n s i n   n θ ) (8) u(r,\theta) =\frac{A_0}{2}+\sum_{n=1}^{+\infty}(\frac{r}{a})^n(A_ncos \space n\theta+B_nsin\space n\theta) \tag{8} u(r,θ)=2A0+n=1+(ar)n(Ancos nθ+Bnsin nθ)(8)
    代入边界条件(5)式得
    u ( a , θ ) = A 0 2 + ∑ n = 1 + ∞ ( A n c o s   n θ + B n s i n   n θ ) = f ( θ ) u(a,\theta)=\frac{A_0}{2}+\sum_{n=1}^{+\infty}(A_ncos\space n\theta+B_nsin\space n\theta)=f(\theta) u(a,θ)=2A0+n=1+(Ancos nθ+Bnsin nθ)=f(θ)
    这是 f ( θ ) f(\theta) f(θ) [ 0 , 2 π ] [0,2\pi] [0,2π]上的Fourier展开式,故对 n = 0 , 1 , 2 , ⋅ ⋅ ⋅ n=0,1,2,··· n=0,1,2,都有
    A n = 1 π ∫ 0 2 π f ( φ ) c o s n φ d φ ( 9 a ) B n = 1 π ∫ 0 2 π f ( φ ) s i n n φ d φ ( 9 b ) A_n=\frac{1}{\pi}\int_0^{2\pi}f(\varphi)cosn\varphi d\varphi \quad (9a)\\ B_n=\frac{1}{\pi}\int_0^{2\pi}f(\varphi)sinn\varphi d\varphi \quad (9b) An=π102πf(φ)cosnφdφ(9a)Bn=π102πf(φ)sinnφdφ(9b)
    将此 A n , B n A_n,B_n An,Bn代入(8)式,便可得二维Laplace方程园内第一边值问题(5)式的形式解

    f ( θ ) f(\theta) f(θ)是圆周上的连续函数时,可直接验证(8)式、(9a)式、(9b)式给出了(5)式的古典解。利用调和函数的极值原理,还可证明Dirichlet问题(5)式的解唯一、稳定。
    u ( r , θ ) = 1 2 π ∫ 0 2 π f ( φ ) d φ + 1 π ∑ n = 1 + ∞ ( r a ) n ∫ 0 2 π f ( φ ) [ c o s   n φ c o s   n θ + s i n   n φ s i n   n θ ] d φ = 1 2 π ∫ 0 2 π f ( φ ) [ 1 + 2 ∑ n = 1 + ∞ ( r a ) n c o s n ( φ − θ ) ] d φ u(r,\theta)=\frac{1}{2\pi}\int_0^{2\pi}f(\varphi)d\varphi+\frac{1}{\pi}\sum_{n=1}^{+\infty}(\frac{r}{a})^n\int_0^{2\pi}f(\varphi)[cos\space n\varphi cos\space n\theta+sin\space n\varphi sin\space n\theta]d\varphi \\ =\frac{1}{2\pi}\int_0^{2\pi}f(\varphi)[1+2\sum_{n=1}^{+\infty}(\frac{r}{a})^ncosn(\varphi-\theta)]d\varphi u(r,θ)=2π102πf(φ)dφ+π1n=1+(ar)n02πf(φ)[cos nφcos nθ+sin nφsin nθ]dφ=2π102πf(φ)[1+2n=1+(ar)ncosn(φθ)]dφ
    若令 z = r a e i ( φ − θ ) z=\frac{r}{a}e^{i(\varphi-\theta)} z=arei(φθ),则 ( r a ) n c o s n ( ξ − θ ) = R e z n (\frac{r}{a})^ncosn(\xi-\theta)=Rez^n (ar)ncosn(ξθ)=Rezn。利用公式 ∑ n = 0 + ∞ z n = 1 1 − z ( ∣ z ∣ < 1 ) \sum_{n=0}^{+\infty}z^n=\frac{1}{1-z}(|z|<1) n=0+zn=1z1(z<1),不难求出上式中
    1 + 2 ∑ n = 1 + ∞ ( r a ) n c o s   n ( φ − θ ) = a 2 − r 2 a 2 − 2 a r   c o s ( φ − θ ) + r 2 1+2\sum_{n=1}^{+\infty}(\frac{r}{a})^ncos\,n(\varphi-\theta)=\frac{a^2-r^2}{a^2-2ar\,cos(\varphi-\theta)+r^2} 1+2n=1+(ar)ncosn(φθ)=a22arcos(φθ)+r2a2r2
    从而
    u ( r , θ ) = a 2 − r 2 2 π ∫ 0 2 π f ( φ ) a 2 − 2 a r c o s ( φ − θ ) + r 2 d φ u(r,\theta)=\frac{a^2-r^2}{2\pi}\int_{0}^{2\pi}\frac{f(\varphi)}{a^2-2arcos(\varphi-\theta)+r^2}d\varphi u(r,θ)=2πa2r202πa22arcos(φθ)+r2f(φ)dφ
    即为Poisson积分公式

    公式(8)可作为园内Laplace方程的一般形式。从以上推导公式(8)的过程可见,圆外 ( r > a ) (r>a) (r>a) L a p l a c e Laplace Laplace方程有界解的一般形式为
    u ( r , θ ) = A 0 2 + ∑ n = 1 + ∞ ( a n ) n ( A n c o s   n θ + B n s i n   n θ ) , (r>a) u(r,\theta)=\frac{A_0}{2}+\sum_{n=1}^{+\infty}(\frac{a}{n})^n(A_ncos\,n\theta+B_nsin\,n\theta),\tag{r>a} u(r,θ)=2A0+n=1+(na)n(Ancosnθ+Bnsinnθ),(r>a)
    而环域 a 1 < r < a 2 a_1<r<a_2 a1<r<a2内Laplace方程解的一般形式则为
    u ( r , θ ) = C 0 2 + D 0 2 l n   r + ∑ n = 1 + ∞ ( C n r n + D n r − n ) ( A n c o s   n θ + B n s i n   n θ ) u(r,\theta)=\frac{C_0}{2}+\frac{D_0}{2}ln\,r+\sum_{n=1}^{+\infty}(C_nr^n+D_nr^{-n})(A_ncos\,n\theta+B_nsin\,n\theta) u(r,θ)=2C0+2D0lnr+n=1+(Cnrn+Dnrn)(Ancosnθ+Bnsinnθ)
    其中,系数 A n , B n , C n , D n A_n,B_n,C_n,D_n An,Bn,Cn,Dn都可根据边界条件,利用Fourier展开的系数公式定出。

    展开全文
  • 二维导热问题的ADI-TDMA算法

    千次阅读 2020-04-29 12:50:15
    交替方向隐式方法是CFD中用于求解对流-扩散问题的一种非常高效的求解方法,它将二维/三维隐式求解问题转化为各方向的一维隐式求解问题处理,并采用TDMA求解,比直接求解多维问题的线性方程组更高效,且时间步长选取...

    二维导热问题的ADI-TDMA算法

    基本原理

    对于多维非稳态导热问题,其时间项的离散格式有两种:显式格式隐式格式。在显式格式中,下一时间步的温度场由上一时间步的已知量求出,不需要求解代数方程组,但其时间步的选取受稳定性条件的限制,故通常只能取很小的时间步。隐式格式不受稳定性限制,但其需要对五对角矩阵(二维)或七对角矩阵(三维)进行求解,故每一时间步的计算量较大。

    相比之下,交替方向隐式方法[1](alternating direction implicit, ADI)是在显式与隐式之间折中的一种方法。对于二维问题,先将 x x x方向的温度值按隐式处理, y y y方向的温度值按显式处理,这样就将二维问题转化为一系列并列的一维问题,可以直接使用三对角阵算法(TDMA)进行计算;然后再倒过来,将 y y y方向隐式, x x x方向显式处理。这种方法被称为Peaceman-Rachford ADI格式。可以证明,该方法对于二维问题是绝对稳定的,故没有稳定性限制,可以取较大的时间步长,同时每一时间步的计算量较小,可大幅加快计算速度。

    下面介绍Peaceman-Rachford ADI格式下二维非稳态导热方程的求解过程。

    基于有限体积法的全隐离散方程可表示为:
    a P 0 ( T P − T P 0 ) = a E ( T E − T P ) − a W ( T P − T W ) + a N ( T N − T P ) − a S ( T P − T S ) + S Δ x Δ y a_P^0(T_P-T_P^0)=a_E(T_E-T_P)-a_W(T_P-T_W)+a_N(T_N-T_P)-a_S(T_P-T_S)+S\Delta x\Delta y aP0(TPTP0)=aE(TETP)aW(TPTW)+aN(TNTP)aS(TPTS)+SΔxΔy

    (1)沿 x x x方向求解:

    先取 Δ t / 2 \Delta t/2 Δt/2的时间步长,对x方向做隐式处理,求解温度场的中间量 U U U,则全隐离散方程可改成:
    a P 0 ( U P − T P 0 ) = a E ( U E − U P ) − a W ( U P − U W ) + a N ( T N 0 − T P 0 ) − a S ( T P 0 − T S 0 ) + S Δ x Δ y a_P^0(U_P-T_P^0)=a_E(U_E-U_P)-a_W(U_P-U_W)+a_N(T_N^0-T_P^0)-a_S(T_P^0-T_S^0)+S\Delta x\Delta y aP0(UPTP0)=aE(UEUP)aW(UPUW)+aN(TN0TP0)aS(TP0TS0)+SΔxΔy

    整理后得到:

    ( a P 0 + a E + a W ) U P = a E U E + a W U W + ( a P 0 − a N − a S ) T P 0 + a N T N 0 + a S T S 0 + S Δ x Δ y (a_P^0+a_E+a_W)U_P=a_EU_E+a_WU_W+(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y (aP0+aE+aW)UP=aEUE+aWUW+(aP0aNaS)TP0+aNTN0+aSTS0+SΔxΔy

    在上式中,除了 U P U_P UP U E U_E UE U W U_W UW,其他量都是已知的。令 a P = a P 0 + a E + a W a_P=a_P^0+a_E+a_W aP=aP0+aE+aW D = ( a P 0 − a N − a S ) T P 0 + a N T N 0 + a S T S 0 + S Δ x Δ y D=(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y D=(aP0aNaS)TP0+aNTN0+aSTS0+SΔxΔy,则上式可以简化为:

    a P U P = a E U E + a W U W + D a_PU_P=a_EU_E+a_WU_W+D aPUP=aEUE+aWUW+D

    通过使用沿 x x x方向的TDMA可以求得温度场中间量 U U U

    (2)沿 y y y方向求解:

    在求得中间量 U U U的基础上,再取 Δ t / 2 \Delta t/2 Δt/2的时间步长,对y方向做隐式处理,求解下一时刻温度场 T 1 T^1 T1,则全隐离散方程可改成:
    a P 0 ( T P − T P 0 ) = a E ( T E 1 − T P 1 ) − a W ( T P 1 − T W 1 ) + a N ( U N − U P ) − a S ( U P − U S ) + S Δ x Δ y a_P^0(T_P-T_P^0)=a_E(T_E^1-T_P^1)-a_W(T_P^1-T_W^1)+a_N(U_N-U_P)-a_S(U_P-U_S)+S\Delta x\Delta y aP0(TPTP0)=aE(TE1TP1)aW(TP1TW1)+aN(UNUP)aS(UPUS)+SΔxΔy

    整理后得到:

    ( a P 0 + a N + a S ) T P 1 = a N T N 1 + a S U S 1 + ( a P 0 − a E − a W ) U P + a E U E + a E U W + S Δ x Δ y (a_P^0+a_N+a_S)T_P^1=a_NT_N^1+a_SU_S^1+(a_P^0-a_E-a_W)U_P+a_EU_E+a_EU_W+S\Delta x\Delta y (aP0+aN+aS)TP1=aNTN1+aSUS1+(aP0aEaW)UP+aEUE+aEUW+SΔxΔy

    在上式中,除了 T p 1 T_p^1 Tp1 T N 1 T_N^1 TN1 T S 1 T_S^1 TS1,其他量都是已知的。令 a P = a P 0 + a N + a S a_P=a_P^0+a_N+a_S aP=aP0+aN+aS D = ( a P 0 − a E − a W ) U P + a E U E + a E U W + S Δ x Δ y D=(a_P^0-a_E-a_W)U_P+a_EU_E+a_EU_W+S\Delta x\Delta y D=(aP0aEaW)UP+aEUE+aEUW+SΔxΔy,则上式可以简化为:

    a P T P 1 = a N T N 1 + a N T N 1 + D a_PT_P^1=a_NT_N^1+a_NT_N^1+D aPTP1=aNTN1+aNTN1+D

    再使用沿 y y y方向的TDMA求得下一时刻的温度场 T 1 T^1 T1

    边界节点的处理方式与上述过程相似,主要区别在于该节点的离散方程中未知项少了一项,例如对于顶部对流边界, a N a_N aN的取值得改为 α Δ x \alpha\Delta x αΔx T N T_N TN由未知的待求解量改为已知量 T f T_f Tf

    求解程序

    下面设计一个算例,并通过Fortran程序给出其求解过程:

    有一个宽400mm,高200mm的矩形计算域。计算域内材料的密度为2500kg/m3,热容为1372J/(kg·K)。含内热源,热生成率恒定为1050J/(m3·s)。导热系数为温度的多项式函数:

    λ = − 1.13588 × 1 0 − 15 T 5 + 3.25358 × 1 0 − 12 T 4 − 3.25305 × 1 0 − 9 T 3 + 1.32926 × 1 0 − 6 T 2 − 9.27637 × 1 0 − 5 T + 1.04478 \lambda=-1.13588\times 10^{-15}T^5+3.25358\times 10^{-12}T^4-3.25305\times 10^{-9}T^3+1.32926\times 10^{-6}T^2-9.27637\times 10^{-5}T+1.04478 λ=1.13588×1015T5+3.25358×1012T43.25305×109T3+1.32926×106T29.27637×105T+1.04478

    计算域的初始温度为1500℃,四周与环境进行对流换热,其中顶部对流给热系数 α 1 \alpha_1 α1为200W/(m2·K),底部和左右侧的对流给热系数 α 2 \alpha_2 α2为80W/(m2·K),环境温度 T f T_f Tf为20℃。模拟的物理时间长度为3h。

    采用均匀结构化网格,每个网格节点的尺寸为2mm×2mm。由于该算例的模型沿 y y y轴是对称的,实际计算域可以只取左半部分,即选取的计算域尺寸为200mm×200mm,网格节点数量为101×101。计算域右侧作为对称边界处理,该边界上 T E = T W T_E=T_W TE=TW a E = a W a_E=a_W aE=aW

    由于Peaceman-Rachford ADI格式的时间步长选取不受稳定性限制,其时间步长的选取比较灵活。由于非稳态导热随着时间的进行,会趋向于稳态,刚开始时温度变化率会很大,随后则会大幅减小,本算例选取时间步长的准则是使选定的时间步长内最大温度变化幅度不超过10℃,这样刚开始时间步长会很小,随后则会变得很大。

    计算过程中选取了6个采样点,并设定每隔5分钟输出一次这6个采样点的温度值。

    求解程序如下,模拟用到的参数在网格材料边界条件初始条件求解控制模块中通过常量定义出来,以便于参数的调整:

    module mesh
        implicit none
        real, parameter :: width = 200e-3, height = 200e-3
        real, parameter :: delta_X = 2e-3, delta_Y = 2e-3
    end module mesh
    
    module material
        implicit none
        real, parameter :: rho = 2500, Cp = 1372, S = 1050
    end module material
    
    module boundary_condition
        implicit none
        real, parameter :: alpha1 = 300, alpha2 = 80
        real, parameter :: Tf = 20
    end module boundary_condition
    
    module initial_condition
        implicit none
        real, parameter :: T0 = 1500
    end module initial_condition
    
    module solution_control
        implicit none
        real, parameter :: time_length = 180 * 60
        real, parameter :: sampling_interval = 5 * 60
        real, parameter :: max_deltaT_per_step = 10
        real, dimension(6) :: sample_points_X = (/30e-3, 100e-3, 200e-3, 30e-3, 100e-3, 200e-3/)
        real, dimension(6) :: sample_points_Y = (/100e-3, 100e-3, 100e-3, 170e-3, 170e-3, 170e-3/)
        integer :: spIndex(6,2)
    end module solution_control
    
    
    program TDMA
        use mesh
        use material
        use solution_control
        implicit none
        integer :: Nx, Ny
        
        !确定X、Y方向的网格节点数量
        Nx = width / delta_X + 1
        Ny = height / delta_Y + 1
    
        !根据采样点位置获得采样点在网格节点中的索引
        spIndex(:,1) = floor(sample_points_X / delta_X + 1)
        spIndex(:,2) = floor(sample_points_Y / delta_Y + 1)
    
        !开始计算
        call calculate(Nx, Ny)
    
    end program TDMA
    
    
    subroutine calculate(Nx, Ny)
        use mesh
        use material
        use boundary_condition
        use initial_condition
        use solution_control
        implicit none
        integer, intent(in) :: Nx, Ny
        integer :: i
        real :: current_time, current_sampling_time, delta_Tau
        real, dimension(0:Nx+1,0:Ny+1) :: T, U, T1, delta_T
        real, dimension(Nx,Ny) :: lambda, aP0, aE, aW, aN, aS, aP, Sv, D, P, Q
        real :: kx(Ny), ky(Nx)
        real :: max_delta_T, max_rate
        real :: tStart, tEnd
    
        !初始化系数矩阵A
        aP0 = 0
        aE = 0
        aW = 0
        aN = 0
        aS = 0
    
        !初始化温度矩阵
        T = Tf
        T(1:Nx,1:Ny) = T0
        U = Tf
        T1 = Tf
    
        !初始化源项矩阵
        Sv = S * delta_X * delta_Y
        Sv(1,:) = Sv(1,:) / 2
        Sv(:,1) = Sv(:,1) / 2
        Sv(:,Ny) = Sv(:,Ny) / 2
    
        !初始化输出
        write(*,"(A4,$)") "Time"
        do i = 1, 5
            write(*,"(A9,I1,$)") "SP", i
        end do
        write(*,"(A10)") "SP6"
    
        !按时间步进行温度场计算
        current_time = 0
        current_sampling_time = sampling_interval
        delta_Tau = 0.0001  !考虑到降温刚开始时的温度变化率很大,初始delta_Tau取一个很小的值
        call CPU_TIME(tStart)
        do while (current_time < time_length)
            !更新当前时间
            current_time = current_time + delta_Tau
    
            !通过多项式函数计算材料导热系数
            lambda = -1.13588E-15 * T(1:Nx,1:Ny)**5 + 3.25358E-12 * T(1:Nx,1:Ny)**4 - 3.25305E-09 * T(1:Nx,1:Ny)**3 &
                   + 1.32926E-06 * T(1:Nx,1:Ny)**2 - 9.27637E-05 * T(1:Nx,1:Ny) + 1.04478
    
            !计算系数矩阵A
            !(1)内部节点的计算
            aP0 = rho * Cp * delta_X * delta_Y / (delta_Tau / 2)
            aE(1:Nx-1,:) = 2 / (1 / lambda(1:Nx-1,:) + 1 / lambda(2:Nx,:)) / delta_X * delta_Y
            aW(2:Nx,:) = 2 / (1 / lambda(2:Nx,:) + 1 / lambda(1:Nx-1,:)) / delta_X * delta_Y
            aN(:,1:Ny-1) = 2 / (1 / lambda(:,1:Ny-1) + 1 / lambda(:,2:Ny)) / delta_Y * delta_X
            aS(:,2:Ny) = 2 / (1 / lambda(:,2:Ny) + 1 / lambda(:,1:Ny-1)) / delta_Y * delta_X
            !(2)左侧、上下侧对流边界节点的处理
            aW(1,:) = alpha2 * delta_Y
            aN(:,Ny) = alpha1 * delta_X
            aS(:,1) = alpha2 * delta_X
            !(3)右侧对称边界的处理
            aW(Nx,:) = aW(Nx,:) * 2
            !(4)边界节点的减半处理
            !(4-1)上侧aP0、aE、aW减半
            aP0(:,Ny) = aP0(:,Ny) / 2
            aE(:,Ny) = aE(:,Ny) / 2
            aW(:,Ny) = aW(:,Ny) / 2
            !(4-2)下侧aP0、aE、aW减半
            aP0(:,1) = aP0(:,1) / 2
            aE(:,1) = aE(:,1) / 2
            aW(:,1) = aW(:,1) / 2
            !(4-3)左侧aP0、aN、aS减半
            aP0(1,:) = aP0(1,:) / 2
            aN(1,:) = aN(1,:) / 2
            aS(1,:) = aS(1,:) / 2
    
            !沿x方向使用TDMA隐式求解
            aP = aP0 + aE + aW
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aN - aS) * T(1:Nx,1:Ny) + aN * T(1:Nx,2:Ny+1) + aS * T(1:Nx,0:Ny-1) + Sv
            !(2)左侧边界的特殊处理
            D(1,:) = D(1,:) + aW(1,:) * Tf
            !!!!!!
            P(1,:) = aE(1,:) / aP(1,:)
            Q(1,:) = D(1,:) / aP(1,:)
            do i = 2, Nx
                kx = aP(i,:) - aW(i,:) * P(i-1,:)
                P(i,:) = aE(i,:) / kx
                Q(i,:) = (D(i,:) + aW(i,:) * Q(i-1,:)) / kx
            end do
            U(Nx,1:Ny) = Q(Nx,:)
            do i = Nx-1, 1, -1
                U(i,1:Ny) = P(i,:) * U(i+1,1:Ny) + Q(i,:)
            end do
    
            !沿y方向使用TDMA隐式求解
            aP = aP0 + aN + aS
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aE - aW) * U(1:Nx,1:Ny) + aE * U(2:Nx+1,1:Ny) + aW * U(0:Nx-1,1:Ny) + Sv
            !(2)上下侧边界的特殊处理
            D(:,1) = D(:,1) + aS(:,1) * Tf
            D(:,Ny) = D(:,Ny) + aN(:,Ny) * Tf
            !!!!!!
            P(:,1) = aN(:,1) / aP(:,1)
            Q(:,1) = D(:,1) / aP(:,1)
            do i = 2, Ny
                ky = aP(:,i) - aS(:,i) * P(:,i-1)
                P(:,i) = aN(:,i) / ky
                Q(:,i) = (D(:,i) + aS(:,i) * Q(:,i-1)) / ky
            end do
            T1(1:Nx,Ny) = Q(:,Ny)
            do i = Ny-1, 1, -1
                T1(1:Nx,i) = P(:,i) * T1(1:Nx,i+1) + Q(:,i)
            end do
    
            !当达到采样时间时,输出采样点温度
            if (current_time >= current_sampling_time) then
                !输出6个采样点的温度值
                write(*,"(I4,$)") floor(current_time / 60)
                do i = 1, 5
                    write(*,"(f10.2,$)") T1(spIndex(i,1), spIndex(i,2))
                end do
                write(*,"(f10.2)") T1(spIndex(6,1), spIndex(6,2))
                !
                !更新下一次需要采样的时间
                current_sampling_time = current_sampling_time + sampling_interval
            end if
    
            !计算最大冷却速率
            delta_T = abs(T1 - T)
            max_delta_T = maxval(delta_T)
            max_rate = max_delta_T / delta_Tau
    
            !根据对每个时间步内的最大温度变化率限制,以及下一次采样时间确定下一个时间步长
            delta_Tau = min(max_deltaT_per_step / max_rate, current_sampling_time - current_time)
    
            !更新温度矩阵T
            T = T1
    
        end do
        call CPU_TIME(tEnd)
        write(*,"('Calculate complete, time consumption: ', f8.4, ' sec')") (tEnd - tStart)
    
    end subroutine calculate
    

    参考文献

    [1] 陶文铨. 数值传热学[M]. 第二版. 西安: 西安交通大学出版社, 2001: 99-104.

    展开全文
  • WindFarmSimulator(WFSim)是基于二维Navier-Stokes方程的中保真,面向控制的风电场模型。 Sjoerd Boersma和Bart Doekemeijer目前正在代尔夫特理工大学积极开发该软件。 可以在网上找到有关WFSim的最新出版物: ...
  • 2. 问题分析一维稳态无源导热控制方程如下: 对于上述方程,利用有限体积法来进行离散求解。离散求解PED方程的步骤为:离散控制域(网格划分)在每一个控制体上离散控制方程插值得到界面值,完成单元离散方程组装...

    1ed71f27d9c31672b07f8cff484ad05c.png

    1. 问题描述

    如图,厚度为

    的无限大平板,导热系数
    ,板内有均匀内热源
    ,表面A温度保持在
    ,表面B温度保持在
    。求板内厚度方向温度分布。

    8cca2e9e64a52184fcd45bdc93cc19a8.png

    2. 问题分析

    一维稳态无源导热控制方程如下:

    对于上述方程,利用有限体积法来进行离散求解。

    离散求解PED方程的步骤为:

    1. 离散控制域(网格划分)
    2. 在每一个控制体上离散控制方程
    3. 插值得到界面值,完成单元离散方程
    4. 组装单元控制方程,形成整体控制方程组(Ax=b)
    5. 求解代数方程组(直接,迭代)
    6. 得到离散场变量

    下面分别按照上面给出的六个步骤进行分析。

    3. 问题求解

    3.1 域离散

    采用均匀网格将棒沿着长度方向离散,如图所示:

    55b0a93ffd3882225614e153765f96c1.png

    为了编写程序,我们需要容器来储存网格数据。对于一维问题,我们需要存储离散节点的位置坐标,单元的体积(在本例中为单元的长度乘以截面积)。

    代码实现

    import 

    3.2 方程离散

    在如图的控制体上,对控制方程进行积分得到:

    上式中间部分继续展开,得到:

    注意,上面已经采用了线性化处理的思路。一般来说,我们不能保证下式绝对成立:

    但是,为了离散的简单性,我们认为上式是成立。并且当k与T不是强耦合的时候,上面的等式带来的误差不是很大。

    在FVM中,材料属性和场变量都是存储在节点上的,但(3)式中需要到界面上的导热系数和温度梯度值。我们需要利用插值函数将节点值插值到界面上。采用不同的插值函数,就可以得到不同的离散方程,不同的离散方程其精度、收敛性和稳定性都不一样。这里不进行深入的讨论(如果讨论可能需要一大本书,并且有些问题也讨论不清楚。)

    3.3 界面值插值

    我们采用如下的差值方式:

    对于导热系数(材料参数),常采用线性加权方法插值。如下:

    其中,为一个几何系数,定义如下:

    对于均匀网格而言,。所以,对于本例中,界面上的导热系数利用如下插值得到:

    对于界面上的梯度值,也有多种插值方法,我们采用简单的线性插值,如下:

    其中,表示C点到E点之间的距离。在本例中,采用均匀网格,

    对于界面面积

    对于方程(3),还存在源项。对于源项,我们一般做如下的处理:

    上式利用了一个近似:认为单元中热源的平均值等于单元中心点处的热源值。这个近似也会产生误差。我们这里利用一点篇幅讨论一下这个误差。如果对这个讨论不感兴趣的读者,可以放心的跳过这部分以继续后面的内容。

    3.4.1 单元均值和中心点处值的差距

    我们以热源为例。假设热源是连续可微的。那么,由下列表达式成立:

    其中定义如下:

    将q在控制体中心展开:

    将(12)带入到(10)中:

    方程(14)右端第一项简化如下:

    方程(14)右端第二项简化如下(利用积分中值定理):

    方程(16)右端项分解为两项:

    对于单元中心点定义如下:

    将(18)代入(17)中,很显然:

    因此我们得到方程(14)右端第二项为0。将此结果和(15)代入到(14)中,得到:

    两边同除以,得到:

    所以,利用控制单元中心值来近似单元平均值,具有二阶精度。

    3.4 单元方程整理

    将导热系数和温度梯度在界面的插值表达式带入到(3)中,得到最后的单元离散方程如下:

    整理得到如下简洁形式:

    显然,对于每一个控制体,都用到了它自己的前后单元值。那么,我们就需要对边界处做特殊的处理。因为边界处的单元只有一个相邻单元。

    代码实现

    a_w 

    3.5 整体矩阵组装

    通过上面的分析,我们成功得到了单元的离散形式的控制方程。因为单元控制方程在离散的时候要用到相邻单元的值,但是相邻单元的值也是未知的。因此这个时候我们可以如下的方法来处理:

    1. 将相邻单元的温度值当做已知(一般给出一个猜测的场变量值),这样就可以显式的解得当前单元处的温度值,但是这样做要对全场进行多次迭代。
    2. 将每个单元的控制方程组装,形成一个对应于全场未知量的系数矩阵,求解这个整体的代数方程组,统一解得场变量。

    这两种方式本质其实是相同的。尤其是当整体系数矩阵较大的时候,CFD中控制方程的非线性导致直接求解方程组很困难,一般常采用迭代求解。那么对于整体方程组采用迭代求解就和第一种思路一样。

    整体矩阵组装的时候,将控制体进行编号,未知温度值当做列向量,形成AT=b的形式。具体的组装方式在后面程序中给出。

    代码实现

    # 方程离散和整体组装
     

    3.6 边界处理

    对于靠近A的边界单元

    从方程开始

    。此时,控制体示意图如下:

    3a2ad5d762f344c6b5da9e5ab6272420.png

    代码实现:

    # 第一个边界单元修正, 此时C对应下标0, k_A = k[0]
     

    对于靠近B的边界

    从方程开始

    。此时,控制体示意图如下:

    d18a8ddb1a529e6d0651765bf34a3d8e.png

    代码实现

    # 第二个边界单元修正,此时,C对应下标N-1, W对应下标N-2, k_B = k[N-1]
     

    至此,我们已经得到了与原来由PDE加边界条件定义的问题近似一致的离散型代数方程组AT=b。对于这类代数方程组的求解,有许许多多的算法可用。我们的问题最终形成的是一个三对角矩阵。可以利用TDM方法求解。这里用python-numpy库自带的函数来求解。

    如果本例中,我们将控制单元设置为5个,则可以得到最终的系数矩阵如下:

    [[ 

    右端列向量为:

    [

    3.7 矩阵求解

    矩阵求解的方法很多。一般采用迭代求解。

    # 代数方程组求解
     

    3.8 结果可视化

    import 

    1562a7f1def38a2cd09e4fef3bfc2d0a.png

    3.9 结果分析

    本问题有理论解,表达式如下:

    将数值解和理论解对比,如下:

    import 

    3b068b9dfbe4991236e58a90a299db59.png

    关于理论解的推导过程如下:

    对方程(27)变形得到:

    对方程(28)两边积分:

    对(29)两边积分:

    所以:

    。代入方程中确定
    两个常数得到:

    得到:

    4. 代码

    完整代码如下:

    import 
    展开全文
  • 问题分析一维稳态无源导热控制方程如下: 对于上述方程,利用有限体积法来进行离散求解。离散求解PED方程的步骤为:离散控制域(网格划分)在每一个控制体上离散控制方程插值得到界面值,完成单元离散方程组装单元...
  • 为了了解矩形腔内非Boussinesq流体自然对流换热的特有现象和规律,利用有限容积法对矩形腔内的冷水自然对流进行了一系列数值模拟,得到了不同宽深比Rayleigh数下的流场和温度场,并对不同条件下的壁面传热特性进行了...
  • 介绍了激光重熔、熔粗温度场的研究现状,叙述了较为典型的二维、三维的稳态和非稳态传热模型,分析比较了从该模型出发所得到的结果与实际情况间的差异。根据已有研究所得到的结论,简要讨论了影响模型建立的几个主要...
  • “Fluid Timescale Control”选择“Physical Timescale”并设置为2s,这个稳态求解和Fluent软件很不一样。残差设置为1E-4。 6. 保存并求解 大概30多步,残差就收敛到1E-4级别。 7. 后处理 CFD-Post后处理模块,由于...
  • 计算了不同负载运行时样机定子的稳态温度场。计算结果与实测值的比较,验证了所采用计算模型及方法的合理性。该电机温度场计算模型可以应用到其他同类电机定子温度场的计算与分析。在该温度场计算模型的基础上,分析...
  • 基于3种不同的正交基函数系(三角函数正交系,Legendre正交多项式和Chebyshev正交多项式),比较研究了一类二维稳态热传导方程反问题中边界温度场的数值反演方法。首先引入泛函,根据线性偏微分方程的叠加原理,将所...
  • 为了研究冻土路基温度场及变形场的动态变化规律,基于伴有相变的路基非稳态温度场控制方程和冻土路基变形场二维数值计算模型,对冬季冻土路基温度场和变形场进行了计算分析,得出路基深层土中的温度变化滞后于表层土...
  • 用有限元法计算了BDFM二维稳态和暂态温度场,揭示了不同结构和运行工况下BDFM的温度分布规律。结果表明,最高温度出现在接线盒附近的定子控制绕组;在控制绕组电流频率一定的情况下,增加散热翅高度、提高散热翅风速...

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 300
精华内容 120
关键字:

二维稳态温度场