python之有限体积法求解一维热传导问题

最新推荐文章于 2024-01-05 17:03:32 发布
最新推荐文章于 2024-01-05 17:03:32 发布
阅读量997 收藏 7
点赞数
文章标签： python numpy matplotlib 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/T20151470/article/details/132890446
版权

1、问题描述

考虑均匀发热无限大平板的稳定导热问题，上图中，A、B两面恒温，控制方程为如下形式：

\frac{\partial }{\partial x}\left(\Gamma {\frac{\partial \phi }{\partial x}} +S_\phi \right )=0

\Gamma为扩散系数，k为材料传热系数，给定厚度L=2cm，k=5W/\left(m^2\cdot K \right )，T_A和T_B分别为100℃和400℃，发热量q为500kW/m^3，确定平板内的稳态温度分布？

2、基于有限体积法的公式推导

假定在y方向和z方向上尺寸无限大，温度梯度只在x方向上有意义。

把整个区域划分为4个控制体，给定\delta x=0.005m

由于传热系数是常数，定义k=k_e=k_w，体机源项S_\phi =q，节点2处的一般离散形式为

a_PT_P=a_ET_E+a_WT_W+b

其中

a_E=\frac{kA}{\delta x^2}

a_W=\frac{kA}{\delta x^2}

a_P=a_E+a_W

b=p

节点3得到控制方程和节点2一样。

对于节点1和4的控制体，在边界点和临近点之间的温度，采用线性近似。

对于恒定传热系数和热流量，方程变为

k\left(\frac{\partial T}{\partial x} \right )_e-k\left(\frac{\partial T}{\partial x} \right )_w+q\delta x=0

在控制体1的左侧和右侧，引入线性梯度近似：

kA\cdot\frac{T_E-T_P}{\delta x}-kA\cdot\frac{T_P-T_A}{\delta x/2}+q\delta x=0

对于节点1，得到离散方程：

a_PT_P=a_ET_E+a_WT_W+b

其中

a_E=\frac{kA}{\delta x^2}

a_W=0

a_P=a_E+a_W+\frac{2kA}{\delta x}

b=q\delta x+\frac{2kA}{\delta x}T_A

同理，节点4的离散方程为：

a_PT_P=a_ET_E+a_WT_W+b

其中

a_E=0

a_W=\frac{kA}{\delta x^2}

a_P=a_E+a_W+\frac{2kA}{\delta x}

b=q\delta x+\frac{2kA}{\delta x}T_B

将传热系数k=5W/\left(m^2\cdot K \right )，热流量q=500kW/m^3，\delta x=0.005m以及单位面积A=1m^2等系数代入离散方程，得到如下方程组：

\begin{bmatrix} 3000&-1000 &0 &0 \\ -1000&2000 &-1000 &0 \\ 0&-1000 &2000 &-1000 \\ 0&0 &-1000 &3000 \end{bmatrix}  \begin{bmatrix} T1\\ T2 \\ T3 \\T4 \end{bmatrix} =\begin{bmatrix} 2000T_A+2500\\ 2500 \\ 2500 \\2000T_B+2500 \end{bmatrix}

3、python求解分析

3.1高斯消去法

高斯消去法是一种线性代数中的算法，用于求解线性方程组。它的优点和缺点如下：

优点：

精度高：高斯消去法通过将系数矩阵化为上三角矩阵，减少了求解过程中的误差，提高了求解的精度。

稳定性好：在求解过程中，高斯消去法通过选取主元来保证计算的稳定性。主元的选取可以避免出现除数为零的情况，从而保证了计算的稳定性。

可扩展性强：高斯消去法可以很容易地扩展到求解更大规模的线性方程组。只需要增加计算机的存储空间和计算能力，就可以求解更大规模的线性方程组。

算法简单易懂：高斯消去法的求解过程可以通过手工计算来理解，也可以通过编程实现来求解。

缺点：

计算量大：高斯消去法的计算量较大，尤其是在求解大规模线性方程组时，计算量更是巨大。这会导致求解时间较长，不适合实时求解。

