体切片图（`slice`）
体切片图（slice）
05_场与矢量 5 张 MATLAB R2023b 100 mm 宽 @ 300 dpi
这个方法适合什么
在体数据上切若干正交或斜切面，看内部结构。适合：三维温度场、变压器内部电位、土壤温度分层，不用做等值面即可看内部。
电气场景
1 铁心温度正交切片
2 土壤电位正交切片
3 开关柜温度正交切片
4 电缆沟温度正交切片
5 变电站电场正交切片
代码
脚本 p5_field_vector.m 中对应 slice 的一段，SEL 用于选择只跑某几个方法。
if any(SEL==10)
% ---- 10.1 铁心温度正交切片 ----
try
xg = linspace(0,1.2,31); yg = linspace(0,0.8,25); zg = linspace(0,1.0,25);
[X,Y,Z] = meshgrid(xg,yg,zg);
T = 34 + 26*exp(-(((X-0.60).^2)/0.09 + ((Y-0.40).^2)/0.05 + ((Z-0.50).^2)/0.09)) ...
+ 8*exp(-(((X-0.60).^2)/0.16 + ((Y-0.40).^2)/0.10 + ((Z-0.10).^2)/0.03));
f = ml_fig(PP,100,90);
ax = ax3(f,true);
axes(ax);
hs = slice(X,Y,Z,T,0.60,0.40,0.50);
set(hs,'EdgeColor','none','FaceAlpha',0.95);
shading(ax,'interp');
colormap(ax,hot); clim(ax,[34 62]);
xlabel(ml_tx('铁心长度 <i>x</i> / m'));
ylabel(ml_tx('铁心宽度 <i>y</i> / m'));
zlabel(ml_tx('铁心高度 <i>z</i> / m'));
xlim(ax,[0 1.2]); ylim(ax,[0 0.8]); zlim(ax,[0 1.0]);
cb3(ax,PP,'铁心温度 <i>T</i> / °C');
save_png(f,FIGD,'铁心温度正交切片');
catch ME
fprintf('[FAIL] 10-1 %s\n',ME.message); close all
end
% ---- 10.2 土壤电位正交切片 ----
try
xg = linspace(0,40,33); yg = linspace(0,40,33); zg = linspace(-4,0,13);
[X,Y,Z] = meshgrid(xg,yg,zg);
V = 95./sqrt((X-10).^2+(Y-10).^2+(Z+0.9).^2+0.20) ...
+ 95./sqrt((X-30).^2+(Y-10).^2+(Z+0.9).^2+0.20) ...
+ 95./sqrt((X-10).^2+(Y-30).^2+(Z+0.9).^2+0.20) ...
+ 95./sqrt((X-30).^2+(Y-30).^2+(Z+0.9).^2+0.20);
f = ml_fig(PP,100,90);
ax = ax3(f,true);
axes(ax);
hs = slice(X,Y,Z,V,20,20,-1.0);
set(hs,'EdgeColor','none','FaceAlpha',0.95);
shading(ax,'interp');
colormap(ax,parula); clim(ax,[20 230]);
xlabel(ml_tx('站区 <i>x</i> / m'));
ylabel(ml_tx('站区 <i>y</i> / m'));
zlabel(ml_tx('土壤深度 <i>z</i> / m'));
xlim(ax,[0 40]); ylim(ax,[0 40]); zlim(ax,[-4 0]);
cb3(ax,PP,'土壤电位 <i>U</i> / V');
save_png(f,FIGD,'土壤电位正交切片');
catch ME
fprintf('[FAIL] 10-2 %s\n',ME.message); close all
end
% ---- 10.3 开关柜温度正交切片 ----
try
xg = linspace(0,1.2,31); yg = linspace(0,0.8,25); zg = linspace(0,1.8,31);
[X,Y,Z] = meshgrid(xg,yg,zg);
T = 30 + 34*exp(-(((X-0.38).^2)/0.045 + ((Y-0.36).^2)/0.045 + ((Z-1.30).^2)/0.10)) ...
+ 16*exp(-(((X-0.90).^2)/0.060 + ((Y-0.50).^2)/0.060 + ((Z-0.55).^2)/0.12));
f = ml_fig(PP,100,90);
ax = ax3(f,true);
axes(ax);
hs = slice(X,Y,Z,T,0.38,0.36,1.30);
set(hs,'EdgeColor','none','FaceAlpha',0.95);
shading(ax,'interp');
colormap(ax,turbo); clim(ax,[30 70]);
xlabel(ml_tx('柜宽 <i>x</i> / m'));
ylabel(ml_tx('柜深 <i>y</i> / m'));
zlabel(ml_tx('柜高 <i>z</i> / m'));
xlim(ax,[0 1.2]); ylim(ax,[0 0.8]); zlim(ax,[0 1.8]);
cb3(ax,PP,'柜内温度 <i>T</i> / °C');
save_png(f,FIGD,'开关柜温度正交切片');
catch ME
fprintf('[FAIL] 10-3 %s\n',ME.message); close all
end
% ---- 10.4 电缆沟温度正交切片 ----
try
xg = linspace(0,12,41); yg = linspace(0,1.2,25); zg = linspace(0,1.0,21);
[X,Y,Z] = meshgrid(xg,yg,zg);
T = 27 + 0.9*(1-Z/1.0) ...
+ 26*exp(-(((Y-0.25).^2)/0.012 + ((Z-0.80).^2)/0.010)).*(0.88+0.12*cos(2*pi*X/12)) ...
+ 26*exp(-(((Y-0.60).^2)/0.012 + ((Z-0.80).^2)/0.010)).*(0.88+0.12*cos(2*pi*X/12)) ...
+ 24*exp(-(((Y-0.95).^2)/0.012 + ((Z-0.80).^2)/0.010)).*(0.88+0.12*cos(2*pi*X/12)) ...
+ 6*exp(-((Z-0.05).^2)/0.020);
f = ml_fig(PP,100,90);
ax = ax3(f,true);
axes(ax);
hs = slice(X,Y,Z,T,6.0,0.60,0.80);
set(hs,'EdgeColor','none','FaceAlpha',0.95);
shading(ax,'interp');
colormap(ax,hot); clim(ax,[27 58]);
xlabel(ml_tx('沿沟长度 <i>x</i> / m'));
ylabel(ml_tx('沟宽 <i>y</i> / m'));
zlabel(ml_tx('敷设深度 <i>z</i> / m'));
xlim(ax,[0 12]); ylim(ax,[0 1.2]); zlim(ax,[0 1.0]);
cb3(ax,PP,'沟内温度 <i>T</i> / °C');
save_png(f,FIGD,'电缆沟温度正交切片');
catch ME
fprintf('[FAIL] 10-4 %s\n',ME.message); close all
end
% ---- 10.5 变电站电场正交切片 ----
try
xg = linspace(-20,20,33); yg = linspace(0,14,25); zg = linspace(0,10,21);
[X,Y,Z] = meshgrid(xg,yg,zg);
E = 0.35 + 2.90*exp(-((Z-9.2).^2)/3.5).*(0.45+0.55*exp(-((Y-7).^2)/45)) ...
+ 0.90*exp(-((X+8).^2)/40);
f = ml_fig(PP,100,90);
ax = ax3(f,true);
axes(ax);
hs = slice(X,Y,Z,E,0,7,4);
set(hs,'EdgeColor','none','FaceAlpha',0.95);
shading(ax,'interp');
colormap(ax,sky); clim(ax,[0.35 3.30]);
xlabel(ml_tx('站区 <i>x</i> / m'));
ylabel(ml_tx('站区 <i>y</i> / m'));
zlabel(ml_tx('高度 <i>z</i> / m'));
xlim(ax,[-20 20]); ylim(ax,[0 14]); zlim(ax,[0 10]);
cb3(ax,PP,'电场强度 <i>E</i> / (kV/m)');
save_png(f,FIGD,'变电站电场正交切片');
catch ME
fprintf('[FAIL] 10-5 %s\n',ME.message); close all
end
end
fprintf('\n=========== 05 场与矢量：全部完成（应 50 张）===========\n');
%% ===================== 局部函数 =====================
function ax = ax2(f)
%AX2 二维坐标区：右侧留白色条位置
ax = axes('Parent',f,'Position',[0.175 0.190 0.700 0.755]);
box(ax,'on'); grid(ax,'off');
set(ax,'LineWidth',0.5,'Layer','top','TickDir','in');
end
function ax = ax3(f,withcb)
%AX3 三维坐标区：左侧留 0.150 给 zlabel；带色条时收窄到 0.605
if nargin < 2, withcb = false; end
if withcb
p = [0.150 0.180 0.605 0.760];
else
p = [0.150 0.180 0.720 0.760];
end
ax = axes('Parent',f,'Position',p);
view(ax,[-37.5 30]); box(ax,'on'); grid(ax,'off');
set(ax,'LineWidth',0.5,'Layer','top','TickDir','in');
end
function cb = cb2(ax,PP,lab)
%CB2 二维图色条：自动避让坐标区
cb = colorbar(ax);
cb.FontName = PP.FontEN; cb.FontSize = PP.FS;
cb.Label.String = ml_tx(lab);
cb.Label.FontName = PP.FontEN; cb.Label.FontSize = PP.FS;
cb.TickDirection = 'in'; cb.LineWidth = 0.5;
end
function cb = cb3(ax,PP,lab)
%CB3 三维图色条：手工定位到坐标区右侧，避免右缘刻度被裁
cb = colorbar(ax);
cb.Position = [0.800 0.180 0.016 0.760];
cb.FontName = PP.FontEN; cb.FontSize = PP.FS;
cb.Label.String = ml_tx(lab);
cb.Label.FontName = PP.FontEN; cb.Label.FontSize = PP.FS;
cb.TickDirection = 'in'; cb.LineWidth = 0.5;
end
function cone_plot3(ax,X,Y,Z,U,V,W,sc,col)
%CONE_PLOT3 圆锥矢量：矢量先归一化，锥体长度 ≈ sc × 网格平均间距
% ★ coneplot 的 scale 会按 max|v| 归一后再乘到矢量上，锥长只由 scale×间距决定，
% 与原始量纲无关；故这里统一归一化，只用相对强弱（0.65~1.0）体现大小差异
MG = sqrt(U.^2 + V.^2 + W.^2);
sm = sort(MG(:));
ref = sm(max(1,round(0.80*numel(sm))));
if ref <= 0, ref = 1; end
LF = 0.65 + 0.35*min(MG/ref,1.2)/1.2; % 0.65 ~ 1.00
U2 = U./max(MG,eps).*LF;
V2 = V./max(MG,eps).*LF;
W2 = W./max(MG,eps).*LF;
% 锥体原点向域内收缩 6%，避免锥尖/锥底被坐标区边框切掉
x0 = min(X(:)); x1 = max(X(:));
y0 = min(Y(:)); y1 = max(Y(:));
z0 = min(Z(:)); z1 = max(Z(:));
cx = x0 + 0.06*(x1-x0) + 0.88*(X-x0);
cy = y0 + 0.06*(y1-y0) + 0.88*(Y-y0);
cz = z0 + 0.06*(z1-z0) + 0.88*(Z-z0);
axes(ax);
h = coneplot(X,Y,Z,U2,V2,W2,cx,cy,cz,sc);
set(h,'FaceColor',col,'EdgeColor','none');
camlight('headlight'); lighting(ax,'gouraud');
end
function qv2(ax,X,Y,U,V,alen,col)
%QV2 二维箭头场：按 80% 分位模长归一后限幅，箭头最大长度固定为 1.35*alen，
% 避免 1/r 型奇异场把少数箭头画爆、越出坐标区
MG = hypot(U,V);
sm = sort(MG(:));
ref = sm(max(1,round(0.80*numel(sm))));
if ref <= 0, ref = 1; end
SC = min(sqrt(MG/ref) , 1.35); % 开方压缩动态范围：弱场也看得见，强场不越界
U2 = U./max(MG,eps).*SC*alen;
V2 = V./max(MG,eps).*SC*alen;
quiver(ax,X,Y,U2,V2,0,'Color',col,'LineWidth',0.5,'MaxHeadSize',0.55);
end
function qv3(ax,X,Y,Z,U,V,W,alen,col)
%QV3 三维箭头场：同 QV2，箭头最大长度 1.35*alen
MG = sqrt(U.^2+V.^2+W.^2);
sm = sort(MG(:));
ref = sm(max(1,round(0.80*numel(sm))));
if ref <= 0, ref = 1; end
SC = min(MG/ref , 1.35);
U2 = U./max(MG,eps).*SC*alen;
V2 = V./max(MG,eps).*SC*alen;
W2 = W./max(MG,eps).*SC*alen;
quiver3(ax,X,Y,Z,U2,V2,W2,0,'Color',col,'LineWidth',0.5);
end
function w = tube_widths(X,Y,Z,U,V,W,sx,sy,sz,wmin,wmax)
%TUBE_WIDTHS 流管直径按种子点处流速线性映射到 [wmin,wmax]，实现「粗细表示强度」
SP = sqrt(U.^2 + V.^2 + W.^2);
sz0 = size(sx);
v = reshape(interp3(X,Y,Z,SP,sx(:),sy(:),sz(:)), sz0); % ★ interp3 对矩阵查询会塌成列，须还原
v(isnan(v)) = min(SP(:));
vlo = min(v(:)); vhi = max(v(:)); % ★ 必须取全元素极值，否则 min/max 按列返回行向量、除法变矩阵除
if vhi - vlo < eps
w = 0.5*(wmin+wmax)*ones(size(v));
else
w = wmin + (wmax-wmin)*(v - vlo)/(vhi - vlo);
end
end
function [X,Y,Z,U,V,W] = gen_rod_E(nx,ny,nz)
%GEN_ROD_E 棒板间隙静电场：棒端 (1.5,1.5,2.1) 对 z=0 接地板
% ★ 电场线从棒端发散后终止在板面，全程留在计算域内，不会在盒子边界被切断
xg = linspace(0,3,nx); yg = linspace(0,3,ny); zg = linspace(0,2.4,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
PH = 26./sqrt((X-1.5).^2+(Y-1.5).^2+(Z-2.1).^2+0.02) + 26*Z;
[gx,gy,gz] = gradient(PH,xg,yg,zg);
U = -gx; V = -gy; W = -gz;
end
function [X,Y,U,V] = gen_oil2(nx,ny)
%GEN_OIL2 变压器油道二维流场（流函数法，天然无散）
xg = linspace(0,4,nx); yg = linspace(0,2.4,ny);
[X,Y] = meshgrid(xg,yg);
PSI = 0.55*sin(pi*X/4).*sin(2*pi*Y/2.4) + 0.10*Y;
[dPdx,dPdy] = gradient(PSI,xg,yg);
U = dPdy; V = -dPdx;
end
function [X,Y,U,V] = gen_wind2(nx,ny)
%GEN_WIND2 变电站风速二维流场：均匀来流 + 绕建筑物旋涡
xg = linspace(0,40,nx); yg = linspace(0,30,ny);
[X,Y] = meshgrid(xg,yg);
PSI = 4.2*Y - 55*(Y-15)./((X-22).^2+(Y-15).^2+30) + 0.60*sin(2*pi*X/40);
[dPdx,dPdy] = gradient(PSI,xg,yg);
U = dPdy; V = -dPdx;
end
function [X,Y,U,V] = gen_cab2(nx,ny)
%GEN_CAB2 开关柜气流二维流场（下进上出，含一次回流）
xg = linspace(0,1.2,nx); yg = linspace(0,1.8,ny);
[X,Y] = meshgrid(xg,yg);
PSI = 0.40*sin(pi*X/1.2).*sin(pi*Y/1.8) + 0.22*(1-X/1.2).*sin(pi*Y/1.8);
[dPdx,dPdy] = gradient(PSI,xg,yg);
U = dPdy; V = -dPdx;
end
function [X,Y,U,V] = gen_duct2(nx,ny)
%GEN_DUCT2 电缆沟自然对流二维流场（双环流）
xg = linspace(0,1.2,nx); yg = linspace(0,1.0,ny);
[X,Y] = meshgrid(xg,yg);
PSI = 0.22*sin(pi*X/1.2).*sin(pi*Y/1.0) ...
+ 0.12*sin(2*pi*X/1.2).*sin(2*pi*Y/1.0);
[dPdx,dPdy] = gradient(PSI,xg,yg);
U = dPdy; V = -dPdx;
end
function [X,Y,U,V] = gen_esp2(nx,ny)
%GEN_ESP2 电除尘器内部气流二维流场（主流 + 极板间回流）
xg = linspace(0,6,nx); yg = linspace(0,2.4,ny);
[X,Y] = meshgrid(xg,yg);
PSI = 1.05*sin(pi*Y/2.4) + 0.30*sin(pi*Y/2.4).*cos(2*pi*X/6);
[dPdx,dPdy] = gradient(PSI,xg,yg);
U = dPdy; V = -dPdx;
end
function [X,Y,Z,U,V,W] = gen_gap_E(nx,ny,nz)
%GEN_GAP_E 球隙静电场：两带电球 + 数值负梯度求 E
xg = linspace(0,4,nx); yg = linspace(0,3,ny); zg = linspace(0,3,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
PH = 7.5./sqrt((X-1.4).^2+(Y-1.5).^2+(Z-1.5).^2+0.03) ...
- 7.5./sqrt((X-2.8).^2+(Y-1.5).^2+(Z-1.5).^2+0.03);
[gx,gy,gz] = gradient(PH,xg,yg,zg);
U = -gx; V = -gy; W = -gz;
end
function [X,Y,Z,U,V,W] = gen_oil(nx,ny,nz)
%GEN_OIL 变压器三维油流流场（轴向主流 + 上下油道环流）
xg = linspace(0,2,nx); yg = linspace(0,1.2,ny); zg = linspace(0,0.8,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
U = 0.50 + 0.10*cos(pi*Y/1.2).*cos(pi*Z/0.8);
V = 0.04*sin(pi*X/2);
W = 0.03*sin(pi*X/2).*cos(pi*Y/1.2);
end
function [X,Y,Z,U,V,W] = gen_wind(nx,ny,nz)
%GEN_WIND 变电站三维风速场（含建筑尾流与上升气流）
xg = linspace(0,40,nx); yg = linspace(0,30,ny); zg = linspace(0,18,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
U = 3.4*(1+0.22*exp(-((Y-15).^2)/60)).*(0.55+0.45*(Z/18).^0.25);
V = 0.55*sin(pi*X/40).*cos(pi*Y/30);
W = 0.30*exp(-((X-22).^2+(Y-15).^2)/70);
end
function [X,Y,Z,U,V,W] = gen_cab(nx,ny,nz)
%GEN_CAB 开关柜三维气流场（下进上出）
xg = linspace(0,1.2,nx); yg = linspace(0,0.8,ny); zg = linspace(0,1.8,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
U = 0.12*sin(pi*Z/1.8);
V = 0.10*sin(pi*X/1.2);
W = 0.55 + 0.30*(X/1.2) + 0.15*cos(pi*Y/0.8);
end
function [X,Y,Z,U,V,W] = gen_tun(nx,ny,nz)
%GEN_TUN 电缆隧道三维通风流场（纵向主流 + 热浮升）
xg = linspace(0,12,nx); yg = linspace(0,1.2,ny); zg = linspace(0,1.0,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
U = 1.60 + 0.50*cos(pi*Y/1.2);
V = 0.10*sin(pi*X/12);
W = 0.18*exp(-((X-6).^2)/18).*(Y/1.2);
end
function [X,Y,Z,QX,QY,QZ] = gen_cab_q(nx,ny,nz)
%GEN_CAB_Q 开关柜三维温度场及其热流密度 q = -k grad T
xg = linspace(0,1.2,nx); yg = linspace(0,0.8,ny); zg = linspace(0,1.8,nz);
[X,Y,Z] = meshgrid(xg,yg,zg);
T = 30 + 34*exp(-(((X-0.38).^2)/0.05 + ((Y-0.36).^2)/0.05 + ((Z-1.30).^2)/0.10));
[gx,gy,gz] = gradient(T,xg,yg,zg);
k = 0.028;
QX = -k*gx; QY = -k*gy; QZ = -k*gz;
end
function p = save_png(f,FIGD,nm)
%SAVE_PNG 统一关网格 → 导出 300 dpi PNG → 边距自检 → 关图
ha = findall(f,'Type','axes');
for k = 1:numel(ha)
set(ha(k),'XGrid','off','YGrid','off','ZGrid','off');
end
p = fullfile(FIGD, ['05_场与矢量_' nm '261107.png']);
print(f, p, '-dpng', '-r300');
ml_chkpng(p, 8);
close(f);
end
浙公网安备 33010602011771号