如何在MATLAB 3D铁相图中实现曲面交线处颜色变化

需求 #

现有MATLAB脚本可生成铁相的3D可视化图，包含3条工况曲线与铁相曲面。需修改脚本，使曲面在与3D曲线的交点区域改变颜色，直观区分铁、wustite、magnetite三种铁相。

效果对比 #

当前效果图 #

目标效果图（手绘示意） #

修改方案 #

核心思路 #

通过给曲面的颜色矩阵C赋值，判断曲面上每个点与三条工况曲线的距离：当点距离某条曲线足够近时，赋予该曲线对应铁相的颜色值；其余点保持原曲面颜色。最终通过surf函数的颜色映射实现分区变色。

步骤1：避免变量覆盖 #

原代码中后续曲面计算会覆盖三条曲线的T、ov_x_co、ov_x_h2变量，需先将三条曲线的数据重命名保存：

% 重命名三条曲线数据，避免被后续曲面计算覆盖
T_curve1 = T1; x_co_curve1 = ov_x_co1; x_h2_curve1 = ov_x_h21;
T_curve2 = T2; x_co_curve2 = ov_x_co2; x_h2_curve2 = ov_x_h22;
T_curve3 = T3; x_co_curve3 = ov_x_co3; x_h2_curve3 = ov_x_h23;

步骤2：构建曲面颜色矩阵 #

在曲面数据计算完成后，创建颜色矩阵C，遍历每个曲面点，计算其与三条曲线的距离，匹配对应相的颜色：

% 初始化颜色矩阵，默认颜色值0.2
C = ones(size(ov_x_co_surf)) * 0.2;

% 设置距离阈值（可根据实际效果调整）
distance_threshold = 0.02;

% 遍历曲面每个点
for i = 1:size(ov_T_surf, 1)
    for j = 1:size(ov_T_surf, 2)
        % 跳过无效值
        if isnan(ov_T_surf(i,j)) || isnan(ov_x_co_surf(i,j)) || isnan(ov_x_h2_surf(i,j))
            continue;
        end
        
        % 计算当前点到三条曲线的最小距离（过滤无效数据）
        valid_idx1 = ~isnan(x_h2_curve1) & ~isnan(x_co_curve1);
        dist1 = any(valid_idx1) ? min(sqrt((T_curve1(valid_idx1)-ov_T_surf(i,j)).^2 + (x_co_curve1(valid_idx1)-ov_x_co_surf(i,j)).^2 + (x_h2_curve1(valid_idx1)-ov_x_h2_surf(i,j)).^2)) : inf;
        
        valid_idx2 = ~isnan(x_h2_curve2) & ~isnan(x_co_curve2);
        dist2 = any(valid_idx2) ? min(sqrt((T_curve2(valid_idx2)-ov_T_surf(i,j)).^2 + (x_co_curve2(valid_idx2)-ov_x_co_surf(i,j)).^2 + (x_h2_curve2(valid_idx2)-ov_x_h2_surf(i,j)).^2)) : inf;
        
        valid_idx3 = ~isnan(x_h2_curve3) & ~isnan(x_co_curve3);
        dist3 = any(valid_idx3) ? min(sqrt((T_curve3(valid_idx3)-ov_T_surf(i,j)).^2 + (x_co_curve3(valid_idx3)-ov_x_co_surf(i,j)).^2 + (x_h2_curve3(valid_idx3)-ov_x_h2_surf(i,j)).^2)) : inf;
        
        % 根据距离设置颜色值（对应parula色阶，可自行调整）
        if dist1 < distance_threshold
            C(i,j) = 0.8; % 铁相：亮色调
        elseif dist2 < distance_threshold
            C(i,j) = 0.5; % wustite相：中间色调
        elseif dist3 < distance_threshold
            C(i,j) = 0.1; % magnetite相：暗色调
        end
    end
end

步骤3：修改曲面绘制代码 #

使用带颜色矩阵的surf函数，关闭边缘线让曲面更平滑：

surf(ov_T_surf, ov_x_co_surf, ov_x_h2_surf, C, 'EdgeColor', 'none');
colormap('parula');

步骤4：修正图例 #

原图例标注不准确，更新为对应铁相：

legend('铁相曲线','Wustite相曲线','Magnetite相曲线','铁相曲面');

完整修改后代码 #

R = 8.314; % J/mol/K
p = 1; % atm

% Script 1（铁相曲线）
T1 = 700:10:1200;
ov_x_co1 = [];
ov_x_h21 = [];