存储空间大：高斯消去法需要存储系数矩阵和常数向量，这会占用较大的存储空间。在求解大规模线性方程组时，存储空间的需求更是巨大。

主元选取困难：高斯消去法的主元选取对计算结果的影响很大。如果选取不当，会导致计算结果的不稳定性和精度下降。因此，主元选取是一个比较困难的问题。

无法处理奇异矩阵：如果系数矩阵是奇异矩阵，那么高斯消去法将无法求解线性方程组。

import numpy as np
def gauss(a, b):
    m, n = a.shape
    c = np.zeros(n)
    for i in range(n):
        if (a[i][i] == 0): # 用高斯消去法解线性方程组时对角线元素不能为0
            print("no answer")

# k表示第一层循环，(0，n-1)行
# i表示第二层循环,(k+1,n)行,计算该行消元的系数
# j表示列
    for k in range(n - 1):
        for i in range(k + 1, n):
            c[i] = a[i][k] / a[k][k] # 计算出系数

            for j in range(k,m): # 从K开始，减少不必要的计算
                a[i][j] = a[i][j] - c[i] * a[k][j] # 对矩阵进行高斯消去
            b[i] = b[i] - c[i] * b[k]

        print('step',k,':n',a,'n')
#print(b)
    x = np.zeros(n)

    x[n - 1] = b[n - 1] / a[n - 1][n - 1] # 解出x[n-1],为回代作准备

# 回代求出方程解
    for i in range(n-2, -1, -1):
        sum= 0.0
        for j in range(n-1, -1, -1):
            sum= sum + a[i][j] * x[j]
        x[i] = (b[i]-sum) / a[i][i]
# 输出计算结果
    for i in range(n):
        print("x" + str(i + 1) + " = ","%.2f" % x[i]) # 输出结果

if __name__ == '__main__':
    A = np.array([
        [3000, -1000, 0, 0],
        [-1000, 2000, -1000, 0],
        [0, -1000, 2000, -1000],
        [0, 0, -1000, 3000],
    ])
    b = np.array([202500, 2500, 2500, 802500])
    gauss(A, b)

step 0 :n [[ 3000 -1000     0     0]
 [    0  1666 -1000     0]
 [    0 -1000  2000 -1000]
 [    0     0 -1000  3000]] n
step 1 :n [[ 3000 -1000     0     0]
 [    0  1666 -1000     0]
 [    0     0  1399 -1000]
 [    0     0 -1000  3000]] n
step 2 :n [[ 3000 -1000     0     0]
 [    0  1666 -1000     0]
 [    0     0  1399 -1000]
 [    0     0     0  2285]] n
x1 =  140.09
x2 =  217.77
x3 =  292.81
x4 =  365.13
 

3.2TDMA算法

TDMA算法，即三对角矩阵算法，是一种用于求解具有三对角线系数矩阵的代数方程组的算法。

其优点有：

高效：TDMA算法的计算复杂度为O(n)，n为方程组的大小。这意味着对于较大规模的方程组，该算法能够快速求解。

稳定性高：该算法在求解过程中，不会出现数值不稳定性或误差累积的问题。

然而，它也有一些缺点：

需要额外的存储空间：对于大规模的问题，需要额外的存储空间来存储中间计算结果。

对初值敏感：TDMA算法的收敛速度可能会受到初值的影响，如果初值选择不当，可能会导致算法不收敛或者收敛速度很慢。

需要精确计算：由于涉及到大量的浮点运算，如果硬件或者软件出现故障，可能会导致计算结果不准确。

需要专业训练：该算法需要专业训练才能理解和使用，对于非专业人员来说可能会有一定的学习难度。

import numpy as np
def tdma(A,b):
    m,n = A.shape
    p = np.zeros(m)
    q = np.zeros(n)
    phi = np.zeros(m)

    p[0] = -A[0,1]/A[0,0]
    q[0] = b[0]/A[0,0]

    for i in range(1,m):
        if i != m-1 :
            p[i] = -A[i,i+1]/(A[i,i]+A[i,i-1]*p[i-1])
        else:
            p[i]= 0
        q[i] = (b[i]-A[i,i-1]*q[i-1])/(A[i,i]+A[i,i-1]*p[i-1])

    # 回代
    phi[m-1] = q[m-1]
    for i in range(m-2,-1,-1):
        phi[i] = p[i]*phi[i+1] +q[i]

    return phi

