对于下底面和侧面绝热，上底面温度与半径平方成正比的柱体，绘制柱体内温度分布图。
这里给出两种尝试：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']);
最后贴上一些效果图哈：
