ODESolver_PID.m

function dx = ODESolver_PID(t, x)
persistent radius_dPrev time_prev 

if isempty(radius_dPrev)
    radius_dPrev = 0;
    time_prev = 0;
end

global error radius_d tex MaxRadExp; 

dx = zeros(12,1);

%Parameters:
sigma = 0.072;                  % Surface tension, kg/s^2
mu = 1e-3;                      % Water kinematic viscosity, m^2/s
rho_L = 1e3;                    % Water density, kg/m^3
k = 1.4;                        % Specific heat ratio
c = 1500;
R_0 = MaxRadExp;
n = 3;
D = 1.76e-09;
K_B = 6.737e09;
alpha_M = .8544;
R_v = 461.52;
Na = 6.02e23;
Rg = 8.314;
Rdot_0 = 0;
mdot_0 = 0;
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;
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;

if (x(6)<258.15)
    P_v = 180.547;
elseif (x(6)>258.15 && x(6)<647.3)
    P_v = 133.322*exp(18.3036-3816.44/(x(6)-46.13));
else
    P_v = 2.075e07;
end
P_v_infty = 133.322*exp(18.3036-3816.44/(T_infty-46.13));

%%%%%%%%%%%%%%%%%%%%%%%%%%% PID IMPLEMENTATION %%%%%%%%%%%%%%%%%%%%%%%%%% 
% radius_d first derivative
radius_desired_current = interp1(tex,radius_d,t);
if t==time_prev
    radius_dDot = 0;
else
    radius_dDot  = ( radius_desired_current - radius_dPrev ) / (t-time_prev);
end
radius_dPrev = radius_desired_current;

% PID tuning
Kp = 6500000;
Kd = 0;%25;%20;
Ki = 0;%250000;

ind = find(tex<t & tex>time_prev);
if ~isempty(ind) 

    % u: feedback signal
    u = Kp*(radius_desired_current - x(11)) + Kd*( radius_dDot - x(12)) + Ki*error;
    error = error + (radius_desired_current - x(1))*(t-time_prev);
else
    u = (Kp*(radius_desired_current - x(11)) + Kd*( radius_dDot - x(12)) + Ki*error)/2.5;
    error = error + (radius_desired_current - x(1))*(t-time_prev);
end


%%%%%%%%%%%%%%%%%%%%%%%%%% Equations to solve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y = x(5)/(Na*x(7)/x(3)) + 0.625*(x(5)/(Na*x(7)/x(3)))^2 + 0.2869*(x(5)/(Na*x(7)/x(3)))^3+0.115*(x(5)/(Na*x(7)/x(3)))^4;
kappa0_water = 9.98e-05 * x(8) - 0.0119;
kappa0_air = 5.39e-05 * x(8) + 0.0108;
kappa0_mix = 0.5*(x(2)/x(3)*kappa0_water + x(1)/x(3)*kappa0_air + 1/(x(2)/(x(3)*kappa0_water)+x(1)/(x(3)*kappa0_air)));
kappa = x(5)/(Na*x(7)/x(3))*(1/Y+1.2+0.755*Y)*kappa0_mix;

if (0.1 * P_infty - x(1) / x(3) * x(9) <0)
	Z = x(8);
else 
	Z = T_0;
end
% keyboard;
G = -1/(2*1.138e-23) * sqrt(pi/1.38e-23) * 1.173 * kappa * 4e-19 / n;
lambda = x(7) / (sqrt(2)* 0.4e-18 * x(3));

E_water_infty = 1.04334e-20 + 6.27210e-23 * (T_infty-250);   
E_water_T_B = 1.04334e-20 + 6.27210e-23 * (x(8)-250);   
E_air_Z = 8.64814e-21 + 4.40808e-23 * (Z-250);   
E_air_T = 8.64814e-21 + 4.40808e-23 * (x(6)-250);   
E_water_T = 1.04334e-20 + 6.27210e-23 * (x(6)-250);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ODEs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


dx(1) = 4*pi*x(11)*D*1000*rho_L*Na/(18*K_B)*(0.1*P_infty-x(1)/x(3)*x(9));

GammaFunc = x(10)*x(3)/(x(2)*x(9))*sqrt(R_v*x(8)/2);
funGamma = @(ttt) exp(-ttt.^2);
IntTerm = integral(@(ttt)funGamma(ttt),0, GammaFunc);
Gamma = exp(-GammaFunc^2)-sqrt(pi)* GammaFunc*(1-2/sqrt(pi)*IntTerm);

dx(2) = 4*pi*x(11)^2*1000*Na/18*alpha_M/sqrt(2*pi*R_v)*(P_v_infty/sqrt(T_infty)- Gamma*x(2)*x(9)/(x(3)*sqrt(x(8))));
dx(3) = dx(1)+ dx(2);

M = 0.001/(Na*x(3))*(18*x(2)+28.97*x(1));
dM = 0.001/(Na*x(3)^2)*(18*(x(3)*dx(2)-x(2)*dx(3))+28.97*(x(3)*dx(1) -x(1)*dx(3)));
X = Na*x(7)/x(3);
dX = 4*pi/3*Na/x(3)^2*(3*x(11)^2*x(12)*x(3)-x(11)^3*dx(3));