# 测试
A = np.array([
    [3000,-1000,0,0],
    [-1000,2000,-1000,0],
    [0,-1000,2000,-1000],
    [0,0,-1000,3000],
])
b = np.array([202500,2500,2500,802500])

phi = tdma(A,b)
print(phi)

[140.  217.5 292.5 365. ]

3.3雅可比迭代法

雅可比迭代法是一种求解线性方程组的迭代方法，其优点和缺点如下：

优点：

计算公式简单：雅可比迭代法的计算公式简单明了，每迭代一次只需要计算一次矩阵和向量的乘法。

并行计算：由于在计算过程中原始矩阵始终不变，因此雅可比迭代法比较容易进行并行计算。

缺点：

收敛速度较慢：雅可比迭代法的收敛速度通常比其他一些迭代方法慢，如高斯-赛德尔迭代法或牛顿法等。

存储空间较大：雅可比迭代法需要存储方程组的系数矩阵和常数向量，因此需要较大的存储空间，这可能在处理大规模问题时成为限制。

import numpy as np
# A为系数矩阵，b为结果矩阵
# k为最大迭代次数，tol为容差限
def Jacobi(A,b,k,tol):
    n = A.shape[1]
    D = np.eye(n)
    D[np.arange(n),np.arange(n)] = A[np.arange(n),np.arange(n)]
    LU = A - D
    X = np.zeros(n)

    # 迭代k次
    for i in range(k):
        X1= X.copy()
        D_inv = np.linalg.inv(D)
        X = -np.dot(np.dot(D_inv,LU),X) + np.dot(D_inv,b)
        if (np.max(np.abs(X-X1)))<=tol:
            print('计算收敛!')
            break
        print('n第',i,'次迭代：',X)
        print('残差值为：',np.abs(X-X1))
        print('当前最大残差为：', np.max(np.abs(X-X1)))
    return X

A = np.array([
    [3000,-1000,0,0],
    [-1000,2000,-1000,0],
    [0,-1000,2000,-1000],
    [0,0,-1000,3000],
])
b = np.array([202500,2500,2500,802500])
X = Jacobi(A, b, 1000,0.005)

n第 32 次迭代： [139.99525386 217.49138441 292.48962427 364.996059  ]
残差值为： [0.00066203 0.00562283 0.00144728 0.00257203]
当前最大残差为： 0.005622829086746606
计算收敛!

3.4高斯赛德尔迭代法

高斯-赛德尔迭代法是一种求解线性方程组的迭代方法，其优点和缺点如下：

优点：

收敛速度快：高斯-赛德尔迭代法通常具有较快的收敛速度，可以在较短的时间内求解出方程组的近似解。

适用范围广：高斯-赛德尔迭代法可以适用于各种线性方程组，无论是大型稀疏矩阵还是小型稠密矩阵都可以得到较好的结果。

计算量较小：相对于其他迭代方法，高斯-赛德尔迭代法的计算量较小，因此在计算时间和内存消耗方面具有优势。

可以处理稀疏矩阵：高斯-赛德尔迭代法可以有效地处理稀疏矩阵，可以减少计算时间和内存消耗。

缺点：

对系数矩阵的条件有要求：高斯-赛德尔迭代法的收敛性与系数矩阵的条件有关，如果系数矩阵不是正定矩阵或其条件数太大，可能会导致迭代不收敛或收敛速度慢。

需要选择适当的迭代参数：高斯-赛德尔迭代法的收敛速度和计算精度与初始迭代参数的选择有关，如果参数选择不当，可能会导致迭代不收敛或收敛速度慢。

对计算机内存要求较高：当处理大型稀疏矩阵时，高斯-赛德尔迭代法需要占用大量的内存空间来存储矩阵和向量，因此对计算机内存要求较高。

