From fa52a961f7964a8cdcb237c88ffd4319eb99655c Mon Sep 17 00:00:00 2001 From: Oliver <oliver-mx@gmx.de> Date: Thu, 8 Aug 2024 13:01:52 +0200 Subject: [PATCH] Vectorized the SWRO/PRO-evaluation --- ERD_boundary1.m => Boundary1.m | 0 ERD_boundary2.m => Boundary2.m | 0 ERD_boundary3.m => Boundary3.m | 0 ERD_boundary4.m => Boundary4.m | 0 DATA_case2.mat | Bin 0 -> 834 bytes ERD_model.m => ODEsystem.m | 0 Sim_0a_data.m | 125 +++++ Sim_0_data.m => Sim_0b_data.m | 8 +- Skript.asv | 166 +++++++ Skript.m | 34 +- fun_1.asv | 836 +++++++++++++++++++++++++++++++++ fun_1.m | 313 ++++++------ 12 files changed, 1340 insertions(+), 142 deletions(-) rename ERD_boundary1.m => Boundary1.m (100%) rename ERD_boundary2.m => Boundary2.m (100%) rename ERD_boundary3.m => Boundary3.m (100%) rename ERD_boundary4.m => Boundary4.m (100%) create mode 100644 DATA_case2.mat rename ERD_model.m => ODEsystem.m (100%) create mode 100644 Sim_0a_data.m rename Sim_0_data.m => Sim_0b_data.m (97%) create mode 100644 Skript.asv create mode 100644 fun_1.asv diff --git a/ERD_boundary1.m b/Boundary1.m similarity index 100% rename from ERD_boundary1.m rename to Boundary1.m diff --git a/ERD_boundary2.m b/Boundary2.m similarity index 100% rename from ERD_boundary2.m rename to Boundary2.m diff --git a/ERD_boundary3.m b/Boundary3.m similarity index 100% rename from ERD_boundary3.m rename to Boundary3.m diff --git a/ERD_boundary4.m b/Boundary4.m similarity index 100% rename from ERD_boundary4.m rename to Boundary4.m diff --git a/DATA_case2.mat b/DATA_case2.mat new file mode 100644 index 0000000000000000000000000000000000000000..f14208f56b7a771b6bfd8753081021c22ffb5451 GIT binary patch literal 834 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQV4Jk_w+L}(NS<NN=+<DO;O0t zvr-7jC{=JQO;=E`P%t#GGB&U>Fj6owFfvgf60kBcKmensFFykVgE|8PL&coO$q5M= z3`u4R&l)Z%9N>Ct%(#(Fz?dmg&`qF@f%z*V!*(ggogn1~aOD;UP8>MEw!~6G!GvRj z(@!6V8m10O$p{Iho=40Me>i)17AH(#WVpGOxd)_A2d-}i!;u_^GmJ|e8$?f~GbZsW zq@^Szb(kbEOyc3;mSD4IVAwI0DF&p+2(E`CAuS<+`H~EefD}u_!WuCJJH{hD4q;0I zYZ!zltUYkXfXnTm#sOxAWs8_LgLL}Cb-n<5iXq97jU$02DNI4{<QZR6rq7WQHFN(q zOD=pB<kc<YZS?u`kA@t6mA`ue{;mnIlfI(kw!)%ApE2=Xs?((AL<Y{O9H&6W2I4Tb z3}S2-Kf@<d4Bp5X9K~4!GT4+676~2+sbHTvoH@zBxF*0(SwMP&Q-(l7(<ey@+gisO zJB4`%nl7Agakb}hJ`#DLs3*X5-?@JW4T^0hoDpnUmELk$Ct<6CoxZ^C)MrcD#8t&a zmWUla?)<e_;idZx_J^An?_^#Q*C@Dszx{%h=}ldXYZDUKlGp=wIDX_(*vD{0e}Tjj zJ#%Aq_02OBmu}J6sA6p2uc*7dg;6}=kY+POPsER&h$AyEoN*36-r>q4|G;^Xu7tqD zqimwbMbFO(C}dyvdVcTwHn}Bj|3rCqJf7akXk0h(?cr&9k`e)zRdv697AVm7nZp=q zuUz}Lx~ZY=OAXKIV1>OGzVP`!dveUY>;J?rn-ebfztLyZuvb2;Wqy6t4v8iEn;GkO ciZw2o&w0vS@4Af#4}Z*4DGmk(iGHJf0M~~n&Hw-a literal 0 HcmV?d00001 diff --git a/ERD_model.m b/ODEsystem.m similarity index 100% rename from ERD_model.m rename to ODEsystem.m diff --git a/Sim_0a_data.m b/Sim_0a_data.m new file mode 100644 index 0000000..bc0f079 --- /dev/null +++ b/Sim_0a_data.m @@ -0,0 +1,125 @@ +function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0a_data(input1, input2) +%% Sim_0a_data(input) +% +% Data for Simulation in chapter 3 +% +% SWRO-PRO hybrid system I (realistic ERD) +% +% option_data = 0 +% +% Input: +% input1 - x +% input2 - x + +%% model versions +version=zeros(1,10); +% version(1)=0 if co-current, 1 otherwise +version(2)=1; % 0 = SWRO beta fixed +version(3)=0; % 0 = PRO beta fixed +version(4)=1; % 0 = ideal SWRO +version(5)=1; % 0 = ideal PRO +% ERD +version(6)=3; +% version(6) = 0 --> only SWRO (no ERD) +% version(6) = 1 --> only SWRO (with ERD) +% version(6) = 2 --> only PRO +% version(6) = 3 --> SWRO-PRO hybrid system (with one ERD) +% version(6) = 4 --> SWRO-PRO hybrid system (with two ERDs) +version(7)=1; % 1 = ICP and ECP for SWRO +version(8)=1; % 1 = ICP and ECP for PRO + +%% Membrane unit properties +H = 1e-3; % height of the membrane [m] +ro_water = 1000; % mass density of water [kg/m^3] +ro_salt = 2165; % mass density of salt [kg/m^3] +Mw = 18; % molecular weight of water [kg/kmol] +Ms = 58.44; % molecular weight of salt [kg/kmol] +Rw = 462; % gas constant of water [J/(kg K)] +T0 = 297; % temperature [K] +eta = 1.3e-3; % seawater viscosity [kg/(m s)] +p_r = 1e5; % pressure [Pa]=[kg/ms^2] +rho_r = 1e3; % density [kg/m^3] +C_r = 1; % salt concentration [%] + +%% SWRO +swro_Z=1; % width of SWRO membrane [m] +swro_L=4; % length of SWRO membrane [m] +if version(6)== 2; swro_L=1; end +swro_alpha = 5.0815e-9; % SWRO water permeablity co-efficient [s/m] +swro_R =0.96; % SWRO salt rejection rate +swro_KK = 1e-2; % SWRO ICP mass transfer coefficient +swro_KD = 1/swro_KK; % SWRO ECP draw side mass transfer coefficient +swro_KF = 1/swro_KK; % SWRO ECP fresh side mass transfer coefficient +swro_x_r=swro_L^2; % x_r=swro_L^2 since x = linspace(0,1,n) (if x = linspace(0,swro_L,n) then x_r=swro_L;) +swro_b1 = H/swro_x_r; % H/swro_L ratio +swro_b2 = swro_Z/swro_x_r; % Z/swro_L ratio +J_r = sqrt(H^3/swro_x_r*p_r*rho_r); % flux [kg/s^2] +swro_gamma= swro_x_r * p_r * swro_alpha /J_r; % SWRO scaling factor - mass balance +swro_gamma2= J_r^2./(swro_x_r^2 * p_r * rho_r); % SWRO scaling factor - momentum balance +swro_W_r=J_r*p_r/rho_r; % net work [W/m^2] +sigma=0.999 ; % Rejection coefficient +swro_beta_fix=4.43e-4/J_r*swro_x_r; % value for fixed SWRO beta [kg/sm^2] + +%% PRO +Z=1; % width of the PRO membrane [m] +L=1.5; % length of the PRO membrane [m] +alpha = 5.47e-9; % water permeablity co-efficient [s/m] +R =0.94; % salt rejection rate [1] +KK=7.13e2; % mass transfer coefficient [sm^2/kg] +KD=1/KK; % PRO ECP draw side mass transfer coefficient [sm^2/kg] +KF=1/KK; % PRO ECP fresh side mass transfer coefficient [sm^2/kg] +x_r=L^2; % x_r=L^2 since x = linspace(0,1,n) (if x = linspace(0,L,n) then x_r=L;) +b1 = H/x_r; % H/L ratio +b2 = Z/x_r; % Z/L ratio +Q_r = sqrt(H^3/x_r*p_r*rho_r); % flux [kg/s^2] +gamma= x_r * p_r * alpha /Q_r ; % PRO scaling factor - mass balance +gamma2= Q_r^2./(x_r^2 * p_r * rho_r); % PRO scaling factor - momentum balance +W_r=Q_r*p_r/rho_r; % net work [W/m^2] +beta_fix =1.71e-4/Q_r*x_r; % value for fixed PRO beta [kg/sm^2] + +%% Sea Water +cE= 35/983/C_r; % salt concentration in seawater +pE= 1e5/p_r; % external pressure +rho_E=(cE + 1)./(ro_water*cE/ro_salt + 1); % density of incomming seawater + +%% SWRO operating conditions +J_sf_0 = 0; % salt flux in fresh side at 0 +J_wf_0 = 0; % water flux in fresh side at L +Pf_L = pE; % pressure fresh side at L + +Pd_0 = 66e5/p_r; % pressure draw side at 0 +Pd_L = 60e5/p_r; % pressure draw side at L (not need in the hybrid system) + +%% PRO operation conditions +Q_sf_0 = 0; % salt flux in fresh side at 0 +pf_L = pE; % pressure of fresh side at L + +pd_0 = 12e5/p_r; % pressure draw side at 0 +pd_L = 14.6e5/p_r; % pressure of fresh side at 0 +pf_0 = 1.3e5/p_r; % pressure draw side at L + +%% ERD/Turbine/Pump parameters +T_eff = .95; % turbine efficiency +HP_eff = .9; % high pressure pump efficiency +LP_eff = .95; % low pressure pump efficiency +V_m=0.052; % Volumetric mixing +ERD_eff=.96; % ERD unit pressure efficiency +ERD_fric=5e-04; % ERD friction coefficient +A_ERD=H*swro_Z; % cross sectional area of ERD inflows/outflows +eta_ERD=0.01; % leak of high pressure brine +mix_density=997/rho_r; % density of mixture in ERD +pw=2; % water price [$/m^3] +pe=0.5; % electricity price [$/kWh] +mixer_ERD=1; % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed) + +%% display figures +fig=[1,0,0,0,0,0,0,0,1,0,0]; % f(i)=1 --> figure i will be displayed +%fig=[0,0,0,0,0,0,0,0,0,0,0]; % f(i)=1 --> figure i will be displayed + +%% model specific changes: +% co-current +if pd_0 > pd_L; version(1)=0; else; version(1)=1; end +% no SWRO needed (--> use trivial data) +if version(6) == 2; version(4)=0; Pd_0=40; Pd_L=35; end + +end \ No newline at end of file diff --git a/Sim_0_data.m b/Sim_0b_data.m similarity index 97% rename from Sim_0_data.m rename to Sim_0b_data.m index 4e8caeb..2cd106e 100644 --- a/Sim_0_data.m +++ b/Sim_0b_data.m @@ -1,11 +1,11 @@ -function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0_data(input1, input2) -%% Sim_32_data(input) +function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0b_data(input1, input2) +%% Sim_0b_data(input) % % Data for Simulation in chapter 3 % % SWRO-PRO hybrid system I (realistic ERD) % -% option_data = 0 +% option_data = 0.1 % % Input: % input1 - x @@ -104,7 +104,7 @@ HP_eff = .9; % high pressure pump efficiency LP_eff = .95; % low pressure pump efficiency V_m=0.052; % Volumetric mixing ERD_eff=.96; % ERD unit pressure efficiency -ERD_fric=5e-06; % ERD friction coefficient +ERD_fric=5e-04; % ERD friction coefficient A_ERD=H*swro_Z; % cross sectional area of ERD inflows/outflows eta_ERD=0.01; % leak of high pressure brine mix_density=997/rho_r; % density of mixture in ERD diff --git a/Skript.asv b/Skript.asv new file mode 100644 index 0000000..e92418f --- /dev/null +++ b/Skript.asv @@ -0,0 +1,166 @@ +%% Matlab skript for the Paper 2024: +% +% - reproduce the simulations +% - reproduce optimisation results +% - reproduce all figures.png (converted to .eps using LibreOfficeDraw) +% - test average computation time +% + +%% Numerical Simulations +% +% Create figures form chapter: Numerical simulations +% +close all; +[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4), +% +%figure(1); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_6.png'); +%figure(2); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_7.png'); +%figure(3); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_8_ERD.png'); +% + +%% Average simulation time: +% +% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data" +% +T=zeros(1,100); +for i=1:length(T) +[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4); +T(i)=time; +end +mean(T) %= 0.9600 +% + +%% Optimisation - case2 (with epsilon constraints) +% +% Create Paretofront for SWRO+ERD +% +load("DATA_case2.mat") +% +% load("DATA_case2.mat") +% save DATA_case2.mat X_2 +% +% X2 - set of operating conditions for epsilon constraint method (2x40) +% Y2 - system output using X2 (4x40) +% +% X2_P - set of operating conditions from paretosearch (2x200) +% Y2_P - system output using X2_P (4x200) +% +% time2 - time for epsilon constraint method (40points) +% time2_P - time for paretosearch solver (200 points) +% + +X2=zeros(2,40);Y2=zeros(4,40);X2_P=zeros(2,200); + + + +load("DATA_T1.mat") +c1=datetime("now"); +x0=[69.9 51.5]; +FWmin=linspace(0.3770, 0.0015, 40); +for i=1:40 +clc, +i, +close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70]; +foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100); +option_mesh = 1e4; option_BVP = 1e-4; option_data = 1; +[x_opt, f_opt]=fmincon(@(x)fun_1(x,1,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',1,option_mesh,option_BVP, FWmin(i)),foptions); +end +c2=datetime("now"); %von 2 bis 39 +save DATA_T1.mat c1 c2 % 00:02:30 + +%% Optimisation - case2 (with paretosearch) +% +% Create Paretofront for SWRO+ERD +% + + + +%% Optimisation-2 (red) +close all;clear all;clc +c3=datetime("now"); +x0=[69.996767760342830 49.753887785404920]; +FWmin=linspace(0.3862, 0.0066, 40); +for i=2:39 +clc, +i, +close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70]; +foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100); +option_mesh = 1e4; option_BVP = 1e-4; option_data = 2; +[x_opt, f_opt]=fmincon(@(x)fun_1(x,2,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',2,option_mesh,option_BVP, FWmin(i)),foptions); +end +c4=datetime("now"); %von 2 bis 39 +save DATA_T2.mat c3 c4 +(c4-c3)/38, %00:01:54 + +%% hybrid paretosearch - old +close all; clear all; clc +rng default % For reproducibility +A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5]; +% initial points: +load('DATA_Approx.mat'); +X0 = X_approx'; %size(X0)= (93,5) +%remove infeasible poitns +z=length(X_approx(1,:)); +for i=1:z + k=0; + for j=1:4 + if X_approx(j,z+1-i) > ub(j); k=1;end + if X_approx(j,z+1-i) < lb(j); k=1;end + end + if k==1;X0(z+1-i,:)=[]; end +end +%After removing points: +%size(X0)= (61,5) +% +options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000); +% x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options) +[X_paretosearch_3, Y_paretosearch_3] = paretosearch(@(x)fun_1(x,3,'Pareto',1e4,1e-4),5,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); + +%% +load DATA_Pareto.mat +scatter(SEC_Pareto_1,FW_Pareto_1);hold on +scatter(SEC_Pareto_2,FW_Pareto_2) + +%% red - Paretosearch +% +close all; clear all; clc +rng default % For reproducibility +% +A= [-1 1]; b= 0; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70]; +%A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5]; +% +% initial points: +load DATA_Pareto.mat +X0 = [X_Pareto_2';X_Pareto_2';X_Pareto_2';X_Pareto_2';X_Pareto_2']; +% +options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000); +c5=datetime("now"); +% x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options) +[X, Y] = paretosearch(@(x)fun_1(x,2,'Pareto',1e4,1e-4),2,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); +c6=datetime("now"); +scatter(-Y(:,1),-Y(:,2));hold on +scatter(SEC_Pareto_2,FW_Pareto_2) +save DATA_T3 c5 c6 +%c6-c5 = 00:56:57 + +%% blue - Paretosearch +% +close all; clear all; clc +rng default % For reproducibility +% +A= [-1 1]; b= 0; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70]; +%A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5]; +% +% initial points: +load DATA_Pareto.mat +X0 = [X_Pareto_1';X_Pareto_1';X_Pareto_1';X_Pareto_1';X_Pareto_1']; +% +options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000); +c7=datetime("now"); +% x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options) +[X, Y] = paretosearch(@(x)fun_1(x,1,'Pareto',1e4,1e-4),2,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); +c8=datetime("now"); +scatter(-Y(:,1),-Y(:,2));hold on +scatter(SEC_Pareto_1,FW_Pareto_1) +save DATA_T4 c7 c8 +%c8-c7 = 01:05:44 diff --git a/Skript.m b/Skript.m index dd69501..2a11bef 100644 --- a/Skript.m +++ b/Skript.m @@ -20,21 +20,40 @@ close all; %% Average simulation time: % -% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data" +% calculate average computation time of previous simulation % T=zeros(1,100); for i=1:length(T) -[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4); +[Output1, Output2, time]=fun_1(1,0.1,'sol',1e4,1e-4); T(i)=time; end mean(T) %= 0.9600 % -%% Optimisation-1(blue) +%% Optimisation - case2 (with epsilon constraints) +% +% Create Paretofront for SWRO+ERD +% +load("DATA_case2.mat") +% +% load("DATA_case2.mat") +% save DATA_case2.mat X2 Y2 X2_P Y2_P time2 time2_P +% +% X2 - set of operating conditions for epsilon constraint method (2x40) +% Y2 - system output using X2 (4x40) +% +% X2_P - set of operating conditions from paretosearch (2x200) +% Y2_P - system output using X2_P (4x200) +% +% time2 - time for epsilon constraint method (40points) +% time2_P - time for paretosearch solver (200 points) +% + +load("DATA_T1.mat") c1=datetime("now"); x0=[69.9 51.5]; FWmin=linspace(0.3770, 0.0015, 40); -for i=2:39 +for i=1:40 clc, i, close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70]; @@ -45,6 +64,13 @@ end c2=datetime("now"); %von 2 bis 39 save DATA_T1.mat c1 c2 % 00:02:30 +%% Optimisation - case2 (with paretosearch) +% +% Create Paretofront for SWRO+ERD +% + + + %% Optimisation-2 (red) close all;clear all;clc c3=datetime("now"); diff --git a/fun_1.asv b/fun_1.asv new file mode 100644 index 0000000..15e4f1a --- /dev/null +++ b/fun_1.asv @@ -0,0 +1,836 @@ +function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,option_BVP) +%% fun Solves BVP for SWRO-PRO hybrid system using bvp5c +% +% fun(input,option_data,obj,option_mesh,option_BVP) will solve the +% nonlinear BVP of the selected SWRO-PRO system. In "option_data" +% all Informations of the system are given with exception of the +% additional "input". +% +% Input: +% input - Input for the Data function +% option_data - Selects Data function used: data(input) +% obj - 'SEC', 'FW', 'Rev', 'PD_net' , 'SE_net', 'approx' or 'sol' ; +% option_mesh - NMax of bvp5c +% option_BVP - RelTol of bvp5c +% +% Output: +% output1 = [SEC_net, FW, Rev, SWRO_Rec, PRO_Rec, PD_net, SE_net]; +% output2 = [W_net, -W_p1, -W_p2, -W_p3, -W_p4, W_t]; +% figures may be displayed if obj='sol' +% +% Example: +% x0=1; +% D=.11; +% obj='sol'; +% mesh=2e4; +% tol=1e-6; +% [Sim_01_output1, Sim_01_output2]=fun(x0,D,obj,mesh,tol), +% +tic; +%% Read data +if option_data == 0; DATA = @(x)Sim_0a_data(input); end +if option_data == 0.1; DATA = @(x)Sim_0b_data(input); end + + + +% Simulations +if option_data == 0.11; DATA = @(x)Sim_11_data(input); end +if option_data == 0.12; DATA = @(x)Sim_12_data(input); end +if option_data == 0.13; DATA = @(x)Sim_13_data(input); end +if option_data == 0.21; DATA = @(x)Sim_21_data(input); end +if option_data == 0.22; DATA = @(x)Sim_22_data(input); end +if option_data == 0.23; DATA = @(x)Sim_23_data(input); end +if option_data == 0.24; DATA = @(x)Sim_24_data(input); end +if option_data == 0.31; DATA = @(x)Sim_31_data(input); end +if option_data == 0.32; DATA = @(x)Sim_32_data(input); end +% PRO Membrane Comparison +if option_data == 0.51; DATA = @(x)Sim_PRO_01_data(input); end +if option_data == 0.52; DATA = @(x)Sim_PRO_02_data(input); end +if option_data == 0.53; DATA = @(x)Sim_PRO_03_data(input); end +% Rough vs. Detailled +if option_data == 0.91; DATA = @(x)Sim_rough_approx_data(input); end +% Optimizations +if option_data==1; DATA = @(x)Opt_1_data(input); end +if option_data==2; DATA = @(x)Opt_2_data(input); end +if option_data==3; DATA = @(x)Opt_3_data(input); end +if option_data==4; DATA = @(x)Opt_4_data(input); end +if option_data==5; DATA = @(x)Opt_L_data(input); end +% for 3d plot +if option_data==11; DATA = @(x)plot3d_data(input); end + +[H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD, KF, KD] ... + = DATA(1); + +%% Set options +ode_options = bvpset('Stats','off','NMax', option_mesh, 'RelTol', option_BVP); + +try +n=25; +x = linspace(0,1,n); + +%% Create initial guess for the BVP +J_wd_0 = (983/1018)/J_r; % guess for J_wd(0) +J_wf_0 = 0.01353/J_r; % should be close to 0 +Pf_0 = 1.1*Pf_L; % guess for P_f(0) +Q_wd_0 = 0.01353*(983/1018)/Q_r; % guess for Q_wd(0) +Q_wf_0 = 0.01353/Q_r; % guess for Q_wf(0) + +if version(1)==0; y_init = [cE; J_wd_0; J_sf_0; J_wf_0; Pd_0; Pf_0; cE; Q_wd_0; Q_sf_0; Q_wf_0; pd_0; pf_0]; end +if version(1)==1; y_init = [cE; J_wd_0; J_sf_0; J_wf_0; Pd_0; Pf_0; cE; -Q_wd_0; Q_sf_0; Q_wf_0; pd_0; pf_0]; end +solinit = bvpinit(x,y_init); + +%% Solve the BVP +if version(6) == 0; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end +if version(6) == 1; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary1(ya,yb, DATA),solinit,ode_options); end +if version(6) == 2; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end +if version(6) == 3; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary3(ya,yb, DATA),solinit,ode_options); end +if version(6) == 4; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary4(ya,yb, DATA),solinit,ode_options); end + +%% Evaluate the solution +n=100; +x=linspace(0,1,n); +y = deval(sol,x); Y = y'; Y = real(Y); + +%% SWRO evaluation +C_d = Y(:,1); +C_f = Y(:,3)./Y(:,4); +J_sd = Y(:,1).*Y(:,2); +J_wd = Y(:,2); +J_d = Y(:,1).*Y(:,2)+Y(:,2); +J_sf = Y(:,3); +J_wf = Y(:,4); +J_f = Y(:,3)+Y(:,4); +P_d = Y(:,5); +P_f = Y(:,6); +% same quantities as in "ODEsystem.m" +swro_p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1)/Ms)/p_r; +swro_p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3)/Ms./Y(:,4))/p_r; +swro_del_c = Y(:,1)./(Y(:,1) + 1) - Y(:,3)./(Y(:,3) + Y(:,4)); +swro_local_ro_d = (Y(:,1) + 1)./(ro_water*Y(:,1)/ro_salt + 1); +swro_local_ro_f = (Y(:,3) + Y(:,4))./(ro_water*Y(:,3)/ro_salt + Y(:,4)); +swro_beta = (1 - swro_R)*((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))/swro_R; +if version(2) == 0 ;swro_beta= swro_beta_fix.*ones(1,n); end +if version(4) == 0 ;swro_beta=zeros(1,n); end +% J_w,in and J_s,in +J_cross = ((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f)); +if version(4) == 0; J_cross = ((Y(:,5) - Y(:,6)) - sigma*(swro_p_osm_d - swro_p_osm_f)); end +J_sin=swro_beta.*swro_del_c; +if version(7)==1 +J_cross = ((Y(:,5)-Y(:,6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta.*(Y(:,5)-Y(:,6)).*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha.*swro_beta.*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d./swro_KD -swro_p_osm_f.*(swro_KK + 1/swro_KF) )); +J_sin = swro_beta.*( swro_del_c-J_cross.*( Y(:,1)./(Y(:,1) + 1)./swro_KD + Y(:,3)./(Y(:,3) + Y(:,4)).*(swro_KK+1/swro_KF)))./(1+swro_alpha.*swro_beta.*(swro_KK+1/swro_KF+1/swro_KD)); +end + +%% PRO evaluation +c_d = Y(:,1+6); % concentration of draw side +c_f = Y(:,3+6)./Y(:,4+6); % concentration of fresh side +Q_sd = Y(:,1+6).*Y(:,2+6); % salt flow rate in draw side +Q_wd = Y(:,2+6); % water flow rate in draw side +Q_d = Y(:,1+6).*Y(:,2+6)+Y(:,2+6); % total flow rate draw side +Q_sf = Y(:,3+6); % salt flow rate in fresh side +Q_wf = Y(:,4+6); % water flow rate in fresh side +Q_f = Y(:,3)+Y(:,4+6); % total flow rate fresh side +p_d = Y(:,5+6); % Pressure draw side +p_f = Y(:,6+6); % Pressure fresh side +% same quantities as in "ODEsystem.m" +p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1+6)/Ms)/p_r; +p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3+6)/Ms./Y(:,4+6))/p_r; +del_c = Y(:,1+6)./(Y(:,1+6) + 1) - Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)); +local_ro_d = (Y(:,1+6) + 1)./(ro_water*Y(:,1+6)/ro_salt + 1); +local_ro_f = (Y(:,3+6) + Y(:,4+6))./(ro_water*Y(:,3+6)/ro_salt + Y(:,4+6)); +beta = (1 - R)*((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)))./R; +if version(3) == 0 ; beta = beta_fix.*ones(1,n); end +if version(5) == 0 ; beta = zeros(1,n); end +% Q_w,in and Q_s,in +Q_cross = (p_osm_d - p_osm_f - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*KK.*alpha.*p_r)./((ones(1,n) + KK*alpha*p_r.*(beta +p_osm_f'))'); +if version(5) == 0; Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6))) ; end +Q_sin = beta'.*del_c; +if version(8)==1 +Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*(KK+1./KF+1./KD).*p_r.*alpha)./(ones(1,n) + alpha.*beta.*(KK+1./KF+1./KD) + alpha.*p_r.*(p_osm_f'.*(KK+1./KF) + p_osm_d'./KD))'; +Q_sin = beta'.*( del_c-Q_cross.*( Y(:,1+6)./(Y(:,1+6) + 1)./KD + Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)).*(KK+1/KF)))./(1+alpha.*beta'); +end + +%% Freswater production and recovery rates +Freshwater = J_wf(end)*J_r; %in [kg/s] +FW = (Freshwater*swro_Z/(swro_local_ro_f(end)*rho_r))*3600; %in [m^3/h] +SWRO_Recovery =(J_f(end)./J_d(1))*100; % SWRO freshwater recovery [%] +PRO_Recovery = (1- Q_f(end)./Q_f(1))*100; % PRO recovery rate [%] +% warning flags: + if FW < 0 + fprintf(2,' \nERROR: Waring: Freshwater production is negative! \n'); + end + if SWRO_Recovery > 100 + fprintf(2,' \nERROR: Waring: SWRO recovery is greater then 100 %% \n'); + end + if SWRO_Recovery < 0 + fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n'); + end + if PRO_Recovery > 100 + fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n'); + end + if PRO_Recovery < 0 + fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n'); + end + +%% Output of System +J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; J_wE = J_E/(cE+1); +W_p1 = 1/HP_eff * (P_d(1)-pE)*(J_E *swro_Z)/rho_E; W_p1 = W_p1*swro_W_r; +W_p2 = 0; W_p2 = W_p2*swro_W_r; +W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r; + + %% Version(6)=0 --> only SWRO (no ERD) + if version(6)==0; W_p3=W_p1; W_p4=0; W_t=0; end + + %% Version(6)=1 --> only SWRO (with ERD) + if version(6)==1 + rho_ERD = V_m*(swro_local_ro_d(end) - rho_E)+rho_E; + C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt); + J_wERD = J_E/(C_ERD+1); + Q_exit = J_d(end)*(1-eta_ERD); + c_exit = (C_d(end)*J_wd(end)*(1 - eta_ERD)+cE*J_wE-C_ERD*J_wERD)/(J_wd(end)*(1 - eta_ERD)+J_wE-J_wERD); + rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE; + f_1= (Q_exit*swro_Z/rho_exit*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2 ); + f_2= (J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0; + if W_p3 < 0 + fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n'); + end + + end + %% Version(6)=2 --> only PRO + if version(6)==2 + if version(1) ==0; W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); end + if version(1) ==1; W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1); end + W_p1=0; W_p3=0; W_t=W_t*W_r; + if version(1) ==0; W_p2= 1/LP_eff * (p_d(1)-pE)*(Q_d(1)*Z)./local_ro_d(1); end + if version(1) ==1; W_p2= 1/LP_eff * (p_d(end)-pE)*(abs(Q_d(end))*Z)./local_ro_d(end); end + W_p2 = W_p2*W_r; + end + + %% Version(6)=3 --> SWRO-PRO hybrid system (with one ERD) + if version(6)==3 + rho_ERD = V_m*(swro_local_ro_d(end) - rho_E)+rho_E; + C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt); + qq=Q_r/J_r; + if version(1) == 0 %co-current + f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2; + f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2; + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; + W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); + end + if version (1) == 1 %counter-current + f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2 ; + f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2 ; + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; + W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1); + end + W_p3=W_p3*swro_W_r; W_t=W_t*W_r; + end + + %% Version(6)=4 --> SWRO-PRO hybrid system (with two ERDs) + if version(6)==4 + if version(1) == 0; rho_ERD2= V_m*(local_ro_d(end) - rho_E)+rho_E; end + if version(1) == 1; rho_ERD2= V_m*(local_ro_d(1) - rho_E)+rho_E; end + C_ERD2= -ro_salt*(rho_ERD2-1)/(rho_ERD2*ro_water-ro_salt); + rho_ERD= V_m*(swro_local_ro_d(end) - rho_ERD2)+rho_ERD2; + C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt); + J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; %J_wE = J_E/(C_ERD2+1); J_wERD= J_ERD/(C_ERD+1); + qq=Q_r/J_r; + J_E_sea = 2*J_E; + J_wE_sea = J_E_sea/(cE+1); + if version(1) == 0 + Q_exit = mixer_ERD*Q_d(end)*(1-eta_ERD); + c_exit = (c_d(end)*mixer_ERD*Q_wd(end)*(1 - eta_ERD)+cE*J_wE_sea-C_ERD2*2*J_E)/(mixer_ERD*Q_wd(end)*(1 - eta_ERD)+J_wE_sea-2*J_E); + end + if version(1) == 1 + Q_exit = mixer_ERD*abs(Q_d(1))*(1-eta_ERD); + c_exit = (c_d(1)*mixer_ERD*abs(Q_wd(1))*(1 - eta_ERD)+cE*J_wE_sea-C_ERD2*2*J_E)/(mixer_ERD*abs(Q_wd(1))*(1 - eta_ERD)+J_wE_sea-2*J_E); + end + + +%% Removed comments + +%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here?? + +% -3.9191 0.3247 0.0131 32.8727 7.1201 NaN NaN +% +% 1.0e+03 * +% -1.2742 -0.9744 0 -0.4873 -0.0010 0.1885 + +% In ideal cases take mean of some quantities. +% In ideal cases salt will fluctuate around 1e-16. We take mean here for +% +% This is just for figures. (), otherwise one could use ylim ) +%if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end +%if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end + + %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + %display(['pERD: ',num2str(pERD)]) % pERD use in the model + %display(['pERD_theoretical: ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency + + %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + + %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + %display(['pERD: ',num2str(pERD)]) % pERD use in the model + %display(['pERD_theoretical: ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + + %J_wERD = J_E/(C_ERD+1); + + + + + + + + + + + + rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE; + %calculation of Pressure P_ERD2 + if version(1) == 0 %co-current + f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2; + f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end)/A_ERD)^2; + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD2 = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; + f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; + if pERD2 > pMax; pERD2=pMax; end + %if pERD2 < 0; pERD2=pMax; end + + end + if version (1) == 1 %counter-current + f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2; + f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1)/A_ERD)^2; + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD2 = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; + f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; + if pERD2 > pMax; pERD2=pMax; end + %if pERD2 < 0; pERD2=pMax; end + end + %calculation of Pressure P_ERD + if version(1) == 0 %co-current + f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2; + f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2; + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; + f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; + if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; + W_t = T_eff * (p_d(end)-pE)*((1-mixer_ERD)*Q_d(end)*Z)/local_ro_d(end); + end + if version (1) == 1 %counter-current + f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2 ; + f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2 ; + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; + f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; + if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; + W_t = T_eff * (p_d(1)-pE)*(abs((1-mixer_ERD)*Q_d(1))*Z)/local_ro_d(1); + end + % Recalculation for Pump 1 + W_p1 = 1/HP_eff * (P_d(1)-pERD2)*(J_E *swro_Z)/rho_ERD2; W_p1 = W_p1*swro_W_r; + end + +%% final output: +W_net= - W_p1 - W_p2 - W_p3 - W_p4 + W_t; % in [W] +PD_net = W_net./(Z*L); %in [W/m^2] +%old SE_net +%if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + Q_d(1)*Q_r*Z/(rho_r*local_ro_d(1)))/1000/3600; end %in [kWh/m^3] +%if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + abs(Q_d(end))*Q_r*Z/(rho_r*local_ro_d(end)))/1000/3600; end %in [kWh/m^3] +%new SE_tot +if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3] +if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3] +SEC_net = W_net*(swro_local_ro_f(end)*rho_r)./(J_f(end)*J_r*swro_Z)/1000/3600; %in [kWh/m^3] +Rev= pw*FW + pe*SEC_net*FW; %in [$/h] +tt=toc; +%% if BVP-solver failed +catch + sol.stats.maxerr =1000; + SEC_net=1000;FW=0;Rev=-20;SWRO_Recovery=NaN;PRO_Recovery=NaN;PD_net=NaN;SE_net=NaN; + W_net=NaN;W_p1=NaN;W_p2=NaN;W_p3=NaN;W_p4=NaN;W_t=NaN; +end +%% Output of objective function +switch (obj) + case 'SEC' + if sol.stats.maxerr > option_BVP + output1 = 20; + else; output1 = -SEC_net; + end + case 'FW' + if sol.stats.maxerr > option_BVP + output1 = 20; + else; output1 = -FW; + end + case 'Rev' + if sol.stats.maxerr > option_BVP + output1 = 0; + else; output1 = -Rev; + end + case 'Pareto' + if sol.stats.maxerr > option_BVP + output1 = [NaN, NaN]; + else; output1 = [-SEC_net, -FW]; + %if SWRO_Recovery < 30;output1 = [NaN, NaN];end + %if SWRO_Recovery > 60;output1 = [NaN, NaN];end + %if PRO_Recovery < 40;output1 = [NaN, NaN];end + %if PRO_Recovery > 95;output1 = [NaN, NaN];end + end + case 'PD_net' + if sol.stats.maxerr > option_BVP + output1 = 0; + else; output1 = -PD_net; + end + case 'SE_net' + if sol.stats.maxerr > option_BVP + output1 = 0; + else; output1 = -SE_net; + end + case 'approx' + if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end + if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end + if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end + if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end + if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end + %sol.stats.nmeshpoints + %sol.stats.maxerr + %sol.stats.nODEevals + %sol.stats.nBCevals + output2=[sol.stats.nmeshpoints sol.stats.maxerr]; + case 'exe' + if sol.stats.maxerr > option_BVP + output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN]; + output2=[NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN]; + else + if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end + if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end + if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end + if version(6) == 3; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end + if version(6) == 4; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end + output2=[x*swro_L; 100*C_d; J_r*J_d; p_r*P_d; p_r*P_f; x*L; 100*c_d; Q_r*Q_d; Q_r*Q_f; p_r*p_d; p_r*p_f]; % 11 + output3=[x(2:end)*swro_L; 100*C_f(2:end); J_r*J_f(2:end); J_r/swro_x_r*J_cross(2:end); J_r/swro_x_r*J_sin(2:end); x(2:end)*L; 100*c_f(2:end); Q_r/x_r*Q_cross(2:end); Q_r/x_r*Q_sin(2:end)]; %9 + % + output1=output1';output2=output2';output3=output3'; + [a,b]=size(output1); + output1=[1:1:b; output1]; + [a,b]=size(output2); + output2=[1:1:b; output2]; + [a,b]=size(output3); + output3=[1:1:b; output3]; + end + case 'sol' + if sol.stats.maxerr > option_BVP + output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN]; + output2=[NaN, NaN, NaN, NaN, NaN, NaN]; + output3=NaN; + else + if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end + if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end + if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end + if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end + if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end + output2=[W_net, -W_p1, -W_p2, -W_p3, -W_p4, W_t]; + output3=tt; + end + + +%% figure 1 +if sol.stats.maxerr ==1000 +fprintf(2,' \nERROR: bvp5c could not satisfy the relative error tolerance ---> no figures could be displayed\n'); +else +lw=1.5; %Linewidth for all figures +if fig(1) == 1 +f=figure(1); % +f.Position = [276.3333 196.3333 1100 1000]; +x2=x(2:end); +% total mass flows +subplot(3,2,1);lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on +yyaxis right +plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% concentrations +subplot(3,2,3);lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on +yyaxis right; C_f=C_f(2:end); +plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Pressures +subplot(3,2,5);lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on +yyaxis right +plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% total mass flows +subplot(3,2,2);lc='#0072bD';rc='#77AC30'; +plot(x*L, Q_d*Q_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on +yyaxis right +plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% concentrations +subplot(3,2,4);lc='#0072bD';rc='#77AC30'; +plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on +yyaxis right; c_f=c_f(2:end); +plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','northwest');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0; +% Pressures +subplot(3,2,6);lc='#0072bD';rc='#77AC30'; +plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on +yyaxis right +plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +f=figure(2); % +f.Position = [199 232.3333 1.0987e+03 639.3333]; +% Permeate flows +subplot(2,2,1); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end); +plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on +yyaxis right +plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Osmotic/ Hydraulic pressure difference +subplot(2,2,3);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]); +plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on +yyaxis right; +plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Permeate flows +subplot(2,2,2); lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end); +plot(x2*L, -Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on +yyaxis right +plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','southwest');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Osmotic/ Hydraulic pressure difference +subplot(2,2,4);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]); +plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); +yyaxis right; +plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','southeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]); + + + +end +%% figure 2 +if fig(2) == 1 +figure(2); % pressure+concentration in both PRO channels +subplot(2,2,1);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_d(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(2,2,2);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(2,2,3);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_d(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +x2=x(2:end);c_f2=100.*c_f(2:end); +subplot(2,2,4);plot(x2*L, c_f2*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +end +%% figure 3 +if fig(3) == 1 +figure(3); % Q_s,in and Q_w,in +x2=x(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end); +subplot(1,2,1);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{w,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(1,2,2);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{s,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +end +%% figure 4 +lw=1.5; %Linewidth for all figures +if fig(4) == 1 +figure(4); % Beta of PRO + if version(2)==0 + plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]);ylim([swro_beta(1)*J_r/swro_x_r*.9 swro_beta(1)*J_r/swro_x_r*1.1]); hold on + else + plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]); hold on + end +end +%% figure 5 +if fig(5) == 1 +figure(5); % pressure+concentration in both SWRO channels +subplot(2,2,1);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5; +x2=x(2:end);C_f2=100.*C_f(2:end); +subplot(2,2,2);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(2,2,4);plot(x2*swro_L, C_f2*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(2,2,3);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +end +%% figure 6 +if fig(6) == 1 +figure(6); % J_s,in and J_w,in +x2=x(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end); +subplot(1,2,2);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(1,2,1);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +end +%% figure 7 +if fig(7) == 1 +f=figure(7); % PRO quantities +f.Position = [100 1 1100 1000]; +subplot(4,3,4);plot(x*L, Q_sd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,1);plot(x*L, Q_wd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,5);plot(x*L, Q_sf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,2);plot(x*L, Q_wf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,7);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(4,3,8);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(4,3,10);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +x2=x(2:end);c_f=100.*c_f(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end); +subplot(4,3,11);plot(x2*L, c_f*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,[3 6]);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,[9 12]);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; +end +%% figure 8 +if fig(8) == 1 +f=figure(8); % SWRO quantities +f.Position = [100 1 1100 1000]; +x2=x(2:end);J_wf2=J_wf(2:end);J_sf2=J_sf(2:end); +subplot(4,3,1);plot(x*swro_L, J_wd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,4);plot(x*swro_L, J_sd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,2);plot(x2*swro_L, J_wf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,5);plot(x2*swro_L, J_sf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,7);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(4,3,8);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5; +subplot(4,3,10);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +C_f=100.*C_f(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end); +subplot(4,3,11);plot(x2*swro_L, C_f*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,[3 6]);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +subplot(4,3,[9 12]);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0; +end +%% figure 9 +if fig(9) == 1 +f=figure(3); % ERD quantities +f.Position = [100 1 900 500]; tiledlayout(1,4); +rot=0; rx=0; %rotation and shift of texts +vals1 = [pE*p_r; pERD*p_r; 0; P_d(end)*p_r; p_d(end)*p_r]; +vals2 = 100*[cE; C_ERD; 0; C_d(end); c_d(1)]; +vals3 = [J_E*J_r; J_ERD*J_r; 0; J_d(end)*J_r; Q_d(1)*Q_r]; +%adjust for counter current / SWRO+ERD / 2ERDs +if version(1)==1; vals2(5) = 100*c_d(end); vals3(5) =-Q_d(end)*Q_r;end +if version(6)==1; vals1(5) = p_exit*p_r; vals2(5) = 100*c_exit; vals3(5) = Q_exit*J_r; end +if version(6)==4; vals1(1)= pERD2*p_r; vals2(1)=100*C_ERD2; end +nexttile +b = bar(1,vals1, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD'; +xticklabels({'[kg/ms^2]'}); ylim([0 1.1*max(vals1)]);ay=gca;ay.YAxis(1).Exponent=5; +t=text(b(1).XEndPoints,b(1).YEndPoints,"P_{s}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); +t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"P_{s}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); +t=text(b(4).XEndPoints,b(4).YEndPoints,"P_{b}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); +t=text(b(5).XEndPoints,b(5).YEndPoints,"P_{b}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); +nexttile +b = bar(1,vals2, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD'; +xticklabels({'[%]'}); ylim([0 1.12*max(vals2)]);ay=gca;ay.YAxis(1).Exponent=0; +t= text(b(1).XEndPoints,b(1).YEndPoints,"C_{s}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"C_{s}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(4).XEndPoints,b(4).YEndPoints,"C_b^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(5).XEndPoints,b(5).YEndPoints,"C_{b}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +nexttile +b = bar(1,vals3, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD'; +xticklabels({'[kg/ms]'}); ylim([0 1.1*max(vals3)]);ay=gca;ay.YAxis(1).Exponent=0; +t=text(b(1).XEndPoints,b(1).YEndPoints,"J_{s}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"J_{s}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(4).XEndPoints ,b(4).YEndPoints,"J_{b}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(5).XEndPoints,b(5).YEndPoints,"J_{b}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot ; +nexttile +if version(6)>3 +lgd = legend([b(1) b(2) b(5) b(4)],'LP ERD 1 Seawater Inlet','HP ERD 1 Seawater Outlet','HP SWRO Brine Outlet','LP PRO Brine Inlet' , 'Location', 'EastOutside'); +else +lgd = legend([b(1) b(2) b(5) b(4)],'Sea water inlet','Sea water outlet','Brine outlet','Brine inlet' , 'Location', 'EastOutside'); +end +lgd.Layout.Tile = 4; axis off ; +end +%% figure 10 +if fig(9) == 1 && version(6) ==4 +close all;f=figure(3); % ERD quantities +f.Position = [950 1 900 500]; tiledlayout(1,4); +rot=0; rx=0; %rotation and shift of texts +vals1 = [pE*p_r; pERD2*p_r; 0; p_d(end)*p_r; p_exit*p_r]; +vals2 = 100*[cE; C_ERD2; 0; c_d(end); c_exit]; +vals3 = [J_E_sea*J_r; 2*J_E*J_r; 0; mixer_ERD*Q_d(end)*Q_r; Q_exit*Q_r]; +%adjust for counter current +if version(1)==1; vals1(4) = p_d(1)*p_r; vals2(4)=c_d(1)*100;vals3(4)=-mixer_ERD*Q_d(1)*Q_r; end +nexttile +b = bar(1,vals1, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD'; +xticklabels({'[kg/ms^2]'}); ylim([0 1.1*max(vals1)]);ay=gca;ay.YAxis(1).Exponent=5; +t=text(b(1).XEndPoints,b(1).YEndPoints,"P_E",'HorizontalAlignment','center','VerticalAlignment','bottom'); +t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"P_{ERD 2}",'HorizontalAlignment','center','VerticalAlignment','bottom'); +if version(1) == 0; t=text(b(4).XEndPoints,b(4).YEndPoints,"p_d^1",'HorizontalAlignment','center','VerticalAlignment','bottom'); end +if version(1) == 1; t=text(b(4).XEndPoints,b(4).YEndPoints,"p_d^0",'HorizontalAlignment','center','VerticalAlignment','bottom'); end +t=text(b(5).XEndPoints,b(5).YEndPoints,"P_E",'HorizontalAlignment','center','VerticalAlignment','bottom'); +nexttile +b = bar(1,vals2, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD'; +xticklabels({'[%]'}); ylim([0 1.12*max(vals2)]);ay=gca;ay.YAxis(1).Exponent=0; +t= text(b(1).XEndPoints,b(1).YEndPoints,"C_E",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(2).XEndPoints+rx ,b(2).YEndPoints," C_{ERD 2}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +if version(1) == 0; t=text(b(4).XEndPoints,b(4).YEndPoints,"c_d^1",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; end +if version(1) == 1; t=text(b(4).XEndPoints,b(4).YEndPoints,"c_d^0",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; end +t=text(b(5).XEndPoints,b(5).YEndPoints,"C_{exit}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +nexttile +b = bar(1,vals3, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD'; +xticklabels({'[kg/ms]'}); ylim([0 1.1*max(vals3)]);ay=gca;ay.YAxis(1).Exponent=0; +t=text(b(1).XEndPoints,b(1).YEndPoints,"J_{sea}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +t=text(b(2).XEndPoints+rx ,b(2).YEndPoints," 2J_{E}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +strg=[num2str(mixer_ERD), '*']; +if version(1) == 0; t=text(b(4).XEndPoints ,b(4).YEndPoints,[strg,"$Q_d^1$"],'Interpreter','latex','HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;end +if version(1) == 1; t=text(b(4).XEndPoints ,b(4).YEndPoints,[strg,"$-Q_d^0$"],'Interpreter','latex','HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;end +t=text(b(5).XEndPoints,b(5).YEndPoints,"Q_{exit}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; +nexttile +lgd = legend([b(1) b(2) b(3) b(4)],'LP ERD 2 Seawater Inlet','HP ERD 2 Seawater Outlet','HP ERD 2 Brine Inlet','LP ERD 2 Brine Outlet' , 'Location', 'EastOutside'); +lgd.Layout.Tile = 4; axis off ; +end +end +%% figures for Paper +if fig(11) == 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f=figure(1); % SWRO quantities +f.Position = [276.3333 196.3333 1100 1000]; +x2=x(2:end); +% total mass flows +subplot(3,2,1);lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on +yyaxis right +plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% concentrations +subplot(3,2,2);lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on +yyaxis right; C_f=C_f(2:end); +plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Pressures +subplot(3,2,3);lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on +yyaxis right +plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Permeate flows +subplot(3,2,5); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end); +plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on +yyaxis right +plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Osmotic/ Hydraulic pressure difference +subplot(3,2,6);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]); +plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on +yyaxis right; +plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f=figure(2); % PRO quantities +f.Position = [936.3333 266.3333 1.0987e+03 632]; +x2=x(2:end); +% total mass flows +subplot(3,2,1);lc='#0072bD';rc='#77AC30'; +plot(x*L, Q_d*Q_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on +yyaxis right +plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% concentrations +subplot(3,2,2);lc='#0072bD';rc='#77AC30'; +plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on +yyaxis right; c_f=c_f(2:end); +plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0; +% Pressures +subplot(3,2,3);lc='#0072bD';rc='#77AC30'; +plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on +yyaxis right +plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Permeate flows +subplot(3,2,5); lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end); +plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on +yyaxis right +plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +% Osmotic/ Hydraulic pressure difference +subplot(3,2,6);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]); +plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); +yyaxis right; +plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]); +end +%% figures for Paper - all in individual plot +if fig(11) == 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f=figure(1);% SWRO quantities +x2=x(2:end); +% total mass flows +lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0; xlim([0 swro_L]); hold on +yyaxis right +plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROa.png'); +% concentrations +close all; f=figure(1); +lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on +yyaxis right; C_f=C_f(2:end); +plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROb.png'); +% Pressures +close all; f=figure(6); +lc='#0072bD';rc='#77AC30'; +plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on +yyaxis right +plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROc.png'); +% Permeate flows +close all;f=figure(1); +lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end); +plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on +yyaxis right +plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROd.png'); +% Osmotic/ Hydraulic pressure difference +close all;f=figure(1); +lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]); +plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on +yyaxis right; +plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROe.png'); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f=figure(2); % PRO quantities +x2=x(2:end); +% total mass flows +lc='#0072bD';rc='#77AC30'; +plot(x*L, Q_d*Q_r,'Color', lc,'LineWidth',lw);ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on +yyaxis right +plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROa.png'); +% concentrations +close all; f=figure(2);lc='#0072bD';rc='#77AC30'; +plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on +yyaxis right; c_f=c_f(2:end); +plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROb.png'); +% Pressures +close all; f=figure(2);lc='#0072bD';rc='#77AC30'; +plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on +yyaxis right +plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROc.png'); +% Permeate flows +close all; f=figure(2);lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end); +plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on +yyaxis right +plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROd.png'); +% Osmotic/ Hydraulic pressure difference +close all; f=figure(2);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]); +plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); +yyaxis right; +plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]); +set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROe.png'); +end + + +%ax = axes; yyaxis('left'); yyaxis('right'); +%ax.YAxis(1).Color = [0 0 0]; ax.YAxis(2).Color = [0 0 0]; + + + + + + + + + + + + + + + + + + + + + + + + + + +end +end \ No newline at end of file diff --git a/fun_1.m b/fun_1.m index 406a9d2..1cfef86 100644 --- a/fun_1.m +++ b/fun_1.m @@ -28,8 +28,8 @@ function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,o % tic; %% Read data -if option_data == 0; DATA = @(x)Sim_0_data(input); end - +if option_data == 0; DATA = @(x)Sim_0a_data(input); end +if option_data == 0.1; DATA = @(x)Sim_0b_data(input); end @@ -80,108 +80,96 @@ if version(1)==1; y_init = [cE; J_wd_0; J_sf_0; J_wf_0; Pd_0; Pf_0; cE; -Q_wd_0; solinit = bvpinit(x,y_init); %% Solve the BVP -if version(6) == 0; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary2(ya,yb, DATA),solinit,ode_options); end -if version(6) == 1; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary1(ya,yb, DATA),solinit,ode_options); end -if version(6) == 2; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary2(ya,yb, DATA),solinit,ode_options); end -if version(6) == 3; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary3(ya,yb, DATA),solinit,ode_options); end -if version(6) == 4; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary4(ya,yb, DATA),solinit,ode_options); end +if version(6) == 0; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end +if version(6) == 1; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary1(ya,yb, DATA),solinit,ode_options); end +if version(6) == 2; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end +if version(6) == 3; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary3(ya,yb, DATA),solinit,ode_options); end +if version(6) == 4; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary4(ya,yb, DATA),solinit,ode_options); end %% Evaluate the solution n=100; x=linspace(0,1,n); y = deval(sol,x); Y = y'; Y = real(Y); -C_d=zeros(1,n);C_f=zeros(1,n);J_sd=zeros(1,n);J_wd=zeros(1,n);J_d=zeros(1,n);J_sf=zeros(1,n);J_wf=zeros(1,n);J_f=zeros(1,n);P_d=zeros(1,n);P_f=zeros(1,n);swro_p_osm_d=zeros(1,n);swro_p_osm_f=zeros(1,n);swro_beta=zeros(1,n);Y(:,4+6)=abs(Y(:,4+6));swro_del_c=zeros(1,n);J_cross=zeros(1,n);J_sin=zeros(1,n); -c_d=zeros(1,n);c_f=zeros(1,n);Q_sd=zeros(1,n);Q_wd=zeros(1,n);Q_d=zeros(1,n);Q_sf=zeros(1,n);Q_wf=zeros(1,n);Q_f=zeros(1,n);p_d=zeros(1,n);p_f=zeros(1,n);p_osm_d=zeros(1,n);p_osm_f=zeros(1,n);beta=zeros(1,n);del_c=zeros(1,n);Q_cross=zeros(1,n);Q_sin=zeros(1,n); local_ro_d=zeros(1,n);local_ro_f=zeros(1,n);swro_local_ro_d=zeros(1,n);swro_local_ro_f=zeros(1,n); - - for i=1:n - %% SWRO evaluation - % Concentrations - C_d(i)= Y(i,1); % concentration of draw side - C_f(i)= Y(i,3)./Y(i,4); % concentration of fresh side - % Flow rates in draw side - J_sd(i)= Y(i,1)*Y(i,2); % salt flow rate in draw side - J_wd(i)= Y(i,2); % water flow rate in draw side - J_d(i) = Y(i,1)*Y(i,2)+Y(i,2); % total flow rate draw side - % Flow rates in fresh side - J_sf(i)= Y(i,3); % salt flow rate in fresh side - J_wf(i)= Y(i,4); % water flow rate in fresh side - J_f(i) = Y(i,3)+Y(i,4); % total flow rate fresh side - % Pressures - P_d(i)=Y(i,5); % Pressure draw side - P_f(i)=Y(i,6); % Pressure fresh side - - % Osmotic pressure - swro_p_osm_d(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,1)/Ms)/p_r; - swro_p_osm_f(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,3)/Ms./Y(i,4))/p_r; - % Concentration difference - swro_del_c(i) = Y(i,1)./(Y(i,1) + 1) - Y(i,3)./(Y(i,3) + Y(i,4)); - % Local density - swro_local_ro_d(i) = (Y(i,1) + 1)./(ro_water*Y(i,1)/ro_salt + 1); % local density of draw side - swro_local_ro_f(i) = (Y(i,3) + Y(i,4))./(ro_water*Y(i,3)/ro_salt + Y(i,4)); % local density of fresh side - % B(x) - swro_beta(i) = (1 - swro_R)*((Y(i,5) - Y(i,6)) - sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)))/swro_R; - if version(2) == 0 ;swro_beta(i)= swro_beta_fix; end - if version(4) == 0 ;swro_beta(i)= 0; end - % J_w,in and J_s,in - - J_cross(i) = ((Y(i,5) - Y(i,6)) - sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i))); - if version(4) == 0; J_cross(i) = ((Y(i,5) - Y(i,6)) - sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i))); end - J_sin(i)=swro_beta(i)*swro_del_c(i); - - % J_w,in and J_s,in with ECP and ICP - if version(7)==1 - J_cross(i) = ((Y(i,5)-Y(i,6))-sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)) + swro_beta(i)*(Y(i,5)-Y(i,6))*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha*swro_beta(i)*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d(i)/swro_KD -swro_p_osm_f(i)*(swro_KK + 1/swro_KF) )); - J_sin(i) = swro_beta(i)*( swro_del_c(i)-J_cross(i)*( Y(i,1)./(Y(i,1) + 1)./swro_KD + Y(i,3)./(Y(i,3) + Y(i,4))*(swro_KK+1/swro_KF)))/(1+swro_alpha*swro_beta(i)*(swro_KK+1/swro_KF+1/swro_KD)); - end - + +%% SWRO evaluation +C_d = Y(:,1); +C_f = Y(:,3)./Y(:,4); +J_sd = Y(:,1).*Y(:,2); +J_wd = Y(:,2); +J_d = Y(:,1).*Y(:,2)+Y(:,2); +J_sf = Y(:,3); +J_wf = Y(:,4); +J_f = Y(:,3)+Y(:,4); +P_d = Y(:,5); +P_f = Y(:,6); +% same quantities as in "ODEsystem.m" +swro_p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1)/Ms)/p_r; +swro_p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3)/Ms./Y(:,4))/p_r; +swro_del_c = Y(:,1)./(Y(:,1) + 1) - Y(:,3)./(Y(:,3) + Y(:,4)); +swro_local_ro_d = (Y(:,1) + 1)./(ro_water*Y(:,1)/ro_salt + 1); +swro_local_ro_f = (Y(:,3) + Y(:,4))./(ro_water*Y(:,3)/ro_salt + Y(:,4)); +swro_beta = (1 - swro_R)*((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))/swro_R; +if version(2) == 0 ;swro_beta= swro_beta_fix.*ones(1,n); end +if version(4) == 0 ;swro_beta=zeros(1,n); end +% J_w,in and J_s,in +J_cross = ((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f)); +if version(4) == 0; J_cross = ((Y(:,5) - Y(:,6)) - sigma*(swro_p_osm_d - swro_p_osm_f)); end +J_sin=swro_beta.*swro_del_c; +if version(7)==1 +J_cross = ((Y(:,5)-Y(:,6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta.*(Y(:,5)-Y(:,6)).*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha.*swro_beta.*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d./swro_KD -swro_p_osm_f.*(swro_KK + 1/swro_KF) )); +J_sin = swro_beta.*( swro_del_c-J_cross.*( Y(:,1)./(Y(:,1) + 1)./swro_KD + Y(:,3)./(Y(:,3) + Y(:,4)).*(swro_KK+1/swro_KF)))./(1+swro_alpha.*swro_beta.*(swro_KK+1/swro_KF+1/swro_KD)); +end + %% PRO evaluation - % Concentrations - c_d(i)= Y(i,1+6); % concentration of draw side - c_f(i)= Y(i,3+6)./Y(i,4+6); % concentration of fresh side - % Flow rates in draw side - Q_sd(i)= Y(i,1+6)*Y(i,2+6); % salt flow rate in draw side - Q_wd(i)= Y(i,2+6); % water flow rate in draw side - Q_d(i) = Y(i,1+6)*Y(i,2+6)+Y(i,2+6); % total flow rate draw side - % Flow rates in fresh side - Q_sf(i)= Y(i,3+6); % salt flow rate in fresh side - Q_wf(i)= Y(i,4+6); % water flow rate in fresh side - Q_f(i) = Y(i,3)+Y(i,4+6); % total flow rate fresh side - % Pressures - p_d(i)=Y(i,5+6); % Pressure draw side - p_f(i)=Y(i,6+6); % Pressure fresh side - - % Osmotic pressure - p_osm_d(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,1+6)/Ms)/p_r; - p_osm_f(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,3+6)/Ms./Y(i,4+6))/p_r; - % Concentration difference - del_c(i) = Y(i,1+6)./(Y(i,1+6) + 1) - Y(i,3+6)./(Y(i,3+6) + Y(i,4+6)); - % local density - local_ro_d(i) = (Y(i,1+6) + 1)./(ro_water*Y(i,1+6)/ro_salt + 1); - local_ro_f(i) = (Y(i,3+6) + Y(i,4+6))./(ro_water*Y(i,3+6)/ro_salt + Y(i,4+6)); - % B(x) - beta(i) = (1 - R)*((p_osm_d(i) - p_osm_f(i)) - (Y(i,5+6) - Y(i,6+6)))/R; - if version(3) == 0 ; beta(i) = beta_fix; end - if version(5) == 0 ; beta(i) = 0; end - % Q_w,in and Q_s,in - Q_cross(i) = (p_osm_d(i) - p_osm_f(i) - (Y(i,5+6) - Y(i,6+6)) - (Y(i,5+6) - Y(i,6+6))*beta(i)*KK*alpha*p_r)./(1 + KK*alpha*p_r*(beta(i) +p_osm_f(i))); - if version(5) == 0; Q_cross(i) = ((p_osm_d(i) - p_osm_f(i)) - (Y(i,5+6) - Y(i,6+6))) ; end - Q_sin(i)=beta(i)*del_c(i); - - if version(8)==1 %% version with ICP and ECP - Q_cross(i)=((p_osm_d(i) - p_osm_f(i)) - (Y(i,5+6) - Y(i,6+6)) - (Y(i,5+6) - Y(i,6+6))*beta(i)*(KK+1/KF+1/KD)*p_r*alpha)./(1 + alpha*beta(i)*(KK+1/KF+1/KD) + alpha*p_r*(p_osm_f(i)*(KK+1/KF) + p_osm_d(i)/KD)); - Q_sin(i)=beta(i)*( del_c(i)-Q_cross(i)*( Y(i,1+6)./(Y(i,1+6) + 1)./KD + Y(i,3+6)./(Y(i,3+6) + Y(i,4+6))*(KK+1/KF)))/(1+alpha*beta(i)*(KK+1/KF+1/KD)); - end - - end - % In ideal cases take mean of some quantities (one could instead use ylim on y-axis) - if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end - if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end +c_d = Y(:,1+6); % concentration of draw side +c_f = Y(:,3+6)./Y(:,4+6); % concentration of fresh side +Q_sd = Y(:,1+6).*Y(:,2+6); % salt flow rate in draw side +Q_wd = Y(:,2+6); % water flow rate in draw side +Q_d = Y(:,1+6).*Y(:,2+6)+Y(:,2+6); % total flow rate draw side +Q_sf = Y(:,3+6); % salt flow rate in fresh side +Q_wf = Y(:,4+6); % water flow rate in fresh side +Q_f = Y(:,3)+Y(:,4+6); % total flow rate fresh side +p_d = Y(:,5+6); % Pressure draw side +p_f = Y(:,6+6); % Pressure fresh side +% same quantities as in "ODEsystem.m" +p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1+6)/Ms)/p_r; +p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3+6)/Ms./Y(:,4+6))/p_r; +del_c = Y(:,1+6)./(Y(:,1+6) + 1) - Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)); +local_ro_d = (Y(:,1+6) + 1)./(ro_water*Y(:,1+6)/ro_salt + 1); +local_ro_f = (Y(:,3+6) + Y(:,4+6))./(ro_water*Y(:,3+6)/ro_salt + Y(:,4+6)); +beta = (1 - R)*((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)))./R; +if version(3) == 0 ; beta = beta_fix.*ones(1,n); end +if version(5) == 0 ; beta = zeros(1,n); end +% Q_w,in and Q_s,in +Q_cross = (p_osm_d - p_osm_f - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*KK.*alpha.*p_r)./((ones(1,n) + KK*alpha*p_r.*(beta +p_osm_f'))'); +if version(5) == 0; Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6))) ; end +Q_sin = beta'.*del_c; +if version(8)==1 +Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*(KK+1./KF+1./KD).*p_r.*alpha)./(ones(1,n) + alpha.*beta.*(KK+1./KF+1./KD) + alpha.*p_r.*(p_osm_f'.*(KK+1./KF) + p_osm_d'./KD))'; +Q_sin = beta'.*( del_c-Q_cross.*( Y(:,1+6)./(Y(:,1+6) + 1)./KD + Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)).*(KK+1/KF)))./(1+alpha.*beta'); +end -%% Analysis of Recovery rates +%% Freswater production and recovery rates Freshwater = J_wf(end)*J_r; %in [kg/s] FW = (Freshwater*swro_Z/(swro_local_ro_f(end)*rho_r))*3600; %in [m^3/h] SWRO_Recovery =(J_f(end)./J_d(1))*100; % SWRO freshwater recovery [%] PRO_Recovery = (1- Q_f(end)./Q_f(1))*100; % PRO recovery rate [%] +% warning flags: + if FW < 0 + fprintf(2,' \nERROR: Waring: Freshwater production is negative! \n'); + end + if SWRO_Recovery > 100 + fprintf(2,' \nERROR: Waring: SWRO recovery is greater then 100 %% \n'); + end + if SWRO_Recovery < 0 + fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n'); + end + if PRO_Recovery > 100 + fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n'); + end + if PRO_Recovery < 0 + fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n'); + end %% Output of System J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; J_wE = J_E/(cE+1); @@ -202,14 +190,13 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r; rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE; f_1= (Q_exit*swro_Z/rho_exit*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2 ); f_2= (J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2); - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; - f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; - if pERD > pMax; pERD=pMax; end - if pERD < 0; pERD=pMax; end - %display(['pERD: ',num2str(pERD)]) % pERD use in the model - %display(['pERD_theoretical: ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0; + if W_p3 < 0 + fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n'); + end + end %% Version(6)=2 --> only PRO if version(6)==2 @@ -220,33 +207,25 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r; if version(1) ==1; W_p2= 1/LP_eff * (p_d(end)-pE)*(abs(Q_d(end))*Z)./local_ro_d(end); end W_p2 = W_p2*W_r; end - + %% Version(6)=3 --> SWRO-PRO hybrid system (with one ERD) if version(6)==3 rho_ERD = V_m*(swro_local_ro_d(end) - rho_E)+rho_E; - C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt); - J_wERD = J_E/(C_ERD+1); qq=Q_r/J_r; + C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt); + qq=Q_r/J_r; if version(1) == 0 %co-current f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2; f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2; - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; - f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; - if pERD > pMax; pERD=pMax; end - %if pERD < 0; pERD=pMax; end W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); end if version (1) == 1 %counter-current f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2 ; f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2 ; - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; - f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; - %display(['pERD: ',num2str(pERD)]) % pERD use in the model - %display(['pERD_theoretical: ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency - if pERD > pMax; pERD=pMax; end - %if pERD < 0; pERD=pMax; end W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1); end @@ -260,8 +239,7 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r; C_ERD2= -ro_salt*(rho_ERD2-1)/(rho_ERD2*ro_water-ro_salt); rho_ERD= V_m*(swro_local_ro_d(end) - rho_ERD2)+rho_ERD2; C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt); - J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; - J_wE = J_E/(C_ERD2+1); J_wERD= J_ERD/(C_ERD+1); + J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; %J_wE = J_E/(C_ERD2+1); J_wERD= J_ERD/(C_ERD+1); qq=Q_r/J_r; J_E_sea = 2*J_E; J_wE_sea = J_E_sea/(cE+1); @@ -278,45 +256,33 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r; if version(1) == 0 %co-current f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2; f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end)/A_ERD)^2; - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; pERD2 = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; - f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; - if pERD2 > pMax; pERD2=pMax; end - %if pERD2 < 0; pERD2=pMax; end - end if version (1) == 1 %counter-current f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2; f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1)/A_ERD)^2; - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; pERD2 = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; - f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; - if pERD2 > pMax; pERD2=pMax; end - %if pERD2 < 0; pERD2=pMax; end + end %calculation of Pressure P_ERD if version(1) == 0 %co-current f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2; f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2; - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; - f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; - if pERD > pMax; pERD=pMax; end - %if pERD < 0; pERD=pMax; end W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_t = T_eff * (p_d(end)-pE)*((1-mixer_ERD)*Q_d(end)*Z)/local_ro_d(end); end if version (1) == 1 %counter-current f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2 ; f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2 ; - f_ERD= ERD_fric*(ERD_eff*f_2 - f_1); - pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; - f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; - if pERD > pMax; pERD=pMax; end - %if pERD < 0; pERD=pMax; end + f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100; + pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_t = T_eff * (p_d(1)-pE)*(abs((1-mixer_ERD)*Q_d(1))*Z)/local_ro_d(1); - end + end % Recalculation for Pump 1 W_p1 = 1/HP_eff * (P_d(1)-pERD2)*(J_E *swro_Z)/rho_ERD2; W_p1 = W_p1*swro_W_r; end @@ -333,6 +299,13 @@ if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600 SEC_net = W_net*(swro_local_ro_f(end)*rho_r)./(J_f(end)*J_r*swro_Z)/1000/3600; %in [kWh/m^3] Rev= pw*FW + pe*SEC_net*FW; %in [$/h] tt=toc; + + + + + + + %% if BVP-solver failed catch sol.stats.maxerr =1000; @@ -808,5 +781,77 @@ end + +end end -end \ No newline at end of file + + + + + + + + + + + + + + + + + + + + + + +%% Removed comments + +%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here?? + +% -3.9191 0.3247 0.0131 32.8727 7.1201 NaN NaN +% +% 1.0e+03 * +% -1.2742 -0.9744 0 -0.4873 -0.0010 0.1885 + +% In ideal cases take mean of some quantities. +% In ideal cases salt will fluctuate around 1e-16. We take mean here for +% +% This is just for figures. (), otherwise one could use ylim ) +%if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end +%if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end + + %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + %display(['pERD: ',num2str(pERD)]) % pERD use in the model + %display(['pERD_theoretical: ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency + + %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + + %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z; + %display(['pERD: ',num2str(pERD)]) % pERD use in the model + %display(['pERD_theoretical: ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + + %J_wERD = J_E/(C_ERD+1); + + %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; + %if pERD2 > pMax; pERD2=pMax; end + %if pERD2 < 0; pERD2=pMax; end + + %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) + pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z; + %if pERD2 > pMax; pERD2=pMax; end + %if pERD2 < 0; pERD2=pMax; end + + %f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end + + %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) + pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z; + %if pERD > pMax; pERD=pMax; end + %if pERD < 0; pERD=pMax; end \ No newline at end of file -- GitLab