用MATLAB实现三元相图的绘制

7月前浏览10853

本文摘要（由AI生成）：

本文主要介绍了相图的绘制方法，以水为案例，讲解了单一组分的相图绘制，并介绍了二组分相图的绘制方法。三元相图的作用和特点，以及如何绘制三元相图。

化学化工生产中对产品进行分离提纯时离不开蒸馏、结晶、萃取等各种单元操作，而这些单元操作过程中的理论知识就要涉及到相平衡的原理。相平衡的重要表现形式就是用相图来表示相平衡系统的状态是如何随着组成、温度、压力等变量而变化，相图就是下图所示的这样的图像

图片1.png

本期过冷水就和大家一起学习一下相图的绘制方式，首先先讲解单一组分的相图绘制一遍大家对相图有更好的了解，以水为案例。

%单组元相图绘制
clear
data1=[7.32E 003.67E 00
1.05E 014.95E 00
1.39E 016.93E 00
1.75E 019.85E 00
2.07E 011.28E 01
2.36E 011.60E 01
2.57E 012.02E 01
2.83E 012.43E 01
];
data2=[2.83E 012.43E 01
3.25E 012.65E 01
4.06E 013.09E 01
4.53E 013.38E 01
5.13E 013.82E 01
5.66E 014.30E 01
6.66E 015.28E 01
7.10E 015.85E 01
7.58E 016.45E 01
7.71E 016.63E 01
8.16E 017.35E 01
];
data3=[2.83E 012.46E 01
2.84E 012.91E 01
2.82E 013.77E 01
2.77E 014.46E 01
2.75E 015.08E 01
2.68E 015.92E 01
2.66E 016.54E 01
2.68E 016.85E 01
2.66E 017.63E 01
];
data4=[8.46E 001.60E 01
1.26E 011.77E 01
1.73E 011.99E 01
2.31E 012.20E 01
2.83E 012.41E 01
];
data5=[8.61E 006.91E 01
2.68E 016.85E 01
5.00E 016.71E 01
7.73E 016.65E 01
8.67E 016.58E 01
];
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(data1(:,1),data1(:,2),'LineWidth',3);
plot(data2(:,1),data2(:,2),'LineWidth',3);
plot(data3(:,1),data3(:,2),'LineWidth',3);
plot(data4(:,1),data4(:,2),'LineWidth',3,'LineStyle','--');
plot(data5(:,1),data5(:,2),'LineWidth',3,'LineStyle','--');
M=[46.23;21.053;41.81;26.7533;7.2732;81.7911; 28.2845;8.4641; 8.6060;27.0622; 51.9470;77.0392;86.3710];
N=[16.55;45.14;74.48;75.8412;3.8563;73.8752;24.4234;15.9546;69.3878;68.4548;67.2886;66.3557;65.6560];
str=['汽';'冰';'水';'A';'B';'C';'D';'E';'e';'d';'c';'b';'a'];
text(axes1,M,N,str,'FontSize',16)
box(axes1,'on');
hold(axes1,'off');
set(axes1,'FontSize',16,'LineWidth',3,'XTick',zeros(1,0),'YTick',zeros(1,0),...
    'YTickLabel','');

相图的特点和绘制方式.图中DA、DB、DC分别表示的是固液平衡线、固气平衡线、气液平衡线；在该线段上的压强和温度条件下水的两相可以共存，不满足温度和压强则不能共存。这样的系统只有一个能独立改变的量。列如水和水蒸气两相平衡就可以用DC线表示，指定了两相平衡的温度，则两相的压力也就确定了。若降低系统的温度，并使其任然保持两相平衡，则水蒸气的压力会沿着DC下降。温度降至0℃系统到达三相点，理论上应该会有冰出现。但对于纯净水来说不会出现冰，这就是水的过冷现象。这种状态下的水就称为过冷水。这就是昵称的由来。相图的应用可以说明系统在外界改变时发生相变化的情况。通过单组分相图的给大家做一点相图作用的简单作用介绍。

二组分相图介绍

