柱体内温度分布图 MATLAB

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

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

1. 散点图仿真

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

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

% function y = cal117(R,h,Nt,N)

close all; clear all

R = ;

h = ;

Nt = ;

%% 求bessel函数零点

pi0 = besal_pi0(,Nt);

%% 计算圆柱体分割点坐标，作为u的前三列坐标

m = ;

cot = ;

for i =  :

    [x_,y_,z_] = cylinder((R*i/),m-);

    for k =  : m

        for j =  :

            cot = cot + ;

            u(cot,) = x_(,k);

            u(cot,) = y_(,k);

            u(cot,) = h*j/();

            u(cot,) = ;       %test

        end

    end

end

%% 计算温度值 u第四列为温度值

for i =  : cot

    for j =  : Nt

        ct(j) = *(R^)*(-/(pi0(j)^))/(pi0(j)*besselj(,pi0(j))*sinh(pi0(j)*h/R));

        u(i,) = u(i,)+ ct(j)*sinh(pi0(j)*u(i,)/R)*besselj(,pi0(j)*sqrt(u(i,)^+u(i,)^)/R);

        if(u(i,)>)

            u(i,)=NaN;

        end

    end

end

%% 图像绘制

[x0,y0,z0] = cylinder(R,); %画圆柱体轮廓

z0 = h * z0;

plot3(x0(,:),y0(,:),z0(,:),'k');hold on

plot3(x0(,:),y0(,:),z0(,:),'k');hold on

for i =  :  :

    plot3([x0(,i),x0(,i)],[y0(,i),y0(,i)],[z0(,i),z0(,i)],'k');hold on

end

x = u(:,);

y = u(:,);

z = u(:,);

c = u(:,);

% [X,Y] = meshgrid(linspace(min(x),max(x),)',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,,'filled','k');

axis([-(R+) (R+) -(R+) (R+)  (h+)]);

% scatter3(x,y,z,,c,'filled');

% axis([-(R+) (R+) -(R+) (R+)  (h+)]);

% colorbar;

2. 切片云图绘制

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

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

clear; clc;

pi0 = besal_pi0(,);

R = ;

h = ;

xa = ;

xb = ;

xc = ;

v = zeros(xa,xb,xc);

for i =  : xa

    for j =  : xb

        for k =  : xc

            for t =  :

                ct(t) = *(R^)*(-/(pi0(t)^))/(pi0(t)*besselj(,pi0(t))*sinh(pi0(t)*h/R));

                v(i,j,k)=v(i,j,k)+ ct(t)*sinh(pi0(t)*k/R)*besselj(,pi0(t)*sqrt((i-)^+(j-)^)/R);

            end

%             if(k>)

%                 v(i,j,k) = NaN;

%             end

            if((i-)^+(j-)^>)

                v(i,j,k) = NaN;

            end

        end

    end

end

[x,y,z]=meshgrid(:xb,:xa,:xc);% 坐标值的范围

h=figure();% 第一幅图

% 各大切片

set(h,'name','取单切片')

subplot()% 切片1

slice(x,y,z,v,[],[],[]);

shading interp

% set(gca,'zdir','reverse');

axis equal

grid on

colorbar

subplot()% 切片2

slice(x,y,z,v,[],[],[]);

shading interp

colormap('jet')

% set(gca,'zdir','reverse');

axis equal

grid on

colorbar

subplot()% 切片3

slice(x,y,z,v,[],[],[]);

shading interp

% set(gca,'zdir','reverse');

axis equal

grid on

colorbar

subplot()% 切片4

slice(x,y,z,v,[],[],[]);

shading interp

% set(gca,'zdir','reverse');

axis equal

grid on

colorbar

%%

h2=figure();

figure('menubar','figure');

% rotate 3D

set(h2,'name','全空间切片','MenuBar','none','ToolBar','none')

slice(x,y,z,v,[::xb],[  ],[   ])

shading interp

colorbar

colormap('jet')

% set(gca,'zdir','reverse');

axis equal

grid on

% box on

%%

h3=figure();

figure('menubar','figure');

set(h3,'name','定值包络立体图','MenuBar','none','ToolBar','none')

set(gcf,'InvertHardcopy','off')

fw=;   %%此值为最外层包络面取值

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([,,])

view()

camlight

camproj perspective

lighting phong % 光照处理

axis equal

grid on

title(['最外层表面的值为： ' , num2str(fw),'^{o}C']);

最后贴上一些效果图哈：