柱体内温度分布图 MATLAB

对于下底面和侧面绝热，上底面温度与半径平方成正比的柱体，绘制柱体内温度分布图。

这里给出两种尝试：1、散点图；2、切片云图

1. 散点图仿真

首先使用解析算法求的场解值的解析表达，其次求解Bessel函数零点，带入场解表达。对于空间点阵划分，采用每一定半径划分圆周上等数量点，在总方向复制累计。

下面给出散点图的仿真代码，Nt为划分精度，Nt越大空间划分点越多。

% function y = cal117(R,h,Nt,N)
close all; clear all
R = 5;
h = 5;
Nt = 10;

%% 求bessel函数零点
pi0 = besal_pi0(0,Nt);


%% 计算圆柱体分割点坐标，作为u的前三列坐标
m = 100;
cot = 0;
for i = 1 : 20
    [x_,y_,z_] = cylinder((R*i/20),m-1);
    for k = 1 : m
        for j = 0 : 19
            cot = cot + 1;
            u(cot,1) = x_(1,k);
            u(cot,2) = y_(1,k);
            u(cot,3) = h*j/(19);
            u(cot,4) = 0;       %test
        end
    end
end


%% 计算温度值 u第四列为温度值
for i = 1 : cot
    for j = 1 : Nt
        ct(j) = 2*(R^2)*(1-4/(pi0(j)^2))/(pi0(j)*besselj(1,pi0(j))*sinh(pi0(j)*h/R));
        u(i,4) = u(i,4)+ ct(j)*sinh(pi0(j)*u(i,3)/R)*besselj(0,pi0(j)*sqrt(u(i,1)^2+u(i,2)^2)/R);
        if(u(i,1)>10)
            u(i,4)=NaN;
        end
    end
end




%% 图像绘制
[x0,y0,z0] = cylinder(R,100); %画圆柱体轮廓
z0 = h * z0;
plot3(x0(1,:),y0(1,:),z0(1,:),'k');hold on
plot3(x0(1,:),y0(1,:),z0(2,:),'k');hold on
for i = 1 : 10 :100
    plot3([x0(1,i),x0(1,i)],[y0(1,i),y0(1,i)],[z0(1,i),z0(2,i)],'k');hold on
end

x = u(:,1);
y = u(:,2);
z = u(:,3);
c = u(:,4);

% [X,Y] = meshgrid(linspace(min(x),max(x),20)',linspace(min(y),max(y),30)');
% Z = griddata(x,y,z,X,Y,'v4');
% C =griddata(x,y,c,X,Y,'v4');

scatter3(x,y,z,24,'filled','k');
axis([-(R+2) (R+2) -(R+2) (R+2) 0 (h+2)]);

% scatter3(x,y,z,24,c,'filled');
% axis([-(R+2) (R+2) -(R+2) (R+2) 0 (h+2)]);
% colorbar;

2. 切片云图绘制

空间点集划分同上，其中将四位数组转化为三维数组，坐标三维即坐标点位置，值即温度值。

可以绘制出纵向切片图、整体剖面图和等温面图三张图：

clear; clc;
pi0 = besal_pi0(0,100);
R = 5;
h = 5;
xa = 10;
xb = 10;
xc = 5;
v = zeros(xa,xb,xc);

for i = 1 : xa
    for j = 1 : xb
        for k = 1 : xc
            for t = 1 : 20
                ct(t) = 2*(R^2)*(1-4/(pi0(t)^2))/(pi0(t)*besselj(1,pi0(t))*sinh(pi0(t)*h/R));
                v(i,j,k)=v(i,j,k)+ ct(t)*sinh(pi0(t)*k/R)*besselj(0,pi0(t)*sqrt((i-5)^2+(j-5)^2)/R);
            end
%             if(k>5)
%                 v(i,j,k) = NaN;
%             end
            if((i-5)^2+(j-5)^2>25)
                v(i,j,k) = NaN;
            end
        end    
    end
end

[x,y,z]=meshgrid(1:xb,1:xa,1:xc);% 坐标值的范围
h=figure(1);% 第一幅图
% 各大切片
set(h,'name','取单切片')
subplot(221)% 切片1
slice(x,y,z,v,[],[1],[]);
shading interp 
% set(gca,'zdir','reverse');
axis equal
grid on
colorbar

subplot(222)% 切片2
slice(x,y,z,v,[],[2],[]);
shading interp 
colormap('jet')
% set(gca,'zdir','reverse');
axis equal
grid on
colorbar

subplot(223)% 切片3
slice(x,y,z,v,[],[3],[]);
shading interp 
% set(gca,'zdir','reverse');
axis equal
grid on
colorbar

subplot(224)% 切片4
slice(x,y,z,v,[],[4],[]);
shading interp 
% set(gca,'zdir','reverse');
axis equal
grid on
colorbar
%%
h2=figure(2);
figure('menubar','figure');
% rotate 3D
set(h2,'name','全空间切片','MenuBar','none','ToolBar','none')
slice(x,y,z,v,[1:2:xb],[2 3 4],[2 3 4 5])
shading interp 
colorbar 
colormap('jet')
% set(gca,'zdir','reverse');
axis equal
grid on
% box on

%%
h3=figure(5);
figure('menubar','figure');
set(h3,'name','定值包络立体图','MenuBar','none','ToolBar','none')
set(gcf,'InvertHardcopy','off')
fw=0;   %%此值为最外层包络面取值
fv=isosurface(x,y,z,v,fw);% 等值面
p=patch(fv);

set(p,'facecolor','b','edgecolor','none');
patch(isocaps(x,y,z,v, fw), 'FaceColor', 'interp', 'EdgeColor', 'none');
colorbar

colormap('jet')
box on
daspect([1,1,1])
view(3)

camlight
camproj perspective
lighting phong % 光照处理
axis equal
grid on
title(['最外层表面的值为： ' , num2str(fw),'^{o}C']);

最后贴上一些效果图哈：