Info

This question is closed. Reopen it to edit or answer.

Can someone help check my codes if I followed the right approach for my problem?

1 view (last 30 days)
Dear Experts,
It will really help me if you can look at the complete code I am trying to solve. The output of the model and experimental results should be similar. But there is a deviation between the two when I check and this make me doubt if my approach is correct. Therefore, your help will mean a lot to me.
The final aim of my code is to converge the difference between the calculated and guess value. But the calculated value should use the guess value first to get the calculated value. My codes are shown below one by one. I will start with the functions I called for the filnal solvers on number 6 and 7.
  1. Function ''Hindrance factor''
function [Kc, Kd, DP, phi] = HindranceFactor(rs, rp, D_inf)
lambda = rs./rp;
phi=(1-lambda);
Kc = (2 - phi).*(1 + 0.054.*lambda - 0.988.*(lambda.^2) + 0.441.*(lambda.^3));
Kd = 1 - 2.3.*lambda + 1.154.*(lambda.^2) + 0.224.*(lambda.^3);
DP = Kd .* D_inf;
end
2. Differential equations to be solved using Ode15i
function dxdz = MS_AA_imp(z,x,xp, x_p, V, Z)
TempK = 25 + 273.15; % Temp in Kelvin
Faraday = 96485.336; % Faraday constant [C/mol]
R_id = 8.314; % R constant [J/(K mol)]
VM_1 = 1.8e-5; % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol
VM_2 = 1.57e-5; % Na molar volume [m3/mol] * Calculated from rS
VM_3 = 4.47e-6; % Cl molar volume [m3/mol] * calcualted form rS
VM_4 = 1.8e-5; % Water molar volume [m3/mol]
D_24 = 1.33E-9; % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9; % D34 value under diluted conditions [m^2/s]
h=7;
hydration= x(1)/((x(1))*(1-x(1)*h)*(R_id*TempK)); % new hydration formula
CT = x(1) / VM_1 + x(2)/ VM_2 + x(3) / VM_3 + (1-sum(x)) / VM_4;
D_14 =10e-9 *(0.665-0.393*x(1))/(1+x(1)* (h/(1-x(1)* h))); % new D_14 formula
IStr = 0.5* (Z(2)^2 * x(2)/(VM_2 *1000) + Z(3)^2 * x(3)/(VM_3 *1000));
D_23 = ((D_24 + D_34)/2 * IStr ^(0.55) / (abs((Z(2)) * (Z(3)))^(2.3)));
J_1 = x_p(1) * V;
J_2 = x_p(2) * V;
J_3 = x_p(3) * V;
J_4 = (1-sum(x_p)) * V;
%equations to be solved
dxdz = [hydration * xp(1) + x(1) * Z(1) * Faraday / (R_id .* TempK)*xp(5) - (J_4 * x(1) - J_1* x(4))/(D_14 * CT);
xp(2) + x(2) * Z(2) * Faraday*xp(5) / ( R_id * TempK) - (J_3 * x(2) - J_2 * x(3))/ D_23 - (J_4 * x(2) - J_2 * x(4))/ (D_24 *CT);
xp(3) + x(3) * Z(3) * Faraday*xp(5) / ( R_id * TempK) - (J_2 * x(3) - J_3 * x(2))/ D_23 - (J_4 * x(3) - J_3 * x(4))/ (D_34* CT); %D_32V replaced with D_
x(1) + x(2) + x(3) + x(4) - 1;
Z(1) * x(1) + Z(2)* x(2) + Z(3) * x(3)];
end
3. function ''partitionFirst". To be called (solved) by fsolve
function F = partitionFirst(C_entrance_MS,C_surf,Z,phi)
k=C_entrance_MS(1); %concentration inside the membrane inlet
l=C_entrance_MS(2);
m=C_entrance_MS(3);
X_mem=C_entrance_MS(4);
VM_1 = 1.8e-5; % lysine molar volume [m3/mol]
VM_2 = 1.57e-5; % Na molar volume [m3/mol] * Calculated from rS
VM_3 = 4.47e-6; % Cl molar volume [m3/mol] * calcualted form rS
F(1)= Z(1)*k + X_mem + Z(2)*phi(2)*C_surf(2)*(k/(phi(1)*C_surf(1)))^(Z(2)/Z(1))...
+ Z(3)*phi(3)*C_surf(3)*(k/(phi(1)*C_surf(1)))^(Z(3)/Z(1));
F(2)= -l + phi(2)*C_surf(2)*(k/(phi(1)*C_surf(1)))^(Z(2)/Z(1));
F(3)= -m + phi(3)*C_surf(3)*(k/(phi(1)*C_surf(1)))^(Z(3)/Z(1));
F(4)= k*Z(1) + Z(2)*l + Z(3)*m + X_mem;
end
4. Function ''NernstPlanck_vol''. To be called (solved) by fsolve
function dydz = NernstPlanck_Vol(Co,C_entracne_MS, Z, Kc, DP) %linearized pore equations( derivated--can be found on other files)
xo_1 = Co(1); %concentration at the outlet of the membrane pore(end of membrane concentration)
xo_2 = Co(2);
xo_3 = Co(3);
TempK = 25 + 273.15; % Temp in Kelvin
Faraday = 96485.336; % Faraday constant [C/mol]
R_id = 8.314; % R constant [J/(K mol)]
VM_1 = 1.8e-5; % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol
VM_2 = 1.57e-5; % Na molar volume [m3/mol] * Calculated from rS
VM_3 = 4.47e-6; % Cl molar volume [m3/mol] * calcualted form rS
VM_4 = 1.8e-5; % Water molar volume [m3/mol]
r_AA = 3.48e-9;
r_Na = 102e-12; % ion size
r_Cl = 181e-12; % ion size
p_7 = 0.6e-3; %Viscosity of pure water at 45C (Pa.s)(engineering toolbox.com)--denoted as p_7
deltaP = 5; %pressure
D_im = 1.2e-9; %diffusivitiy not yet worked out for all the solutes
Le = 1.65e-4; % effective membrane thickness
C(1:3)=C_entracne_MS;
x_1 = C_entracne_MS(1); %correct by rreplaceing C() with C_entrnce()
x_2 = C_entracne_MS(2);
x_3 = C_entracne_MS(3);% here was s added by abel and I think it is an error
x_ave=(Co(1:3)+C(1:3))/2;
%corrected equtions
eq4 = ((x_ave(1)*Z(1)*r_AA^2*deltaP/(8*p_7*Le)*(Kc(1)-1)*R_id*TempK)/DP(1) - x_ave(1)*Z(1)*VM_1*deltaP/Le)/(Faraday*(Z(1)^2*x_ave(1)+Z(2)^2*x_ave(2)+ Z(3)^2*x_ave(3))) +...
((x_ave(2)*Z(2)*r_Na^2*deltaP/(8*p_7*Le)*(Kc(2)-1)*R_id*TempK)/DP(2) - x_ave(2)*Z(2)*VM_2*deltaP/Le)/(Faraday*(Z(1)^2*x_ave(1)+Z(2)^2*x_ave(2)+ Z(3)^2*x_ave(3))) + ...
((x_ave(3)*Z(3)*r_Cl^2*deltaP/(8*p_7*Le)*(Kc(3)-1)*R_id*TempK)/DP(3) - x_ave(3)*Z(3)*VM_3*deltaP/Le)/(Faraday*(Z(1)^2*x_ave(1)+Z(2)^2*x_ave(2)+ Z(3)^2*x_ave(3)));
dydz(1) =x_1 - xo_1 + (x_ave(1)*(r_AA^2*deltaP/(8*p_7*Le)*(Kc(1)-1)*R_id*TempK)/DP(1) - x_ave(1)*VM_1*deltaP/Le - x_ave(1)*Z(1).*eq4)*Le; %assigning Elctric potential with number x(7)1
dydz(2) =x_2 - xo_2 + (x_ave(2)*(r_Na^2*deltaP/(8*p_7*Le)*(Kc(2)-1)*R_id*TempK)/DP(2) - x_ave(2)*VM_2*deltaP/Le - x_ave(2)*Z(2).*eq4)*Le; %assigning Elctric potential with number x(7)1
dydz(3) =x_3 - xo_3 + (x_ave(3)*(r_Cl^2*deltaP/(8*p_7*Le)*(Kc(3)-1)*R_id*TempK)/DP(3) - x_ave(3)*VM_2*deltaP/Le - x_ave(3)*Z(3).*eq4)*Le; %assigning Elctric potential with number x(7)1
end
5. function ''Partionsecond''. To be solved (called) by fsolve
function F=partitionsecond(Cp_calc,C_outlet,Z,phi)
kk=Cp_calc(1);
ll=Cp_calc(2);
mm=Cp_calc(3);
%bsed on electro neutrality in the permeate
F(1)= kk*Z(1) + Z(2)*C_outlet(2)/phi(2)*1/((C_outlet(1)/(phi(1)*kk))^(Z(2)/Z(1)))...
+ Z(3)*C_outlet(3)/phi(3)*1/((C_outlet(1)/(phi(1)*kk))^(Z(3)/Z(1))); %%% pationing equation for rearranging equating for Cp_cal
F(2)= -ll + C_outlet(2)/phi(2)*1/((C_outlet(1)/(phi(1)*kk))^(Z(2)/Z(1))); %partionig equation 2
F(3)= -mm + C_outlet(3)/phi(3)*1/((C_outlet(1)/(phi(1)*kk))^(Z(3)/Z(1))); %partionig equation 3
end
6. The above functions will be used to solve my problem. To do this, I used anonymous functions to use the output of one function as an input variable for the next function and so on. Below is the code for the cascaded functions.
function [sol] = find_CP_MSNP(Cp_guess, phi_start, Jv,Z, z_pol,phi,Kc,DP)
global VM_1 VM_2 VM_3 VM_4
global TempK Faraday R_id D_14_inf D_24 D_34
%ode15i for MS_AA_imp
Cp_guess_MS = [Cp_guess 1- sum(Cp_guess)];
[y0, yp0] = decic(@MS_AA_imp, 0, [phi_start 0], [0 0 0 0 0], [1 1 0 0 0], [],[],Cp_guess_MS, Jv, Z);
[zsol_MS, ysol_MS] = ode15i(@(z,x,xp)MS_AA_imp(z,x,xp ,Cp_guess_MS, Jv, Z), [0 z_pol], y0, yp0);
C_surf = ysol_MS(end,:)./[VM_1 VM_2 VM_3 VM_4 1];
%fsolve for partitionFirst
% define anonymous function 'f' for the partition function which accepts C_surf as an input parameter
zz=[0.00001 0.00001 0.00001 0.000001];
f = @(zz)partitionFirst(zz,C_surf(1:3),Z,phi);
C_entrance_pore = fsolve(f,zz);
C_entrance = C_entrance_pore.*[VM_1 VM_2 VM_3 1]; %to conver to mole fraction
%fsolve for NernstPlanck_Vol
% define anonymous function 'nernstPlanckAnonymous' for the NernstPlanck_Vol function which accepts C_entrance as a parameter
z=[0.00001 0.00001 0.00001];
nernstPlanckAnonymous = @(z)NernstPlanck_Vol(z,C_entrance(1:3),Z,Kc,DP);
C_end = fsolve(nernstPlanckAnonymous,z); %concentration on the end of membrane
% fsolve for partitionsecond uses the output of the above equation(C_end) as an
% input
secondpartitionAnonymous = @(z)partitionsecond(z, C_end(1:3),Z,phi);
Cp_calc = fsolve(secondpartitionAnonymous,z);
sol = Cp_calc(1:3) - Cp_guess(1:3);
end
7. The last solver to solve the function on number 6 is shown below. Here the solver solves the function on number 6 (''find_CP_MSNP'') sothat the difference between calculated value (Cp_calc(1:3)) and guess value (Cp_guess(1:3)) to converge or approaches zero.
N.B. the guess value Cp_guess(1:3) is the starting point for the functions on number 6 (''find_CP_MSNP'') and used on the last point.
global TempK Faraday R_id D_14_inf D_24 D_34
global VM_1 VM_2 VM_3 VM_4
TempK = 25 + 273.15; % Temp in Kelvin
VM_1 = 1.8e-5; % lysine molar volume [m3/mol]
VM_2 = 1.57e-5; % Na molar volume [m3/mol] * Calculated from rS
VM_3 = 4.47e-6; % Cl molar volume [m3/mol] * calcualted form rS
VM_4 = 1.8e-5; % Water molar volume [m3/mol]
D_14_inf = 6.1e-11; % D14 value under diluted conditions and IEP [m^2/s]
D_24 = 1.33E-9; % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9; % D34 value under diluted conditions [m^2/s]
Faraday = 96485.336; % Faraday constant [C/mol]
R_id = 8.314; % R constant [J/(K mol)]
z_pol = 1.23e-5;
r_AA = 3.48e-9;
r_Na = 102e-12; % ion size
r_Cl = 181e-12; % ion size
rs = [r_AA r_Na r_Cl];
D = [ D_14_inf D_24 D_34];
Vs = [VM_1 VM_2 VM_3];
rp = 9.75e-9;
Le = 1.65e-4; % effective membrane thickness
[Kc, Kd, DP, phi] = HindranceFactor(rs, rp, D)
C_ret_arb1 = [0.02 0.02 0.02]; % using average value of pH 7 0.15 M
C_ret_arb=[C_ret_arb1 1-sum(C_ret_arb1)];
Z_arb = [-1 1 -1];
Cp_guess_arb = [0.025 0.033 0.023]; % I changed the above line into the current line by removing C_ret
%% for arbitrary flux value
options = optimoptions('fsolve', 'Display','iter','algorithm', 'levenberg-marquardt','FunctionTolerance', 1e-4,'StepTolerance', 1e-4,'MaxFunctionEvaluation', 300);
Jv_arb = 5e-6;
[Cp_calc_arb, fval, exitflag] = fsolve(@(y)find_CP_MSNP( y, C_ret_arb,Jv_arb,Z_arb, z_pol,phi,Kc,DP), Cp_guess_arb(1:3),options); %from the line above
%have to find out why.
if exitflag < 1
iter = 0;
mod = 0.01;
while xor(exitflag < 1, iter >50)
Cp_guess_arb = C_ret_arb .* [mod 1 1 1]; % can change all the value, this is just for demonstration purposes
[Cp_calc_arb,fval, exitflag] = fsolve(@(y)find_CP_MSNP( y, C_ret_arb, Jv_arb, Z_arb, ...
z_pol,phi,Kc,DP), Cp_guess_arb(1:3), options);
iter = iter +1
mod = mod + 0.01;
end
end
Cp_calc_arb
I hope the codes are clear and you will be able to help me.
Thank you very much in advance.
  2 Comments
Jan
Jan on 7 Feb 2023
Edited: Jan on 7 Feb 2023
I do not see any chance, that a reader could find a bug in the code. The question is far too complicated to be answered seriously.
John D'Errico
John D'Errico on 7 Feb 2023
Exactly as Jan says. You cannot expect someone to spend the many hours of time to read that lengthy code and understand what it is that you WANT it to do, as opposed to what you may have WRITTEN it to do. Is there a bug in the code? How can we know, since we cannot possibly know if you wrote exactly what should be there? We have not even been told what the code is supposed to do.
So, learn to use the debugger. Spend the time, and trace thorough your code. Verify on each line that it did what you would expect it to do. Plot things if necessary. Think about what you see. Does it make sense? If so, then go on to the next line.
Seriously. What you are asking is far beyond the scope of Answers.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!