%%二组分相图
clear
data1=[0.00E 000.00E 00
9.79E-015.59E 01];
data2=[1.66E-032.08E 01
9.79E-015.58E 01];
data3=[1.66E-032.08E 01
1.00E 000.00E 00
];
data4=[0.00E 000.00E 00
6.18E-022.09E 00
8.86E-023.13E 00
1.50E-015.09E 00
2.18E-017.18E 00
2.71E-019.51E 00
3.28E-011.20E 01
3.62E-011.34E 01
4.00E-011.57E 01
4.48E-011.84E 01
4.83E-012.08E 01
5.15E-012.31E 01
5.44E-012.53E 01
5.71E-012.74E 01
5.97E-012.92E 01
6.13E-013.07E 01
6.51E-013.33E 01
6.75E-013.56E 01
7.07E-013.82E 01
7.47E-014.08E 01
7.79E-014.31E 01
8.13E-014.55E 01
8.59E-014.86E 01
9.79E-015.58E 01
];
data5=[1.68E-032.07E 01
5.83E-022.12E 01
1.34E-012.16E 01
1.99E-012.19E 01
2.50E-012.25E 01
3.07E-012.33E 01
3.60E-012.46E 01
4.12E-012.59E 01
4.71E-012.79E 01
5.08E-012.93E 01
5.65E-013.19E 01
5.94E-013.36E 01
6.42E-013.60E 01
6.98E-013.89E 01
7.47E-014.18E 01
7.84E-014.43E 01
9.79E-015.58E 01
];
data6=[1.66E-032.08E 01
6.41E-021.92E 01
1.32E-011.73E 01
1.84E-011.58E 01
2.46E-011.43E 01
2.87E-011.31E 01
3.63E-011.10E 01
4.28E-019.18E 00
5.12E-017.01E 00
5.77E-015.44E 00
6.53E-014.14E 00
7.29E-012.84E 00
8.40E-011.46E 00
9.21E-017.93E-01
1.00E 00-1.23E-01
];


figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(data1(:,1),data1(:,2),'LineWidth',3,'LineStyle','--');
plot(data2(:,1),data2(:,2),'LineWidth',3,'LineStyle','--');
plot(data3(:,1),data3(:,2),'LineWidth',3,'LineStyle','--');
plot(data4(:,1),data4(:,2),'LineWidth',3);
plot(data5(:,1),data5(:,2),'LineWidth',3);
plot(data6(:,1),data6(:,2),'LineWidth',3);
%[m,n]=ginput(1)
M=[0;0.9750;0;1;0.3303];
N=[0;55.8318;20.7656;0;12.2968];
str=['A';'B';'C';'D';'E'];
text(axes1,M,N,str,'FontSize',16)
ylabel('$P$','Interpreter','latex');
xlabel('$x_B$','FontSize',20,'Interpreter','latex');
title('氯仿-乙醚蒸气压与液相组成关系');
box(axes1,'on');
hold(axes1,'off');
set(axes1,'FontSize',15,'LineWidth',3);

该图表示氯仿-乙醚体系在温度恒定的情况下,组分数-压强-相组成的变化情况.BC以上区域为液相互溶区间，此时压力条件下不同组分都变成一相；AEDBC区间为气液共存区;AED为气相区。过冷水给大家分析量一组分相图和二个组分的相图的作用以及绘制案例，是为了引出三组分的相图绘制。本期过冷水主要是想和告诉大家三元相图应该怎么绘制。首先先带大家了解一下三元相图的作用和特点，然后在根据要求实现三元相图的绘制。

Matlab-模型.png

三个组分的组成变换可以用质量分数x表示。以质量分数为例。如果在A、B、C三个组分中选两个独立变量，则另外一个就是因变量物化中常用等边三角形来表示三组分变化，组成变换有下列关系：

Matlab-模型1.png

过冷水和大家一起学习一下等腰三角形上的相图点是怎么读的。过状态点D想三角形的底边AB做平行线，于AC变的交点可以确定组分C的含量为c;过状态点D做AC的平行线交BC边于b点可以确定B的含量；过状态点D做BC的平行线交AB与a点，可以确定A的含量，所以D的组成为(a,b,c);这是过冷水在实践中摸索出来第一种方法，书上给的方法是：

Matlab-模型2.png

过状态点D向底边AB作平行于其它两边的直线，其分别于底边交于a,b两点。此两点将底边分为三个线段，左边Aa表示B的含量，ab表示C的含量；bB表示B的含量。讲完了三元相图的基本知识点，现在进入正题和说一说应该怎么程序化实现三元相图的绘制

需要申明的一点是三角形图不就是在直角坐标系中绘制了一个三角形吗。将直角坐标系改成了三角形坐标系。自然坐标就需要从直角坐标改成三角形坐标下。对应的转换公式为：

利用转换完后的的坐标进行绘图即可得到三角形对应的图。代码可以分为：

1 进行坐标转换

