Download presentation
Presentation is loading. Please wait.
1
Phase Diagram (Subregular solution) 신소재공학과 20041008 송양희 Liquid 의 subregular solution 에서의 Gibbs energy
2
BCC 와 FCC solution
3
Regression subroutine mid_1(al1,bl1) double precision origin(39,3) double precision A,B,C,D,E,al1,bl1,r1 integer i open(unit=10,file='abc1.txt',status='old') do 100 i=1,39 read(10,*) origin(i,1), origin(i,2) 100 continue do 110 i=1,39 origin(i,3)=origin(i,2)/origin(i,1)/(1-origin(i,1)) 110 continue A=0 B=0 C=0 D=0 E=0 do 120 i=1,39 A=A+origin(i,1)*origin(i,3) B=B+origin(i,1) C=C+origin(i,3) D=D+origin(i,1)*origin(i,1) E=E+origin(i,3)*origin(i,3) 120 continue bl1=(39*A-B*C)/(39*D-B*B) al1=C/39-bl1*B/39 return end
4
Regression 결과
5
Eutectic Point 찾기
6
각 chemical potential
7
각 평형함수 subroutine F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl,Xb,R1) double precision al1,al2,bl1,bl2,ab1,ab2,bb1,bb2 double precision T,Xl,Xb,R1 double precision ex,ex1 ex=16250-12.5*T+8.3144*T*log(Xl/Xb)+(al1+al2*T)*(1-Xl)*(1-Xl) ex1=ex+2*(bl1+bl2*T)*(1-Xl)*(1-Xl)*Xl R1=ex1-(ab1+ab2*T)*(1-Xb)*(1-Xb)-2*(bb1+bb2*T)*(1-Xb)*(1- Xb)*Xb return end
8
Jacobian matrix call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T+0.0001,Xl,Xb,R1) call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T-0.0001,Xl,Xb,R2) J(1,1)=(R1-R2)/2/0.0001 call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl+0.0001,Xb,R1) call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl-0.0001,Xb,R2) J(1,2)=(R1-R2)/2/0.0001 call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl,Xb,R1) call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl,Xb,R2) J(1,3)=(R1-R2)/2/0.0001 callF_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl,Xb+0.0001,R1) call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl,Xb-0.0001,R2) J(1,4)=(R1-R2)/2/0.0001
9
Eutectic Point ** 처음 임의의 값 지정 T=1000 Xl=0.5 Xb=0.8 Xf=0.2 number=0 ** jacobian matrix 입력 ** 역함수 구하기 ** 새로운 값 구하기 call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T,Xl,Xb,R1) call F_2(af1,af2,bf1,bf2,ab1,ab2,bb1,bb2,T,Xf,Xb,R2) call F_3(al1,al2,bl1,bl2,af1,af2,bf1,bf2,T,Xl,Xf,R3) call F_4(ab1,ab2,bb1,bb2,af1,af2,bf1,bf2,T,Xb,Xf,R4) T1=T-RJ(1,1)*R1-RJ(1,2)*R2-RJ(1,3)*R3-RJ(1,4)*R4 Xl1=Xl-RJ(2,1)*R1-RJ(2,2)*R2-RJ(2,3)*R3-RJ(2,4)*R4 Xf1=Xf-RJ(3,1)*R1-RJ(3,2)*R2-RJ(3,3)*R3-RJ(3,4)*R4 Xb1=Xb-RJ(4,1)*R1-RJ(4,2)*R2-RJ(4,3)*R3-RJ(4,4)*R4 ** 처음에 구한값이 0~1 를 벗어났을때 임의로 값 보정 100 if(Xl1.GE.1) then Xl1=Xl1-1 go to 100 elseif(Xl1.LE.0) then Xl1=-Xl1 go to 100 endif 200 if(Xf1.GE.1) then Xf1=Xf1-1 go to 200 elseif(Xf1.LE.0) then Xf1=-Xf1 go to 200 endif 300 if(Xb1.GE.1) then Xb1=Xb1-1 go to 300 elseif(Xb1.LE.0) then Xb1=-Xb1 go to 300 endif
10
Eutectic Point ** 식에 넣고 0 이 나오는지를 통해 error 계산 call F_1(al1,al2,bl1,bl2,ab1,ab2,bb1,bb2,T1,Xl1,Xb1,R1) error=dabs(R1) number=number+1 T=T1 Xl=Xl1 Xf=Xf1 Xb=Xb1 ** 일정횟수가 지나거나 error 가 작을때에만 다시 돌아가 지 않음 if(number.LT.100000.AND. error.Gt.0.00000001) then go to 20 endif ** 100000 이 되어도 값을 못찾으면 이전 온도에서 약간의 변화를 주어 보정후 다시 계산 if(number.EQ.100000) then print *, "I can't solution" else print *, "========RESULT============" print *, T,Xl,Xf,Xb print *, "error",error endif
11
다른 온도에서 equilibrium FCC_liquid equilibriumBCC_liquid equilibriumFCC_BCC equilibrium
12
소스 파일 ** FCC_liquid equilibrium ** 처음의 임의의 mole fraction xe1=Xf1 xe2=Xl1 ** 각 온도에서 loop 를 통해 mole fration 구하기 do 11 Te=T1,1500 number=0 ** jacobian matrix 직접입력 ** 역함수 구하기 ** 새로운 값 구하기 call F_3(al1,al2,bl1,bl2,af1,af2,bf1,bf2,Te,Xe2,Xe1,R1) call F_6(al1,al2,bl1,bl2,af1,af2,bf1,bf2,Te,Xe2,Xe1,R2) xn1=xe1-RJ(1,1)*R1-RJ(1,2)*R2 xn2=xe2-RJ(2,1)*R1-RJ(2,2)*R2 ** 처음에 구한값이 0~1 를 벗어났을때 임의로 값 보정 101 if(xn1.GE.1) then xn1=xn1-1 go to 101 elseif(xn1.LE.0) then xn1=-xn1 go to 101 endif 201 if(xn2.GE.1) then xn2=xn2-1 go to 201 elseif(xn2.LE.0) then xn2=-xn2 go to 201 endif ** 식에 넣고 0 이 나오는지를 통해 error 계산 call F_3(al1,al2,bl1,bl2,af1,af2,bf1,bf2,Te,Xn2,Xn1,R1) error=dabs(R1) number=number+1 xe1=xn1 xe2=xn2 ** 일정횟수가 지나거나 error 가 작을때에만 다시 돌아가지 않음 if(number.LT.100000.AND. error.Gt.0.0000000001) then go to 21 endif ** 결과 제출 open(unit=71,file="result1",status="new") if(number.NE.100000) then write(71,*) Te,xe1,xe2 else write(71,*) Te,"no solution" endif 11 continue
13
RESULT
Similar presentations
