% main file for void growth simulations during solid-state nanofoaming 
% see PRSA paper for theoretical framework (and assumptions)
% suggestions and questions: frederikvanloock05@outlook.com

% cd('C:\Users\fvanloock\Dropbox\Data PRSA')
	
	
	
clear all	
close all	
clc	
	
global T_g_0 C1 C2 eps_ref R_g gamma  C_0 c_0 R_0 b_0 x_mesh counter pg_0 pa pressure_vec ...	
    counter_smooth rho_L t_end temp_up temp_low c_vec  eta_g tau t_vec...	
    rho_0 M_w rho_vec temp_vec tg_vec ...	
    v Q eps_dot_0 E_R  eps_dot_R alpha_R n sigma_plastic_vec stress_vec sigma_rubber  mechanism_vec K_H eos_index KH_index ...	
    thickness_vec t_min thickness_check ductility_check T_Tg_vec T_g_0 stress_profile_vec f_int_vec	
	
tic % start timer	

	
% Simulation parameters and initial conditions	
N_0 = 2E20; % Nucleation density in cells/m3	
pg_0 = 31*10^6; % Initial gass pressure in Pa	
pa = 0.1013*10^6;% 1*1013.25*10^2; % Ambient pressure in Pa	
	
% Material propertiesc_0	
M_w = 44*10^-3; % Molecular weight CO2 kg/mol	
R_g = 8.31441; % Universal gas constant J/kg*K	
	
gamma =  0.028; %Surface tension N/m	
rho_L = 1190; %Density solid in kg/m3	

	
R_0_th = 1*(2*gamma./(pg_0-pa)); % min. initial cell radius to check	
R_0 = 10.5*10^-9; % initial cell radius in m  (high Mw 5nm with 20E20 Nd  low w 10.5nm with Nd 2E20)	
K_H = 2.12*10^-4; 

K_H2 = K_H*M_w/rho_L  % -/Pa	
	
% Approx initial concentration in weight percent 
c_0 = pg_0*K_H2 % 	
	
por_0 = (3/(R_0^3*4*3.1415*N_0)+1)^-1	
	
b_0 = (R_0^3/por_0)^(1/3);	

K_H_0 = K_H; % mol/(m3Pa) 


% Compute initial concentration
C_0 = pg_0*K_H_0; %(mol/m3)

% gas free Tg	
T_g_0 = 114.5+273; % gas-free Tg in K	

	

% simulation time period	
t_0 =   0;	
t_end = 600;	
t_n = 1e2;	
	
t_span = linspace(t_0,t_end,t_n);	
delta_t_global = t_span(2)-t_span(1);	
	
% define mesh 	
x_0 = R_0;	
x_end = b_0;	
x_n = 1e2;	
x_mesh = logspace(log10(x_0),log10(x_end),x_n); 	
	
a_0  = R_0;	

t_par_vec = [27 39 60 73]; % deg	(match with measured foaming bath temperatures as void growth is highly sensitive to Tf)

tau_vec = ones(1,length(t_par_vec))*20; % see temp profile

t_par_vec =    (t_par_vec + 273)./(T_g_0)	

for i = 1:length(t_par_vec)
  
  
tau = tau_vec(i);
temp_low = 273-15; % measured sample temperature prior to foaming (match with measured temperature, can be depending on pol grade)
temp_up = t_par_vec(i)*T_g_0;


opts = odeset('RelTol',1e-11,'AbsTol',1e-16,'Stats','on');; %might need to play with this to get converging solution at higher temperatures
counter = 0;
[t,x]=ode15s(@odesolver_uniformconcentration_new, t_span,[a_0],opts);

% extract solutions to plot, for example:
a_vec = x(:,1);
b_vec = (a_vec.^3+b_0^3-R_0^3).^(1/3);
f_vec = a_vec.^3./b_vec.^3;
plot(t_span,f_vec)
hold on
xlim([0 700])
ylim([0 1])
end