有限差分法实现三维物体的传热过程【MATLAB仿真】

最新推荐文章于 2024-02-13 22:49:18 发布

攻城狮不暴燥

最新推荐文章于 2024-02-13 22:49:18 发布

阅读量3.1k

点赞数 14

分类专栏：计算电磁学文章标签： matlab

本文链接：https://blog.csdn.net/weixin_46096864/article/details/130117403

版权

计算电磁学专栏收录该内容

1 篇文章 2 订阅

订阅专栏

有限差分法实现三维物体的传热过程

理论推导
代码实现
结果展示

理论推导

从热传导方程出发，有：
$\rho C \frac{\partial T}{\partial t} - \nabla \cdot (\kappa \nabla T) = Q$
其中， $\rho$ 为材料的密度； $C$ 为材料的比热容； $\kappa$ 为材料的热导率(或导热系数)； $Q$ 为焦耳热源。
若不考虑热源项 $Q$ ，则上述的传热方程可写为：
$\begin{cases}\frac{\partial T}{\partial t} = \alpha(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2}) \\\\ \alpha = \frac{\kappa}{\rho C} \end{cases}$
时间变量采用前向差分，空间变量先采用前向差分再采用后向差分，那么有：
$\frac{T(x,y,z,t+\Delta t)-T(x,y,z,t)}{\Delta t} = \alpha[\frac{T(x+\Delta x,y,z,t)+T(x-\Delta x,y,z,t)-2T(x,y,z,t)}{(\Delta x)^2} + \frac{T(x,y+\Delta y,z,t)+T(x,y-\Delta y,z,t)-2T(x,y,z,t)}{(\Delta y)^2} + \frac{T(x,y,z+\Delta z,t)+T(x,y,z-\Delta z,t)-2T(x,y,z,t)}{(\Delta z)^2}]$
假定各个方向上的空间步长均相等： $\Delta x = \Delta y = \Delta z = \Delta$ ，那么材料在 $t+\Delta t$ 时刻的温度分布为：
$T(x,y,z,t+\Delta t) = T(x,y,z,t) + \frac{\Delta t \cdot \alpha}{\Delta^2}[T(x+\Delta x,y,z,t)+T(x-\Delta x,y,z,t)-2T(x,y,z,t) + T(x,y+\Delta y,z,t)+T(x,y-\Delta y,z,t)-2T(x,y,z,t) + T(x,y,z+\Delta z,t)+T(x,y,z-\Delta z,t)-2T(x,y,z,t)]$

代码实现

主函数包括【几何建模】、【参数设置】、【热计算】三大部分。

close all;clear;
%% 几何部分：建模与网格剖分
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%建模部分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = (1:0.25:18)';
y = (1:0.25:8)';
z = (1:0.25:2)';
%%%%%%%%%%%%%%%%%%%%%%%%%%构造网格坐标点信息%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l = length(x)*length(y)*length(z);
OutPut = zeros(l,3); % 重要参数一：用于存储网格坐标信息
% 网格x轴坐标点信息
countx = 1;
for i = 1:length(x)
    for j=1:length(y)*length(z)
        OutPut(j+length(y)*length(z)*(i-1),1) = x(countx);
    end
    countx = countx+1;
end
% 网格y轴坐标点信息
for i=1:length(x)
    county = 1;
    for j=1:length(y)
        for k=1:length(z)
            OutPut(k+length(z)*(j-1)+length(y)*length(z)*(i-1),2) = y(county);
        end
        county = county+1;
    end
end
% 网格z轴坐标点信息
OutPut(:,3) = repmat(z,length(x)*length(y),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 参数部分：为模型赋热传导所需数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%输入参数数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输入材料的热传导参数，Param=[介质的热导率 介质的密度 介质的比热容 材料的热导率 材料的密度 材料的比热容]
Param = [130 2330 700 398 8930 386];%介质是硅，材料是铜
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%处理网格信息%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MeshTria = delaunayTriangulation(OutPut);
MeshEdges = meshLenth(MeshTria);
MeshEdgesSelected = MeshEdges(MeshEdges(:,4)==1,:);
% MeshAttach找到每个点的所有连接点
MeshEdges2 = [MeshEdgesSelected(:,1:3);MeshEdgesSelected(:,2),MeshEdgesSelected(:,1),MeshEdgesSelected(:,3)];
[~,MeshEdgesLabel] = sort(MeshEdges2(:,1));
MeshEdgesAll = MeshEdges2(MeshEdgesLabel,:);
% EdgesCenter所有网格边的中心点
EdgesCenter = (OutPut(MeshEdgesAll(:,1),:) + OutPut(MeshEdgesAll(:,2),:))./2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%显示网格剖分的结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MeshEdgesSelectedSide = MeshEdgesSelected;
% figure;
% T1=datetime;
% for i = 1:size(MeshEdgesSelectedSide,1)
%     Temp = [OutPut(MeshEdgesSelectedSide(i,1),:);OutPut(MeshEdgesSelectedSide(i,2),:)]';
%     plot3(Temp(1,:),Temp(2,:),Temp(3,:),'k.-')
%     hold on
% end
% T2=datetime;
% disp(['显示用时:' num2str(seconds(T2-T1)) 's' ])
% axis equal
% title({[ '点数=' num2str(size(OutPut,1)) '  边数=' num2str(size(MeshEdges,1)) ]} )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算alpha%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 设置填充介质的位置
AirSurf11 = logical((EdgesCenter(:,1) == 1)+(EdgesCenter(:,2) == 4.5)+(EdgesCenter(:,3) == 1));
AirSurf12 = logical((EdgesCenter(:,1) == 4)+(EdgesCenter(:,2) == 8)+(EdgesCenter(:,3) == 2));
AirSurf21 = logical((EdgesCenter(:,1) == 4)+(EdgesCenter(:,2) == 5.5)+(EdgesCenter(:,3) == 1));
AirSurf22 = logical((EdgesCenter(:,1) == 14)+(EdgesCenter(:,2) == 8)+(EdgesCenter(:,3) == 2));
AirSurf31 = logical((EdgesCenter(:,1) == 14)+(EdgesCenter(:,2) == 4.5)+(EdgesCenter(:,3) == 1));
AirSurf32 = logical((EdgesCenter(:,1) == 18)+(EdgesCenter(:,2) == 8)+(EdgesCenter(:,3) == 2));
AirSurf41 = logical((EdgesCenter(:,1) == 1)+(EdgesCenter(:,2) == 1)+(EdgesCenter(:,3) == 1));
AirSurf42 = logical((EdgesCenter(:,1) == 4)+(EdgesCenter(:,2) == 3.5)+(EdgesCenter(:,3) == 2));
AirSurf51 = logical((EdgesCenter(:,1) == 4)+(EdgesCenter(:,2) == 1)+(EdgesCenter(:,3) == 1));
AirSurf52 = logical((EdgesCenter(:,1) == 14)+(EdgesCenter(:,2) == 2.5)+(EdgesCenter(:,3) == 2));
AirSurf61 = logical((EdgesCenter(:,1) == 14)+(EdgesCenter(:,2) == 1)+(EdgesCenter(:,3) == 1));
AirSurf62 = logical((EdgesCenter(:,1) == 18)+(EdgesCenter(:,2) == 3.5)+(EdgesCenter(:,3) == 2));
AirSurf = logical((AirSurf11.*AirSurf12) + (AirSurf21.*AirSurf22) + (AirSurf31.*AirSurf32) + (AirSurf41.*AirSurf42) + (AirSurf51.*AirSurf52) + (AirSurf61.*AirSurf62));
% 设置介质与材料连接边的位置
Surf11 = logical((EdgesCenter(:,1) == 1.25)+(EdgesCenter(:,2) == 4.75)+(EdgesCenter(:,3) == 1.25));
Surf12 = logical((EdgesCenter(:,1) == 3.75)+(EdgesCenter(:,2) == 7.75)+(EdgesCenter(:,3) == 1.75));
Surf21 = logical((EdgesCenter(:,1) == 4.25)+(EdgesCenter(:,2) == 5.75)+(EdgesCenter(:,3) == 1.25));
Surf22 = logical((EdgesCenter(:,1) == 13.75)+(EdgesCenter(:,2) == 7.75)+(EdgesCenter(:,3) == 1.75));
Surf31 = logical((EdgesCenter(:,1) == 14.25)+(EdgesCenter(:,2) == 4.75)+(EdgesCenter(:,3) == 1.25));
Surf32 = logical((EdgesCenter(:,1) == 17.75)+(EdgesCenter(:,2) == 7.75)+(EdgesCenter(:,3) == 1.75));
Surf41 = logical((EdgesCenter(:,1) == 1.25)+(EdgesCenter(:,2) == 1.25)+(EdgesCenter(:,3) == 1.25));
Surf42 = logical((EdgesCenter(:,1) == 3.75)+(EdgesCenter(:,2) == 3.25)+(EdgesCenter(:,3) == 1.75));
Surf51 = logical((EdgesCenter(:,1) == 4.25)+(EdgesCenter(:,2) == 1.25)+(EdgesCenter(:,3) == 1.25));
Surf52 = logical((EdgesCenter(:,1) == 13.75)+(EdgesCenter(:,2) == 2.25)+(EdgesCenter(:,3) == 1.75));
Surf61 = logical((EdgesCenter(:,1) == 14.25)+(EdgesCenter(:,2) == 1.25)+(EdgesCenter(:,3) == 1.25));
Surf62 = logical((EdgesCenter(:,1) == 17.75)+(EdgesCenter(:,2) == 3.25)+(EdgesCenter(:,3) == 1.75));
% 介质的Alpha
MeshEdgesAll(:,4) = Param(1)/(Param(2)*Param(3));
% 材料的Alpha
MeshEdgesAll(AirSurf,4) = Param(4)/(Param(5)*Param(6));
% 边界的Alpha
MeshEdgesAll(Surf11,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(2)*Param(3));
MeshEdgesAll(Surf12,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(5)*Param(6));
MeshEdgesAll(Surf21,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(2)*Param(3));
MeshEdgesAll(Surf22,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(5)*Param(6));
MeshEdgesAll(Surf31,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(2)*Param(3));
MeshEdgesAll(Surf32,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(5)*Param(6));
MeshEdgesAll(Surf41,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(2)*Param(3));
MeshEdgesAll(Surf42,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(5)*Param(6));
MeshEdgesAll(Surf51,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(2)*Param(3));
MeshEdgesAll(Surf52,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(5)*Param(6));
MeshEdgesAll(Surf61,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(2)*Param(3));
MeshEdgesAll(Surf62,4) = Param(1)*Param(4)/(Param(1)+Param(4))/(Param(5)*Param(6));
MeshAttach=MeshEdgesAll;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 计算部分：热传导过程仿真
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%输入仿真时空步数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = 300;% 总仿真时间（单位：秒）
deltat = 0.01;deltax = 0.01;% 时空步长
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设置初始温度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PointNum = size(OutPut,1);
T = 80.*ones(PointNum,f);
for i=1:10005
    if(OutPut(i,1)>=1 && OutPut(i,1)<=4)
        if(OutPut(i,2)>=4.5 && OutPut(i,2)<=8)
            if(OutPut(i,3)>=1 && OutPut(i,3)<=2)
                T(i,1) = 20;
            end
        end
    end
end
for i=1:10005
    if(OutPut(i,1)>=4 && OutPut(i,1)<=14)
        if(OutPut(i,2)>=5.5 && OutPut(i,2)<=8)
            if(OutPut(i,3)>=1 && OutPut(i,3)<=2)
                T(i,1) = 20;
            end
        end
    end
end
for i=1:10005
    if(OutPut(i,1)>=14 && OutPut(i,1)<=18)
        if(OutPut(i,2)>=4.5 && OutPut(i,2)<=8)
            if(OutPut(i,3)>=1 && OutPut(i,3)<=2)
                T(i,1) = 20;
            end
        end
    end
end
for i=1:10005
    if(OutPut(i,1)>=1 && OutPut(i,1)<=4)
        if(OutPut(i,2)>=1 && OutPut(i,2)<=3.5)
            if(OutPut(i,3)>=1 && OutPut(i,3)<=2)
                T(i,1) = 20;
            end
        end
    end
end
for i=1:10005
    if(OutPut(i,1)>=4 && OutPut(i,1)<=14)
        if(OutPut(i,2)>=1 && OutPut(i,2)<=2.5)
            if(OutPut(i,3)>=1 && OutPut(i,3)<=2)
                T(i,1) = 20;
            end
        end
    end
end
for i=1:10005
    if(OutPut(i,1)>=14 && OutPut(i,1)<=18)
        if(OutPut(i,2)>=1 && OutPut(i,2)<=3.5)
            if(OutPut(i,3)>=1 && OutPut(i,3)<=2)
                T(i,1) = 20;
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%热传导更新方程迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MeshAttachTemp = [diff(MeshAttach(:,1));1];
MeshAttachTemp1 = find(MeshAttachTemp ~= 0);
MeshAttachTemp1 = [0;MeshAttachTemp1];
for t = 1:f-1
    for n = 1:size(MeshAttachTemp1,1)-1
        nn = MeshAttach(MeshAttachTemp1(n+1),1);
        Ttemp = 0;
        for m = MeshAttachTemp1(n)+1:MeshAttachTemp1(n+1)
            Ttemp = Ttemp + MeshAttach(m,4)*deltat*((T(MeshAttach(m,2),t) - T(nn,t))/ deltax^2);
        end
        T(nn,t+1) = T(nn,t) + Ttemp;
        if Ttemp > 10
            error('计算结果震荡，请减小导热系数Alpha!   Results is volatile!  Please reduce the parameter Alpha!')
        end
        if isnan(T(n,t+1))
            error('计算结果不存在，请检查网格间距是否为0!   Results is NaN!  Please Check the length of mesh!')
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%显示热传导结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Section = logical((OutPut(:,3) == 1.5));% 显示的xoy截面Section
DisplayTime = 1:60:300;% 显示的时间序列
for i = DisplayTime
    figure('color','w');
    surf(reshape(T(Section,i),length(y),length(x)),'LineStyle','none')
    colormap jet
    title('截面温度')
    colorbar;
%     caxis([20,40]);
    axis equal
    set(0,'defaultfigurecolor','w')
    view(0,90)
    pause(1e-4)
end

子函数为处理模型的剖分部分

function MeshEdges = meshLenth(MeshTria)

% Compute the length of mesh side

MeshPoint = MeshTria.Points;
MeshEdges = edges(MeshTria);

for i = 1:size(MeshEdges,1)
    Temp = MeshPoint(MeshEdges(i,1),:) - MeshPoint(MeshEdges(i,2),:);
    MeshEdges(i,3) = sqrt( sum(MeshPoint(MeshEdges(i,1),:) - MeshPoint(MeshEdges(i,2),:)).^2 );
    MeshEdges(i,4) = sum(abs(Temp)<0.15)==2;
end

end