A=[0.7;0.3];B=[0;0.4];C=[0.3;0.3];
Total = (A B C);
fA = A./Total;fB = B./Total;fC = 1-(fA fB);
y = fC*sin(deg2rad(60))
x = 1 - fA - y*cot(deg2rad(60))
[x, i] = sort(x);
y = y(i);

2 设置三角坐标系

direction = 'clockwise';
percentage = false;
xoffset = 0.25;
yoffset = 0.01;
 
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;
tc = get(cax,'xcolor');
ls = get(cax,'gridlinestyle');
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');
 
set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
         'DefaultTextFontName',   get(cax, 'FontName'), ...
         'DefaultTextFontSize',   get(cax, 'FontSize'), ...
         'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
         'DefaultTextUnits','data')
hold on;
plot ([0 1 0.5 0],[0 0 sin(1/3*pi) 0], 'color', tc, 'linewidth',1,'handlevisibility','off');
set(gca, 'visible', 'off');
patch('xdata', [0 1 0.5 0], 'ydata', [0 0 sin(1/3*pi) 0],'edgecolor',tc,'facecolor',get(gca,'color'),'handlevisibility','off');
majorticks = linspace(0, 1, 10   1);
majorticks = majorticks(1:end-1);
multiplier = 1;
labels = num2str(majorticks(end:-1:1)'*multiplier);
zerocomp = zeros(size(majorticks)); % represents zero composition
[lxc, lyc] = terncoords(1-majorticks, majorticks, zerocomp);
text(lxc 0.05, lyc-0.025, [repmat('  ', length(labels), 1) labels]);
[lxb, lyb] = terncoords(majorticks, zerocomp, 1-majorticks); % fB = 1-fA
text(lxb-0.115, lyb-0.065, labels, 'VerticalAlignment', 'Top');
[lxa, lya] = terncoords(zerocomp, 1-majorticks, majorticks);
text(lxa-0.035, lya 0.09, labels);
nlabels = length(labels)-1;
for i = 1:nlabels
        plot([lxa(i 1) lxb(nlabels - i   2)], [lya(i 1) lyb(nlabels - i   2)], ls, 'color', tc, 'linewidth',0.25,...
           'handlevisibility','off');
        plot([lxb(i 1) lxc(nlabels - i   2)], [lyb(i 1) lyc(nlabels - i   2)], ls, 'color', tc, 'linewidth',0.25,...
           'handlevisibility','off');
        plot([lxc(i 1) lxa(nlabels - i   2)], [lyc(i 1) lya(nlabels - i   2)], ls, 'color', tc, 'linewidth',0.25,...
           'handlevisibility','off');
end;
 
% Reset defaults
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits', fUnits );

3 使用plot 绘图

q = plot(x, y);
hold on
A=[0.3;0.2];B=[0.4;0.7];C=[0.3;0.1];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);

4 在原图像上进行叠加绘图

A=[0.2;0.15];B=[0.7;0.75];C=[0.1;0.1];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
A=[0.2;0];B=[0.7;0.7];C=[0.1;0.3];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
A=[0.15;0];B=[0.75;0.7];C=[0.1;0.3];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
A=[0.15;0];B=[0.75;0.95];C=[0.1;0.05];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
A=[0.7;0.35];B=[0;0.45];C=[0.3;0.2];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
 
A=[0.35;0.24];B=[0.45;0.68];C=[0.2;0.08];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
A=[0.24;0.5];B=[0.68;0.5];C=[0.08;0];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
A=[0.2;0.24];B=[0.7;0.68];C=[0.1;0.08];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
A=[0.15 0.35];B=[0.75 0.75];C=[0.1;0];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
 
 
A=[0.15 0.35];B=[0.75  0.65];C=[0.1 0];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
ternplot([0.15 0.35],[0.75  0.65],[0.1 0])
 
A=[0.15 0.1];B=[0.75  0.9];C=[0.1 0];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);
ternplot([0.15 0.35],[0.75  0.65],[0.1 0])
 
A=[0.2;0];B=[0.7;0.5];C=[0.1;0.5];
[fA, fB, fC] = fractions(A, B, C);
[x, y] = terncoords(fA, fB, fC);
[x, i] = sort(x);
y = y(i);
q = plot(x, y);

最终我们就得到了上图的三元相图，实际三元相图的核心是坐标系的转换，除此之外没有特别神秘的地方，同理可的常见的：半等高图、三维轮廓图，直角立体图经过坐标转换都可以实现三角坐标系上的绘制。由于做起来比较麻烦，放在本期不太合适，本次过冷水只是和大家分享最基础的三元相图的绘制，在本期的基础上会陆续实现上述提及的在直角转三角坐标的具体实现案例。