# How to draw binary phase diagram ??

19 views (last 30 days)
관 윤 on 17 Jan 2021
Answered: SaiDileep Kola on 20 Jan 2021
clear all;
clc;
%%
G1_s=300; % 1 element pure gibbs energy
G2_s=500; % 2 element pure gibbs energy
Om_s=5000; % 1,2 element mixing enthalpy
Om_l=3000; % 1,2 element mixing enthalpy
R=1.987; % gas constant
Tm_1=600; % element1 melting temparture
Tm_2=900; % element2 melting temparture
Hm_1=2000; % element1 Heat of fusion
Hm_2=1300; % element2 Heat of fusion
Tset=[0:700:1400]; %Temparture Setting
%%
for i =1:length(Tset)
T=Tset(i);
syms x a b G_l(x) G_s(x) dG_s(x) dG_l(x) Sola Solb
G1_l=G1_s+Hm_1/Tm_1*(Tm_1-T); % element1 liquid
G2_l=G2_s+Hm_2/Tm_2*(Tm_2-T); % element2 liquid
G_s(x) = G1_s*x+G2_s*(1-x)+Om_s*x*(1-x)+R*T*(x*log(x)+(1-x)*log(1-x)); % x-G function of Gibbs energy of solid
G_l(x) = G1_l*x+G2_l*(1-x)+Om_l*x*(1-x)+R*T*(x*log(x)+(1-x)*log(1-x)); % x-G function of Gibbs energy of liquid
dG_s(x)=diff(G_s(x),x);
dG_l(x)=diff(G_l(x),x);
denk1=(((G_s(b)-G_l(a))/(b-a))==(dG_s(a))); %eq 1
denk2=(dG_s(b)==dG_l(a)); %eq 2
[Sola,Solb] = vpasolve([denk1 denk2],[a b],[0.1;0.9]);
a=Sola;
b=Solb;
disp([T,a,b])
end
i think there is something wrong in my code and i want to ask you what is wrong
i want to plot T versus a graph and T versus b

SaiDileep Kola on 20 Jan 2021
Since the solution from vpasolve is complex in this case, try to plot as real vs complex values as described here.