用Matlab求解一维非稳态周期性导热问题（有限单元法+隐式离散+高斯赛德尔迭代法）

Silver_777

已于 2023-08-03 16:50:49 修改

阅读量6k

点赞数 11

分类专栏： Matlab&CFD 文章标签： matlab 算法

于 2022-05-28 12:48:45 首次发布

本文链接：https://blog.csdn.net/Silver_777/article/details/125016866

版权

Matlab&CFD 专栏收录该内容

4 篇文章 4 订阅

订阅专栏

一、问题描述与分析

本次问题条件如下：

计算模拟如下一维常物性无内热源非稳态导热的温度场，以及内外壁面的热流密度，并进行温度场和热流的特点分析，相关参数如下。

室内温度恒定为20℃，无限大平壁的比热为1.0kJ/kg℃，室外气温周期性变化，周期为1天（即24小时）。

二、网格划分与节点方程式

由问题所给条件可知，两侧边界条件均为对流边界条件，室内侧温度恒定，室外侧温度周期性变化，无限大平壁内部没有热源，将该无限大平壁沿厚度方向均匀的划分为300个单元，左侧室内右侧室外，节点编号依次为1-2-3-...-300-301：

则该一维常物性无内热源非稳态导热问题的有限单元节点方程如下（注这里不是有限差分法，是有限单元法，离散格式为隐式）：

左端节点方程

右端节点方程

中间节点方程

其中，i表示空间维度节点，t表示时间维度节点。上述公司中，左端是t时刻的节点温度，是未知量需要求解的，右端是t-1时刻的节点温度，是已知的。这部分的公式要注意和后面高斯赛德尔迭代进行区分，下一节的迭代次数专门用了不同字母表示，就是为了防止非稳态传热的时间递进与高斯赛德尔本身的迭代递进弄混。

三、高斯赛德尔迭代法

上述节点方程式可以通过高斯赛德尔迭代法求解。节点方程式本质是一个n+1元的线性方程组（假定划分单元数为n）：

$\begin{cases} a_{1,1}x_1 + a_{1,2}x_2 + \cdots + a_{1,n}x_{n+1} = b_1 \\ a_{2,1}x_1 + a_{2,2}x_2 + \cdots + a_{2,n}x_{n+1} = b_2 \\ \cdots \\ a_{n+1,1}x_1 + a_{n+1,2}x_2 + \cdots + a_{n+1,n+1}x_{n+1} = b_{n+1} \end{cases}$

其矩阵形式如下：

$\begin{pmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n+1} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n+1} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n+1,1} & a_{n+1,2} & \cdots & a_{n+1,n+1} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n+1} \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_{n+1} \end{pmatrix}$

根据二中总结的节点方程式，大部分系数项是0，要求解k时刻下n+1个节点处的温度，则有：

$\begin{pmatrix} a_{1,1} & a_{1,2} & 0 & 0 & \cdots & 0 & 0 \\ a_{2,1} & a_{2,2} & a_{2,3} & 0 &\cdots & 0 & 0 \\ 0 & a_{3,2} & a_{3,3} & a_{3,4} &\cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & a_{n,n+1} & a_{n+1,n+1} \end{pmatrix} \begin{pmatrix} T_1 \\ T_2 \\ \vdots \\ T_{n+1} \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_{n+1} \end{pmatrix}$

系数的表示由下图可得，左右两边的系数均是已知的，且矩阵对角线系数部位0，因此n+1个方程求解n+1个节点的温度是可求解的。

高斯赛德尔迭代法的迭代公式为：

$T_i^{(k+1)}=\frac{1}{a_{i,i}}\left(b_i-\sum_{j=1}^{i-1}a_{i,j}T_j^{(k+1)}-\sum_{j=i+1}^{n+1}a_{i,j}T_j^{(k)}\right) \\ \quad i=1,2,\cdots,n+1$

注意，这里T的上标k不是时间，是迭代次数，k表示第k次迭代，k+1表示第k+1次迭代。由公式可以看到，高斯赛德尔迭代法同时利用了刚迭代出来k+1次数据以及上一次迭代的k次数据，因此收敛得更快。和大部分迭代法一样，它需要使用一个初始解去尝试，求解出来的T会一步一步逼近正确解，因此要设置一个误差范围，当相邻两次解差值小于设定值时，则认为已经收敛到足够精度的解。
每个t 时刻下，进行多轮迭代求出该时刻下的温度场近似解，然后根据节点方程式进行t+1时刻的求解，这样逐时把所有时刻的温度场求解出来。因此在代码编写中，有两个循环，需要注意区分，外层循环是时间循环，是对非稳态进程的求解，内层循环是高斯赛德尔迭代，是对具体一个时刻温度场的求解。

