CavitationPID.m

global error radius_d dt time_array tex MaxRadExp;

error = 0;
        
% Reading experimental data
% data2 = bubbleRadiuscell{i,j}';
% data1 = bubbleTimeVectorcell{i,j}';


% Find the maximum radius and the corresponding time of maximum occurence
[MaxRadExp,I_MaxRadExp] = max(data2);


tex = data1(I_MaxRadExp:end)-data1(I_MaxRadExp);
rex = data2(I_MaxRadExp:end);

radius_d = rex;
dt = 0.00001;
time_array = 0:dt:dt*(numel(tex)-1);

%% CONSTANTS %%
R_0 = MaxRadExp;
% R_0 = 1.03813e-03;
mu = 1e-3;
n = 3;
D = 1.76e-09;
K_B = 6.737e09;
alpha_M = .8544;
R_v = 461.52;
Na = 6.02e23;
Rg = 8.314;
T_infty = 293.15; 
P_infty = 1.01e05; 
Pinf = 1.01e05; 
T_0 = T_infty;
TB_0 = T_infty;
Vol_0 = 4/3*pi*R_0^3;
P_0 = 50000; %%%%%%%%%%%%%%%%%%%
P_water = 0.83775*55000; %%%%%%%%%%%%%%%%%%%
n_t_0 = P_0*Vol_0*Na/(Rg*T_0);
n_water_0 = P_water / P_0 * n_t_0;
n_air_0 = n_t_0 - n_water_0;

% Van der Waals Constants
a_air_VdW = 0.1402;
a_air_water_VdW = 0.278594;
a_water_VdW = 0.5536;
b_air_VdW = 3.753e-05;
b_air_water_VdW = 3.38882e-05;
b_water_VdW = 3.049e-05;

% Initial Conditions
R_0 = R_0;
Rdot_0 = 0;
mdot_0 = 0;
n_t_0 = P_0*Vol_0*Na/(Rg*T_0);
n_water_0 = P_water / P_0 * n_t_0;
n_air_0 = n_t_0 - n_water_0;
a_0_VdW = a_air_VdW*(n_air_0/n_t_0)^2 + 2*a_air_water_VdW*( n_air_0* n_water_0/n_t_0^2) + a_water_VdW*(n_water_0/n_t_0)^2;
b_0_VdW = b_air_VdW*(n_air_0/n_t_0)^2 + 2*b_air_water_VdW*( n_air_0* n_water_0/n_t_0^2) + b_water_VdW*(n_water_0/n_t_0)^2;
T_0 = T_infty;
P_0 = 50000; 
Vol_0 = 4/3*pi*R_0^3;
TB_0 = T_infty;

% Initialization
x0 = [n_air_0; n_water_0; n_t_0; a_0_VdW; b_0_VdW; T_0; Vol_0; TB_0; P_0; mdot_0; R_0; Rdot_0];


%% Main Section: Solving the Equations
tic
options = odeset('AbsTol',1e-3,'Stats','on');%'Reltol',1e-9,
[t_ode, x] = ode45(@ODESolver_PID, [time_array(1),time_array(end)], x0, options);

% Error calculation
radius_desired_ode = interp1(tex,radius_d,t_ode);
e = x(:,11) - radius_desired_ode; % Error radius
toc

figure(205)
plot(t_ode,x(:,11),'-r','LineWidth',2.75)
hold on 
h=plot(tex,rex,'o', 'MarkerSize',8,'Color',[0, 0.4470, 0.7410]),set(h, 'markerfacecolor', get(h, 'color'));
title('Bubble Radius','fontsize',18,'FontName','Garamond','FontWeight','bold');
xlabel('time (sec)', 'fontsize',28,'FontName','Garamond','FontWeight','bold');
ylabel('Radius (m)', 'fontsize',28,'FontName','Garamond','FontWeight','bold');
legend('Predicted Radius (MODEL)', 'Desired Radius (EXP)')
grid on
set(gcf,'Color',[1 1 1]);
set(gca,'fontsize',28,'FontName','Garamond','FontWeight','bold','Color','w');