柱體內溫度分布圖 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']);

 

最后貼上一些效果圖哈：