四、在Matlab上求解

1. 前处理

这部分主要是前面的参数设置，同时定义了该无限大平壁的初始温度（初始温度场会影响影响初期的图像，所以为了好看设置了一个接近室内外环境的梯度温度场）。计算周期day设置成7天，时间不长为60s，这两个参数对程序运行时间的影响比较大，我这里是为了能计算更准以及看出周期性变化，计算的时间也相当长，实际运行的时候可以进行适当调整。

在本问题中，室外温度是周期性变化的，提前计算好以便后面循环调用。已知的A系数矩阵也进行了计算。求解对象T被设置成了一个维度为（n+1,Tn）的二维数组，列向量代表时间维度，行向量表示空间维度。

本模块相关参数定义如下：

命名	意义	值	单位
dx	空间步长	0.001	m
i	空间下标
n	厚度划分数目	300
dt	时间步长	60	s
t	时间下标
day	计算天数	7	天
Tn	时间划分数目	602460*day/dt
lambda	平壁导热系数	2	W/(m*℃)
rou	平壁密度	2000	kg/m^3
c	平壁比热容	1000	J/kg*K
h1	室内侧对流换热系数	10	W/(m^2*℃)
t1	室内侧空气温度	20	℃
h2	室外侧对流换热系数	20	W/(m^2*℃)
t2_0	室外侧空气温度定值部分	30	℃
Aw	室外侧空气温度变化系数	15	℃
t2	室外侧空气温度
ep	求解精度	0.000001

%% 前处理

% 室内侧
h1 = 10; % 对流换热系数 W/(m^2*K)
t1 = 20; % 空气温度 

% 室外侧
h2 = 20; % 对流换热系数 W/(m^2*K)
t2_0 = 30; % 空气温度定值部分 K
Aw = 15; % 空气温度变化系数

% 无限大平壁
D = 0.3; % 厚度 m
c = 1000; % 比热容 J/kg*K
lambda = 2; % 导热系数 W/(m*K)
rou = 2000; % 密度 kg/m^3

% 网格划分
dx = 0.001; % 空间步长 m
n = D/dx; % 厚度划分数目
dt = 60; % 时间步长（10分钟） s
day = 7; % 所计算的天数
Tn = 60*24*60*day/dt; % 所计算时间节点的数目
ep = 0.000001; % 求解精度
MAX = 20000; % 最大求解次数

% 室外侧空气温度周期性变化设定
for i = 1:Tn
    t2(i) = t2_0 + Aw * cos(2 * pi * i * dt / 24 / 60 / 60 );
end

A = zeros(n+1,n+1); % 系数矩阵
B = zeros(n+1,1); % 系数矩阵
T = zeros(n+1,Tn); % 所有时刻的平壁温度场

for i = 1:n+1
    T(i,1) = 20 + i*10/n; % 初始温度
end

% 定义线性方程组系数矩阵
A(1,1) = lambda/dx + h1 + rou*c*dx/3/dt;
A(1,2) = rou*c*dx/6/dt - lambda/dx;
A(n+1,n) = rou*c*dx/6/dt - lambda/dx;
A(n+1,n+1) = lambda/dx + h2 + rou*c*dx/3/dt;
for i = 2:n
    A(i,i-1) = rou*c*dx/6/dt - lambda/dx;
    A(i,i) = 2*lambda/dx + 2*rou*c*dx/3/dt;
    A(i,i+1) = rou*c*dx/6/dt - lambda/dx;
end

2. 高斯赛德尔迭代法求解

从K=2时刻开始，采用高斯赛德尔迭代法求解上述矩阵形式的线性代数方程组。这里需要注意的地方是，高斯赛德尔迭代法能相对更快速收敛的地方就在于它每次计算都使用了最新解，因此计算时需要按顺序进行；同时，在高斯赛德尔迭代求解中，一开始也需要提供一个任意初试解T_prev，以便推动方程组的迭代求解，这个地方要注意同初始时刻温度场区分，因为高斯赛德尔的初始解在每个时刻下的循环中都会用到，而初始温度场只有在最开始时刻的计算中会用到。

对于内层循环，每个Iteration计算完后，T_prev存储该次迭代结果，进入下一次iteration后，新的数据实时更新在T数组中，同时T_prev上次迭代的数据也被用到。高斯赛德尔迭代方程式的正确表达非常重要，这个地方很容易出错，之前的代码也是因为这部分混淆所以温度场算的有些问题，在本次修改中已经全部更正。

