solving set of equations using iteration

clear all
err=10000; iter =1;
T_mf = 60;
T_mp= 30;
while (err>0.0001)
    iter=iter+1;
    % disp(err);
    "+++++++++++Tempretures +++++++++++";
    T_bf =70; % feed temperature
    T_bp = 20; % permete temperature
    Pi =22/7;
    R =8314; %Ideal gas constant
    T_m = (T_mf + T_mp)/2; % Mean temperature across the memmbrane (C)
    "++++++++++++Flowrates ++++++++++++";
    V_dot_Feed=4.5; %volume flow rate for feed (L/min)
    V_dot_permete=3.65;%volume flow rate for permete (L/min)
    "============================ Properties of The hot Side =============================" ;
    "++++++++++++Dimensions ++++++++++++";
    L_Feed = 0.066; % channel's length
    W_Feed = 0.024; % channel's wedth
    Depth_Feed = 0.005; % channel's depth 		                       
    A_Feed_Cross = W_Feed*Depth_Feed; 		                        
    D_h_Feed = (4*A_Feed_Cross)/(2*(W_Feed+Depth_Feed)); % channel's hydrolic diameter	                         
    A_Feed = Depth_Feed*W_Feed;% channel's area (three channels)
    "============================ Properties of The cold Side =============================" ;
    "++++++++++++Dimensions ++++++++++++";
    L_permete = 0.066; 
    W_permete = 0.024;
    Depth_permete = 0.005; 		                       
    A_permete_Cross = W_permete*Depth_permete; 		                        
    D_h_permete = (4*A_permete_Cross)/(2*(W_permete+Depth_permete)); 	                         
    A_permete = Depth_permete*W_permete;
    "+++++++++++Concentrations ++++++++++";
    Mass_NaCl = 100; %mass concentration of NaCl (g/L)		                          
    %Mass_NaCl_Permeate = 0; 		                          
    Mass_Water = 1000  - Mass_NaCl; %mass concentration of feed water (g/L)
    Mol_w = 18; %molocular weight of feed water
    Mol_NaCl = 58.8; %molocular weight of Nacl g/mol
    n_Water = Mass_Water/Mol_w; %Number of moles in feed water
    n_NaCl = Mass_NaCl/Mol_NaCl;%Number of moles in Nacl 
    n_Total = n_NaCl + n_Water; %Total number of moles in feed		                         
    x_NaCl = Mass_NaCl/1000  %Dr.Atia code                                            
    %x_NaCl = n_NaCl/n_Total; %Nacl mole fraction 		                        
    Gamma_Feed = 1 - (0.5*x_NaCl) - (10*(x_NaCl)^2);
    x_Feed = 1- x_NaCl; %feed water mole fraction 
    C_Feed=0.14 % concentration of salt in the feed stream (g/L)
    "++++++++++ feed Water Velocity +++++++++++";
    pr= [3.321 2.826  2.442  2.14  1.898];%Prandtle number of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    Pr_Feed =interp1(T2,pr,T_bf);
    k2= [0.6879 0.6989  0.7086 0.7171  0.7244];%Thermal conductivity of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    K_Feed =interp1(T2,k2,T_bf);%interpolation to get to get the thermal conductivity between a bove temperature(T2) 
    mu2= [0.000547 0.000466  0.000404 0.000354  0.000314];%dynamic viscosity of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    mu_Feed =interp1(T2,mu2,T_bf);%interpolation to get to get the dynamic viscosity between a bove temperature(T2)
    mu_mean_Feed=1/(557.82468+19.408782*T_mf +0.1360459*T_mf^2 +((-3.1160832*10^-4)*T_mf^3));%dynamic viscosity at the surface of membrane in the hot side
    
    rho2= [988.05 983.21  977.78 971.80 965.32];%density of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    rho_Feed =interp1(T2,rho2,T_bf);%interpolation to get to get the density between a bove temperature(T2)
    v_Feed = (V_dot_Feed*0.001*1/60)/(3*A_Feed_Cross);
    "++++++++++ permete Water Velocity +++++++++++";
    pr= [11.19 9.45  8.09  7  6.13];
    T1 = [5 10  15  20  25];
    Pr_permete =interp1(T1,pr,T_bp);
    k1= [0.571 0.580  0.589 0.598 0.607];
    T1 = [5 10  15  20  25];
    K_permete =interp1(T1,k1,T_bp);
    mu1= [0.001518 0.001306  0.001138 0.001002  0.000890];
    T1 = [5 10  15  20  25];
    mu_permete =interp1(T1,mu1,T_bp);
    mu_mean_permete=1/(557.82468+19.408782*T_mp +0.1360459*T_mp^2 +((-3.1160832*10^-4)*T_mp^3));
    
    rho1= [999.97 999.70  999.10 998.21 997.05];
    T1 = [5 10  15  20  25];
    rho_permete =interp1(T1,rho1,T_bp);
    v_permete = (V_dot_permete*0.001*1/60)/(3*A_permete_Cross);
    "++++++++++ Renolds number and nusselt for feed +++++++++++";
    Re_Feed  = rho_Feed*v_Feed*D_h_Feed/mu_Feed; 
    Nusselt_Laminar_Feed = 1.86*(Re_Feed*Pr_Feed*D_h_Feed/L_Feed)^(0.33);  	                      
    Nusselt_Turblent_Feed= 0.023*(1+(6*D_h_Feed/L_Feed))*((Re_Feed)^0.8)*((Pr_Feed)^0.33)  %Dr. atia code
    %Nusselt_Turblent_Feed = 0.027*((Re_Feed)^(0.8))*(Pr_Feed)^(0.4)*((mu_Feed/mu_mean_Feed))^(0.14);                 
    if Re_Feed <2300
        Nusselt_Feed = Nusselt_Laminar_Feed ;
    else
        Nusselt_Feed =Nusselt_Turblent_Feed; 
    end
    h_Feed= (Nusselt_Feed*K_Feed)/D_h_Feed;%Convective heat transfer coefficient W/m2-K
    h_fg = (1.7535*(T_m+273)+2024.3)*1000; %Enthalpy or Latent Heat of vaporization of water, J/kg		                        
    
    "++++++++++ Renolds number and nusselt for permete +++++++++++";
    Re_permete  = (rho_permete*v_permete*D_h_permete)/mu_permete; 
    Nusselt_Laminar_permete = 1.86*(Re_permete*Pr_permete*D_h_permete/L_permete)^0.33;  	                      
    Nusselt_Turblent_permete = 0.023*(1+(6*D_h_permete/L_permete))*((Re_permete)^0.8)*((Pr_permete)^0.33)  %Dr. atia code
    %Nusselt_Turblent_permete = 0.027*(Re_permete)^(0.8)*(Pr_permete)^(0.33)*((mu_mean_permete/mu_permete))^(0.14);	                     
    if Re_permete<2300
        Nusselt_permete = Nusselt_Laminar_permete ;
    else
        Nusselt_permete =Nusselt_Turblent_permete; 
    end
    h_permete= (Nusselt_permete*K_permete)/D_h_permete;                         	                         
    "============================== Properties of Membrane =============================="; 
    alpha = 0.9;%a factor α is the ratio of Knudsen diffusion to mass diffusion(can vary between 0 and 1)
    delta_m = 153.9*(10^-6);%membrane thickness 		                       
    Epsilon = 0.8;%membrane porosity
    d_pore = 0.45*(10^-6);%mean pore size
    A_membrane = 0.013395;%membrane area		                       
    K_gas = (0.027 + 0.0184)/2;%gas conductivity
    K_mem = 0.25;%polymer membrane conductivity 
    Tau = 1/Epsilon;% membrane tortuosity 		                        
    K_m = ((Epsilon/K_gas)+((1-Epsilon)/K_mem))^(-1);%membrane conductivity
    "================================= Mass Transfer ================================="; 
    P_owf = exp((23.1964)-((3816.44)/((T_mf+273)-46.13)));%vapor pressures of feed sides at the membrane surfaces (pa)
    P_owp =exp((23.1964)-((3816.44)/((T_mp+273)-46.13)));%vapor pressures of permeate sides at the membrane surfaces(pa)
    P_f=101325; %(pa)
    P_p=101325; %(pa)
    
    P_wvp  = exp((23.1964)-((3816.44)/((T_m+273)-46.13))); %partial pressure of water vapors inside the pores
    P_pore = (P_f +P_p)/2;%The total pressure inside the pores and it is assumed to be the average of feed and permeate side pressures(pa)
    P_airpore = P_pore - P_wvp;	%the partial pressure of air inside the membrane pores 	                                      		                      
    Dwa  = 1.895*10^(-5)*(T_m+273)^(2.072);%the diffusivity of water vapors in air."Product of Pressue and Diffusion of water vapour in air [Pa.m^2/s]"                                                             
    D_k = (((3*delta_m*Tau)/(2*Epsilon*d_pore))*((Pi*R*(T_m+273))/(8*Mol_w))^(0.5))^(-1);% Knudsen diffusion coefficients	                   
    D_m  = ((R*(T_m+273)*delta_m*Tau*P_airpore)/(Epsilon*Mol_w*Dwa ))^(-1);%molecular diffusion coefficients
    D_e  = ((alpha/D_k )+((1-alpha)/D_m ))^(-1);%equivalent diffusio coefficient
    %D_e = ((2*Epsilon*d_pore)/(3*delta_m))*(((2*Mol_w)/(pi*R*(T_m+273)))^0.5)%            "M. Essalhi, M. Khayet, Dr atia code"
    J_w = D_e *((P_owf*Gamma_Feed*x_Feed)-P_owp);%is the mass flux of permeate
    Flux_MD=J_w*3600;% is the mass flux of permeate in one 
    Sc=mu_mean_Feed/(rho_Feed*D_e)%the Schmidt number
    Sh_Laminer_Feed=1.86*(Re_Feed*Sc*(D_h_Feed/L_Feed)) %the Sherwood number for Laminer flow
    Sh_Turblent_Feed=0.023*((Re_Feed)^0.8)*(Sc)^0.33 %the Sherwood number for Turblent flow
    if Re_Feed <2300
        Sh = Sh_Laminer_Feed;
    else
        Sh =Sh_Turblent_Feed; 
    end
    k_s=Sh*(D_e/D_h_Feed)% the solute mass transfer coefficient for the diffusive mass transfer through the concentration boundary layer in the feed side
    C_mf=C_Feed*exp(J_w/(k_s*rho_Feed))%concentration at membrane surface
    T_mf2 =(K_m*(T_bp+(h_Feed/h_permete)*T_bf)+((delta_m*(h_Feed*T_bf-J_w*h_fg))))/((K_m)+(h_Feed*(delta_m+(K_m/h_permete))));%The temperatures of the membrane surfaces at feed side 
    T_mp2 =(K_m*(T_bf+(h_permete/h_Feed)*T_bp)+((delta_m*(h_permete*T_bp+J_w*h_fg))))/((K_m)+(h_permete*(delta_m+(K_m/h_Feed))));%The temperatures of the membrane surfaces at permeate side 
    %T_mf2 = ((A_Feed *h_Feed*T_bf)+((K_m/delta_m)*(A_membrane*T_mp2))-(J_w*h_fg*A_membrane))/((((K_m/delta_m)*(A_membrane))+(A_Feed*h_Feed)));
    %T_mp2 = ((A_permete*h_permete*T_bp)+((K_m/delta_m)*(A_membrane*T_mf2))+(J_w*h_fg*A_membrane))/((((K_m/delta_m)*(A_membrane))+(A_permete*h_permete)));
    U=((1/h_Feed)+(1/((K_m/delta_m)+((J_w*h_fg)/(T_mf-T_mp))))+(1/h_permete))^(-1)%The overall heat transfer coefficient between feed and permete
    Q_Feed = h_Feed*(T_bf-T_mf2);%Convective heat transfer (Watt) from the feed stream to membrane surface,
    Q_permete =h_permete*(T_mp2-T_bp);%Convective heat transfer (Watt) from the feed stream to membrane surface,
    Q_c = (K_m/delta_m)*(T_mf2- T_mp2)%Conductive heat transfer through membrane (Watt)  	                      
    Q_v = J_w*h_fg;%Heat carried by the vapors (Evaporative heat transfer) (Watt)
    Q_Membrane = Q_c+ Q_v;%The total heat transfer across the membrane
    EE=(Q_v/Q_Membrane)*100 %The evaporative (thermal) efficiency, EE,
    err1=(abs(T_mf-T_mf2 ));
    err2=(abs(T_mp-T_mp2 ));
    err= err1+err2
    T_mf=T_mf2;
    T_mp=T_mp2;
    if iter>500
        break;
    end
    
