Since the solution from vpasolve is complex in this case, try to plot as real vs complex values as described here.
How to draw binary phase diagram ??
17 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
采纳的回答
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Thermal Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!