"The V_2 plot should give a zero velocity when imaginary number outcome occurs."
You can replace complex elements of V_2 with 0 as soon as they are calculated:
if ~isreal(V_2(i,ii))
V_2(i,ii) = 0;
end
See below:
clc;
clear;
close all;
%basic assumption/variables
g = 9.81; %Gravity constant
d_50 = 0.018; %Rock/particle diameter(50% of the rocks) [m]
rho_w = 1000; %Water density [kg/m^3]
rho_s = 2600; %Stone density [kg/m^3]
rho_r = (rho_s-rho_w)/rho_w; %relative density [kg/m^3] also desribed as R_sd
vis_l = 0.0000013 ; %Kinematic viscosity Miedema model excel[m^2/s])
C_fr = 0.416; %assumption friction factor bed and pipe, miedema book p451 regime 3
C_rgh = 1.5*10^-3 ; %Roughness steel--> OMAE2014-23437 epsilon eq:17 --> warn cast iron: https://www.engineeringtoolbox.com/surface-roughness-ventilation-ducts-d_209.html
w_a = 0; %added water to pipe[m^3/s]
Cd = 0.4; %Form coefficient of rock, Notes rhee p.12
%situation dependent variable/inputs
L = 2 ; %Length of pipe
D_p = 0.094 ; %Pipe Diameter
lampda_w= d_50/D_p ; %Ratio particle size and pipe diameter--> used for hindered settling
D_H = D_p ; %Hydraulic diameter, assumde to be te same as Pipe diamter in this case --> cirular pipe
R_e = (3*D_H)/vis_l; %Reynolds number
C_vb = 0.20 ; %Volumetric bed concentration
A_p = (pi/4)*D_p^2; %Surface cross section pipe, described in figure 1, paper OMAW2014
a_wilson= 2.75; %Wilson factor p3 at eq 18 paper OMOW2014
n = 2.4 ; %The value chosen for the hindered sttling exponent n = 2.4, which is a value for high particle Reynolds numbers. This value is chosen because relative large particles are used for subsea rock installation (D = 0.02 - 0.1 [m]).
v_ls_1dv= 1.34*sqrt(2*g*D_p*rho_r)*cosd(80)^(1/3); %5.2.10 The Limit Deposit Velocity, from miedema
% D_H=sqrt((4*(A_p-(0.25*pi*D_p^2*(beta-sin(beta)*cos(beta))))/pi));
% alpha_list = 55:1:90;
alpha_list = 45:1:90;
beta_list = linspace(0.01*pi,0.8*pi);
%this loop is with the original bed density [C_vb]
for i = 1:length(alpha_list)
for ii = 1:length(beta_list)
beta = beta_list(ii);
alpha = alpha_list(i);
A_2=0.25*D_p^2*(beta-sin(beta)*cos(beta));
A_1=A_p-A_2;
F_gx=g*sind(alpha)*rho_r*C_vb*A_2;
F_12fl=(a_wilson*1.325)/(log(0.27*d_50/D_H+5.75/R_e^0.9))^2*0.125*(((A_2)*(1-C_vb)/(A_1))+1)^2*sin(beta)*D_p;
F_2fr=((C_fr*g*cosd(alpha)*rho_r*C_vb*A_p*(beta-sin(beta)*cos(beta))*beta)/(beta*D_p*pi))*D_p*beta;
F_1fl=0.125*1.325/(log(0.27*C_rgh/D_H+5.75/R_e^0.9))^2*((A_2*(1-C_vb))/A_1)^2*(pi-beta)*D_p;
F_2fl= 1.325/(log(0.27*C_rgh/d_50+5.75/R_e^0.9))^2*0.125*beta*D_p*(1-C_vb);
V_2(i, ii) =sqrt((F_gx-F_2fr)/(F_12fl+F_1fl+F_2fl));
if ~isreal(V_2(i,ii))
V_2(i,ii) = 0;
end
P_o(i,ii)=V_2(i,ii)*A_2*C_vb*rho_s;
end
end
plot(beta_list, V_2(i,:),'LineWidth',3)
figure(1)
surf(beta_list, alpha_list, V_2);
title('Velocity sliding bed')
xlabel('Beta size of sliding bed')
ylabel('Alpha, angle pipe[degrees]')
zlabel('Velocity [m/s]')
grid on
%Draw horizontal lines for a certain production:
hold on
contour3(beta_list,alpha_list,V_2,[0.3 0.3], '-r', 'LineWidth',2) % Draw some velocity contour lines
hold on
contour3(beta_list,alpha_list,V_2,[0.45 0.45], '-r', 'LineWidth',2) % Draw some velocity contour lines
hold on
contour3(beta_list,alpha_list,V_2,[0.63 0.63], '-r', 'LineWidth',2) % Draw some velocity contour lines
hold off
figure(2)
surf(beta_list, alpha_list, P_o);
title('Production')
xlabel('Beta size of sliding bed')
ylabel('Alpha, angle pipe[degrees]')
zlabel('Production [kg/s]')
grid on
%Draw horizontal lines for a certain production:
hold on
contour3(beta_list,alpha_list,P_o,[0.47 0.47], '-r', 'LineWidth',2) % Draw Contour At 1500 tons/H
hold on
contour3(beta_list,alpha_list,P_o,[0.235 0.235], '-r', 'LineWidth',2) % Draw Contour At 750 tons/H
hold off