end

[Mohammad Abu abbas]

the erroe always increase , I do not know why , i tried many times but with bad results . I am sure that the equation is right but the erroe is increasing , it should be decreasing. any one can help, I will be thankdul

[Torsten]

Don't use iteration - try "fsolve" to solve your equations.

Answer 1

[Torsten]

% Initial guess
T_mf0 = 60;
T_mp0 = 30;
% Call solver
x = fsolve(@fun,[T_mf0,T_mp0]);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

%Output result
T_mf = x(1)

T_mf = 59.4094

T_mp = x(2)

T_mp = 47.5374

function res = fun(x)
T_mf = x(1);
T_mp= x(2);
    "+++++++++++Tempretures +++++++++++";
    T_bf =70; % feed temperature
    T_bp = 20; % permete temperature
    Pi =22/7;
    R =8314; %Ideal gas constant
    T_m = (T_mf + T_mp)/2; % Mean temperature across the memmbrane (C)
    "++++++++++++Flowrates ++++++++++++";
    V_dot_Feed=4.5; %volume flow rate for feed (L/min)
    V_dot_permete=3.65;%volume flow rate for permete (L/min)
    "============================ Properties of The hot Side =============================" ;
    "++++++++++++Dimensions ++++++++++++";
    L_Feed = 0.066; % channel's length
    W_Feed = 0.024; % channel's wedth
    Depth_Feed = 0.005; % channel's depth 		                       
    A_Feed_Cross = W_Feed*Depth_Feed; 		                        
    D_h_Feed = (4*A_Feed_Cross)/(2*(W_Feed+Depth_Feed)); % channel's hydrolic diameter	                         
    A_Feed = Depth_Feed*W_Feed;% channel's area (three channels)
    "============================ Properties of The cold Side =============================" ;
    "++++++++++++Dimensions ++++++++++++";
    L_permete = 0.066; 
    W_permete = 0.024;
    Depth_permete = 0.005; 		                       
    A_permete_Cross = W_permete*Depth_permete; 		                        
    D_h_permete = (4*A_permete_Cross)/(2*(W_permete+Depth_permete)); 	                         
    A_permete = Depth_permete*W_permete;
    "+++++++++++Concentrations ++++++++++";
    Mass_NaCl = 100; %mass concentration of NaCl (g/L)		                          
    %Mass_NaCl_Permeate = 0; 		                          
    Mass_Water = 1000  - Mass_NaCl; %mass concentration of feed water (g/L)
    Mol_w = 18; %molocular weight of feed water
    Mol_NaCl = 58.8; %molocular weight of Nacl g/mol
    n_Water = Mass_Water/Mol_w; %Number of moles in feed water
    n_NaCl = Mass_NaCl/Mol_NaCl;%Number of moles in Nacl 
    n_Total = n_NaCl + n_Water; %Total number of moles in feed		                         
    x_NaCl = Mass_NaCl/1000;  %Dr.Atia code                                            
    %x_NaCl = n_NaCl/n_Total; %Nacl mole fraction 		                        
    Gamma_Feed = 1 - (0.5*x_NaCl) - (10*(x_NaCl)^2);
    x_Feed = 1- x_NaCl; %feed water mole fraction 
    C_Feed=0.14; % concentration of salt in the feed stream (g/L)
    "++++++++++ feed Water Velocity +++++++++++";
    pr= [3.321 2.826  2.442  2.14  1.898];%Prandtle number of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    Pr_Feed =interp1(T2,pr,T_bf);
    k2= [0.6879 0.6989  0.7086 0.7171  0.7244];%Thermal conductivity of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    K_Feed =interp1(T2,k2,T_bf);%interpolation to get to get the thermal conductivity between a bove temperature(T2) 
    mu2= [0.000547 0.000466  0.000404 0.000354  0.000314];%dynamic viscosity of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    mu_Feed =interp1(T2,mu2,T_bf);%interpolation to get to get the dynamic viscosity between a bove temperature(T2)
    mu_mean_Feed=1/(557.82468+19.408782*T_mf +0.1360459*T_mf^2 +((-3.1160832*10^-4)*T_mf^3));%dynamic viscosity at the surface of membrane in the hot side
    
    rho2= [988.05 983.21  977.78 971.80 965.32];%density of feed water at different temperatures, W/m-K
    T2 = [50 60  70  80  90];
    rho_Feed =interp1(T2,rho2,T_bf);%interpolation to get to get the density between a bove temperature(T2)
    v_Feed = (V_dot_Feed*0.001*1/60)/(3*A_Feed_Cross);
    "++++++++++ permete Water Velocity +++++++++++";
    pr= [11.19 9.45  8.09  7  6.13];
    T1 = [5 10  15  20  25];
    Pr_permete =interp1(T1,pr,T_bp);
    k1= [0.571 0.580  0.589 0.598 0.607];
    T1 = [5 10  15  20  25];
    K_permete =interp1(T1,k1,T_bp);
    mu1= [0.001518 0.001306  0.001138 0.001002  0.000890];
    T1 = [5 10  15  20  25];
    mu_permete =interp1(T1,mu1,T_bp);
    mu_mean_permete=1/(557.82468+19.408782*T_mp +0.1360459*T_mp^2 +((-3.1160832*10^-4)*T_mp^3));
    
    rho1= [999.97 999.70  999.10 998.21 997.05];
    T1 = [5 10  15  20  25];
    rho_permete =interp1(T1,rho1,T_bp);
    v_permete = (V_dot_permete*0.001*1/60)/(3*A_permete_Cross);
    "++++++++++ Renolds number and nusselt for feed +++++++++++";
    Re_Feed  = rho_Feed*v_Feed*D_h_Feed/mu_Feed; 
    Nusselt_Laminar_Feed = 1.86*(Re_Feed*Pr_Feed*D_h_Feed/L_Feed)^(0.33);  	                      
    Nusselt_Turblent_Feed= 0.023*(1+(6*D_h_Feed/L_Feed))*((Re_Feed)^0.8)*((Pr_Feed)^0.33);  %Dr. atia code
    %Nusselt_Turblent_Feed = 0.027*((Re_Feed)^(0.8))*(Pr_Feed)^(0.4)*((mu_Feed/mu_mean_Feed))^(0.14);                 
    if Re_Feed <2300
        Nusselt_Feed = Nusselt_Laminar_Feed ;
    else
        Nusselt_Feed =Nusselt_Turblent_Feed; 
    end
    h_Feed= (Nusselt_Feed*K_Feed)/D_h_Feed;%Convective heat transfer coefficient W/m2-K
    h_fg = (1.7535*(T_m+273)+2024.3)*1000; %Enthalpy or Latent Heat of vaporization of water, J/kg		                        
    
    "++++++++++ Renolds number and nusselt for permete +++++++++++";
    Re_permete  = (rho_permete*v_permete*D_h_permete)/mu_permete; 
    Nusselt_Laminar_permete = 1.86*(Re_permete*Pr_permete*D_h_permete/L_permete)^0.33;  	                      
    Nusselt_Turblent_permete = 0.023*(1+(6*D_h_permete/L_permete))*((Re_permete)^0.8)*((Pr_permete)^0.33);  %Dr. atia code
    %Nusselt_Turblent_permete = 0.027*(Re_permete)^(0.8)*(Pr_permete)^(0.33)*((mu_mean_permete/mu_permete))^(0.14);	                     
    if Re_permete<2300
        Nusselt_permete = Nusselt_Laminar_permete ;
    else
        Nusselt_permete =Nusselt_Turblent_permete; 
    end
    h_permete= (Nusselt_permete*K_permete)/D_h_permete;                         	                         
    "============================== Properties of Membrane =============================="; 
    alpha = 0.9;%a factor α is the ratio of Knudsen diffusion to mass diffusion(can vary between 0 and 1)
    delta_m = 153.9*(10^-6);%membrane thickness 		                       
    Epsilon = 0.8;%membrane porosity
    d_pore = 0.45*(10^-6);%mean pore size
    A_membrane = 0.013395;%membrane area		                       
    K_gas = (0.027 + 0.0184)/2;%gas conductivity
    K_mem = 0.25;%polymer membrane conductivity 
    Tau = 1/Epsilon;% membrane tortuosity 		                        
    K_m = ((Epsilon/K_gas)+((1-Epsilon)/K_mem))^(-1);%membrane conductivity
    "================================= Mass Transfer ================================="; 
    P_owf = exp((23.1964)-((3816.44)/((T_mf+273)-46.13)));%vapor pressures of feed sides at the membrane surfaces (pa)
    P_owp =exp((23.1964)-((3816.44)/((T_mp+273)-46.13)));%vapor pressures of permeate sides at the membrane surfaces(pa)
    P_f=101325; %(pa)
    P_p=101325; %(pa)
    
    P_wvp  = exp((23.1964)-((3816.44)/((T_m+273)-46.13))); %partial pressure of water vapors inside the pores
    P_pore = (P_f +P_p)/2;%The total pressure inside the pores and it is assumed to be the average of feed and permeate side pressures(pa)
    P_airpore = P_pore - P_wvp;	%the partial pressure of air inside the membrane pores 	                                      		                      
    Dwa  = 1.895*10^(-5)*(T_m+273)^(2.072);%the diffusivity of water vapors in air."Product of Pressue and Diffusion of water vapour in air [Pa.m^2/s]"                                                             
    D_k = (((3*delta_m*Tau)/(2*Epsilon*d_pore))*((Pi*R*(T_m+273))/(8*Mol_w))^(0.5))^(-1);% Knudsen diffusion coefficients	                   
    D_m  = ((R*(T_m+273)*delta_m*Tau*P_airpore)/(Epsilon*Mol_w*Dwa ))^(-1);%molecular diffusion coefficients
    D_e  = ((alpha/D_k )+((1-alpha)/D_m ))^(-1);%equivalent diffusio coefficient
    %D_e = ((2*Epsilon*d_pore)/(3*delta_m))*(((2*Mol_w)/(pi*R*(T_m+273)))^0.5)%            "M. Essalhi, M. Khayet, Dr atia code"
    J_w = D_e *((P_owf*Gamma_Feed*x_Feed)-P_owp);%is the mass flux of permeate
    Flux_MD=J_w*3600;% is the mass flux of permeate in one 
    Sc=mu_mean_Feed/(rho_Feed*D_e);%the Schmidt number
    Sh_Laminer_Feed=1.86*(Re_Feed*Sc*(D_h_Feed/L_Feed)); %the Sherwood number for Laminer flow
    Sh_Turblent_Feed=0.023*((Re_Feed)^0.8)*(Sc)^0.33; %the Sherwood number for Turblent flow
    if Re_Feed <2300
        Sh = Sh_Laminer_Feed;
    else
        Sh =Sh_Turblent_Feed; 
    end
    k_s=Sh*(D_e/D_h_Feed);% the solute mass transfer coefficient for the diffusive mass transfer through the concentration boundary layer in the feed side
    C_mf=C_Feed*exp(J_w/(k_s*rho_Feed));%concentration at membrane surface
    T_mf2 =(K_m*(T_bp+(h_Feed/h_permete)*T_bf)+((delta_m*(h_Feed*T_bf-J_w*h_fg))))/((K_m)+(h_Feed*(delta_m+(K_m/h_permete))));%The temperatures of the membrane surfaces at feed side 
    T_mp2 =(K_m*(T_bf+(h_permete/h_Feed)*T_bp)+((delta_m*(h_permete*T_bp+J_w*h_fg))))/((K_m)+(h_permete*(delta_m+(K_m/h_Feed))));%The temperatures of the membrane surfaces at permeate side
    res(1) = T_mf - T_mf2;
    res(2) = T_mp - T_mp2;
end