Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@jeshragh
jeshragh authored Mar 2, 2021
1 parent b197976 commit fde2fe8
Show file tree
Hide file tree
Showing 2 changed files with 236 additions and 0 deletions.
90 changes: 90 additions & 0 deletions CavitationPID.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
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');
146 changes: 146 additions & 0 deletions ODESolver_PID.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
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

0 comments on commit fde2fe8

Please sign in to comment.