I'm trying to solve the following system of nonlinear equations with fsolve:
function F=equations(x)F(1)=h_hot*(Tav_hot-Thot_membrane)-((Thot_in-Thot_out)*Mav_hot*CP_w)/Am; F(2)=h_hot*(Tav_hot-Thot_membrane)-(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane));F(3)=(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane))-Kw/dhelta3*(Tcold_membrane-Thot_coolingplate);F(4)=Kw/dhelta3*(Tcold_membrane-Thot_coolingplate)-Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate);F(5)=Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate)-h_cold*(Tcold_coolingplate-Tav_cold);F(6)=h_cold*(Tcold_coolingplate-Tav_cold)-CP_w*Mcold_in*(Tcold_out-Tcold_in)/Am;x0=[310 300 315 312 311 295] %[310 300 315 312 311 295];
options=optimset('Display','iter')[x,fval,exit,output]=fsolve(@equations,x0,options)equations_guess=equations(x0)And it returns:
Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 7 1.81204e+12 1.19e+10 1 1 14 1.7734e+12 1 1.55e+10 1 2 21 1.72466e+12 2.5 1.67e+10 2.5 3 22 1.72466e+12 2.5 1.67e+10 2.5 4 29 1.68914e+12 0.625 3.03e+10 0.625 5 30 1.68914e+12 1.5625 3.03e+10 1.56 6 37 1.64038e+12 0.390625 1.27e+11 0.391 7 38 1.64038e+12 0.976562 1.27e+11 0.977 8 39 1.64038e+12 0.244141 1.27e+11 0.244 9 46 1.60886e+12 0.0610352 5.91e+11 0.061 10 47 1.60886e+12 0.152588 5.91e+11 0.153 11 48 1.60886e+12 0.038147 5.91e+11 0.0381 12 55 1.5969e+12 0.00953674 1.43e+12 0.00954 13 56 1.5969e+12 0.00953674 1.43e+12 0.00954 14 63 1.57561e+12 0.00238419 1.59e+13 0.00238 15 64 1.57561e+12 0.00596046 1.59e+13 0.00596 16 65 1.57561e+12 0.00149012 1.59e+13 0.00149 17 66 1.57561e+12 0.000372529 1.59e+13 0.000373 18 73 1.56808e+12 9.31323e-05 7.44e+13 9.31e-05 19 74 1.56808e+12 0.000232831 7.44e+13 0.000233 20 75 1.56808e+12 5.82077e-05 7.44e+13 5.82e-05 21 82 1.56316e+12 1.45519e-05 1.02e+15 1.46e-05 22 83 1.56316e+12 3.63798e-05 1.02e+15 3.64e-05 23 84 1.56316e+12 9.09495e-06 1.02e+15 9.09e-06 24 91 1.56124e+12 2.27374e-06 7.19e+14 2.27e-06 25 92 1.56124e+12 2.27374e-06 7.19e+14 2.27e-06 26 99 1.56106e+12 5.68434e-07 6.46e+14 5.68e-07 27 106 1.56073e+12 5.68434e-07 7.09e+14 5.68e-07 28 113 1.56057e+12 5.68434e-07 6.43e+14 5.68e-07 29 120 1.56023e+12 5.68434e-07 7.04e+14 5.68e-07 30 127 1.5601e+12 5.68434e-07 6.43e+14 5.68e-07 31 134 1.55975e+12 5.68434e-07 7.01e+14 5.68e-07 32 141 1.55965e+12 5.68434e-07 6.44e+14 5.68e-07 33 148 1.55927e+12 5.68434e-07 6.99e+14 5.68e-07 34 155 1.5592e+12 5.68434e-07 6.46e+14 5.68e-07 35 162 1.55881e+12 5.68434e-07 6.96e+14 5.68e-07 36 163 1.55881e+12 5.68434e-07 6.96e+14 5.68e-07 37 170 1.55868e+12 1.42109e-07 7.35e+14 1.42e-07 38 177 1.55862e+12 1.42109e-07 7.28e+14 1.42e-07 39 184 1.55855e+12 1.42109e-07 7.36e+14 1.42e-07 40 191 1.55849e+12 1.42109e-07 7.3e+14 1.42e-07 41 198 1.55841e+12 1.42109e-07 7.38e+14 1.42e-07 42 205 1.55836e+12 1.42109e-07 7.32e+14 1.42e-07 43 212 1.55828e+12 1.42109e-07 7.4e+14 1.42e-07 44 219 1.55822e+12 1.42109e-07 7.34e+14 1.42e-07 45 226 1.55814e+12 1.42109e-07 7.42e+14 1.42e-07 46 233 1.55808e+12 1.42109e-07 7.36e+14 1.42e-07 47 240 1.558e+12 1.42109e-07 7.44e+14 1.42e-07 48 247 1.55794e+12 1.42109e-07 7.38e+14 1.42e-07 49 254 1.55785e+12 1.42109e-07 7.46e+14 1.42e-07 50 261 1.55779e+12 1.42109e-07 7.4e+14 1.42e-07 51 268 1.5577e+12 1.42109e-07 7.48e+14 1.42e-07 52 275 1.55764e+12 1.42109e-07 7.42e+14 1.42e-07 53 282 1.55755e+12 1.42109e-07 7.5e+14 1.42e-07 54 289 1.55749e+12 1.42109e-07 7.44e+14 1.42e-07 55 296 1.5574e+12 1.42109e-07 7.51e+14 1.42e-07 56 303 1.55733e+12 1.42109e-07 7.46e+14 1.42e-07 57 310 1.55724e+12 1.42109e-07 7.53e+14 1.42e-07 58 317 1.55717e+12 1.42109e-07 7.48e+14 1.42e-07 59 324 1.55707e+12 1.42109e-07 7.55e+14 1.42e-07 60 331 1.55701e+12 1.42109e-07 7.5e+14 1.42e-07 61 338 1.5569e+12 1.42109e-07 7.57e+14 1.42e-07 62 345 1.55683e+12 1.42109e-07 7.52e+14 1.42e-07 63 352 1.55673e+12 1.42109e-07 7.58e+14 1.42e-07 64 359 1.55666e+12 1.42109e-07 7.53e+14 1.42e-07 65 366 1.55655e+12 1.42109e-07 7.59e+14 1.42e-07 66 373 1.55648e+12 1.42109e-07 7.55e+14 1.42e-07 67 380 1.55637e+12 1.42109e-07 7.6e+14 1.42e-07 68 387 1.5563e+12 1.42109e-07 7.56e+14 1.42e-07 69 394 1.55618e+12 1.42109e-07 7.6e+14 1.42e-07 70 401 1.55612e+12 1.42109e-07 7.56e+14 1.42e-07 71 408 1.556e+12 1.42109e-07 7.6e+14 1.42e-07 72 415 1.55594e+12 1.42109e-07 7.56e+14 1.42e-07 73 422 1.55583e+12 1.42109e-07 7.59e+14 1.42e-07 74 429 1.55577e+12 1.42109e-07 7.55e+14 1.42e-07 75 436 1.55566e+12 1.42109e-07 7.57e+14 1.42e-07 76 443 1.55561e+12 1.42109e-07 7.54e+14 1.42e-07 77 450 1.55552e+12 1.42109e-07 7.55e+14 1.42e-07 78 457 1.55548e+12 1.42109e-07 7.52e+14 1.42e-07 79 464 1.55541e+12 1.42109e-07 7.52e+14 1.42e-07 80 471 1.55538e+12 1.42109e-07 7.49e+14 1.42e-07 81 478 1.55532e+12 1.42109e-07 7.49e+14 1.42e-07 82 485 1.5553e+12 1.42109e-07 7.47e+14 1.42e-07 83 492 1.55526e+12 1.42109e-07 7.46e+14 1.42e-07 84 499 1.55525e+12 1.42109e-07 7.44e+14 1.42e-07 85 506 1.55522e+12 1.42109e-07 7.44e+14 1.42e-07 86 507 1.55522e+12 1.42109e-07 7.44e+14 1.42e-07 87 514 1.55469e+12 3.55271e-08 7.77e+14 3.55e-08 88 515 1.55469e+12 8.88178e-08 7.77e+14 8.88e-08 89 522 1.55435e+12 2.22045e-08 8.05e+14 2.22e-08 90 523 1.55435e+12 5.55112e-08 8.05e+14 5.55e-08 91 530 1.55421e+12 1.38778e-08 8.28e+14 1.39e-08 92 531 1.55421e+12 3.46945e-08 8.28e+14 3.47e-08 93 538 1.5542e+12 8.67362e-09 8.24e+14 8.67e-09 94 545 1.55416e+12 8.67362e-09 8.28e+14 8.67e-09 95 546 1.55416e+12 2.1684e-08 8.28e+14 2.17e-08 96 553 1.55415e+12 5.42101e-09 8.3e+14 5.42e-09 97 554 1.55415e+12 1.35525e-08 8.3e+14 1.36e-08 98 561 1.55413e+12 3.38813e-09 8.32e+14 3.39e-09 99 562 1.55413e+12 8.47033e-09 8.32e+14 8.47e-09 100 569 1.55413e+12 2.11758e-09 8.32e+14 2.12e-09 101 576 1.55412e+12 5.29396e-09 8.27e+14 5.29e-09 102 583 1.55409e+12 5.29396e-09 8.33e+14 5.29e-09 103 584 1.55409e+12 1.32349e-08 8.33e+14 1.32e-08 104 591 1.55408e+12 3.30872e-09 8.3e+14 3.31e-09 105 598 1.55407e+12 8.27181e-09 8.23e+14 8.27e-09 106 605 1.55403e+12 8.27181e-09 8.31e+14 8.27e-09 Solver stopped prematurely. fsolve stopped because it exceeded the function evaluation limit, options.MaxFunEvals = 600 (the default value). x = 1.0e+02 * 3.0994 - 0.0000i 3.0095 + 0.0000i 3.1319 + 0.0000i 3.1319 + 0.0000i 3.1148 - 0.0000i 2.9335 - 0.0000i fval = 1.0e+05 * -0.4949 + 0.1306i -0.4949 + 0.1306i -0.0036 - 0.0000i -0.5144 + 0.0000i -8.6309 - 0.0000i 8.9511 - 0.0000i exit = 0 output = iterations: 106 funcCount: 605 algorithm: 'trust-region-dogleg' firstorderopt: 8.3052e+14 message: 'Solver stopped prematurely. fsolve stopped because it exceeded the function evaluation limit, options.MaxFunEvals = 60...' equations_guess= 1.0e+05 * -3.2452 -3.4664 0.2218 -0.4550 -8.7591 9.0376Could anyone tell me what I have to do to get through the problem? Moreover I don't know why the solution x and fval become complex numbers. It's clear that increasing the number of iterations won't help!!
Thank you!!!!
The complete code is copied below, if anyone finds it necessary:
function solve()clc%Parameters
%modulo
%L=1; %lunghezza m
%W=0.5; %larghezza m
%spacer
%df= %filament diameter m
%deltha1=3; %thickness m
%es= %porosity
%THspacer= %angle porosity
%hot channel
deltha1=0.003; %thickness m
%dh= %equivalent diameter m
Apax=deltha1*0.03; %trasverse surface m^2
Deq=2*deltha1; %equivalent diameter m
%membrane
epsilon=0.8; %porosity
Am=0.0417; %surface m^2
dp=2.00E-07; %pore diameter m
dhelta2 =2.50E-04; %thickness mKm= 0.25; %conductivity W/mK
tau= 1.5; %tortuosity
K1= epsilon/tau;K0= (2/3)*K1*(dp/2);%airgap
dhelta3= 0.003; %thickness m%Kag= %conductivity water-air W/mK
%coolant plate
dhelta4= 7.00E-05; %thickness mKp=0.2; %conductivity W/mK%Antoine law log(10)P[mmHg]=A-B/(C+T[°C])
A=8.07131; B=1730.63; C=233.426;%other parameters
lambda=2200000; %latent heat of vapor j/kg
Patm= 101325; %Pa
CP_w= 4186; %Cp water j/kgK
CP_vap= 1900; %Cp vap j/kgK
MU_cold= 8.55E-04; %viscosity cold Pa*s
MU_hot= 3.65E-04; %viscosity hot Pa*s
PM= 0.018; % KG/mol water
R= 8.314; %universal constant Pa*m^3/(mol*k)
Kw= 0.63; %conductivity water W/mK
Kair= 0.03; %conductivity air W/mK
D12= 2.20E-05; %diffusivity water-air m^2/s
PR_hot= CP_w*MU_hot/Kw; %prandtl hot channel
PR_cold= CP_w*MU_cold/Kw; %prandtl cold channel
%boundary condition
Mhot_in= 0.02 ; %Feed mass in kg/s
Thot_in= 323.2 ; %Feed temperature in K
Mcold_in= 0.02 ; %Coolant mass in kg/s
Tcold_in= 291.1 ; %Coolant temperature in K
%v_hot=Mhot_in/(Apax*1000); %speed hot feed
%resolution
x0=[310 300 315 312 311 295]; %[310 300 315 312 311 295]
options=optimset('Display','iter')[x,fval,exit,output]=fsolve(@equations,x0,options)equations_guess=equations(x0)function F=equations(x)%J=x(1); %heat flux W/m^2
Thot_out= x(1); %Feed temperature out K
Tcold_out= x(2); %Coolant temperature out K
Thot_membrane= x(3); %temperature on the membrane hot side K
Tcold_membrane= x(4); %temperature on the membrane permeate side K
Thot_coolingplate= x(5); %temperature on cooling plate hot K
Tcold_coolingplate= x(6); %temperature on cooling plate cold K
Tav_hot= (Thot_in-Thot_out)/2; %temperature avarage in-out hot K
Tav_cold= (Tcold_out-Tcold_in)/2; %temperature avarage in-out cold K
Tav_m= (Thot_membrane-Tcold_membrane)/2; %temperature avarage membrane K
Psat_hot=133.3*10^(A-(B/(C+Thot_membrane-273))); %saturation pressure
Psat_cold=133.3*10^(A-(B/(C+Tcold_membrane-273))); %saturation pressure hmembrane=1/dhelta2*(epsilon*Kair+(1-epsilon)*Km); %membrane transfer heat coefficent
Pair_Hot_membrane= Patm-Psat_hot;Pair_Cold_membrane= Patm-Psat_cold;Pair_ml=(Pair_Hot_membrane-Pair_Cold_membrane)/log(Pair_Hot_membrane/Pair_Cold_membrane);Dknudsen=K0*(8*R*Tav_m/(pi*PM))^(1/2); %knudsen diffusivity
Dstagnant=D12*Patm/Pair_ml; %stagnant diffusivity
Koverall=1/(R*Tav_m)*(1/((1/Dknudsen+1/Dstagnant))); %membrane transfer mass coefficent
Nd=Koverall/dhelta2*(Psat_hot-Psat_cold); %distillate flux mol/m^2*s
Nd_l=Nd*PM; %distillate flux l/m^2*s o kg/m^2*s
Mhot_out= Mhot_in-Nd_l*Am; %mass hot out kg/s
Mav_hot= (Mhot_in-Mhot_out)/2; %mass avarage hot kg/s
Re_hot=Mav_hot*Deq/(Apax*MU_hot); %Reinolds hot channel
Re_cold=Mcold_in*Deq/(Apax*MU_cold); %Reinolds cold channel
h_hot=((8.4577*(Re_hot^0.1652)+0.7214*(Re_hot^0.5155))/2)*Kw/Deq; %transfer heat coeffcient hot W/m^2K
h_cold=((8.4577*(Re_cold^0.1652)+0.7214*(Re_cold^0.5155))/2)*Kw/Deq; %transfer heat coeffcient cold W/m^2K
F(1)=h_hot*(Tav_hot-Thot_membrane)-((Thot_in-Thot_out)*Mav_hot*CP_w)/Am; F(2)=h_hot*(Tav_hot-Thot_membrane)-(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane));F(3)=(lambda*Nd_l+hmembrane*(Thot_membrane-Tcold_membrane))-Kw/dhelta3*(Tcold_membrane-Thot_coolingplate);F(4)=Kw/dhelta3*(Tcold_membrane-Thot_coolingplate)-Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate);F(5)=Kp/dhelta4*(Thot_coolingplate-Tcold_coolingplate)-h_cold*(Tcold_coolingplate-Tav_cold);F(6)=h_cold*(Tcold_coolingplate-Tav_cold)-CP_w*Mcold_in*(Tcold_out-Tcold_in)/Am;endsave end
Best Answer