profit.m

These are matlab and simulink files to model the membrane crystallization system, including the matl

function f=profit(y)
global sigma Tmem Mentercry P rhoLac kV kB kG V Tcry Wentercry Wsatcry Wleavecry Mleavecry Lm epsilon

Tmem = y(1);
Mentercry = y(2);

P = 1;
rhoLac = 1550;
rhoWat = 1000;
Mlac=342.30;
kV=0.43;
kB =1.8e11/60;
kG = 2.27e-8;
V=20;
Tcry=25;
Wentercry=( 0.0053*(Tmem+273)^2-2.825*(Tmem+273)+390.06 )/100;
Wsatcry = ( 0.0053*(Tcry+273)^2-2.825*(Tcry+273)+390.06 )/100;
    
Wleavecry = (Mentercry*Wentercry - P)/(Mentercry-P);
Mleavecry = Mentercry - P;

%initial x0,x1,x2,x3,x4,Fout,epsilon
x0=[11.637,2.694,1.29,0.9272,0.888,0.001618,0.75];
options=optimset('Display','iter');   % Option to display output
[x,fval] = fsolve(@Calc_epsilon,x0,options);
m0=x(1)/sigma^3;
m1=x(2)/sigma^2;
m2=x(3)/sigma;
m3=x(4);
m4=x(5)*sigma;
epsilon=x(7);

    if (epsilon>=0.6)
    Lm=m4/m3;
    display(epsilon);
    Lmf=Lm*1e6;
    display(Lmf);
    LmL = Lm/0.0002;
    Value=0.8;
        if(LmL<0.2)
            Cpro=Value*(1-0.3);
        elseif((LmL>=0.2)&&(LmL<0.6))
            Cpro=Value*(1-0.2);
        elseif((LmL>=0.6)&&(LmL<0.85))
            Cpro=Value*(1-0.1);
        elseif((LmL>=0.85)&&(LmL<1.15))
            Cpro=Value;
        elseif((LmL>=1.15)&&(LmL<1.4))
            Cpro=Value*(1+0.1);
        elseif((LmL>=1.4)&&(LmL<1.8))
            Cpro=Value*(1+0.2);
        elseif(LmL>=1.8)
            Cpro=Value*(1+0.3);
        end
                
    Cequ = ((20.875*(0.75*V^(1/3))^2-0.776*0.75*V^(1/3)+38.25)*1000)/10;

    Wfeed=0.05;
    Feed = P/Wfeed;
    Pmem=60;
    Pcry=1;
    Wmix= (Mleavecry*Wleavecry+Feed*Wfeed)/(Mleavecry+Feed);
    rhomix = 1/(Wmix/rhoLac+(1-Wmix)/rhoWat);
    VFpump =(Feed+Mleavecry)/rhomix;
    Pumppower =VFpump*(Pmem-Pcry)*1e5/(0.6*60*1000);%kW
    Cpump = 0.08 * Pumppower*8000;
    Tfeed=40;
    Cproduct = Cpro*1*8000*60; 
    MixT = (Tcry*Mleavecry+Tfeed*Feed)/(Feed+Mleavecry);
    Cfeed = Feed*0.023*8000*60;   
    Tbh=(MixT+Tmem)/2+273.15;
    Cp= ((1-Wmix)*4.18+Wmix*(0.0076*Tbh-0.2328));
    DeltaT=Tmem-MixT; % K;
    Heat = (Feed+Mleavecry)*Cp*DeltaT/(60);%kW
                
    Cheat = 0.01*Heat*8000;
 
    f = -(Cproduct-Cfeed-Cheat-Cequ-Cpump);
    pro=-f;
    display(pro);
    
    else
    f=0;
    %error('epsilon is out of the range 0.6~1');
    end