import numpy as np
# A为系数矩阵，b为结果矩阵
# k为最大迭代次数，tol为容差限
def Guauss_seidel(A,b,k,tol):
    n = A.shape[1]
    D = np.eye(n)
    D[np.arange(n),np.arange(n)] = A[np.arange(n),np.arange(n)]
    LU = A - D
    # 得到LU的上、下三角矩阵
    L = np.tril(LU)
    U = np.triu(LU)
    D_L = D + L
    X = np.zeros(n)

    # 迭代k次
    for i in range(k):
        X1= X.copy()
        D_L_inv = np.linalg.inv(D_L)
        X = -np.dot(np.dot(D_L_inv,U),X) + np.dot(D_L_inv,b)

        if (np.max(np.abs(X-X1)))<=tol:
            print('计算收敛!')
            break
        print('n第',i,'次迭代：',X)
        print('残差值为：',np.abs(X-X1))
        print('当前最大残差为：', np.max(np.abs(X-X1)))
    return X

A = np.array([
    [3000,-1000,0,0],
    [-1000,2000,-1000,0],
    [0,-1000,2000,-1000],
    [0,0,-1000,3000],
])
b = np.array([202500,2500,2500,802500])
X = Guauss_seidel(A, b, 1000,0.005)

n第 17 次迭代： [139.99560885 217.49300459 292.49490235 364.99830078]
残差值为： [0.00387807 0.00617804 0.00450202 0.00150067]
当前最大残差为： 0.006178036008577692
计算收敛!