for i = 1:numel(T1)
    % Calculation for Script 1
    k6 = exp((-10034 -38.635*T1(i)*log(T1(i))+271.78*T1(i))./(-R.*T1(i)));
    k7 = exp((26546-38.635*T1(i)*log(T1(i))+238.315*T1(i))./(-R.*T1(i)));
    k2 = exp((172140 - 177.7.*T1(i))./(-R.*T1(i)));

    x_co = k6 * k2 / p;
    x_h2 = (p-k6*k2*(k6+1))/(p*(k7+1));

    ov_x_co1 = [ov_x_co1; x_co];
    ov_x_h21 = [ov_x_h21; x_h2];
end

% Script 2（Wustite相曲线）
T2 = 700:10:1200;
ov_x_co2 = [];
ov_x_h22 = [];

for i = 1:numel(T2)
    % Calculation for Script 2
    k6 = exp((-21785+25*T2(i))./(-R.*T2(i)));
    k7 = exp((14799-8.465*T2(i))./(-R.*T2(i)));
    k2 = exp((172140 - 177.7.*T2(i))./(-R.*T2(i)));

    x_co = k6 * k2 / p;
    x_h2 = (p-k6*k2*(k6+1))/(p*(k7+1));

    ov_x_co2 = [ov_x_co2; x_co];
    ov_x_h22 = [ov_x_h22; x_h2];
end

% Script 3（Magnetite相曲线）
T3 = 700:10:1200;
ov_x_co3 = [];
ov_x_h23 = [];

for i = 1:numel(T3)
    % Calculation for Script 3
    k6 = exp((-18844 -9.66*T3(i)*log(T3(i))+86.695*T3(i))./(-R.*T3(i)));
    k7 = exp((17736-9.66*T3(i)*log(T3(i))+52.23*T3(i))./(-R.*T3(i)));
    k2 = exp((172140 - 177.7.*T3(i))./(-R.*T3(i)));

    x_co = k6 * k2 / p;
    x_h2 = (p-k6*k2*(k6+1))/(p*(k7+1));

    ov_x_co3 = [ov_x_co3; x_co];
    ov_x_h23 = [ov_x_h23; x_h2];
end

% 重命名三条曲线数据，避免被后续曲面计算覆盖
T_curve1 = T1; x_co_curve1 = ov_x_co1; x_h2_curve1 = ov_x_h21;
T_curve2 = T2; x_co_curve2 = ov_x_co2; x_h2_curve2 = ov_x_h22;
T_curve3 = T3; x_co_curve3 = ov_x_co3; x_h2_curve3 = ov_x_h23;

% 处理曲线数据的无效值
for i = 1:size(x_h2_curve1, 1)
    if x_h2_curve1(i) < 0 || x_h2_curve1(i) > 1
        x_h2_curve1(i) = NaN;
    end
end
for i = 1:size(x_co_curve1, 1)
    if x_co_curve1(i) < 0 || x_co_curve1(i) > 1
        x_co_curve1(i) = NaN;
    end
end
% 重复处理曲线2、3的无效值
for i = 1:size(x_h2_curve2, 1)
    if x_h2_curve2(i) < 0 || x_h2_curve2(i) > 1
        x_h2_curve2(i) = NaN;
    end
end
for i = 1:size(x_co_curve2, 1)
    if x_co_curve2(i) < 0 || x_co_curve2(i) > 1
        x_co_curve2(i) = NaN;
    end
end
for i = 1:size(x_h2_curve3, 1)
    if x_h2_curve3(i) < 0 || x_h2_curve3(i) > 1
        x_h2_curve3(i) = NaN;
    end
end
for i = 1:size(x_co_curve3, 1)
    if x_co_curve3(i) < 0 || x_co_curve3(i) > 1
        x_co_curve3(i) = NaN;
    end
end

% 绘制三条工况曲线
figure;
plot3(T_curve1, x_co_curve1, x_h2_curve1, '-', 'LineWidth', 2.5, 'Color', 'red');
hold on;
plot3(T_curve2, x_co_curve2, x_h2_curve2, '-', 'LineWidth', 2.5, 'Color', 'green');
plot3(T_curve3, x_co_curve3, x_h2_curve3, '-', 'LineWidth', 2.5, 'Color', 'blue');

% 计算铁相曲面数据
T_surf = 700:10:1200;
ov_x_co_surf = [];
ov_x_h2_surf = [];
ov_T_surf = [];