当error满足误差要求时，结果基本收敛，t 时刻下的温度场求解完成，跳出内层循环。在外层循环中，从头进行t+1时刻的计算，以此完成所有时刻的温度场求解。

%% 求解温度场
for t = 2:Tn

    T_prev = 20*ones(n+1,1); % 每个k时刻下的试错初始解

    % 高斯赛德尔迭代
    for iter = 1:MAX
        error = 0;
        for i = 1:n+1
            if i == 1
                sum_A = A(i,i+1:n+1) * T_prev(i+1:n+1);
                B(i) = h1*t1 + (rou*c*dx/dt/6) * (2*T(1,t-1)+T(2,t-1));
                T(i,t) = (B(i) - sum_A) / A(i,i);
            elseif i == n+1
                sum_A = A(i,1:i-1) * T(1:i-1,t);
                B(i) = h2*t2(t) + (rou*c*dx/dt/6) * (2*T(n+1,t-1)+T(n,t-1));
                T(i,t) = (B(i) - sum_A) / A(i,i);
            else
                sum_A = A(i,1:i-1) * T(1:i-1,t) + A(i,i+1:end) * T_prev(i+1:end);
                B(i) = (rou*c*dx/dt/6) * (T(i-1,t-1)+4*T(i,t-1)+T(i+1,t-1));
                T(i,t) = (B(i) - sum_A) / A(i,i);
            end
            error = error + abs(T(i,t) - T_prev(i));
        end
        if error < ep
            break; % 当解收敛时退出迭代
        end
        T_prev = T(:,t);
    end
end

3. 作图

在画图的时候有个小tips，因为时间步长小，时间节点非常多，会导致网格显示不清，因此这里是选取每个小时的数据显示。

%% 绘制温度场随时间的变化图像

%作图
X = zeros(1,n+1);
Y = zeros(1,24*day);
TT = zeros(24*day,n+1);
 
for i = 1:24*day
    Y(i) = 60*(i-1)+2; % 60分钟一个数
    for j = 1:n+1
        X(j) = j;
        TT(i,j) = T(X(j),Y(i));
    end
end
 
figure(1)
hold on;
[XX,YY] = meshgrid(X,Y);
surfc(YY,XX,TT)
xlabel('Time - AXIS')
ylabel('Thickness - AXIS')
zlabel('T (^{\circ}C)')
title(colorbar,'^{\circ}C')

4. 运行时间

如上参数设置， CPU是AMD Ryzen 7 4800H，跑了三个多小时。每一个时刻都需要成百上千次的高斯赛德尔迭代计算，因此减少时间节点数目是有效的方法。

四、结果分析

1. 温度场可视化

2. 表面热流密度分析

根据牛顿冷却公式

由于已经求解出温度场，可根据上式计算本问题中的无限大平壁各个时刻下的内外表面热流密度，在Excel中对内外壁面的热流密度进行分析，下图为内外壁面热流密度、室内外空气温度、平壁内外表面壁温随时间的变化曲线。

室内外壁温、室外侧空气温度的关系：

室外侧壁温、室外侧热流密度、空气温度的关系：

室内侧壁温、室内侧热流密度、空气温度的关系

五、问题说明

由上述求解结果可以推断得到

室外温度周期性变化时，无限大平壁外表面壁温和热流密度的变化也呈现出周期性的特点，三者震荡频率一致；
室外温度周期性变化时，无限大平壁外表面壁温波比空气温度波落后一些，具有一定的延迟性，表面热流密度的震荡是壁温波和空气温度波综合作用的结果；
无限大平壁外表面的温度震荡最明显，越靠近内侧，震荡幅度越小，可见温度波自外向内发生了衰减，这是由于平壁材料对温度波的阻尼作用。

Silver_777

关注

11
点赞
踩
83

收藏

觉得还不错? 一键收藏
3
评论
用Matlab求解一维非稳态周期性导热问题（有限单元法+隐式离散+高斯赛德尔迭代法）

使用有限单元法计算一维常物性无内热源非稳态导热的温度场，以及内外壁面的热流密度，并进行温度场和热流的特点分析。室内温度恒定为20℃，无限大平壁的比热为1.0kJ/kg℃，室外气温周期性变化，周期为1天（即24小时）。二、网格划分与节点方程式由问题所给条件可知，边界条件为对流边界条件，则该一维常物性无内热源非稳态导热问题的有限单元节点方程如下（注这里不是有限差分法，是有限单元法，离散格式为隐式）：.
复制链接

扫一扫

专栏目录