老歌老听老掉牙
关注
  • 0
    点赞
  • 踩
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
fipy:FiPy 是一个用 Python 编写的有限体积 PDE 求解器
07-24
概述 .. only:: latex :term:`FiPy` is an object oriented, partial differential equation (PDE) solver, written in :term:`Python`, based on a standard finite volume (FV) approach. The framework has been developed in the Materials Science and Engineering Division (MSED_) and Center for Theoretical and Computational Materials Science (CTCMS_), in the Material Measurement Laboratory (MML
GS迭代_高斯赛德尔迭代法_稳态迭代_
09-29
一维稳态导热问题的高斯赛德尔迭代法求解方法
python matplot包绘制有限体积解
Tankrun1997的博客
03-03 397
python matplotlib包绘制非结构(三角)网格有限体积解
三维导热问题的ADI-TDMA算法
weixin_41654618的博客
04-30 3224
本文在上一篇文章中关于二维非稳态导热问题的交替方向隐式方法的基础上，设计了一个三维算例，并通过Fortran程序，演示适用于三维模型的Brian ADI格式的交替方向隐式方法的求解过程。
Python-FiPy有限体积安装
liuyhoo
08-11 1770
FiPy是基于标准有限体积（FV）方法，用Python编写的面向对象的偏微分方程（PDE）求解器。通过Pygmsh或gmsh对求解域进行网格划分，然后基于Fipy求解偏微分方程。
利用有限元法求解一维热传导问题
TATANICYELLOW的博客
04-01 1456
首先是进行网格划分，网格个数当然是越多越好，但是求解的时间也会相应大大增加，比如在下方程序中因为我把网格数设置为了500，相应的刚度矩阵就是501行501列，再加上我将时间步长设置较小，往往运行一次需要几分钟时间（不过我也有在考虑利用pytorch进行GPU加速来减小运行时间，希望可行吧）。实际上，运用有限元更简单的例子可能是稳态情况下的一维热传导问题，但是，笔者受外界不可抗力的影响，先上手的是非稳态的问题（不过稳态也没什么难度）。先写这些吧，接下来的一周要着手完成二维热传导的问题，祝我好运！
二维泊松_泊松方程_poisson_二维泊松_
10-03
二维泊松方程是偏微分方程领域中一个重要的基本方程，广泛应用于电磁学、流体力学、热传导、量子力学等多个物理领域。它描述了在一个区域内，某个物理量（如电势、温度或密度）的梯度与该物理量的负二阶导数之间的...
pde显示求解_偏微分方程显式求解_
09-29
在计算机科学和数值分析领域，偏微分方程（PDE）的显式求解是一种常用的数值...在实际应用中，可能需要结合隐式方法或半隐式方法以克服稳定性限制，或者采用更高级的算法如有限元法或有限体积法来提高解决方案的质量。
PyFluids:用于在各种情况下模拟和分析流体流动的Python软件包
03-25
2. **数值方法**：软件包采用了常见的数值求解技术，如有限差分法、有限体积法和有限元法，保证了模拟的精度和稳定性。 3. **边界条件处理**：PyFluids允许用户设定各种边界条件，如无滑移壁、自由表面、对流边界...
Diffusion_计算物理_二维扩散_
09-29
这一现象广泛存在于自然界，如热传导、化学反应、粒子迁移等。"Diffusion_计算物理_二维扩散_"这个主题着重讨论如何利用数学模型和数值方法解决二维扩散方程。 二维扩散方程是偏微分方程的一种，它描述的是在二维...
实战十五：差分法解热传导方程
03-15
使用Python写用差分法解热传导方程的程序
有限体积法源代码
02-28
有限体积法 源代码 包括算例 网格刨分 各种方法
热传导各向异性扩散降噪（python代码实现）
04-24
各向异性扩散滤波主要是用来平滑图像的，克服了高斯模糊的缺陷，各向异性扩散在平滑图像时是保留图像边缘的，和双边滤波很像。 通常我们有将图像看作矩阵的，看作图的，看作随机过程的，记得过去还有看作力场的。这次新鲜，将图像看作热量场了。每个像素看作热流，根据当前像素和周围像素的关系，来确定是否要向周围扩散。比如某个邻域像素和当前像素差别较大，则代表这个邻域像素很可能是个边界，那么当前像素就不向这个方向扩散了，这个边界也就得到保留了。
非稳态导热问题有限体积法
05-12
为了研究非稳态导热问题中的常用有限体积方案(FVM1)和新提出的有限体积方案(FVM2)的精度差异问题,采用理论推导和实例验证的方法,建立基于有限体积法原理的热平衡积分方程,利用矩形网格进行离散化,建立了平面非稳态热传导问题有限体积法数学模型,推导出与有限单元法(FEM)方案相类似的线性方程组,对比分析所建立的各个方程组中参数的差异.研究结果表明:方案FVM2与FEM精度相同,为较优方案.最后结合一个典型实例,用VB编制了数值计算软件,对比分析这两种方案,验证了上述结果.该结果对工程计算数值模拟有一定的参考价值.
manapy：并行二维非结构化浅层方程
02-20
- 解决方案可能涉及数值方法，如有限元法（FEM）、有限体积法（FVM）或有限差分法（FDM），这些方法都是非结构化网格上的常用求解策略。 - 用户可以自定义边界条件、初始条件，并对结果进行后处理，如可视化。 7....
有限体积PDE求解器--基于Python的FiPy
stereohomology
10-22 1万+
NIST Home > MML > CTCMS >  FiPy > FiPy: A Finite Volume PDE Solver Using Python Version 3.1 Cahn-Hilliard Phase Field Level Set Questions? Suggestions? J
python+julia+matlab：有限体积法计算二维流场温度场
最新发布
走向CTO的路上...
01-05 487
它将流场划分为离散的控制体积，通过对控制体积进行质量、动量和能量守恒方程的积分，得到离散的代数方程组。这些方程组描述了控制体积内部的物理过程，通过求解这些方程组，可以得到流场的各种属性，如流速、压力和温度分布。它可以用于解决各种流动和传热问题，例如流体流动、热传导、对流传热、辐射传热等。具体的应用场景包括但不限于：空气动力学研究、工程流体力学、燃烧模拟、热交换器设计、风扇和涡轮机设计等。给定入口的初始温度和流速，固体为热源，经过流体冷却后，出口为压力出口边界，采用SIMPLE算法求解出口的温度和流速。
python：有限体积法温度求解器
走向CTO的路上...
12-26 984
它将求解区域分解成离散的控制体（体积单元），通过对控制体进行质量、能量和动量的守恒方程进行积分和离散化，得到代表守恒量平衡的方程组，进而通过迭代求解得到数值解。1.如何用有限体积法求解非稳态问题，求解一个长方体高温物体在空气中冷却的三维的温度场，物体上表面和侧面的边界条件是对流换热和辐射，物体底面的边界条件是和地面的导热，假设散热时长是一小时，并用python编写一个可以运行的实例，最后能输出这块物体的温度分布彩虹图。对于给定的问题，如长方体高温物体在空气中冷却的三维温度场，可以采用有限体积法进行求解。
matlab求解一维热传导问题
热门推荐
qq_38870183的博客
01-18 1万+
Matlab 求解一维热传导方程 %%传热模型 clear ; g = ‘squareg’; %定义以为正方形 b = ‘squareb1’; %定义边界条件为零 c = 1; a = 0; f = 0; d = 1; %产生初始三角网格 [p,e,t]=initmesh(g); %定义初始条件 u0 = ze...
python 使用有限差分求解一维热传导方程
05-13
一维热传导方程可以表示为： $$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$$ 其中，$u(x,t)$ 表示温度分布，$\alpha$ 为热扩散系数。 使用有限差分方法求解一维热传导方程可以分为以下几个步骤： 1. 离散化：将时间和空间分别离散化为 $t_i = i \Delta t$ 和 $x_j = j \Delta x$，其中 $\Delta t$ 和 $\Delta x$ 分别表示时间和空间的步长。 2. 近似：将偏导数用差分来近似，例如： $$\frac{\partial u}{\partial t} \approx \frac{u_{j,i+1} - u_{j,i}}{\Delta t}$$ $$\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{j+1,i} - 2u_{j,i} + u_{j-1,i}}{\Delta x^2}$$ 3. 替换：将近似后的偏导数代入原方程，得到差分方程： $$u_{j,i+1} = u_{j,i} + \alpha \frac{\Delta t}{\Delta x^2} (u_{j+1,i} - 2u_{j,i} + u_{j-1,i})$$ 4. 边界条件：确定边界条件，例如固定边界条件 $u(0,t) = u(L,t) = 0$。 5. 时间推进：使用差分方程逐步推进时间，例如使用显式欧拉法： $$u_{j,i+1} = u_{j,i} + \alpha \frac{\Delta t}{\Delta x^2} (u_{j+1,i} - 2u_{j,i} + u_{j-1,i})$$ 6. 初始条件：确定初始条件，例如 $u(x,0) = f(x)$。 下面是一个 Python 实现的示例代码： ```python import numpy as np import matplotlib.pyplot as plt # 设置参数 L = 1.0 # 空间长度 T = 0.1 # 时间长度 alpha = 1.0 # 热扩散系数 N = 100 # 空间网格数 M = 100 # 时间网格数 dx = L / N # 空间步长 dt = T / M # 时间步长 # 初始化温度分布 u = np.zeros((N+1, M+1)) u[0,:] = u[N,:] = 0.0 u[:,0] = np.sin(np.pi*np.arange(N+1)/(N+1)) # u[:,0] = np.exp(-100*(np.arange(N+1)/(N+1) - 0.5)**2) # 高斯分布 # 使用有限差分方法求解热传导方程 for i in range(M): for j in range(1, N): u[j,i+1] = u[j,i] + alpha*dt/dx**2 * (u[j+1,i] - 2*u[j,i] + u[j-1,i]) # 绘制温度分布曲线 x = np.linspace(0, L, N+1) t = np.linspace(0, T, M+1) X, T = np.meshgrid(x, t) fig = plt.figure(figsize=(8,6)) ax = fig.add_subplot(1,1,1,projection='3d') ax.plot_surface(X, T, u.T, cmap='coolwarm') ax.set_xlabel('x') ax.set_ylabel('t') ax.set_zlabel('u') plt.show() ``` 其中，初始化温度分布的代码可以根据实际情况进行修改。绘制温度分布曲线的代码使用了 `matplotlib` 库中的 `plot_surface` 函数。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值