MATLAB-输出在目标温度、目标压力下平衡时的气相摩尔分量及气液组成

MATLAB-输出在目标温度、目标压力下平衡时的气相摩尔分量及气液组成

MATLAB-输出在目标温度、目标压力下平衡时的气相摩尔分量及气液组成

MATLAB-输出在目标温度、目标压力下平衡时的气相摩尔分量及气液组成

MATLAB-输出在目标温度、目标压力下平衡时的气相摩尔分量及气液组成

MATLAB-输出在目标温度、目标压力下平衡时的气相摩尔分量及气液组成

第三题
%第三题 求气相摩尔分量
function [V]=GPMF(T,P)
A=[ 1  1.297  33.2    -0.22  65    0.637503;
    2  4.600  190.6   0.008  99    0.093293;
    3  4.884  305.4   0.098  148   0.017343;
    4  22.05  647.3   0.344  56    0.168547;
    
    5  11.28  405.6   0.250  72.5  0.002714;
    6  8.937  373.2   0.100  98.5  0.012105;
    7  8.309  324.6   0.12   81    0.000000283;
    8  3.384  460.4   0.227  306   0.068494;];
    Bi=A(:,2);%Pci
    Ci=A(:,3);%Tci
    D=A(:,4);%Wi
    Z=A(:,6);%Zi
    B=P./Bi;%公式2
    C=T./Ci;%公式3
    %公式1
    K=exp(5.37*(1+D).*(1-1./C))./B;%Ki
    %V的初值
    V0=0.9;
    while(1)
        %公式4
        F=sum((Z.*(K-1))./(1+((K-1)*V0)));
        %公式5
        F1=sum(-(Z.*(K-1).*(K-1))./((1+(K-1)*V0).*(1+(K-1)*V0)));
        %公式6
        V=V0-F/F1;
        %公式7
        dif=abs((V-V0)/V);
        if (dif<0.001)
            break
        else
            V0=V;
        end
    end
end
第四题
%第四题 求气相组成yV、液相组成xL(摩尔分数)
function [YV,XL]=MF(T,P,V)
A=[ 1  1.297  33.2    -0.22  65    0.637503;
    2  4.600  190.6   0.008  99    0.093293;
    3  4.884  305.4   0.098  148   0.017343;
    4  22.05  647.3   0.344  56    0.168547;
    5  11.28  405.6   0.250  72.5  0.002714;
    6  8.937  373.2   0.100  98.5  0.012105;
    7  8.309  324.6   0.12   81    0.000000283;
    8  3.384  460.4   0.227  306   0.068494;];
    Bi=A(:,2);%Pci
    Ci=A(:,3);%Tci
    D=A(:,4);%Wi
    Z=A(:,6);%Zi
    B=P./Bi;%公式2
    C=T./Ci;%公式3
    %公式1
    K=exp(5.37*(1+D).*(1-1./C))./B;%Ki
    %公式8
    Yi=Z.*K./(1+(K-1).*V);
    %公式9
    YV=Yi./sum(Yi);
    %公式10
    Xi=Z./(1+(K-1).*V);
    XL=Xi./sum(Xi);