dx(4) = 2*a_air_VdW*x(1)/x(3)^3*(x(3)*dx(1)-x(1)*dx(3)) + 2*a_air_water_VdW*x(2)/x(3)^3*(x(3)*dx(1)-x(1)*dx(3)) + 2*a_air_water_VdW*x(1)/x(3)^3*(x(3)*dx(2)-x(2)*dx(3))+ 2*a_water_VdW*x(2)/x(3)^3*(x(3)*dx(2)-x(2)*dx(3));
dx(5) = 2/x(3)^3*((b_air_VdW*x(1)+b_air_water_VdW*x(2))*(x(3)*dx(1)-x(1)*dx(3))+( b_water_VdW*x(2)+b_air_water_VdW*x(1))*(x(3)*dx(2)-x(2)*dx(3)));
dx(6) = 1/(4.40808e-23*x(1)+6.2721e-23*x(2))*(-x(9)*4*pi*x(11)^2*x(12)+kappa*4*pi*x(11)^2*(x(8)-x(6))/(n*lambda)+4*pi*x(11)^2*1000/18*Na*alpha_M/(sqrt(2*pi*R_v))*(P_v_infty/sqrt(T_infty)*E_water_infty- Gamma*x(2)*x(9)/(x(3)*sqrt(x(8)))*E_water_T_B)+4*pi*x(11)*D*1000*rho_L*Na/(18*K_B)*(0.1*P_infty-x(1)/x(3)*x(9))*E_air_Z-x(3)^2*x(4)/Na^2/x(7)^2*4*pi*x(11)^2*x(12)+1/(Na^2*x(7))*(2*x(3)*dx(3)*x(4)+x(3)^2*dx(4))-dx(1)*E_air_T-dx(2)*E_water_T);
dx(7) = 4*pi*x(11)^2*x(12);
dx(8) = G*(-dx(6)* sqrt(M/x(8))+(x(8)-x(6))/(2*x(8))* sqrt(x(8)/M)*dM)/(1-G*sqrt(M/x(8))+M/(2*x(8)^2)*(x(8)-x(6))*sqrt(x(8)/M));
dx(9) = (Rg*dx(6)-(x(9)+x(4)/X^2)*(dX-dx(5)))/(X-x(5))+2*x(4)/X^3*dX-dx(4)/X^2;
dx(10) = -alpha_M/sqrt(2*pi*R_v)*Gamma*(x(9)/sqrt(x(8))/x(3)^2*(x(3)*dx(2)-x(2)*dx(3))+x(2)/x(3)*(dx(9)/sqrt(x(8))-0.5*x(9)/x(8)^1.5*dx(8)));

Rho_g = 0.001/(Na*x(7))*(18*x(2)+28.97*x(1));
dRho_g = 0.001/(Na*x(7)^2)*(x(7)*(18*dx(2)+28.97*dx(1))-dx(7)*(18*x(2)+28.97*x(1)));


dx(11) = x(12);

% dx(12) = ((1+x(12)/c)/rho_L*(x(9)-2*sigma/x(11)-4*mu/x(11)*(x(12)-x(10)/rho_L)-x(10)^2*(1/rho_L-1/Rho_g)-Pinf)+dx(10)*x(11)/rho_L*(1-x(12)/c+x(10)/(rho_L*c))+ x(10)/rho_L*(x(12)+ x(10)/(2*rho_L)+x(12)* x(10)/(2*rho_L*c))-1.5*x(12)^2*(1-x(12)/(3*c)+2* x(10)/(3*rho_L*c))+x(11)/(rho_L*c)*(dx(9)+2*sigma*x(12)/x(11)^2+4*mu*(x(12)^2/x(11)^2+(x(11)*dx(10)-x(12)* x(10))/(rho_L*x(11)^2))-2*x(10)*dx(10)*(1/rho_L-1/Rho_g)-x(10)^2/Rho_g^2*dRho_g))/(x(11)*(1-x(12)/c+x(10)/(c*rho_L))+4*mu/(c*rho_L));
%%%%%%%%%%%%%%%%%%%%%%%%%%% ADDING FEEDBACK SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%
dx(12) = (u+(1+x(12)/c)/rho_L*(x(9)-2*sigma/x(11)-4*mu/x(11)*(x(12)-x(10)/rho_L)-x(10)^2*(1/rho_L-1/Rho_g)-Pinf)+dx(10)*x(11)/rho_L*(1-x(12)/c+x(10)/(rho_L*c))+ x(10)/rho_L*(x(12)+ x(10)/(2*rho_L)+x(12)* x(10)/(2*rho_L*c))-1.5*x(12)^2*(1-x(12)/(3*c)+2* x(10)/(3*rho_L*c))+x(11)/(rho_L*c)*(dx(9)+2*sigma*x(12)/x(11)^2+4*mu*(x(12)^2/x(11)^2+(x(11)*dx(10)-x(12)* x(10))/(rho_L*x(11)^2))-2*x(10)*dx(10)*(1/rho_L-1/Rho_g)-x(10)^2/Rho_g^2*dRho_g))/(x(11)*(1-x(12)/c+x(10)/(c*rho_L))+4*mu/(c*rho_L));

time_prev = t;
end