for i = 1:numel(T_surf)
    k1 = exp((36580 - 33.465.*T_surf(i))./(-R.*T_surf(i)));
    k2 = exp((172140 - 177.7.*T_surf(i))./(-R.*T_surf(i)));

    x_co = 0.0:0.05:1.0;
    x_h2 = zeros(1, numel(x_co));
    var_temp = T_surf(i)*ones(1, numel(x_co));
    ov_T_surf = [ov_T_surf; var_temp];

    for j = 1:numel(x_co)    
        x_co2 = x_co(j).^2 ./ k2;
        x_h2(j) = ((x_co(j)./(x_co2 .* k1)).*(1-x_co(j)-x_co2))./(1+(x_co(j)./(x_co2.*k1)));
    end

    ov_x_co_surf = [ov_x_co_surf;x_co];
    ov_x_h2_surf = [ov_x_h2_surf; x_h2];
end

% 处理曲面数据的无效值
for i = 1:size(ov_x_h2_surf, 1)
    for j = 1:size(ov_x_h2_surf, 2)
        if ov_x_h2_surf(i, j) < 0 || ov_x_h2_surf(i, j) > 1
            ov_x_h2_surf(i, j) = NaN;
        end
    end
end

% 构建曲面颜色矩阵
C = ones(size(ov_x_co_surf)) * 0.2; % 默认颜色
distance_threshold = 0.02; % 距离阈值，可调整

for i = 1:size(ov_T_surf, 1)
    for j = 1:size(ov_T_surf, 2)
        if isnan(ov_T_surf(i,j)) || isnan(ov_x_co_surf(i,j)) || isnan(ov_x_h2_surf(i,j))
            continue;
        end
        
        % 计算到三条曲线的最小距离
        valid_idx1 = ~isnan(x_h2_curve1) & ~isnan(x_co_curve1);
        if any(valid_idx1)
            dist1 = min(sqrt((T_curve1(valid_idx1) - ov_T_surf(i,j)).^2 + ...
                (x_co_curve1(valid_idx1) - ov_x_co_surf(i,j)).^2 + ...
                (x_h2_curve1(valid_idx1) - ov_x_h2_surf(i,j)).^2));
        else
            dist1 = inf;
        end
        
        valid_idx2 = ~isnan(x_h2_curve2) & ~isnan(x_co_curve2);
        if any(valid_idx2)
            dist2 = min(sqrt((T_curve2(valid_idx2) - ov_T_surf(i,j)).^2 + ...
                (x_co_curve2(valid_idx2) - ov_x_co_surf(i,j)).^2 + ...
                (x_h2_curve2(valid_idx2) - ov_x_h2_surf(i,j)).^2));
        else
            dist2 = inf;
        end
        
        valid_idx3 = ~isnan(x_h2_curve3) & ~isnan(x_co_curve3);
        if any(valid_idx3)
            dist3 = min(sqrt((T_curve3(valid_idx3) - ov_T_surf(i,j)).^2 + ...
                (x_co_curve3(valid_idx3) - ov_x_co_surf(i,j)).^2 + ...
                (x_h2_curve3(valid_idx3) - ov_x_h2_surf(i,j)).^2));
        else
            dist3 = inf;
        end
        
        % 分配颜色
        if dist1 < distance_threshold
            C(i,j) = 0.8; % 铁相区域颜色
        elseif dist2 < distance_threshold
            C(i,j) = 0.5; % Wustite相区域颜色
        elseif dist3 < distance_threshold
            C(i,j) = 0.1; % Magnetite相区域颜色
        end
    end
end

% 绘制带颜色的曲面
surf(ov_T_surf, ov_x_co_surf, ov_x_h2_surf, C, 'EdgeColor', 'none');
colormap('parula');

% 设置坐标轴与图例
ylim([0, 1]);
zlim([0, 1]);
xlabel('T (K)');
ylabel('x_{CO}');
zlabel('x_{H2}');
title('铁相3D分区可视化');
legend('铁相曲线','Wustite相曲线','Magnetite相曲线','铁相曲面');
grid on;
hold off;

注意事项 #

• 距离阈值调整：distance_threshold的值需根据曲线与曲面的贴合度调整，过大会导致颜色区域过宽，过小则可能无法识别。
• 颜色值调整：颜色值对应parula色阶的0-1范围，可根据需求更换colormap（如jet、cool）或直接使用RGB颜色（需修改颜色矩阵为三维数组）。
• 无效值处理：代码中保留了NaN处理逻辑，避免无效数据干扰可视化效果。

内容的提问来源于stack exchange，提问作者Bhavya Nagda