精华内容
下载资源
问答
  • 一维非稳态导热热传导Matlab程序,通过此可以解决偏微方程
  • 精品文档 第五次作业前三题写在作业纸上 一用有限差分方法求解一维非定常热传导方程 初始条件和边界条件见说明 .pdf 文件 热扩散系数 =const 2 T T 2 t x 1. 用 Tylaor 展开法推导出 FTCS格式的差分方程 2. 讨论该...
  • 采用理论推导和实例验证的方法,建立基于有限体积法原理的热平衡积分方程,利用矩形网格进行离散化,建立了平面非稳态热传导问题有限体积法数学模型,推导出与有限单元法(FEM)方案相类似的线性方程组,对比分析所建立的...
  • 传热学第四版 高等教育出版社 杨世铭 陶文铨 编著 第四章非稳态热传导例题4-6的数值解法的C++源文件,老师要做的大作业,结果可以做出例题里面的等温线图,给以后要做大作业的同学提供个参考
  • 个围绕偏微分方程的设立求解的数学模型。赛题连接 声明:本文所有图片均出自论文、未经允许,谢绝使用。 思路 如何建立个室内温度的控制模型呢?当时刚看到这个赛题的时候,我首先想到的是各种物理公式,传热...

    前言

    本文为东三省数学建模B赛题,室温调控模型的建模以及仿真思路分享。是一个围绕偏微分方程的设立求解的数学模型。赛题连接
    声明:本文所有图片均出自论文、未经允许,谢绝使用。

    思路

    如何建立一个室内温度的控制模型呢?当时刚看到这个赛题的时候,我首先想到的是各种物理公式,传热学,对流换热,牛顿冷却定律什么的,这些传热学的定律公式确实是非常有用的,但是还需要温度分布图,这就没办法单单只靠这些公式完成。那么存在热源的一个温室,他的温度分布用什么样的模型来搭建,偏微分方程解决了这个问题。偏微分方程在这个模型中简单来说就是是用来求热量传递的时候,在一段时间后热量如何分布的。一维的偏微分方程就是求一条线上的温度分布,二维就是求一个面上的温度分布。当然三维的偏微分方程就是求室内的温度分布。但是这里面临一个问题,三维的偏微分方程的数值解是非常难解的,为了计算的方便,可以采用二维的偏微分方程求解南北墙面的温度分布,再利用南北墙面的温度作为底面的一个边界条件来求XOY面的一个温度分布。这个就能够大概的描述整个室内的一个温度分布情况。作为有热源的情况还得单独设立热源。并且要做出许多假设,各个边界是绝热条件还是不绝热,这都是问题,具体情况具体分析。

    建模假设

    针对模型,我做了如下的假设:
    1:暖气片可以完全转化热量,进水口温度就是暖气片温度;
    2:暖气片散热方式以对流传热为主;
    3:东西墙绝热,没有热量散失;
    4:不考虑热辐射的影响;
    5:房间初始温度与室外一致;
    6:假设各个墙体内不含内热源。

    偏微分方程的matlab求解

    这部分我当时也有很大的困惑,不明白怎么将这个和温度联系起来,怎么和室外温度联系起来。后来我发现,是我对偏微分方程的理解不够全面,可能现在也还是。因为我不是数学专业的,这些数学原理我也不太懂。我说下我的理解:求解温度分布,简单来说就是把二维墙面的温度分布设置为一个矩阵,矩阵上的每一个点就代表墙面上一个点的温度,然后建立二维偏微分方程,求解它的数值解,就能求出在一段时间后这个墙面的温度分布矩阵。二维偏微分方程的求解用的是有限差分法,网上也有很多教程,这边对于求解就不具体说了。说下简单的步骤:

    • 设出三个方向传递的A矩阵
    • 设出墙面的初始温度分布矩阵
    • 设定热源矩阵
    • 设置边界条件
    • 通过有限差分法,迭代求解偏微分方程。
    • 画出图像

    matlab仿真画图

    给出第一问的matlab程序,可以复制粘贴运行试试。但是求解偏微分方程是给很大的计算量,通常要计算非常的久(看时间间隔怎么取的),通常我为了画出一个图,跑程序就得跑个把小时,甚至更久。

    % 定义变量
    C_w = 4200;             % 定义水的比热容
    lou_w = 1000;           % 定义水的密度
    V_w = 0.5;              % 定义水流量
    m_w = lou_w*V_w;        % 定义水的质量
    K_s = 50;               % 暖气片传热系数
    F_s = 2;                % 暖气片对流换热面积
    T_h = 40;               % 暖气水温
    T_out = 0;              % 室外温度
    sigma1 = 0.25;          % 墙体厚度
    lamda1 = 1.2;           % 墙体传热系数
    lamda2 = 5.7;           % 玻璃传热系数
    h1 = 0.9;               % 内表面传热系数
    h2 = 1.6;               % 外表面传热系数
    dx = 0.05;              % 给定差分变量取值范围
    Lx = 5;                 % 定义x轴距离
    x = 0:dx:Lx;            % 差分密度
    Lxx = 5;                % 同上
    xx = 0:dx:Lxx;
    dz = 0.05;
    Lz = 2.8;               % 同上
    z = 0:dz:Lz;                
    dy = 0.05;
    Ly = 10;                % 同上
    y = 0:dy:Ly;
    time = 360;             % *10即为实际时间
    dt = 0.0001;
    t = 0:dt:time;
    
    Ax = (-2*eye(length(x))+diag(ones(1, length(x)-1), 1)+diag(ones(1, length(x)-1), -1));  % 定义X方向的A矩阵
    Ay = (-2*eye(length(y))+diag(ones(1, length(y)-1), 1)+diag(ones(1, length(y)-1), -1));  % 定义Y方向上的A矩阵
    Az = (-2*eye(length(z))+diag(ones(1, length(z)-1), 1)+diag(ones(1, length(z)-1), -1));  % 定义Z方向上的A矩阵
    
    [Z, X] = meshgrid(z, x);  
    [Y, XX] = meshgrid(y, xx);
    
    U0xoz = ones(length(Ax), length(Az))*0;  % 定义南北墙初始温度为0℃
    U0xoy = ones(length(Ax), length(Ay))*0;  % 定义xoy面初始温度0℃
    
    fxozs = ones(length(Ax), length(Az));
    fxozs(20:90, 1:25) = 40;                 % 定义热源(仅南墙暖气作用)
    fxozn = ones(length(Ax), length(Az));
    fxozn(20:90, 1:25) = 30;                 % 定义热源(仅北墙暖气作用)
    fxoy1 = 5*exp(-1/2*(((XX-2.5).^2)/1+((Y-8.75).^2)/1));
    fxoy2 = 5*exp(-1/2*(((XX-2.5).^2)/1+((Y-1.25).^2)/1));
    
    border_n = 8*exp(-1/2*(((z-1).^2)/1));       % 定义边界条件
    border_nx = 10*exp(-1/2*(((x-2.5).^2)/1));   % 定义边界条件
    Uxozs = U0xoz;
    Uxozn = U0xoz;
    Uxoy = U0xoy;
    a = 0.22;                                    % 定义散热系数
    
    xm = zeros(1, time*10);                      % 控制变量的定义
    xm_in_s = zeros(1, time*10);
    xm_in_n = zeros(1, time*10);
    xm_in_xoy = zeros(1, time*10);
    
    temp = 1;                                    % 控制变量的定义 
    temp_in_s = 1;
    temp_in_n = 1;
    temp_in_xoy = 1;
    
    Q_w = 0;                                     %  墙体导热初值设置
    Q_s = 0;                                     % 对流换热初值设置
    
    for n = 1:length(t)-1                        % 迭代时间,计算热传导偏微分方程
        Uxozs = Uxozs+a^2*(1/dx^2*Ax*Uxozs+1/dz^2*Uxozs*Az+fxozs )*dt;   % 南墙的温度分布
        min_xoz_s = mean(Uxozs,2);      		 % 计算南墙平均温度分布,压缩为一维向量作为xoy面的边界条件
        Uxozs(:, 1) = border_nx;             	 % 边界设置
        Uxozs(:, end) = Uxozs(:,end-1);
        Uxozs(1, :) =  border_n;
        Uxozs(end, :) = border_n;
        
        Uxozn = Uxozn+a^2*(1/dx^2*Ax*Uxozn+1/dz^2*Uxozn*Az+fxozn)*dt;   % 北墙的温度分布
        min_xoz_n = mean(Uxozn,2);      		% 计算北墙平均温度分布,压缩为一维向量
        Uxozn(:, 1) = border_nx;                % 边界设置
        Uxozn(:, end) = Uxozn(:,end-1);
        Uxozn(1, :) =  border_n;
        Uxozn(end, :) = border_n;
        
        Uxoy = Uxoy+a^2*(1/dx^2*Ax*Uxoy+1/dy^2*Uxoy*Ay+fxoy1+fxoy2)*dt;  % xoy墙的温度分布
        Uxoy(:, 1) =min_xoz_s;                	% 绝缘设置
        Uxoy(:, end) = min_xoz_n;
        Uxoy(1, :) =  Uxoy(2, :);
        Uxoy(end, :) = Uxoy(end-1, :);
        
        min_in_s = mean(Uxozs(:));     % 南墙的均值
        min_in_n = mean(Uxozn(:));     % 北墙的均值
        min_in_xoy = mean(Uxoy(:));    % xoy的均值
        %------------------------------------------回水温度的计算--
        Q_w = ((min_in_s'-T_out)/((1/h1)+(sigma1/lamda1)+(1/h2)))*5*2.8*2*3600;  % 墙体换热
        Q_s = K_s*F_s*(T_h-min_in_s')*3600;   % 对流换热第一个暖气片
        Q_t = Q_w+Q_s;                        % 计算总能量
        T_x = T_h-(Q_t/(C_w*m_w));            % 计算回水温度
        fxozn(20:90, 1:25) = T_x;             % 将回水温度设置为北墙的热源
        xm(1, temp) = T_x;                    % 将回水温度存入向量,之后绘图
        xm_in_s(1, temp_in_s) = min_in_s;     % 保存数据画图
        xm_in_n(1, temp_in_n) = min_in_n;
        xm_in_xoy(1, temp_in_xoy) = min_in_xoy;
        
        if mod(n,1000) == 1                   % 加速画图
            temp = temp + 1;
            temp_in_s = temp_in_s + 1;
            temp_in_n = temp_in_n + 1;
            temp_in_xoy = temp_in_xoy + 1;
            subplot(311);
            surf(X,Z,Uxozs);                  % 绘制南墙温度分布
            shading interp;
            axis([x(1) x(end) z(1) z(end) 0 40]);
            title('南墙温度分布')
            xlabel('x/m');
            ylabel('z/m');
            zlabel('T/℃');
            subplot(312);
            surf(X,Z,Uxozn);                  % 绘制北墙温度分布
            shading interp;
            axis([x(1) x(end) z(1) z(end) 0 40]);
            title('北墙温度分布');
            xlabel('x/m');
            ylabel('z/m');
            zlabel('T/℃');
            subplot(313);
            surf(XX,Y,Uxoy);                  % 绘制xoy墙温度分布
            shading interp;
            axis([xx(1) xx(end) y(1) y(end) 0 40]);
            title('xoy面温度分布');
            xlabel('x/m');
            ylabel('y/m');
            zlabel('T/℃');
    %     view(2);  % 改变视角方向
            getframe;
        end
    end
    % 绘图
    x = 1:1:time*10;
    figure(4);
    plot(x, xm(1, 1:time*10));
    legend('出水口温度')
    title('暖气回水温度与时间的关系');
    xlabel('时间/s');
    ylabel('回水温度/℃');
    axis([x(1) x(time*10) 0 40]);
    figure(5);
    plot(x, xm_in_s(1, 1:time*10));
    hold on
    plot(x, xm_in_n(1, 1:time*10));
    hold on
    plot(x, xm_in_xoy(1, 1:time*10));
    legend('南墙平均温度', '北墙平均温度',  '底面xoy面平均温度')
    axis([x(1) x(time*10) 0 30]);
    title('各个面温度随时间变化曲线');
    xlabel('时间/s');
    ylabel('温度/℃');
    
    

    在这里插入图片描述
    左边为十分钟左右后的三个面的温度分布,右边为一小时后的三个面的温度分布
    在这里插入图片描述
    另一个角度
    在这里插入图片描述
    在这里插入图片描述

    心得分享

    对于时间还算比较充分的东北三省联赛数学建模,时间上是完全够用的。关于分工,打数学建模一定要明确自己的定位,分工明确。在刚开始的时候可以不用那么分工明确,因为三个人都没有思路,就是一起讨论,用什么模型,用什么方法。之后有了一定的方向后再开始分工。对于这个赛题,我们最初也走了一些弯路,没有想到偏微分方程的方法,第一眼咋一看以为是道物理题,就全往物理方面想了,什么热传导,牛顿冷却定律,传热学,对流换热什么的都看了下,后来明白数学建模怎么可能没有数学方法。但是这也不算弯路,这些东西对后面建模也有重要作用。后来就开始研究偏微分方程,主要是研究怎么用上偏微分方程,然后一步步的把模型建起来。关于时间安排,留给论文排版的时间是绝对要有了,如果没有写论文的经历的话那更要有了,通常我们写的论文是非常不合格的,就单单改论文就要改很长时间。其他时间安排就看自己的进度了,论文这里是肯定要留时间的,不能写个初稿就完事了。如果到了后期,有些分工任务弄完的,比如建模代码写完的,图画完的,仿真弄完的,就可以帮着写论文的同学,帮忙把公式用word写出来,好让他直接能复制粘贴,帮忙看看语言逻辑通不通顺,错别字有没有等等。就是不能闲着,不是自己的任务做完就完事了。
    如果你觉得文章对你有一点点帮助或启示的话,不妨点个赞噢😊

    展开全文
  • 维非稳态导热

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

    计算题目

    计算求解在有内热源和无内热源情况下,固体材料分别为铜体和保温材料( λ = 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()
    
    展开全文
  • 关于假人皮肤外侧热传导问题的差分法求解

    千次阅读 多人点赞 2018-10-29 00:17:23
    本人还只是名编程小白,在本问题中也只是用了自己现有的这点水平,此外,这还是我第次在此发表我的文章,对于文章中出现的一些错误希望读者能恰当的指出,欢迎在评论区交流。 问题的提出 今年的数学建模国赛A...

    说明

    本人还只是一名编程小白,在本问题中也只是用了自己现有的这点水平,此外,这还是我第一次在此发表我的文章,对于文章中出现的一些错误希望读者能恰当的指出,欢迎在评论区交流。

    问题的提出

    今年的数学建模国赛A题就是我问题的来源,不过限于我的数学能力,对于问题的求解我只会给出第一问的求解方案。不过从最终的实际效果来看,我在第一问的处理上还是比较合适的,毕竟好歹得了个省奖。
    在高温环境下工作时,人们需要穿着专用服装以避免灼伤。专用服装通常由三层织物材料构成,记为I、II、III层,其中I层与外界环境接触,III层与皮肤之间还存在空隙,将此空隙记为IV层。
    为设计专用服装,将体内温度控制在37ºC的假人放置在实验室的高温环境中,测量假人皮肤外侧的温度。为了降低研发成本、缩短研发周期,请你们利用数学模型来确定假人皮肤外侧的温度变化情况,并解决以下问题:
    (1) 专用服装材料的某些参数值由附件1给出,对环境温度为75ºC、II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟的情形开展实验,测量得到假人皮肤外侧的温度。建立数学模型,计算温度分布。

    问题的分析

    这是个比较常见的热力学问题,只不过这个问题中,将介质分成了四层而已,如果你之前有学习过简单的热传导方程的相关知识,不难得出的是:在问题(Ⅰ)中,你需要把已知的皮肤外界温度分布用上,加入已知条件中,作为一个边界条件,这样,对这个问题,我们就有了两个边界条件:

    1. 第一层与外界环境接触边界温度始终假设为75℃不变;
    2. 假人皮肤外侧的温度分布。

    以及一个初始条件,在这里我是这样处理的:假定假人放入75℃环境中的时刻,四层介质温度都为37℃,这样就有了一个初始条件。
    为了得到温度的分布,我会使用有限差分的方法来解决这个问题,毕竟对于连续的热传导偏微分方程想要用连续的思想来求解还是挺有难度的。

    问题的求解

    热传导偏微分方程的建立

    对于这个偏微分方程的建立,会需要用到非稳态传热方程:在这里插入图片描述
    由于隔层介质是没有热源的,所以qv=0,在这里插入图片描述

    以及我们刚开始提到的边界条件和初始条件,另外的连续性条件则说明了热传导的过程中,热量的传导是连续的,(3)中的qi(x,t)代表热流密度。
    在这里插入图片描述

    有限差分方法

    好了,通过上面的内容我们已经得到了非稳态热传导的偏微分方程以及必要的几个条件,现在我们需要做的是如何把这样一个连续的方程变成一个离散的方程,所谓的离散,就是让时间和空间位置都离散为一个个的点,我们只需要得出在这个点上的温度。例如用一阶向前差分代替一阶偏导:在这里插入图片描述
    按照这样的思路,我们可以得到(2)式的差分形式:在这里插入图片描述
    除此之外,需要得到边界条和初始条件的差分方程,这并不难,需要处理的是要将假人皮肤外表面温度分布拟合出来,通过使用matlab自带的NNtool工具箱,选择一种拟合函数,我在这里用的是Gaussian拟合,最终得出的近似边界方程为:在这里插入图片描述
    综合前面的所有差分方程,我们得到一个方程组:在这里插入图片描述
    有了这个方程组之后,剩下的工作就是编程序了,不过对于这种差分方法,有一个弊端就是:差分会带来一定得误差,并且会传递,导致误差的一个积累,为了使得结果在一定的误差允许范围内,差分时需要满足一个条件:网格傅里叶准则,在本问题当中即进行以下处理:在这里插入图片描述
    我这里之所以将空间上的划分为两个部分,目的在于减小运算次数,减少划分出的网格空间。
    最后献上我的matlab代码:

    %%非稳态热传导差分
    clc;clear;
    p(1:4)=[ 300	862	74.2	1.18];       %介质密度向量
    c(1:4)=[ 1377	2100	1726	1005];         %介质比热容向量
    s(1:4)=[ 0.082	0.37	0.045	0.028];   %介质热传导率向量
    d(1:4)=[ 0.6	6	3.6	5];     %介质厚度向量
    X=15.2*0.001;                             %介质总长度M(1:4)=[7,67,103,108];                  %差分后各个界面的位置
    N=540000;                                  %对时间差分的总间隔
    dt=0.01;                                   %时间间隔为0.01s
    M(1:4)=[7,67,103,108];   
    dx(1:4)=[0.1*0.001 0.1*0.001 0.1*0.001 1*0.001];%不同介质对x的差分间隔向量
    a(1:4)=[1.98499E-07 ,2.04397E-07 ,3.51373E-07 ,2.36108E-05];  %常系数a=s/(c*p)
    alpha(1:4)=a*dt/dx.^2;
    for i=1:M(4)
        U(i,1)=37;                                              %初始条件
    end
    for j=2:N+1
        U(1,j)=75;                                              %边界条件
    end
    for i=1:N+1                                                 %Ⅳ层右侧温度分布
           a1 =        47.44 ;
           b1 =        5212 ;
           c1 =        4664  ;
           a2 =        5.791  ;
           b2 =        250.6  ;
           c2 =        428.2  ;
           a3 =        9.624  ;
           b3 =        617.6  ;
           c3 =        899.2 ;
           a4 =        19.28 ;
           b4 =        1377 ;
           c4 =        2036  ;
           x=i*dt;
    U(108,i)= a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2) + a4*exp(-((x-b4)/c4)^2)+0.9;             %拟合函数
    end 
     for j=1:N
        for i=2:M(4)-1
            if(i<=M(1))
                alpha_use=alpha(1);
            else if(i>M(1) && i<=M(2))
                  alpha_use=alpha(2);
                else if(i>M(2) && i<=M(3))
                        alpha_use=alpha(3);
                    else if(i>M(3) && i<=M(4))
                            alpha_use=alpha(4);               %判断不同介质以区分
                          end 
                      end 
                    end 
             end              
        U(i,j+1)=U(i,j)+alpha_use*(U(i-1,j)+U(i+1,j)-2*U(i,j));  %差分迭代公式
        end
    end
    h(:,:)=U(1:108,1:100:N+1);                                 %存储温度分布数据
    i=1:5401;
    %做出介质分界面处温度分布
    figure(1);
    plot(i,h(7,:),'-r',i,h(67,:),'--b',i,h(103,:));
    legend('Ⅰ、Ⅱ层分界面温度分布℃','Ⅱ、Ⅲ层分界面温度分布℃','Ⅲ、Ⅳ层分界面温度分布℃');
    %做出三维温度分布图
    figure(2);
    i=1:108;
    j=1:5401;
    surf(i*0.1,j,h(i,j)');
    colorbar;
    shading interp;         %取消图像网格线
    xlabel('位置(mm)');
    ylabel('时间(s)');
    zlabel('温度℃');
    %做出不同时刻温度随位置变化曲线图
    figure(3);
    plot(i*0.1,h(:,10),'-r',i*0.1,h(:,50),'-b',i*0.1,h(:,100),'-g');
    legend(' 10s时刻不同位置温度分布℃','50s时刻不同位置温度分布℃','100s时刻不同位置温度分布℃');
    
    
    

    运行结果图为:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    展开全文
  • 基于热传导方程的高温作业专用服装设计(

    千次阅读 多人点赞 2020-08-03 01:42:38
    基于热传导方程的高温作业专用服装设计(一) 此文为不才为备赛准备而已,若有先生雅正,不胜荣幸. 感谢,感谢…… 摘要 本文以傅里叶定律、牛顿冷却定律为理论依据,主要研究高温作业专用服装的...求解一维非稳态热

    基于热传导方程的高温作业专用服装设计(一)

    此文为不才为备赛准备而已,若有先生雅正,不胜荣幸.
    感谢,感谢……

    摘要

    本文以傅里叶定律、牛顿冷却定律为理论依据,主要研究高温作业专用服装的设计,建立了热传导模型,使用追赶法求解结果。
    在问题一中,因为在传导过程中可以不计空气对流,将模型简化为一维热传导模型,并借助拉依达准则剔除异常数据,确保数据得精确性,并观察到题目分为稳态与非稳态两种情况。
    对于一维稳态热传导模型,由热流密度函数求解出热流密度具体值,由温度分布函数求解出各介质的温度如下:
    在这里插入图片描述
    求解一维非稳态热传导问题模型,借助牛顿冷却定律、枚举法,确定空气与皮肤表面的转化系数,并将初值温度37度、左边界Dirichlet边值条件,右边界Robin边值条件,及基于临界面热流量密度和温度相等作为边界值条件使用,建立一维非稳态热传导问题模型。将连续定解区域作网格剖分,用隐式向后差分格式对原微分方程组离散化,得到三对角线性方程组,借助追赶法求解,得到时间与空间维度下的温度分布,见problem1.xlsx。通过查阅文献对差分方法进行对比分析,确定所选模型的精确性,对模型进行误差分析,定义偏差指数 f并求得其值为0.4593,最大误差1.99。
    在问题二中,确定II介质的最优厚度,建立关于厚度的优化模型。在实际生产生活中,降低产品成本是首要考虑因素,制定II介质的最优厚度,实质上就是确定符合生产生活条件的最小厚度。使用问题一所建立的模型,根据问题二的要求对目标进行优化,利用循环遍历的枚举法,得出II介质的最小厚度为19.3mm,即II介质的最优厚度。
    在问题三中,求Ⅱ、IV介质厚度的问题中,建立多目标优化模型。首先,从成本与穿着舒适度两方面,制定“最优”准则,并确定两个不同的优化目标;确定不同的权重之后,将双目标问题转化为单目标问题,可用问题二的循环遍历法继续求解,搜索出II介质的最优厚度为21.6mm,IV层介质的最优厚度为6.4mm。

    关键字:傅里叶定律 隐式差分 追赶法 枚举法 多目标优化。

    一、问题重述

    在高温环境工作时,人需要穿着专用服装以避免灼伤。专用服装通常由三层织物材料构成,记为I、II、III层,其中I层与外界环境接触,III层与皮肤之间还存有空隙,将此空隙记为IV层。
    为设计专用服装,将体内温度控制在37ºC的假人放置在实验室的高温环境中,测量假人皮肤外侧的温度。为了降低研发成本、缩短研发周期,请你们利用数学模型来确定假人皮肤外侧的温度变化情况,并解决以下问题:
    (1)专用服装材料的某些参数值由附件1给出,对环境温度为75ºC、II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟的情形开展实验,测量得到假人皮肤外侧的温度(见附件2)。建立数学模型,计算温度分布,并生成温度分布的Excel文件(文件名为problem1.xlsx)。
    (2) 当环境温度为65ºC、IV层的厚度为5.5 mm时,确定II层的最优厚度,确保工作60分钟时,假人皮肤外侧温度不超过47ºC,且超过44ºC的时间不超过5分钟。
    (3) 当环境温度为80 时,确定II层和IV层的最优厚度,确保工作30分钟时,假人皮肤外侧温度不超过47ºC,且超过44ºC的时间不超过5分钟。

    二、问题分析

    2.1问题一分析
    问题一,根据假设建立一维热传导方程,建立基于热传导方程的温度分布模型,确定在一维空间中介质在不 同时刻,不同厚度下的温度。在模型建立时本文首先借助导热基本定律——傅里叶定律和能量守恒定律推导热传导方程。
    数据处理之后,将题目分解为稳态与非稳态两种不同的形式。然后,求解一维稳态热传导问题模型,之后借助牛顿冷却定律,借助枚举法,确定空气与皮肤表面的转化系数,并将初值温度、左边界Dirichlet边值条件,右边界Robin边值条件,临界面热流量密度和温度相等作为边界值条件,建立维非稳态热传导问题模型,采用有限差分法中隐式差分方法进行模型求解,最后,对模型进行误差分析,定义偏差指数 f并求得其值为0.4593,最大误差1.99。
    2.2问题二分析
    问题二,确定II介质最优厚度是一个最优化问题,首先从服装成本与穿着舒适度两个方面讨论“最优"标准的制定,确定优化问题的目标为II层介质厚度最小。
    其次,考虑问题二关于假人皮肤外侧温度的两个要求,同时结合问题–建立的基于热传导方程的温度分布模型,确定最优化问题的约束条件,从而建立II层最优厚度的单目标优化模型。
    问题二模型的求解利用循环遍历的变步长枚举法,对II层介质的所有可能厚度进行遍历,求出满足约束条件的最小厚度。
    2.3问题三分析
    问题三,有两个变量,分别为Ⅱ和Ⅳ的厚度,相比于问题二增加了一个变量,所以求解Ⅱ,IV两层的最优厚度是一个多目标的优化问题。
    通过查阅相关文献[1]得知,人体外表皮在温度大于44℃时开始发生热损伤,但是在此题中给出30分钟内不超过47℃,并且由于外表温度是单调非递减,所以必定在25分钟之后升至44℃。这可以作为问题三中的约束条件。
    从服装成本与穿着舒适度两个方面考虑,分别制定出不同方面下的“最优”厚度标准,确定多目标优化问题的两个不同的目标。确定不同的权重之后,将双目标问题转化为单目标问题,可用问题二的循环遍历法继续求解。

    三、模型假设

    1、假设各层的介质都是均匀的,且保持各向同性;
    2、假设热防护服的材质均匀无空隙,防护服与假人的形状抽象为嵌套合柱体;
    3、假设各层介质的热传导率在各个方向一致;
    4、假设外界环境无辐射;
    5、忽略材料内部热源对模型的影响;
    6、热传导沿垂直于皮肤方向进行,故系统可假设为一维模型。

    四、符号说明

    在这里插入图片描述
    在这里插入图片描述

    五、问题一的模型建立与求解

    5.1模型准备

    5.1.1简化模型分析

    高温作业专用服装的设计问题,实质上是将模糊现实问题转化为严谨数学模型,对参数进行数字化处理,以数学模型的优化设计,达到实际问题求解的目的。
    图5.1.1.1高温防护服三维立体图
    如图5.1.1所示,I层为材料层且与外界环境接触层,厚度0.6mm;II层为材料层,厚度6mm;III层为材料层,厚度为3.6mm;IV为与皮肤之间还存在空隙层,厚度为5mm;
    易知,在高温作业服工作过程中,实际为空间温度的扩散过程,根据空间三维热传导方程[2]:

    在这里插入图片描述
    其中系数v > 0为介质的热扩散系数, 为已知热源函数。
    可知,本题扩散过程与厚度、高度(防护服底部为基准面)、宽度有关,但因为在传导过程中可以不计空气对流,且不考虑人本身体积对模型影响,因此高度与宽度两维数将作为忽略量处理。即本题模型可近似看为只与厚度有关的一维平面模型。
    一维热传导方程,即抛物型方程[2]为:

    在这里插入图片描述
    一维简化平面图如下:
    图5.1.1.2一维简化平面图
    根据题意,环境温度为75ºC,体内温度37ºC。

    5.1.2名词解释

    温度分布[3]:
    在给定时间的空间某区域内温度随空间位置的变化,一般地讲,物体的温度分布是坐标与时间的函数,即
    在这里插入图片描述
    在这里插入图片描述
    温度梯度:
    沿等温线、等温面法线方向温度的增量与方向距离比值的极限称之为温 度梯度。
    温度梯度为矢量,方向垂直于等温线,且指向温度增加的方向。
    在这里插入图片描述
    傅里叶定律:
    单位时间通过一定截面的导热量,正比于垂直于截面的温度梯度和截面 面积。
    热流量:
    在这里插入图片描述
    热流密度:在这里插入图片描述
    牛顿冷却定律:
    牛顿冷却定律是温度高于周围环境的物体向周围媒质传递热量逐渐冷却 时所遵循的规律。具体的表述为:当物体表面与周围存在温差时,单位时间 从单位面积散失的热量与温度差成正比,比例系数为热传导率。
    计算公式为:
    在这里插入图片描述

    5.1.3数据处理与分析

    拉依达准则剔除异常数据
    附件二所示为皮肤外侧的测量温度随时间变化关系,为确保数据中无 可疑数据影响模型结果,对其进行拉依达准则剔除异常数据处理。
    拉依达准则是指先假设一组检测数据只含有随机误差,对其进行计算处 理得到标准偏差,按一定概率确定一个区间,认为凡超过这个区间的误差, 就不属于随机误差而是粗大误差,含有该误差的数据应予以剔除。
    在这里插入图片描述
    时间与温度关系图分析
    借助Matlab软件,将剔除异常数据过后的数据导入绘制时间-温度关系 图,如图5.1.3.1所示。
    通过对温度随时间变化图像的分析,我们从中可以看出,温度在 t=0 时刻从接近人体体温的位置开始缓慢上升,直到达到48℃左右开始接近平 缓至稳定不变,那么我们可将温度场在时间维度上的变化分为两个阶段。 分析温度分布时将两种状态分开,即非稳态分析和稳态分析。
    稳态:即是一个传热达到动态平衡,从宏观上讲温度在时间上基本达到 恒定的过程,那么我们可以初步考虑,试图用空间维度的角度,描述温度 场在 x 轴方向上的分布特点,其中可以采用逐层差分的方式构建方程,基 于稳态传热学相关理论进行合理的分析和求解。
    非稳态:是一个变化的过程,因此可以 基于传热学的一些基本定律在时 间和空间两个维度上取微元 构建一维偏微分方程,运用数学物理方法中的 一些手段以及 Matlab相关化简及求数值解的算法,期望得到一个可描述非 稳态阶段的温度场函数T(x,t),并绘制出其空间曲面图。
    图5.1.3.1时间温度关系图

    5.1.4模型流程图

    根据假设建立一维温度分度模型,即建立根据距离x的变化,温度 的方程。
    具体模型建立思路如下图5.2.1所示:
    图5.2.1温度分布模型流程图

    5.1.5导热问题的边界条件说明

    在这里插入图片描述

    5.2一维稳态热传导问题及其分析

    在这里插入图片描述
    图5.2.1.1一维高温服装坐标图
    在这里插入图片描述
    在这里插入图片描述

    5.3一维非稳态热传导问题及其分析

    5.3.1热传导方程的推导

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    5.3.2定解问题条件建立一维非稳态热传导模型

    图5.3.2.1一维详细图解
    首先,建立所需的一维图解如图5.3.2.1,确定不同介质本身的热传导方程;
    其次,根据在同一临界面具有相同的热流量密度与温度作为第二边界条件;
    最后,根据牛顿冷却定律,给出Robin右边界条件,即可进行模型建立求解,如下图5.3.2.2所示。
    图5.3.2.2定界问题分析简图
    Step1:确定介质a,b,c,d的热传导方程
    在这里插入图片描述
    Step2:第一、二边界条件定解条件运用
    在这里插入图片描述
    在这里插入图片描述
    Step3:第三边界条件定解条件运用
    已知外界温度为75°,并将热量传递给介质,又附件可得外界温度始终高于人体温度37°,进而在外界与人体热量交换过程中,符合牛顿冷却定律定义式,即是温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律。
    依据牛顿冷却定律,给出 Robin 边界条件:
    在这里插入图片描述
    其中, 表示为介质d在贴近假人皮肤处的温度值; 表示为假人内部恒定温度,即 37°C,h为介质与皮肤之间的转化系数,即冷却系数。
    在 t=0 时刻,介质 a,b,c,d的温度均与假人的恒定皮肤外侧温度相 同,即均为 37°,由此确定初值条件为:
    在这里插入图片描述
    在x=0时,第I层与外界直接接触,温度为T0=75°
    在这里插入图片描述
    步骤3完整定解条件分析结果为:
    在这里插入图片描述
    Step4:转化系数h的确定
    通过查阅资料可知,空气对流换热系数取值范围为
    其次,通过附件 2 给出的不同时刻下假人皮肤表面的温度值,借助变步 长多次枚举法确定最佳转化系数的值。
    图5.3.2.2在不同转化系数h下的皮肤温度变化图
    在这里插入图片描述
    Step5:最终基于热传导方程的温度分布模型
    在这里插入图片描述
    在这里插入图片描述

    5.3.3基于有限差分法的一维非稳态热传导模型求解[4]

    一维非稳态热传导模型属于抛物型方程,因边值条件条件复杂难以求得 解析解,故本文采用有限差分法。
    差分方法又称为有限差分方法或网格法,是求偏微分方程定解问题的数 值解中应用最广泛的方法之一。基本思想:将连续的定解区域用有限个离 散点构成的网格来代替;把定解区域上的连续变量的函数用网格上定义的 离散变量的函数来近似;把原方程和定解条件中的微商用差商来近似。最 终,把原微分方程和定解条件用代数方程组来代替,即有限差分方程组。 解此方程组就可以得到原问题在离散点上的近似值。
    Step1:隐式差分格式方法的确定
    在显示差分格式中,需满足求时间步长与空间步长的比 r 0.5 ,即时 间步长要比空间步长小得多。若不满足此条件,极易发生解的爆破。在本 文温度分布模型中,时间步长与空间步长达到 50,远远大于0.5,因此不适 合用显示差分格式,故本文选用隐式差分格式。
    Step2:对求解区域进行网格剖分
    在这里插入图片描述
    图5.3.3.1矩形网格图
    Step3:建立隐式差分格式
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    Step4:求隐式差分矩阵
    在这里插入图片描述
    Step5: 利用追赶法解三对角线性方程组[5]
    在这里插入图片描述
    在这里插入图片描述

    5.3.4求解结果分析

    四层介质的四个临界面相同时间间隔 下的温度值。绘制成温度关于时间的二维曲线,如下图5.2.1.1所示.
    图5.3.4.1在不同时刻下临界面处的温度变化情况
    在这里插入图片描述
    图5.3.4.2不同时刻不同厚度的温度变化图
    在这里插入图片描述

    5.4差分格式精度分析[4]

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    将 IV 层介质右侧与皮肤直接接触的临界面的温度与题目中附件2假人 皮肤外侧的测量温度进行对比,可检验求解结果的正确性。
    图5.4.5临界面 IV 处误差分析图
    借用matalb编程得标准偏差f=0.4593,基本符合题意。

    参考文献

    [1] 潘斌.热防护服装热传递数学建模及参数决定反问题[D].浙江理工大学.2016.12.28
    [2]司守奎.数学建模方法与应用[M].北京:国防大学出版社,1994年:360页~370页.
    [3]靳涛.高温作业专用服装的温度分布模型[J].消防界(电子版),2020, 6(10):
    24-28.
    [4]zzlghust.第二章-稳态热传导[Z].@百度文库,2018.
    [5]zhangchao3322218.追赶法解三对角线性方程组(Matlab)[Z]. @CSDN,2011.

    clear;%清除工作区变量
    clc;%清屏
    close all;%关闭所有图形窗口
    
    z=[];
    h=8.6125; %空气交换系数
    %% 材料参数输入
    m1=6;m2=60;m3=36;m4=50;% 分别对四种介质分割
    m=m1+m2+m3+m4;% 介质分割和
    n=5400;% 对时间分割
    t=5400;% 总时长
    l1=0.6/1000;l2=6/1000;l3=3.6/1000;l4=5/1000;% 四种材料厚度
    lam_1=0.082;lam_2=0.37;lam_3=0.045;lam_4=0.028;% 四种材料的热传导率
    de_1=300;de_2=862;de_3=74.2;de_4=1.18;% 四种材料的密度
    c1=1377;c2=2100;c3=1726;c4=1005;% 四种材料的比热容
    
    %% 计算热扩散率
    a1=lam_1/(c1*de_1);% I层材料的热扩散率
    a2=lam_2/(c2*de_2);% II层材料的热扩散率
    a3=lam_3/(c3*de_3);% III层材料的热扩散率
    a4=lam_4/(c4*de_4);% IV层材料的热扩散率
    
    %% 材料长度分割和时间步长分割求解
    derta_x1=l1/m1;% I层材料的分割长度
    derta_x2=l2/m2;% II层材料的分割长度
    derta_x3=l3/m3;% III层材料的分割长度
    derta_x4=l4/m4;% IV层材料的分割长度
    derta_t=t/n;% 时间步长分割
    
    %% 计算各层介质剖分的步长比
    r1=derta_t/derta_x1^2*a1;% 第I层介质剖分的步长比
    r2=derta_t/derta_x2^2*a2;% 第II层介质剖分的步长比
    r3=derta_t/derta_x3^2*a3;% 第III层介质剖分的步长比
    r4=derta_t/derta_x4^2*a4;% 第IV层介质剖分的步长比
    
    u=zeros(m+1,n+1);% 定义四层耦合介质温度分布矩阵
    %% 初始条件和边界条件
    u(:,1)=37;%初始条件
    u(1,:)=75;%边界条件
    
    %% 差分格式的系数矩阵的构造
    A=zeros(m,m);
    for i=1:m1-1
     A(i,i)=1+2*r1;
     A(i,i+1)=-r1;
    if i>=2
    A(i,i-1)=-r1;
    end
    end
    A(m1,m1)=(lam_1/derta_x1+lam_2/derta_x2);
    A(m1,m1-1)=-lam_1/derta_x1;
    A(m1,m1+1)=-lam_2/derta_x2;
    
    for i=m1+1:m1+m2-1
    A(i,i)=1+2*r2;
     A(i,i+1)=-r2;    
    A(i,i-1)=-r2;
    end
    A(m1+m2,m1+m2)=(lam_2/derta_x2+lam_3/derta_x3);
    A(m1+m2,m1+m2-1)=-lam_2/derta_x2;
    A(m1+m2,m1+m2+1)=-lam_3/derta_x3;
    
    for i=m1+m2+1:m1+m2+m3-1
    A(i,i)=1+2*r3;
     A(i,i+1)=-r3;    
    A(i,i-1)=-r3;
    end
    A(m1+m2+m3,m1+m2+m3)=(lam_3/derta_x3+lam_4/derta_x4);
    A(m1+m2+m3,m1+m2+m3-1)=-lam_3/derta_x3;
    A(m1+m2+m3,m1+m2+m3+1)=-lam_4/derta_x4;
    
    for i=m1+m2+m3+1:m1+m2+m3+m4-1
    A(i,i)=1+2*r4;
    A(i,i-1)=-r4;
     A(i,i+1)=-r4;    
    end
    A(m,m)=h+lam_4/derta_x4;
    A(m,m-1)=-lam_4/derta_x4;
    
    %% 构造右端项
    for k=2:n+1 
    b=zeros(m,1);
    for i=2:m-1
    b(i,1)=u(i+1,k-1);
    end
    b(1,1)=u(2,k-1)+r1*u(1,k);
    b(m1,1)=0;
    b(m1+m2,1)=0;
    b(m1+m2+m3,1)=0;
    b(m,1)=37*h;
    
    %% 追赶法求解
    bb=diag(A)';
    aa=[0,diag(A,-1)'];
    c=diag(A,1)';
    N=length(bb);
    L=zeros(N);
    uu0=0;y0=0;aa(1)=0;
    L(1)=bb(1)-aa(1)*uu0;
    y(1)=(b(1)-y0*aa(1))/L(1);
    uu(1)=c(1)/L(1);
    for i=2:(N-1)
    L(i)=bb(i)-aa(i)*uu(i-1);
    y(i)=(b(i)-y(i-1)*aa(i))/L(i);
    uu(i)=c(i)/L(i);
    end
    L(N)=bb(N)-aa(N)*uu(N-1);
    y(N)=(b(N)-y(N-1)*aa(N))/L(N);
    x(N)=y(N);
    for i=(N-1):-1:1
    x(i)=y(i)-uu(i)*x(i+1);
    end
    u(2:m+1,k)=x';
    end
    
    %% 绘制不同时刻不同厚度温度分布图
    x=1:1:m+1;
    t=1:1:t+1;
    surf(t,x,u)
    shading interp
    xlabel('t')
    ylabel('x')
    zlabel('u')
    
    u=round(u,2);
    xlswrite('不同时间不同厚度下的温度分布.xlsx',u)% 生成温度分布的excle文件
    
    %% 四个临界面下的部分温度分布表
    U=zeros(5401,4);
    U(:,1)=u(m1+1,:)';
    U(:,2)=u(m1+m2+1,:)';
    U(:,3)=u(m1+m2+m3+1,:)';
    U(:,4)=u(m1+m2+m3+m4+1,:)';
    
    xlswrite('problem1.xlsx',U)% 存储生成四个临界面下的部分温度分布表于problem1
    
    figure
    subplot(2,2,1)
    plot(U(:,1),'r')
    xlabel('t');ylabel('u');
    title('临界面I')
    axis([0 5400 30 80])
    subplot(2,2,2)
    plot(U(:,2),'r')
    xlabel('t');ylabel('u');
    title('临界面II')
    axis([0 5400 30 80])
    subplot(2,2,3)
    plot(U(:,3),'r')
    title('临界面III')
    axis([0 5400 30 80])
    xlabel('t');ylabel('u');
    subplot(2,2,4)
    plot(U(:,4),'r')
    title('临界面IV')
    axis([0 5400 30 80])
    xlabel('t');ylabel('u');
    

    求解空气系数范围

    clear;%清除工作区变量
    clc;%清屏
    close all;%关闭所有图形窗口
    
    z=[];
    for h=8:0.01:9 %确定空气交换系数范围
    %% 材料参数输入
    m1=6;m2=60;m3=36;m4=50;% 分别对四种介质分割
    m=m1+m2+m3+m4;% 介质分割和
    n=5400;% 对时间分割
    t=5400;% 总时长
    l1=0.6/1000;l2=6/1000;l3=3.6/1000;l4=5/1000;% 四种材料厚度
    lam_1=0.082;lam_2=0.37;lam_3=0.045;lam_4=0.028;% 四种材料的热传导率
    de_1=300;de_2=862;de_3=74.2;de_4=1.18;% 四种材料的密度
    c1=1377;c2=2100;c3=1726;c4=1005;% 四种材料的比热容
    
    %% 计算热扩散率
    a1=lam_1/(c1*de_1);% I层材料的热扩散率
    a2=lam_2/(c2*de_2);% II层材料的热扩散率
    a3=lam_3/(c3*de_3);% III层材料的热扩散率
    a4=lam_4/(c4*de_4);% IV层材料的热扩散率
    
    %% 材料长度分割和时间步长分割求解
    derta_x1=l1/m1;% I层材料的分割长度
    derta_x2=l2/m2;% II层材料的分割长度
    derta_x3=l3/m3;% III层材料的分割长度
    derta_x4=l4/m4;% IV层材料的分割长度
    derta_t=t/n;% 时间步长分割
    
    %% 计算各层介质剖分的步长比
    r1=derta_t/derta_x1^2*a1;% 第I层介质剖分的步长比
    r2=derta_t/derta_x2^2*a2;% 第II层介质剖分的步长比
    r3=derta_t/derta_x3^2*a3;% 第III层介质剖分的步长比
    r4=derta_t/derta_x4^2*a4;% 第IV层介质剖分的步长比
    
    u=zeros(m+1,n+1);% 定义四层耦合介质温度分布矩阵
    %% 初始条件和边界条件
    u(:,1)=37;%初始条件
    u(1,:)=75;%边界条件
    
    %% 差分格式的系数矩阵的构造
    A=zeros(m,m);
    for i=1:m1-1
     A(i,i)=1+2*r1;
     A(i,i+1)=-r1;
    if i>=2
    A(i,i-1)=-r1;
    end
    end
    A(m1,m1)=(lam_1/derta_x1+lam_2/derta_x2);
    A(m1,m1-1)=-lam_1/derta_x1;
    A(m1,m1+1)=-lam_2/derta_x2;
    
    for i=m1+1:m1+m2-1
    A(i,i)=1+2*r2;
     A(i,i+1)=-r2;    
    A(i,i-1)=-r2;
    end
    A(m1+m2,m1+m2)=(lam_2/derta_x2+lam_3/derta_x3);
    A(m1+m2,m1+m2-1)=-lam_2/derta_x2;
    A(m1+m2,m1+m2+1)=-lam_3/derta_x3;
    
    for i=m1+m2+1:m1+m2+m3-1
    A(i,i)=1+2*r3;
     A(i,i+1)=-r3;    
    A(i,i-1)=-r3;
    end
    A(m1+m2+m3,m1+m2+m3)=(lam_3/derta_x3+lam_4/derta_x4);
    A(m1+m2+m3,m1+m2+m3-1)=-lam_3/derta_x3;
    A(m1+m2+m3,m1+m2+m3+1)=-lam_4/derta_x4;
    
    for i=m1+m2+m3+1:m1+m2+m3+m4-1
    A(i,i)=1+2*r4;
    A(i,i-1)=-r4;
     A(i,i+1)=-r4;    
    end
    A(m,m)=h+lam_4/derta_x4;
    A(m,m-1)=-lam_4/derta_x4;
    
    %% 构造右端项
    for k=2:n+1 
    b=zeros(m,1);
    for i=2:m-1
    b(i,1)=u(i+1,k-1);
    end
    b(1,1)=u(2,k-1)+r1*u(1,k);
    b(m1,1)=0;
    b(m1+m2,1)=0;
    b(m1+m2+m3,1)=0;
    b(m,1)=37*h;
    
    %% 追赶法求解
    bb=diag(A)';
    aa=[0,diag(A,-1)'];
    c=diag(A,1)';
    N=length(bb);
    L=zeros(N);
    uu0=0;y0=0;aa(1)=0;
    L(1)=bb(1)-aa(1)*uu0;
    y(1)=(b(1)-y0*aa(1))/L(1);
    uu(1)=c(1)/L(1);
    for i=2:(N-1)
    L(i)=bb(i)-aa(i)*uu(i-1);
    y(i)=(b(i)-y(i-1)*aa(i))/L(i);
    uu(i)=c(i)/L(i);
    end
    L(N)=bb(N)-aa(N)*uu(N-1);
    y(N)=(b(N)-y(N-1)*aa(N))/L(N);
    x(N)=y(N);
    for i=(N-1):-1:1
    x(i)=y(i)-uu(i)*x(i+1);
    end
    u(2:m+1,k)=x';
    end
    
    %% 绘制h=8和h=9的皮肤温度变换图
    if h==8
    B1=u(m+1,:);
    plot(B1,'color','b');
    hold on  
    end
    if h==9
    C1=u(m+1,:);
    plot(C1,'color','y') ;
    end
    q=u(m+1,t+1)-48.08;
    z=[z q];
    [d p]=min(abs(z));
    end
    
    fprintf('空气系数范围:\n')
    fprintf('%.2f---%.2f\n',8+(p-1)*0.01,8+(p+1)*0.01)
    
    A1=xlsread('CUMCM-2018-Problem-A-Chinese-Appendix.xlsx',2,'B3:B5403');
    plot(A1,'color','r')
    
    title('不同转化系数h下的皮肤温度变化图')
    xlabel('时间')
    ylabel('温度') 
    
    legend('h=9.0000','h=8.0000','实际值')
    x=0:1:6000;y1=48.08*ones(6001,1);
    plot(x,y1,'--','color','c')
    text(150.6685,48.6752,'48.08')
    
    

    求解空气系数

    clear;%清除工作区变量
    clc;%清屏
    close all;%关闭所有图形窗口
    
    z=[];
    for h=8.61:0.0001:8.63 %确定空气交换系数
    %% 材料参数输入
    m1=6;m2=60;m3=36;m4=50;% 分别对四种介质分割
    m=m1+m2+m3+m4;% 介质分割和
    n=5400;% 对时间分割
    t=5400;% 总时长
    l1=0.6/1000;l2=6/1000;l3=3.6/1000;l4=5/1000;% 四种材料厚度
    lam_1=0.082;lam_2=0.37;lam_3=0.045;lam_4=0.028;% 四种材料的热传导率
    de_1=300;de_2=862;de_3=74.2;de_4=1.18;% 四种材料的密度
    c1=1377;c2=2100;c3=1726;c4=1005;% 四种材料的比热容
    
    %% 计算热扩散率
    a1=lam_1/(c1*de_1);% I层材料的热扩散率
    a2=lam_2/(c2*de_2);% II层材料的热扩散率
    a3=lam_3/(c3*de_3);% III层材料的热扩散率
    a4=lam_4/(c4*de_4);% IV层材料的热扩散率
    
    %% 材料长度分割和时间步长分割求解
    derta_x1=l1/m1;% I层材料的分割长度
    derta_x2=l2/m2;% II层材料的分割长度
    derta_x3=l3/m3;% III层材料的分割长度
    derta_x4=l4/m4;% IV层材料的分割长度
    derta_t=t/n;% 时间步长分割
    
    %% 计算各层介质剖分的步长比
    r1=derta_t/derta_x1^2*a1;% 第I层介质剖分的步长比
    r2=derta_t/derta_x2^2*a2;% 第II层介质剖分的步长比
    r3=derta_t/derta_x3^2*a3;% 第III层介质剖分的步长比
    r4=derta_t/derta_x4^2*a4;% 第IV层介质剖分的步长比
    
    u=zeros(m+1,n+1);% 定义四层耦合介质温度分布矩阵
    %% 初始条件和边界条件
    u(:,1)=37;%初始条件
    u(1,:)=75;%边界条件
    
    %% 差分格式的系数矩阵的构造
    A=zeros(m,m);
    for i=1:m1-1
     A(i,i)=1+2*r1;
     A(i,i+1)=-r1;
    if i>=2
    A(i,i-1)=-r1;
    end
    end
    A(m1,m1)=(lam_1/derta_x1+lam_2/derta_x2);
    A(m1,m1-1)=-lam_1/derta_x1;
    A(m1,m1+1)=-lam_2/derta_x2;
    
    for i=m1+1:m1+m2-1
    A(i,i)=1+2*r2;
     A(i,i+1)=-r2;    
    A(i,i-1)=-r2;
    end
    A(m1+m2,m1+m2)=(lam_2/derta_x2+lam_3/derta_x3);
    A(m1+m2,m1+m2-1)=-lam_2/derta_x2;
    A(m1+m2,m1+m2+1)=-lam_3/derta_x3;
    
    for i=m1+m2+1:m1+m2+m3-1
    A(i,i)=1+2*r3;
     A(i,i+1)=-r3;    
    A(i,i-1)=-r3;
    end
    A(m1+m2+m3,m1+m2+m3)=(lam_3/derta_x3+lam_4/derta_x4);
    A(m1+m2+m3,m1+m2+m3-1)=-lam_3/derta_x3;
    A(m1+m2+m3,m1+m2+m3+1)=-lam_4/derta_x4;
    
    for i=m1+m2+m3+1:m1+m2+m3+m4-1
    A(i,i)=1+2*r4;
    A(i,i-1)=-r4;
    A(i,i+1)=-r4;    
    end
    A(m,m)=h+lam_4/derta_x4;
    A(m,m-1)=-lam_4/derta_x4;
    
    %% 构造右端项
    for k=2:n+1 
    b=zeros(m,1);
    for i=2:m-1
    b(i,1)=u(i+1,k-1);
    end
    b(1,1)=u(2,k-1)+r1*u(1,k);
    b(m1,1)=0;
    b(m1+m2,1)=0;
    b(m1+m2+m3,1)=0;
    b(m,1)=37*h;
    
    %% 追赶法求解
    bb=diag(A)';
    aa=[0,diag(A,-1)'];
    c=diag(A,1)';
    N=length(bb);
    L=zeros(N);
    uu0=0;y0=0;aa(1)=0;
    L(1)=bb(1)-aa(1)*uu0;
    y(1)=(b(1)-y0*aa(1))/L(1);
    uu(1)=c(1)/L(1);
    for i=2:(N-1)
    L(i)=bb(i)-aa(i)*uu(i-1);
    y(i)=(b(i)-y(i-1)*aa(i))/L(i);
    uu(i)=c(i)/L(i);
    end
    L(N)=bb(N)-aa(N)*uu(N-1);
    y(N)=(b(N)-y(N-1)*aa(N))/L(N);
    x(N)=y(N);
    for i=(N-1):-1:1
    x(i)=y(i)-uu(i)*x(i+1);
    end
    u(2:m+1,k)=x';
    end
    q=u(m+1,t+1)-48.08;
    z=[z q];
    [d p]=min(abs(z));
    end
    
    fprintf('空气交换系数:\n')
    fprintf('    %.4f\n',8.61+0.0001*p)
    
    

    链接:https://pan.baidu.com/s/1DhKwv8OECyCAWAbbubkGYQ
    提取码:qelh

    展开全文
  • 一维非稳态导热求解,绝热和热流边界条件,附Matlab程序,共同学习提高function oned% 根据网上一维算例改了下边界条件,边界条件为绝热和恒定热流 % by hxgclear all;clc;%%%%%%%%%%%%%%需要输入的物性参数Lambda=10;%...
  • 多维稳态含内热源导热问题的数值分析,张敏,John C. Chai,用分离变量法和变量变换原理,求解内热源稳态热传导问题。通过数值计算和与精确解相比较,展示和证明这种方法的实用性。
  • Workbench中提供了两种热传导分析模块:Steady-State Thermal ...两种热传导计算模块本文主要介绍固体的热传导的分析方法(包括稳态热传导和瞬态热传导),边界条件的解释(包括热载荷,接触传热的参数解释),...
  • 为了研究高温服装的导热性,开展假人温度分布实验,将问题转化为一维复合介质传热问题,建立数学模型,模拟假人皮肤外侧温度变化情况 给定I层和IV层厚度及各层参数,求解问题 问题描述 在高温环境下工作时,人们...
  • 连铸结晶器内钢液的传热可视为个稳态过程,可用依赖于拉坯速度的三维稳态热传导方程描述。本文针对该稳态模型,采用Galerkin加权余量法推导考虑第类边界条件、第二类边界条件和第三类边界条件的有限元方程,得到...
  • 导热问题的ADI-TDMA算法

    千次阅读 2020-04-30 17:06:03
    本文在上篇文章中关于二维非稳态导热问题的交替方向隐式方法的基础上,设计了个三维算例,并通过Fortran程序,演示适用于三维模型的Brian ADI格式的交替方向隐式方法的求解过程。
  • 通过非稳态热流法和PCK实验装置测定秸秆纤维复合墙体的基本热物理量,研究了温度对秸秆纤维复合墙体的热传导和湿迁移特性的影响,并得到了导热系数、湿迁移系数与温度的简化关系式.考虑秸秆纤维增强水泥非均匀的材料...
  • 为探求该混合物蒸 发过程中丙酮浓度变化与温度变化的关系,根据一维非稳态热传导模型设计出可以实时测量温度数据的试 验装置。运用 Levenberg- Marquardt(L- M)方法以及气液相平衡原理,结合试验结果,求得丙酮-水溶液...
  • 前言:随着电子技术的不断发展,高能量密度,小型化,快速迭代已经成为电子产品的趋势,电子产品中的管理也日益成为个挑战。ICEPAK作为个主流的仿真软件,日益受到设计者的关注。在面向电子产品的仿真,...
  • 导热问题的ADI-TDMA算法

    千次阅读 2020-04-29 12:50:15
    交替方向隐式方法是CFD中用于求解对流-扩散问题的一种非常高效的求解方法,它将二维/三维隐式求解问题转化为各方向的一维隐式求解问题处理,并采用TDMA求解,比直接求解多维问题的线性方程组更高效,且时间步长选取...
  • 摘要:本文主要针对各向同性、有限尺寸、高导热材料样品因闪光加热所引起的非一维传热过程,建立与实际测试更接近的传热模型,采用数值计算方法分析闪光法测试中的背温曲线测量误差,由此明确闪光光斑尺寸、样品截...
  • 假人处于恒高温环境中,不考虑防护服织物的边缘热量损失,且人体和防护服的空气间隔很小,忽略空气的自然对流,只考虑热传导;故可以把织物视为导热多层平面,且属于非稳态导热过程。建立“高温环境-防护服-假人体表...
  • 基于FBG传感模型,利用线性薛定谔方程和一维稳态热传导方程分析了线性效应和光热效应各自对于光谱变化的贡献比例。通过分析波长漂移量与相位改变的关系,建立了基于FBG的线性折射率系数计算模型;通过理论仿真...
  • T_rect_rod_isotup 基于通过变量分离方法获得的二维热传导方程的解析解。 T [K] - 多层矩形棒横截面的相对温度(超过环境温度的温度)。 边界条件假设: - 棒底部的等温条件; - 隔热侧壁; - 上表面的等温条件...
  • 今天我们就来看篇跟热成像相关的文章,这篇文章是由海康、华师大、上交大以及教育部人工智能重点实验室联合发布。 背景 当前科学技术的进步在解决疑难案件中发挥着越来越重要的作用。像仪可以捕获肇事者在...
  • 一维情况下热传导的傅立叶定律有: q = -k∂u/∂x;(q为热流强度) 在x方向上(微元法): dt时间流入左表面的热量为:q|x Adt; dt时间流出右表面的热量为:q|x+dx Adt; 所以净流入为:q|x Adt - q|x+dx Adt = -∂q...
  • 作动器,或者叫螺线管,有时候也称为制动器、电磁阀、电磁铁等,是种通过通电产生磁场来控制衔铁实现理想力矩和位移的设备。衔铁为铁磁物质,受到磁场作用后,产生吸力并把电能转化成机械能,用于对负载的速度、...
  • 目前已经建立起各种各样的理论模型(如热传导模型、流体动力学模型,分子动力学模型、双温模型等) 大多数只是在特定的范围内有用。 经典的热传导模型主要用在长脉冲烧蚀的领域( 如纳秒级脉冲激光) ,当激光脉冲达到...
  • % GeometricThermalICs - 区域或区域边界上的初始温度 % NodalThermalICs - 在网格节点指定的初始温度 % ThermalResults - 热解及其派生量 % SteadyStateThermalResults - 稳态热模型解及其派生量 % ...
  • 如何使用 COMSOL 进行电热分析?

    千次阅读 2020-06-28 16:17:09
    、电磁损耗的热源有多种类型。我们可以使用 COMSOL软件的内置功能计算所有的电磁热源(准静态或高频状态);软件中预定义的接口包括焦耳,感应加热,微波加热 和激光加热。 1、焦耳多物理接口耦合了固体传热 ...
  •   、线性代数基本方程组   基本方程组:   矩阵表示: ...要点:矩阵的每行代表个方程,m行代表m个线性联立方程。 n列代表n个变量。如果m是独立方程数,根据m&lt;n、m=n、...
  • 状态转移概率:在事件的发展变化过程中,从有种状态出发,下时刻转移到其他状态的可能性,记为 P { E i ∣ E j } = ⋅ p i j \mathbb{P}\{E_i|E_j\}\overset{\cdot}= p_{ij} P{Ei​∣Ej​}=⋅pij​ 。...

空空如也

空空如也

1 2 3 4 5 ... 8
收藏数 157
精华内容 62
关键字:

一维非稳态热传导模型