end
第五题
%第五题 求气相混合物平均引力系数amV、平均斥力系数bmV,液相混合物平均引力系数amL、平均斥力系数bmL
function [amV,bmV,amL,bmL]=AC(T,P,YV,XL)
A=[ 1  1.297  33.2    -0.22  65    0.637503;
    2  4.600  190.6   0.008  99    0.093293;
    3  4.884  305.4   0.098  148   0.017343;
    4  22.05  647.3   0.344  56    0.168547;
    5  11.28  405.6   0.250  72.5  0.002714;
    6  8.937  373.2   0.100  98.5  0.012105;
    7  8.309  324.6   0.12   81    0.000000283;
    8  3.384  460.4   0.227  306   0.068494;];
    Bi=A(:,2);%Pci
    Ci=A(:,3);%Tci
    D=A(:,4);%Wi
    Vi=A(:,5);%Vi
    B=P./Bi;%公式2
    C=T./Ci;%公式3
    %公式12 
    %各组分引力系数
    a=0.45724*82.06*82.06*C.*C./B;
    %公式13
    %各组分斥力系数
    b=0.07780*82.06*C./B;
    %公式15
    m=0.37464+1.5422*D-0.26992*D.*D;
    %公式14
    %各组分温度可调函数
    ai=power((1+m.*(1-power(C,0.5))),2);
    %公式17和18
    k=zeros(8,8);
    for i=1:8
        for j=1:8
            if i==j
                k(i,j)=0;
            else
                k(i,j)=1-power((2*power(Vi(i),1/6).*power(Vi(j),1/6))./(power(Vi(i),1/3)+power(Vi(j),1/3)),exp(1));
            end
        end
    end
    bmV=sum(YV.*b);
    bmL=sum(XL.*b);
    amV=0;
    amL=0;
    for i=1:8
        for j=1:8
            amV=amV+YV(i).*YV(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
            amL=amL+XL(i).*XL(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
        end
    end
end
第六题
%第六题 气相混合物偏差系数ZmV,液相混合物偏差系数ZmL 
function [ZmV,ZmL]=MDC(T,P,amV,bmV,amL,bmL)
    %ZMV
    %公式23
    AmV=amV.*P./(82.06*82.06*T.*T);
    %公式24
    BmV=bmV.*P./(82.06*T);
    ZmV0=1;
    while(1)
        %公式25
        F_ZmV=power(ZmV0,3)-(1-BmV).*power(ZmV0,2)+(AmV-2*BmV-3*power(BmV,2)).*ZmV0-(AmV.*BmV-power(BmV,2)-power(BmV,3));
        %公式26
        F1_ZmV=3*power(ZmV0,2)-2*(1-BmV).*ZmV0+(AmV-2*BmV-3*power(BmV,2)).*ZmV0;
        %公式27
        ZmV=ZmV0-F_ZmV./F1_ZmV;
        %公式28
        dif1=abs((ZmV-ZmV0)/ZmV);
        if (dif1<0.0005)
            break
        else
            ZmV0=ZmV;
        end
    end
    
    %ZmL
    %公式29
    AmL=amL.*P./(82.06*82.06*T.*T);
    %公式30
    BmL=bmL.*P./(82.06*T);
    ZmL0=1;
    while(1)
        %公式31
        F_ZmL=power(ZmL0,3)-(1-BmL).*power(ZmL0,2)+(AmL-2*BmL-3*power(BmL,2)).*ZmL0-(AmL.*BmL-power(BmL,2)-power(BmL,3));
        %公式32
        F1_ZmL=3*power(ZmL0,2)-2*(1-BmL).*ZmL0+(AmL-2*BmL-3*power(BmL,2)).*ZmL0;
        %公式33
        ZmL=ZmL0-F_ZmL./F1_ZmL;
        %公式34
        dif2=abs((ZmL-ZmL0)/ZmL);
        if (dif2<0.0005)
            break
        else
            ZmL0=ZmL;
        end
    end
end
第七题
%第七题 气相中各组分逸度fiV,液相中各组分逸度fiL
function [fiV,fiL]=FOEC(T,P,ZmV,ZmL,a,ai,k,YV,XL,bmV,bmL,AmV,AmL,BmV,BmL,amV,amL)
    A=[ 1  1.297  33.2    -0.22  65    0.637503;
    2  4.600  190.6   0.008  99    0.093293;
    3  4.884  305.4   0.098  148   0.017343;
    4  22.05  647.3   0.344  56    0.168547;
    5  11.28  405.6   0.250  72.5  0.002714;
    6  8.937  373.2   0.100  98.5  0.012105;
    7  8.309  324.6   0.12   81    0.000000283;
    8  3.384  460.4   0.227  306   0.068494;];
    Bi=A(:,2);%Pci
    Ci=A(:,3);%Tci
    D=A(:,4);%Wi
    Vi=A(:,5);%Vi
    B=P./Bi;%公式2
    C=T./Ci;%公式3
    %各组分斥力系数
    b=0.07780*82.06*C./B;
    %fiV
    %公式36
    HiV=0;
    for j=1:8
        HiV=HiV+YV(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
    end
    fiV=YV.*P.*exp(b./bmV.*(ZmV-1)-log(ZmV-BmV)-AmV./(2.*sqrt(2).*BmV).*(2.*HiV./amV-b./bmV).*(log((ZmV+2.414.*BmV)./(ZmV-0.414.*BmV))));
    %fiL
    HiL=0;
    for j=1:8
        HiL=HiL+XL(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
    end
    fiL=XL.*P.*exp(b./bmL.*(ZmL-1)-log(ZmL-BmL)-AmL./(2.*sqrt(2).*BmL).*(2.*HiL./amL-b./bmL).*(log((ZmL+2.414.*BmL)./(ZmL-0.414.*BmL))));
end
第八题
%第八题 热力学平衡检验
function [V,YV,XL]=TMET(T,P)
    A=[ 1  1.297  33.2    -0.22  65    0.637503;
    2  4.600  190.6   0.008  99    0.093293;
    3  4.884  305.4   0.098  148   0.017343;
    4  22.05  647.3   0.344  56    0.168547;
    5  11.28  405.6   0.250  72.5  0.002714;
    6  8.937  373.2   0.100  98.5  0.012105;
    7  8.309  324.6   0.12   81    0.000000283;
    8  3.384  460.4   0.227  306   0.068494;];
    Bi=A(:,2);%Pci
    Ci=A(:,3);%Tci
    D=A(:,4);%Wi
    Vi=A(:,5);%Vi
    Z=A(:,6);%Zi
    B=P./Bi;%公式2
    C=T./Ci;%公式3
    while(1)
        %公式1
        K=exp(5.37*(1+D).*(1-1./C))./B;%Ki
        %V的初值
        V0=0.9;
        while(1)
            %公式4
            F=sum((Z.*(K-1))./(1+((K-1)*V0)));
            %公式5
            F1=sum(-(Z.*(K-1).*(K-1))./((1+(K-1)*V0).*(1+(K-1)*V0)));
            %公式6
            V=V0-F/F1;
            %公式7
            dif=abs((V-V0)/V);
            if (dif<0.001)
                break
            else
                V0=V;
            end
        end
        %公式8
        Yi=Z.*K./(1+(K-1).*V);
        %公式9
        YV=Yi./sum(Yi);
        %公式10
        Xi=Z./(1+(K-1).*V);
        XL=Xi./sum(Xi);
        %公式12 
        %各组分引力系数
        a=0.45724*82.06*82.06*C.*C./B;
        %公式13
        %各组分斥力系数
        b=0.07780*82.06*C./B;
        %公式15
        m=0.37464+1.5422*D-0.26992*D.*D;
        %公式14
        %各组分温度可调函数
        ai=power((1+m.*(1-power(C,0.5))),2);
        %公式17和18
        k=zeros(8,8);
        for i=1:8
            for j=1:8
                if i==j
                    k(i,j)=0;
                else
                    k(i,j)=1-power((2*power(Vi(i),1/6).*power(Vi(j),1/6))./(power(Vi(i),1/3)+power(Vi(j),1/3)),exp(1));
                end
            end
        end
        bmV=sum(YV.*b);
        bmL=sum(XL.*b);
        amV=0;
        amL=0;
        for i=1:8
            for j=1:8
                amV=amV+YV(i).*YV(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
                amL=amL+XL(i).*XL(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
            end
        end
         %ZMV
        %公式23
        AmV=amV.*P./(82.06*82.06*T.*T);
        %公式24
        BmV=bmV.*P./(82.06*T);
        ZmV0=1;
        while(1)
            %公式25
            F_ZmV=power(ZmV0,3)-(1-BmV).*power(ZmV0,2)+(AmV-2*BmV-3*power(BmV,2)).*ZmV0-(AmV.*BmV-power(BmV,2)-power(BmV,3));
            %公式26
            F1_ZmV=3*power(ZmV0,2)-2*(1-BmV).*ZmV0+(AmV-2*BmV-3*power(BmV,2)).*ZmV0;
            %公式27
            ZmV=ZmV0-F_ZmV./F1_ZmV;
            %公式28
            dif1=abs((ZmV-ZmV0)/ZmV);
            if (dif1<0.0005)
                break
            else
                ZmV0=ZmV;
            end
        end
         %ZmL
        %公式29
        AmL=amL.*P./(82.06*82.06*T.*T);
        %公式30
        BmL=bmL.*P./(82.06*T);
        ZmL0=1;
        while(1)
            %公式31
            F_ZmL=power(ZmL0,3)-(1-BmL).*power(ZmL0,2)+(AmL-2*BmL-3*power(BmL,2)).*ZmL0-(AmL.*BmL-power(BmL,2)-power(BmL,3));
            %公式32
            F1_ZmL=3*power(ZmL0,2)-2*(1-BmL).*ZmL0+(AmL-2*BmL-3*power(BmL,2)).*ZmL0;
            %公式33
            ZmL=ZmL0-F_ZmL./F1_ZmL;
            %公式34
            dif2=abs((ZmL-ZmL0)/ZmL);
            if (dif2<0.0005)
                break
            else
                ZmL0=ZmL;
            end
        end

        %fiV
        %公式36
        HiV=0;
        for j=1:8
            HiV=HiV+YV(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
        end
        fiV=YV.*P.*exp(b./bmV.*(ZmV-1)-log(ZmV-BmV)-AmV./(2.*sqrt(2).*BmV).*(2.*HiV./amV-b./bmV).*(log((ZmV+2.414.*BmV)./(ZmV-0.414.*BmV))));

        HiL=0;
        for j=1:8
            HiL=HiL+XL(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
        end
        fiL=XL.*P.*exp(b./bmL.*(ZmL-1)-log(ZmL-BmL)-AmL./(2.*sqrt(2).*BmL).*(2.*HiL./amL-b./bmL).*(log((ZmL+2.414.*BmL)./(ZmL-0.414.*BmL))));
        dif2=abs(fiV-fiL);
        
        if (dif2<0.0001)
            break;
        else
            K=YV.*fiL./(XL.*fiV);
        end
    end
end
上一篇:Codeforces Round #735 (Div. 2) C. Mikasa


下一篇:Codeforces Round #735 (Div. 2)