From 4792f3f5bb0185dd079824f3828eddff76554432 Mon Sep 17 00:00:00 2001 From: is05 Date: Thu, 3 Apr 2025 17:16:09 +0000 Subject: [PATCH] Create function Meidanshahi_Complete_gaoptimset() --- function Meidanshahi_Complete_gaoptimset() | 476 +++++++++++++++++++++ 1 file changed, 476 insertions(+) create mode 100644 function Meidanshahi_Complete_gaoptimset() diff --git a/function Meidanshahi_Complete_gaoptimset() b/function Meidanshahi_Complete_gaoptimset() new file mode 100644 index 0000000..762f565 --- /dev/null +++ b/function Meidanshahi_Complete_gaoptimset() @@ -0,0 +1,476 @@ +function Meidanshahi_Complete_gaoptimset() +%-------------------------------------------------------------------------- +% MEIDANSHAHI_COMPLETE_GAOPTIMSET : +% Script complet en un seul fichier pour l'optimisation du dessalage +% électrostatique (PBE). Identique à la version 'optimoptions', sauf +% qu'on utilise 'gaoptimset' pour configurer ga(). +% +% Optimisation de 4 variables (X, Qd, T, X_recycle) : +% 0) X ∈ [0..10] % (taux eau de lavage) +% 1) Qd ∈ [0..5] L/h (désémulsifiant) +% 2) T ∈ [50..120] °C (température) +% 3) X_recycle ∈ [0..100] % (recyclage purge) +% +% Objectif : minimiser f = (salinité / 40) - (efficacité / 100). +%-------------------------------------------------------------------------- + + %% 1. PARAMÈTRES GLOBAUX, PHYSIQUES ET DISCRÉTISATION + + % Densités (kg/m³) + rho_h2o = 1000; + rho_oil = 850; + + % Viscosités + mu_water = 1e-3; % eau + A_mu = 11.98; + B_mu = 0.0198; % => mu_oil = A_mu * exp(-B_mu*(T en K)) + + % Tension interfaciale + sigma_base = 0.03; % ~30 mN/m + + % Gravité, permittivité + g = 9.81; + epsilon0 = 8.854e-12; + + % Vanne de mélange + DeltaP_valve = 1e5; % Pa + t_res_valve = 0.1; % s + v_char_valve = sqrt(2*DeltaP_valve / rho_oil); + + % Cassure / coalescence + We_crit = 1.0; + Ca_crit = 0.5; + K_break = 1.0; + K_coal = 1.0; + + % Brute initiale + initial_water_frac = 0.03; % ~3 % + initial_salt_conc = 100; % mg sel / L brut + oil_flow_Lh = 1e5; % 100 m³/h => 100000 L/h + + % Temps de séjour dessaleurs + t_res_desalter = 60; % s + dt_desalter = 1; % s (pas) + + % Séparation gravitaire + d_sep = 3e-4; % 300 µm + + % Recyclage + recycle_tol = 1e-3; + + % Désémulsifiant + a_sigma = 1; + b_sigma = 0.1; + sigma_min_factor = 0.01; + + % Champ électrique (AC) + E_field = 3e5; % 3 kV/cm => 300000 V/m + + % Discrétisation PBE (méthode des classes) + N_classes = 20; + d_min = 1e-6; % 1 µm + d_max = 1e-3; % 1 mm + d_classes = logspace(log10(d_min), log10(d_max), N_classes); + v_classes = (pi/6)*d_classes.^3; + lambda = d_min; % échelle turbulence + + %% 2. OPTIMISATION AVEC gaoptimset (pour versions MATLAB plus anciennes) + + % Bornes : [X, Qd, T, X_recycle] + lb = [0, 0, 50, 0 ]; + ub = [10, 5, 120, 100 ]; + + % Définition des options via gaoptimset + opts = gaoptimset('Display','iter', ... + 'PopulationSize',20, ... + 'Generations',20); + % Remarque : si vous aviez un MATLAB plus récent ET le Global Optimization Toolbox, + % vous pourriez mettre 'MaxGenerations' => 'Generations'. + % Mais sur de vieilles versions, le param s'appelle 'Generations'. + + % Appel de ga (4 variables) + [opt_vars, fval] = ga(@(vars) localObj(vars, ... + rho_h2o, rho_oil, mu_water, A_mu, B_mu, sigma_base, g, epsilon0, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, recycle_tol, oil_flow_Lh, ... + a_sigma, b_sigma, sigma_min_factor, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field), ... + 4, [], [], [], [], lb, ub, [], opts); + + X_opt = opt_vars(1); + Qd_opt = opt_vars(2); + T_opt = opt_vars(3); + Xrec_opt = opt_vars(4); + + % Simulation finale + [sal_opt, eff_opt, prof_opt] = simulateDesalting(X_opt, Qd_opt, T_opt, Xrec_opt, ... + rho_h2o, rho_oil, mu_water, A_mu, B_mu, sigma_base, g, epsilon0, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, recycle_tol, oil_flow_Lh, ... + a_sigma, b_sigma, sigma_min_factor, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field); + + % Affichage + fprintf('\n=== Résultats Optimaux (gaoptimset) ===\n'); + fprintf('X = %.2f %%\n', X_opt); + fprintf('Qd = %.2f L/h\n', Qd_opt); + fprintf('T = %.1f °C\n', T_opt); + fprintf('X_recycle = %.1f %%\n', Xrec_opt); + fprintf('Salinité = %.2f mg/L\n', sal_opt); + fprintf('Efficacité = %.1f %%\n', eff_opt); + fprintf('fval = %.2f\n', fval); + + %% 3. FIGURES : DISTRIBUTION & ANALYSE DE SENSIBILITÉ + + % 3.1 Distribution granulométrique + N_init = prof_opt.initial.N; + N_mix = prof_opt.afterMix.N; + N_d1 = prof_opt.afterDesal1.N; + N_d2 = prof_opt.afterDesal2.N; + + vol_init = N_init.*v_classes; + vol_mix = N_mix .*v_classes; + vol_d1 = N_d1 .*v_classes; + vol_d2 = N_d2 .*v_classes; + + sum_init = sum(vol_init); + pct_init = 100*(vol_init / sum_init); + pct_mix = 100*(vol_mix / sum_init); + pct_d1 = 100*(vol_d1 / sum_init); + pct_d2 = 100*(vol_d2 / sum_init); + + figure; + semilogx(d_classes*1e6, pct_init,'-k','LineWidth',1.5); hold on; + semilogx(d_classes*1e6, pct_mix ,'--b','LineWidth',1.5); + semilogx(d_classes*1e6, pct_d1 ,'-.r','LineWidth',1.5); + semilogx(d_classes*1e6, pct_d2 ,':m','LineWidth',2); + hold off; + xlabel('Diamètre (µm)'); + ylabel('Volume eau (%) relatif à l''initial'); + legend('Initial','Après vanne','Après 1er dessal','Après 2e dessal','Location','best'); + title('Distribution de taille des gouttelettes (gaoptimset)'); + grid on; + + % 3.2 Analyse de sensibilité + vars_base = [X_opt, Qd_opt, T_opt, Xrec_opt]; + var_names = {'X(%)','Qd(L/h)','T(°C)','X_{recycle}(%)'}; + var_ranges = { + linspace(0,10,6), + linspace(0,5,6), + linspace(50,120,6), + linspace(0,100,6) + }; + + for v=1:4 + vals = var_ranges{v}; + sal_vals = zeros(size(vals)); + eff_vals = zeros(size(vals)); + + for i=1:length(vals) + testVars = vars_base; + testVars(v) = vals(i); + + [sal_vals(i), eff_vals(i)] = simulateDesalting(testVars(1), testVars(2), testVars(3), testVars(4), ... + rho_h2o, rho_oil, mu_water, A_mu, B_mu, sigma_base, g, epsilon0, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, recycle_tol, oil_flow_Lh, ... + a_sigma, b_sigma, sigma_min_factor, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field); + end + + figure; + yyaxis left; + plot(vals, sal_vals,'-o','LineWidth',1.5); ylabel('Salinité (mg/L)'); + yyaxis right; + plot(vals, eff_vals,'-s','LineWidth',1.5); ylabel('Efficacité (%)'); + xlabel(var_names{v}); + legend('Salinité','Efficacité','Location','best'); grid on; + title(['Sensibilité : ' var_names{v} ' (gaoptimset)']); + end +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% FONCTIONS LOCALES (mêmes que précédemment) %%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function f = localObj(vars, ... + rho_h2o, rho_oil, mu_water, A_mu, B_mu, sigma_base, g, epsilon0, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, recycle_tol, oil_flow_Lh, ... + a_sigma, b_sigma, sigma_min_factor, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field) + + [sal, eff] = simulateDesalting(vars(1), vars(2), vars(3), vars(4), ... + rho_h2o, rho_oil, mu_water, A_mu, B_mu, sigma_base, g, epsilon0, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, recycle_tol, oil_flow_Lh, ... + a_sigma, b_sigma, sigma_min_factor, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field); + + % objectif = (sal / 40) - (eff / 100) + f = (sal/40) - (eff/100); +end + + +function [sal_out, eff_out, dist_profiles] = simulateDesalting( ... + X, Qd, T, X_recycle, ... + rho_h2o, rho_oil, mu_water, A_mu, B_mu, sigma_base, g, epsilon0, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, recycle_tol, oil_flow_Lh, ... + a_sigma, b_sigma, sigma_min_factor, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field) + + % Viscosité de l'huile en fonction de T + T_K = T + 273.15; + mu_oil_local = A_mu * exp(-B_mu * T_K); + + % Désémulsifiant => sigma local + conc_demul = Qd / oil_flow_Lh; + sigma_local = sigma_base*(1 - a_sigma*(conc_demul)/(b_sigma + conc_demul)); + if sigma_local < sigma_min_factor*sigma_base + sigma_local = sigma_min_factor*sigma_base; + end + + % Eau de lavage + recycle + X_frac = X/100; + recycle_frac = X_recycle/100; + + water_feed_guess = initial_water_frac + X_frac; + s_out=0; e_out=0; p_out=struct(); w_left=0; + + for iter=1:20 + water_feed_frac = initial_water_frac + X_frac + ... + recycle_frac*(water_feed_guess - initial_water_frac - X_frac); + + [s_out,e_out,p_out,w_left] = oneSimulationStep(water_feed_frac, mu_oil_local, sigma_local, ... + rho_h2o, rho_oil, mu_water, ... + initial_water_frac, initial_salt_conc, ... + t_res_desalter, dt_desalter, d_sep, ... + DeltaP_valve, t_res_valve, v_char_valve, ... + We_crit, Ca_crit, K_break, K_coal, lambda, ... + d_classes, v_classes, N_classes, E_field, ... + g, epsilon0); + + water_purged = water_feed_frac - w_left; + new_feed = initial_water_frac + X_frac + recycle_frac*water_purged; + if abs(new_feed - water_feed_frac) coalescence + N_afterDesal1 = solveCoalescence(N_afterMix, mu_oil_local, sigma_local, ... + rho_h2o, rho_oil, mu_water, t_res_desalter, dt_desalter, ... + d_sep, K_coal, E_field, g, epsilon0, d_classes, v_classes, N_classes); + + for i=1:N_classes + if d_classes(i)>= d_sep + N_afterDesal1(i)=0; + end + end + + %% 2e dessaleur => coalescence + N_afterDesal2 = solveCoalescence(N_afterDesal1, mu_oil_local, sigma_local, ... + rho_h2o, rho_oil, mu_water, t_res_desalter, dt_desalter, ... + d_sep, K_coal, E_field, g, epsilon0, d_classes, v_classes, N_classes); + + for i=1:N_classes + if d_classes(i)>= d_sep + N_afterDesal2(i)=0; + end + end + + water_left_frac = sum(N_afterDesal2.*v_classes); + + salt_conc_eau = initial_salt_conc / water_feed_frac; + sal_out = salt_conc_eau * water_left_frac; + eff_out = (initial_salt_conc - sal_out)/initial_salt_conc*100; + if eff_out>100 + eff_out=100; + end + + profiles.initial.N = N0; + profiles.initial.water_frac = tot_water_vol; + profiles.afterMix.N = N_afterMix; + profiles.afterMix.water_frac = sum(N_afterMix.*v_classes); + profiles.afterDesal1.N = N_afterDesal1; + profiles.afterDesal1.water_frac = sum(N_afterDesal1.*v_classes); + profiles.afterDesal2.N = N_afterDesal2; + profiles.afterDesal2.water_frac = water_left_frac; +end + + +function N_out = solveBreakagePBE(N_in, Ca_valve, We_valve, mu_oil_local, ... + DeltaP_valve, t_res_valve, rho_oil, K_break, We_crit, Ca_crit, lambda, ... + d_classes, v_classes, N_classes) + + dtv = t_res_valve/10; + N_current = N_in; + + for step=1:10 + births = zeros(1,N_classes); + deaths = zeros(1,N_classes); + for k=1:N_classes + if N_current(k)<=0, continue; end + g_k = local_g(d_classes(k), Ca_valve, We_valve, mu_oil_local, ... + DeltaP_valve, t_res_valve, rho_oil, K_break, We_crit, Ca_crit, lambda); + if g_k<=0, continue; end + + nb_breaks = g_k*N_current(k)*dtv; + nb_breaks = min(nb_breaks, N_current(k)); + if nb_breaks<=0, continue; end + + v0 = v_classes(k); + for i=1:(k-1) + v_mid = v_classes(i); + f_val = breakage_distribution(v_mid, v0); + vol_to_i = f_val*v0; + num_frag = vol_to_i / v_mid; + births(i) = births(i) + nb_breaks*num_frag; + end + deaths(k) = deaths(k) + nb_breaks; + end + + N_current = N_current + births - deaths; + N_current(N_current<0)=0; + end + + N_out = N_current; +end + + +function val_g = local_g(dk, Ca_valve, We_valve, mu_oil_local, ... + DeltaP_valve, t_res_valve, rho_oil, K_break, We_crit, Ca_crit, lambda) + + if Ca_valve<=Ca_crit + val_g=0; + return; + end + + val_g = K_break * 63.927 * ((We_valve^(11/5))/We_crit) * ... + sqrt((DeltaP_valve/(t_res_valve*rho_oil))/mu_oil_local) * ... + (Ca_valve^2.2) * (dk/(2*lambda))^(4/5); +end + + +function f_val = breakage_distribution(v, v0) + % Distribution gaussienne (Coulaloglou-Tavlarides) + f_val = 2.4./v0 .* exp(-4.5.*( (2.*v - v0).^2 ./ v0.^2 )); +end + + +function N_out = solveCoalescence(N_in, mu_oil_local, sigma_local, ... + rho_h2o, rho_oil, mu_water, t_res_desalter, dt_desalter, ... + d_sep, K_coal, E_field, g, epsilon0, d_classes, v_classes, N_classes) + + N_current = N_in; + steps = floor(t_res_desalter/dt_desalter); + + for s=1:steps + births = zeros(1,N_classes); + deaths = zeros(1,N_classes); + vol_water = sum(N_current.*v_classes); + delta = min(vol_water,0.999); + mu_ratio = mu_water/mu_oil_local; + + for i=1:N_classes + for j=i:N_classes + if N_current(i)<=0 || N_current(j)<=0, continue; end + + beta_ij = K_coal * pi*(d_classes(i)+d_classes(j))^2 * ... + ((mu_ratio+1)/((3*mu_ratio+2)*mu_oil_local)*abs(rho_h2o-rho_oil)*(d_classes(i)/2)^2*(1-delta^2)*g); + + eff_coal = collision_eff_withEF(delta, d_classes(i), E_field, ... + epsilon0, rho_h2o, rho_oil, g); + beta_ij = beta_ij * eff_coal; + + col = beta_ij * N_current(i)*N_current(j)*dt_desalter; + if i==j + col=0.5*col; + end + col=min(col, min(N_current(i),N_current(j))); + if col>0 + v_new = v_classes(i)+v_classes(j); + [~, k] = min(abs(v_classes - v_new)); + births(k)= births(k)+ col; + deaths(i)= deaths(i)+ col; + deaths(j)= deaths(j)+ col; + end + end + end + + N_current = N_current + births - deaths; + N_current(N_current<0)=0; + end + + N_out = N_current; +end + + +function eff_val = collision_eff_withEF(delta, d_i, E, epsilon0, rho_h2o, rho_oil, g) + eff_val = 0.45 * ( ... + (2.*delta.*abs(rho_h2o - rho_oil).*(1-delta).*g.*(d_i/2).^3) ./ ... + (3.*epsilon0.*(1+delta).^2.*E.^2) ... + ).^(-0.55); +end