I need help writing a code that can find the common tangent between two gibbs free energy curves at a specific temperature. As of right now, I am able to make the inidvidual curves, but I have no idea how I would go about finding their common tangents and the points those tangents intersect. Below is the code and the figures produced.
#the goal is to plot: 1. G vs composition of B, then determine the minimas when they intersect.
#when the meet at a point on either 100% A or 100% B those are the melting points.
#In the ideal solution W_H=0, but for our regular solutino models this value will be negative.
s_b=28.93 #in J/mol*K This is for ammonia, calculated by H/T so for ammonia: 5653/195.4, for iso: 5280/184 = 11.64
s_a=22.043
H_a=6020 #6020 water, 5653 for ammonia 5280 for iso
H_b=5653 #ammonia
X_b=np.linspace(0.0000000000001,1.0)
X_a=np.linspace(1.0,0.0000000000001)
#W_h=np.linspace(-2.262145855e-12,-2.1700789464999998e-11)
N_a=6.022e23
z=4
R=8.3144 #J/molK
def G_a(H_a, T, s_a):
return H_a-T*s_a
def G_b(H_b, T, s_b):
return H_b-T*s_b
def A(G_a, X_a):
return X_a*G_a(H_a, T, s_a)
def B(G_b, X_b):
return X_b*G_b(H_b, T, s_b)
def H_mix(X_a, X_b, W_h):
#This is right DO NOT edit eqn
#This value will either be 0 (ideal) or negative (mixing favored)
#def W_h(N_a, z, w_aa, w_bb, w_ab):
#return 0.5*N_a*z*(2*w_ab-w_aa-w_bb)
return X_a*X_b*W_h
def S_mix(R, X_a, X_b):
#This is Right DO NOT edit eqn
return -R*(X_a*np.log(X_a)+X_b*np.log(X_b))
def G_ideal(R,T,X_a,X_b):
#Correct, DO NOT edit eqn
return R*T*(X_a*np.log(X_a)+X_b*np.log(X_b))
def G_mix(S_mix, H_mix, T):
return H_mix(X_a, X_b, W_h)-T*S_mix(R, X_a, X_b)
def G_ex(G_mix,G_ideal):
return G_mix(S_mix, H_mix, T)-G_ideal(R,T,X_a,X_b)
def G_AB(G_a, G_b, X_a, X_b):
return (X_a*G_a(H_a, T, s_a)) + (X_b*G_b(H_b, T, s_b))
def Ga(A, G_mix): #Water
W_h=-2.1700789464999998e-11
return A(G_a, X_a) + G_mix(S_mix, H_mix, T)
def Gb(B, G_mix): #ammonia
W_h=-2.262145855e-12
return B(G_b, X_b) + G_mix(S_mix, H_mix, T)
def G(X_a, X_b, G_mix, W_h, B, A):
return (X_a*A(G_a, X_a))+(X_b*B(G_b, X_b))+G_mix(S_mix,H_mix, T)
T=1000
W_h=-2
plt.plot(X_b,Ga(A, G_mix),'b-', label='1000 K')
plt.plot(X_b,Gb(B, G_mix),'b-.', label='1000 K')
plt.title('1000 K')
plt.show()
Produces @1000 K: Curves at 1000 K
Can you help?