Skip to content
Snippets Groups Projects
Commit b12b97a7 authored by Halfwerk, Ruben's avatar Halfwerk, Ruben
Browse files

Upload New File

parent c41a8237
No related branches found
No related tags found
No related merge requests found
clear
clc
close all
%% loading data
files = ["Lactose 25% Exp. 3.mat";
"Lactose 25% Exp. 4.mat";...
"Lactose 25% Exp. 5.mat";...
"Lactose 30% Exp. 36.mat";...
"Lactose 30% Exp. 37.mat";...
"Lactose 40% Exp. 1.mat";...
"Lactose 40% Exp. 3.mat";...
"Lactose 40% Exp. 8.mat";...
"Lactose 40% Exp. 9.mat";...
"Lactose 40% Exp. 10.mat";...
"Lactose 40% Exp. 11.mat";...
"Lactose 40% Exp. 12.mat"]
for iii=1:length(files)
seed = [250 250 250 250 250 250 250 250 250 250 250 250];
concentration = [25.7 25.73 25.6 30.31 30.4 40.27 40 40 40 40 40 41.11]
nummer = [3 4 5 36 37 1 3 8 9 10 11 12];
SeedUsed = [1 1 1 1 1 1 1 1 1 1 1 1];
load(files(iii));
if Time1(1) == 0
con3 = w_lactose(1)*100;
con3 = concentration(iii)
else
con3 = concentration(iii);
end
num3 = seed(iii)
num5 = 1;
Seeded = SeedUsed(iii);
%% Setting up Constants
global M_lactose M_H2O M_lactose_H2O rho kv ks w_eq m_sol
%load('Lactose 25% Exp. 3.mat') % Loading CSD and concentration data, slightly modified
M_lactose = 342.3; % mol mass lactose in gram/mol
M_H2O = 18.01; % mol mass H2O in gram/mol
M_lactose_H2O = M_lactose+M_H2O; % mol mass Lactose.H2O in gram/mol
rho = 1590; % Densicleraty of crystal in kg/m3
kv = pi/6; % volume shape factor
ks = pi; % Surface shape factor
w_eq = 10.75/100; % Equilibrium concentration
m_sol = 1; % mass of total solution
%% Combining CSD and desaturation data/time
if Time1(1)==0
else
Time1=[0;Time1]
dCC =[(con3-10.75)/100;dCC]
end
if Time1(end)>=t_CSD(end)
[dt_lastCSD,i_lastCSD] = min(abs(t_CSD-Time1(end)));
Time0=Time1(1:find(Time1>=t_CSD(i_lastCSD),1));
else
i_lastCSD = find(t_CSD>Time1(end),1)
Time0=Time1;
end
t_lastCSD = t_CSD(i_lastCSD);
[dt_lastCSD,i_lastCSD] = min(abs(t_CSD-24))
% Reduce amount of data points (some measurements were taken every few seconds, overloading the model)
if length(Time1) < 200
dCCL = dCC;
else
Time1=Time1(1:10:end)
dCC =dCC(1:10:end)
dCCL = dCC;
end
TimeL = 0:0.01:Time1(end)
%% Time scale not the same functions to find TimeL overlapping with the emausrements
for i=1:length(t_CSD)
[val(i),vec_CSD(i)]=min(abs(TimeL-t_CSD(i)));
end
for i=1:length(Time1)
[va(i),vec(i)]=min(abs(TimeL-Time1(i)));
end
%% Determening moments and average sizes
dat1 =MeaninBandt;
dat2 = Fraction_t(1:100,num5);
M_4 = sum(Fraction_t(1:100,:).*(MeaninBandt));
M_3 = sum(Fraction_t(1:100,:));
M_2 = sum(Fraction_t(1:100,:)./(MeaninBandt));
M_1 = sum(Fraction_t(1:100,:)./(MeaninBandt.^2));
M_0 = sum(Fraction_t(1:100,:)./(MeaninBandt.^3));
d43_meas = M_4./M_3;
d32_meas = M_3./M_2;
d10_meas = M_1./M_0;
%% Creating new size classes removing zero values
global UpperBoundN datAll
datAll =10.^linspace(log10(dat1(18)),log10((dat1(95))),101)'; % 100 size bins
LowerboundN = datAll(1:(end-1));
UpperBoundN = datAll(2:(end));
BandwithN = (UpperBoundN - LowerboundN); % Band size
MeaninBandtN = (UpperBoundN + LowerboundN)/2; % Average size
LogMeaninBandt = log10(MeaninBandtN);
datF = MeaninBandtN;
%% For calculating the start number distribution, can be seed or the second point
dCC_point2 = interp1(Time1,dCC,t_CSD(num5));
mseed = num3*10^-6;
W_initial = con3/100;
w_point2=dCC_point2+w_eq;
m_con_point2 =(m_sol*(1-W_initial)-m_sol*W_initial*M_H2O/M_lactose)./(1-w_point2-w_point2*M_H2O/M_lactose); % Mass of concentrate, acounting for hydrate (in kg)
if Seeded == 1
mt_point2 = mseed
else
mt_point2 = (m_sol-m_con_point2)+mseed; % Mass of lactose monohydrate crystals (in kg)
end
%% Creating cummalative distribution rahter then normal CSD (data from laser difraction measurements)
cumCSD = cumsum(Fraction_t./sum(Fraction_t));
dat4 =cumCSD(1:end,num5);
dat4F = zeros(length(dat1),1);
dat4F(1) = dat4(1);
dat4F(2:end,1) = dat4(2:end)-dat4(1:end-1);
Cum_start = interp1(UpperBound,cumCSD(1:100,num5),datAll,'spline');
Cum_end = interp1(UpperBound,cumCSD(1:100,i_lastCSD),datAll,'spline');
Cum_total = interp1(UpperBound,cumCSD(1:100,:),datAll,'spline');
%% Creating new CSD values at every point of new size classes
NewdisFraction= Cum_start(2:end)-Cum_start(1:end-1);
NewdisFractionend = Cum_end(2:end)-Cum_end(1:end-1);
for i = 1:length(NewdisFraction)
if NewdisFraction(i)<10^-8
NewdisFraction(i)=0;
end
if NewdisFractionend(i)<10^-8
NewdisFractionend(i)=0;
end
end
% Plotting CSD
figure
semilogx(dat1,Fraction_t(1:100,num5)./(log10(UpperBound)-log10(Lowerbound)),'o','DisplayName','Start measured')
hold on
semilogx(dat1,Fraction_t(1:100,i_lastCSD)./(log10(UpperBound)-log10(Lowerbound)),'o','DisplayName','final measured')
semilogx(MeaninBandtN,NewdisFraction./(log10(UpperBoundN)-log10(LowerboundN)),'+' ,'DisplayName','Gaussian fit start')
semilogx(MeaninBandtN,NewdisFractionend./(log10(UpperBoundN)-log10(LowerboundN)),'+' ,'DisplayName','Gaussian fit start')
semilogx(MeaninBandtN,NewdisFractionend./(log10(UpperBoundN)-log10(LowerboundN)) , 'DisplayName','Gaussian fit start')
semilogx(MeaninBandtN,NewdisFraction./(log10(UpperBoundN)-log10(LowerboundN)),'DisplayName','Gaussian fit start')
ylabel('Volume density(in %)')
xlabel('Size (in um)')
legend()
hold off
%%
figure
hold on
plot(log10(dat1),Fraction_t(1:100,num5),'*','DisplayName','Laser diffraction')
bar(log10(MeaninBandtN),NewdisFraction, 'DisplayName','Cummalative transfered')
ylabel('Volume density(in %)')
xlabel('Size (in um)')
legend()
hold off
%% Determening inital number distribution
N_t_00 = (mt_point2*NewdisFraction)./(rho.*kv*((MeaninBandtN).^3)); % number denstity per size class
N_total = sum(N_t_00);
A0 =ks*N_t_00.*MeaninBandtN.^2; % Surface area at the start
% conversion to water based concentration (in kg/kg water rather then kg/kg
% solution)
dCCL_H2O =(dCCL+w_eq)./(1-(dCCL+w_eq))-w_eq/(1-w_eq);
%% Determening the number dsitribution at the each time the CSD was measured
% Determening mass crystallized at each point
dCC_point = interp1(Time1,dCC,t_CSD);
w_point=dCC_point+w_eq;
m_con_point =(m_sol*(1-W_initial)-m_sol*W_initial*M_H2O/M_lactose)./(1-w_point-w_point*M_H2O/M_lactose); % Mass of concentrate, acounting for hydrate (in kg)
mt_point = (m_sol-m_con_point)+mseed; % Mass of lactose monohydrate crystals (in kg)
N_t = zeros(length(MeaninBandt),length(t_CSD));
% Determening number distributionat each point
for i =1:length(t_CSD)
N_t(:,i) = (mt_point(i).*Fraction_t(:,i))./(rho.*kv*((MeaninBandt).^3));
end
% Final number distribution
N_end = mt_point(i_lastCSD).*NewdisFractionend./(rho.*kv*((MeaninBandtN).^3));
N_end1 = N_t(:,i_lastCSD)
N_end2 = (mt_point(i_lastCSD).*Fraction_t(:,i_lastCSD))./(rho.*kv*((MeaninBandt).^3));
%% Parameter estimation by minimization of error
% Initial gues
theta0 = [ 2.51 0.7 14.6]
Fraction_final =1;
Fraction_t_int = 1;
f = @(theta)findmin(theta,TimeL(1:end),dCCL,N_t_00,A0,W_initial,Fraction_t_int,MeaninBandtN,BandwithN,dCCL_H2O,NewdisFractionend,N_end,Cum_total,vec,vec_CSD,mseed); % Function to calculate the difference between measured and moddeled
% boundary conditions
A = [];
b = [];
Aeq = [];
beq = [];
lb = [2,0,10];
ub = [10,1,15];
%options0 = optimoptions('fmincon','DiffMinChange',0.1,'OptimalityTolerance',0,'StepTolerance', 0,'StepTolerance', 0,'MaxIterations',20000,'MaxFunctionEvaluations',50000);
options = optimoptions('fmincon','MaxFunctionEvaluations',1000);
theta = fmincon(f,theta0,A,b,Aeq,beq,lb,ub,[],options);
%% Running the model with the minimized value of theta
[n,t,dCt_model,mt_model,mt_t_model,B,G,Final_fraction_model,d43,d10,d50,d90,ninit1,dCt_model_H2O,M_4,M_3,w_alpha,w_beta,dC_alpha,Final_fraction_model_cum]=Mass_balance_FVM(theta,TimeL,N_t_00,A0,W_initial,MeaninBandtN,BandwithN,mseed);
% Determening percentiles of the model, moments and average sizes
for i= 1:length(TimeL)
d10(i) = UpperBoundN(find(Final_fraction_model_cum(i,:)>0.1,1));
d50(i) = UpperBoundN(find(Final_fraction_model_cum(i,:)>0.5,1));
d90(i) = UpperBoundN(find(Final_fraction_model_cum(i,:)>0.9,1));
end
%% Statistic al analysis
%determening +- 10% theta
frac =2;
for i = 1:length(theta)
theta_90(i,:) = theta;
theta_110(i,:) = theta;
theta_90(i,i)=theta(i)*(1-frac/100);
theta_110(i,i)=theta(i)*(1+frac/100);
[n90,~,~,~,~,~,~,~,~,~,~,~,~,dCt_model_H2O90,~,~,~,~,~,Final_fraction_model_cum98]=Mass_balance_FVM(theta_90(i,:),TimeL,N_t_00,A0,W_initial,MeaninBandtN,BandwithN,mseed);
[n110,~,~,~,~,~,~,~,~,~,~,~,~,dCt_model_H2O110,~,~,~,~,~,Final_fraction_model_cum102]=Mass_balance_FVM(theta_110(i,:),TimeL,N_t_00,A0,W_initial,MeaninBandtN,BandwithN,mseed);
n_90(:,i)= n90(end,:)';
n_110(:,i)= n110(end,:)';
dCt_model_H2O_90(:,i) =dCt_model_H2O90(vec);
dCt_model_H2O_110(:,i) =dCt_model_H2O110(vec);
for j=1:length(vec_CSD)
if j == 1
Final_fraction_model_cum_98(j:100,i) = Final_fraction_model_cum98(vec(j),:)';
Final_fraction_model_cum_102(j:100,i) = Final_fraction_model_cum102(vec(j),:)';
Final_fraction_model_cum_98_norm(j:100,i) = Final_fraction_model_cum98(vec(j),:)'./mean(Cum_total(1:100,j));
Final_fraction_model_cum_102_norm(j:100,i) = Final_fraction_model_cum102(vec(j),:)'./mean(Cum_total(1:100,j));
else
Final_fraction_model_cum_98((j-1)*100+1:(j)*100,i) = Final_fraction_model_cum98(vec(j),:)';
Final_fraction_model_cum_102((j-1)*100+1:(j)*100,i) = Final_fraction_model_cum102(vec(j),:)';
Final_fraction_model_cum_98_norm((j-1)*100+1:(j)*100,i) = Final_fraction_model_cum98(vec(j),:)'./mean(Cum_total(1:100,j));
Final_fraction_model_cum_102_norm((j-1)*100+1:(j)*100,i) = Final_fraction_model_cum102(vec(j),:)'./mean(Cum_total(1:100,j));
end
end
clear n90 n110 dCt_model_H2O90 dCt_model_H2O110 Final_fraction_model_cum98 Final_fraction_model_cum102
end
%%
J1 = [((dCt_model_H2O_110-dCt_model_H2O_90)./mean(dCCL_H2O));(Final_fraction_model_cum_102_norm-Final_fraction_model_cum_98_norm)]
% Determening Jacobian matrix: Y110-Y90/(0.2*Y)
for i =1:length(theta)
J(:,i) = J1(:,i)./((frac/100*2)*theta(i));
end
%
figure
plot(TimeL,dCt_model_H2O)
hold on
%plot(Time0,dCt_model_H2O_110,'--')
%plot(Time0,dCt_model_H2O_90,':')
hold off
figure
semilogx(MeaninBandtN,n(end,:))
hold on
semilogx(MeaninBandtN,n_90,':')
semilogx(MeaninBandtN,n_110,'--')
hold off
%tsteps
figure
semilogx(MeaninBandtN,Final_fraction_model_cum(end,:))
hold on
semilogx(MeaninBandtN,Final_fraction_model_cum_98((end-99):end),':')
semilogx(MeaninBandtN,Final_fraction_model_cum_102((end-99):end),'--')
hold off
%
epsilon = [(dCCL_H2O-dCt_model_H2O(vec))./(mean(dCCL_H2O));(Cum_end(1:end-1)-Final_fraction_model_cum(end,1:end)')./(mean(Cum_end(1:end-1)))];
sigma2= (1/(length(epsilon)-length(theta)))*transpose(epsilon)*epsilon;
variance = sigma2*((transpose(J)*J)^-1);
sigma = sqrt(diag(variance))
tvalue = tinv(0.975,length(epsilon));
confidence95 = sigma.*tvalue
correlation=variance./(sigma*sigma');
[v,d] = eig(variance)
%% Calculating final number and size distributions at t = tend
Final=ninit1(vec_CSD(i_lastCSD),:)'.*(rho.*kv*(MeaninBandtN(1:end).^3)); % Calculate volume per class at t = tend
Final_fraction=Final/sum(Final); % Calculate volume fraction per class at t= tend
q3_n_final = ninit1(vec_CSD(i_lastCSD),:)./sum(ninit1(vec_CSD(i_lastCSD),:))./(log10(UpperBoundN(1:end)')-log10(LowerboundN(1:end)')); % Calculate number distributionat t= tend
Q3_start=NewdisFraction./(log10(UpperBoundN(1:end))-log10(LowerboundN(1:end))); % Calculate volume distribution at t = tend
Q3_final=Final_fraction./(log10(UpperBoundN(1:end))-log10(LowerboundN(1:end))); % Calculate volume distribution at t = tend
%% Drawing new figures
% Input for saving
con3 = concentration(iii);
% experimetn number
num4 = nummer(iii);
%% Alpha and beta lactose concentration over time
figure
plot(Time1,(dCC+10.75/100)./(1-(dCC+10.75/100)),'o','DisplayName','Lactose concentration Measured')
hold on
%plot(TimeL,(dCt_model+10.75/100)./(1-(dCt_model+10.75/100)),'DisplayName','Lactose concentration Model')
plot(TimeL,w_alpha,'DisplayName','Alpha lactose Concentration')
plot(TimeL,w_beta,'DisplayName','Beta lactose Concentration')
plot(TimeL,w_beta+w_alpha,'DisplayName','Alpha+Beta lactose Concentration')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' Desaturation with alpha_beta'+'.fig')
hold off
%% Different figures
figure
hold on
plot(TimeL,dCt_model_H2O)
%plot(Time0,dCCL_H2O,'o')
legend('model','Measured')
xlabel('Time elapsed (in hours)')
ylabel('Supersaturation (in kg/kg solution)')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' Desaturation'+'.fig')
hold off
%
% Final number distrbution
figure
semilogx(MeaninBandtN*10^6,N_t_00(1:end),'*','DisplayName','CSD start')
hold on
semilogx(MeaninBandtN*10^6,n(end,1:end),'o','DisplayName','Number distribution model')
semilogx(MeaninBandtN*10^6,N_end,'o','DisplayName','Number distribution meas')
xlabel('Crystal size in µm')
ylabel('frequency')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' Number Distribution'+'.fig')
hold off
% Final number distrbution
figure
semilogx(MeaninBandtN*10^6,N_t_00/sum(N_t_00),'*','DisplayName','CSD start')
hold on
semilogx(MeaninBandtN*10^6,n(end,1:end)/sum(n(end,1:end)),'o','DisplayName','Number distribution model')
semilogx(MeaninBandtN*10^6,N_end/sum(N_end),'o','DisplayName','Number distribution meas')
xlabel('Crystal size in µm')
ylabel('frequency')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' Number Distribution'+'.fig')
hold off
% Final volume distrbution
figure
semilogx(MeaninBandt*10^6,Fraction_t(1:100,num5)./(log10(UpperBound(1:end))-log10(Lowerbound(1:end))),'b*','DisplayName','start measured')
hold on
semilogx(MeaninBandtN*10^6,Q3_final,'DisplayName','CSD model')
semilogx(MeaninBandt*10^6,Fraction_t(1:100,i_lastCSD)./(log10(UpperBound(1:end))-log10(Lowerbound(1:end))),'r*','DisplayName','Final measured')
semilogx(MeaninBandtN*10^6,NewdisFractionend./(log10(UpperBoundN(1:end))-log10(LowerboundN(1:end))),'r','DisplayName','CSD end gaussian')
xlabel('Crystal size in µm')
ylabel('Fraction')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' Volume distribution'+'.fig')
hold off
% Final cummalative distrbution
figure
semilogx(UpperBoundN*10^6,Cum_start(2:end),'b*','DisplayName','start measured')
hold on
semilogx(UpperBoundN*10^6,Cum_end(2:end),'r*','DisplayName','Final measured')
semilogx(UpperBoundN*10^6,Final_fraction_model_cum(end,:),'DisplayName','CSD model')
xlabel('Crystal size in µm')
ylabel('Fraction')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' Cummulative volume distribution'+'.fig')
hold off
% d43 over time
figure
plot(t_CSD,d43_meas,'*','DisplayName','d43 measured')
hold on
plot(TimeL,d43,'DisplayName','d43 model')
%plot(TimeL,d43_measIntp,'DisplayName','d43 interpolated')
xlabel('time elapsed (in hours)')
ylabel('Volume averaged particle size (in um)')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' d43'+'.fig')
hold off
% d10/d50/d90 over time
figure
plot(TimeL,d10*10^6,'b','DisplayName','d10 model')
hold on
plot(t_CSD,d10_meas*10^6,'b*','DisplayName','d10 measured')
plot(TimeL,d50*10^6,'r','DisplayName','d50 model')
plot(t_CSD,d50_meas*10^6,'r*','DisplayName','d50 measured')
plot(TimeL,d90*10^6,'c','DisplayName','d90 model')
plot(t_CSD,d90_meas*10^6,'c*','DisplayName','d90 measured')
legend()
savefig('Lactose '+string(con3)+'% Exp. '+string(num4)+' d10_90_d50'+'.fig')
hold off
% % Exporting data to excel and saving it to a mat file
% creating large matrix PSD
%% Saving to table
filename = 'Data model '+string(con3)+'% Exp. '+string(num4)+'.xlsx'
z1 = zeros(length(TimeL),1);
Parameter_T1 = repmat({'Lactose concentration model'},length(TimeL),1)
Parameter_T2 = repmat({date},length(TimeL),1)
T1 = table(W_initial*100+z1,con3*100+num4+z1,Parameter_T1,theta(1)+z1,confidence95(1)+z1,2.162+z1,z1,theta(2)+z1,confidence95(2)+z1,0+z1,0+z1,theta(3)+z1,confidence95(3)+z1,1.5+z1,0+z1,TimeL',dCt_model,dC_alpha,mt_model,G,B,Parameter_T2);
T1.Properties.VariableNames = {'Start Concentration','Exp.','Parameter','theta 1','95% confidence 1','theta 2','95% confidence 2','theta 3','95% confidence 3','theta 4','95% confidence 4','theta 5','95% confidence 5','theta 6','95% confidence 6','Time','Supersaturation model','Supersaturation alpha lactose','suspension density','Growth rate model','Nucleation rate model','Model version'};
writetable(T1,filename,'Sheet','Concentration data model');
%% extracting psd data
Q_meas = NewdisFractionend./(log10(UpperBoundN(1:end))-log10(LowerboundN(1:end)));
Q3_final;
% Making table
filename = 'CSD model '+string(con3)+'% Exp. '+string(num4)+'.xlsx'
z1 = zeros(length(Q_meas),1);
Parameter_T3 = repmat({'CSD'},length(z1),1);
Parameter_T4 = repmat({date},length(z1),1);
T2 = table(W_initial*100+z1,con3*100+num4+z1,Parameter_T3, t_lastCSD+z1,MeaninBandtN,Q_meas,Q3_final,Parameter_T4);
T2.Properties.VariableNames = {'Start Concentration','Exp.','Parameter','Time','Average size','Measured','Fraction','Model version'};
writetable(T2,filename,'Sheet','CSD data model');
%% saving all data to matfile
save('Results model '+string(con3)+'% Exp. '+string(num4)+'.mat')
%%
clearvars -except iii files nummer seed concentration
close all
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment