% Experimental determination of maximum slip velocities clear; clc; close all %% 1. General parameters and variables % Flow velocity profile Deltap = 100; % Pa eta = 1e-3; % m^2/s L = 1e-2; % m -- Length of the channel W = 51.5e-6; % m -- width of the channel h = 60e-6; % m -- height of the channel U = UvsP(Deltap, L,h,W); epsilon = 80*8.85e-12; D = 2e-9; % m^2/s phith = 25e-3; % V sigma = 2.7e-3; % S/m Ks = 1e-9; % S colorR = [212 56 37]/255; colorO = [239 135 51]/255; colorY = [248 206 70]/255; colorG = [70 156 118]/255; colorB = [56 117 176]/255; %% 2. Particles of 1 um in size a = 0.5*1*1e-6; % m Particle radius zeta0Dim = -71.4*1e-3; % V Zeta potential y0 = 0.5e-6*ones(1,3); % 2.1. Frequency Sweep f = [50 90 167]; % Hz -- Frequency of the field E0 = 60e3; % V/m Eth = phith./a; % V/m zeta0 = zeta0Dim/phith; alpha1 = epsilon*Eth.^2.*a.^2/eta/D; D1 = 1./alpha1; v0 = epsilon*a.*Eth.^2/eta; % Velocity Scale Du = Ks/sigma./a; omega0 = D./a.^2; % rad/s f0 = omega0/2/pi; % Hz omega=2*pi*f./omega0; k = sqrt(1i*omega); % 2.1.1. Predicted slip velocities from CPEO theory N=length(f); U0_th = zeros(1,N); Up=Uslip(k,Du,zeta0)+Uslip_conv(k,Du,zeta0,D1); for j=1:N U0_th(j) = v0*(E0/Eth)^2*(Uslip(k(j),Du,zeta0)+... Uslip_conv(k(j),Du,zeta0,D1)); end % 2.1.2. Experimental Results % 2.1.3. Wall separation observed at the end of the channel y_Exp = [7.1415 6.8309 6.8309]*1e-6; % 2.1.4. Calculation of maximum slip velocity u0_WR = zeros(size(f)); for ii = 1:N u0_WR(ii) = U0vsY_v4(a,y_Exp(ii),y0(ii),L,W,U); end % 2.1.5. Representation figure hold on box on ax=gca; ax.FontSize = 20; ax.FontName = 'Arial'; h3=plot(f,U0_th*1e6,'-k','Linewidth',1.5,'color',colorG); h4 = plot(f, u0_WR*1e6,'^r','Linewidth',1.25,'Color',colorR,... 'MarkerFaceColor',colorR); legend([h3 h4],'Theoretical','Experiments','Location','SE') xlabel('Frequency (Hz)') ylabel('Slip Velocity (\mum)') title('1 \mum Particles. Frequency Sweep 60 kV/m') set(gca,'XScale','log') % 2.2. Voltage Sweep f = 50; % Hz -- Frequency of the field E0 = [60 70 80]*1e3; % V/m Eth = phith./a; % V/m zeta0 = zeta0Dim/phith; alpha1 = epsilon*Eth.^2.*a.^2/eta/D; D1 = 1./alpha1; v0 = epsilon*a.*Eth.^2/eta; % Velocity scale omega0 = D./a.^2; % rad/s f0 = omega0/2/pi; % Hz omega=2*pi*f./omega0; k = sqrt(1i*omega); % 2.2.1. Predicted slip velocities from CPEO theory N=length(E0); U0_th = zeros(1,N); Up=Uslip(k,Du,zeta0)+Uslip_conv(k,Du,zeta0,D1); for j=1:N U0_th(j) = v0*(E0(j)/Eth)^2*(Uslip(k,Du,zeta0)+... Uslip_conv(k,Du,zeta0,D1)); end % 2.2.2. Experimental Results % 2.2.3. Wall separation observed at the end of the channel y_Exp = [7.1415 8.0731 8.6943]*1e-6; % 2.2.4. Calculation of maximum slip velocity u0_WR = zeros(size(E0)); for ii = 1:N u0_WR(ii) = U0vsY_v4(a,y_Exp(ii),y0(ii),L,W,U); end % 2.2.5. Representation figure hold on box on ax=gca; ax.FontSize = 20; ax.FontName = 'Arial'; h3=plot(E0,U0_th*1e6,'-k','Linewidth',1.5,'color',colorG); h4 = plot(E0, u0_WR*1e6,'^r','Linewidth',1.25,'Color',colorR,... 'MarkerFaceColor',colorR); legend([h3 h4],'Theoretical','Experiments','Location','SE') xlabel('Electric field magnitude (V/m)') ylabel('Slip Velocity (\mum)') title('1 \mum Particles. Voltage Sweep 50 Hz') set(gca,'XScale','log') %% 3. Particles of 3 um in size a = 0.5*3*1e-6; % m Particle radius zeta0Dim = -78.6*1e-3; % V Zeta potential y0 = 1.5e-6*ones(1,3); % 3.1. Frequency Sweep f = [50 90 167]; % Hz -- Frequency of the field E0 = 60e3; % V/m Eth = phith./a; % V/m zeta0 = zeta0Dim/phith; alpha1 = epsilon*Eth.^2.*a.^2/eta/D; D1 = 1./alpha1; v0 = epsilon*a.*Eth.^2/eta; % Velocity Scale Du = Ks/sigma./a; omega0 = D./a.^2; % rad/s f0 = omega0/2/pi; % Hz omega=2*pi*f./omega0; k = sqrt(1i*omega); % 3.1.1. Predicted slip velocities from CPEO theory N=length(f); U0_th = zeros(1,N); Up=Uslip(k,Du,zeta0)+Uslip_conv(k,Du,zeta0,D1); for j=1:N U0_th(j) = v0*(E0/Eth)^2*(Uslip(k(j),Du,zeta0)+... Uslip_conv(k(j),Du,zeta0,D1)); end % 3.1.2. Experimental Results % 3.1.3. Wall separation observed at the end of the channel y_Exp = [16.4582 16.1477 15.2160]*1e-6; % 3.1.4. Calculation of maximum slip velocity u0_WR = zeros(size(f)); for ii = 1:N u0_WR(ii) = U0vsY_v4(a,y_Exp(ii),y0(ii),L,W,U); end % 3.1.5. Representation figure hold on box on ax=gca; ax.FontSize = 20; ax.FontName = 'Arial'; h3=plot(f,U0_th*1e6,'-k','Linewidth',1.5,'color',colorG); h4 = plot(f, u0_WR*1e6,'^r','Linewidth',1.25,'Color',colorR,... 'MarkerFaceColor',colorR); legend([h3 h4],'Theoretical','Experiments','Location','SE') xlabel('Frequency (Hz)') ylabel('Slip Velocity (\mum)') title('3 \mum Particles. Frequency Sweep 60 kV/m') set(gca,'XScale','log') % 3.2. Voltage Sweep f = 50; % Hz -- Frequency of the field E0 = [40 50 60]*1e3; % V/m Eth = phith./a; % V/m zeta0 = zeta0Dim/phith; alpha1 = epsilon*Eth.^2.*a.^2/eta/D; D1 = 1./alpha1; v0 = epsilon*a.*Eth.^2/eta; % Velocity scale omega0 = D./a.^2; % rad/s f0 = omega0/2/pi; % Hz omega=2*pi*f./omega0; k = sqrt(1i*omega); % 3.2.1. Predicted slip velocities from CPEO theory N=length(E0); U0_th = zeros(1,N); Up=Uslip(k,Du,zeta0)+Uslip_conv(k,Du,zeta0,D1); for j=1:N U0_th(j) = v0*(E0(j)/Eth)^2*(Uslip(k,Du,zeta0)+... Uslip_conv(k,Du,zeta0,D1)); end % 3.2.2. Experimental Results % 3.2.3. Wall separation observed at the end of the channel y_Exp = [13.3526 14.9054 16.4582]*1e-6; % 3.2.4. Calculation of maximum slip velocity u0_WR = zeros(size(E0)); for ii = 1:N u0_WR(ii) = U0vsY_v4(a,y_Exp(ii),y0(ii),L,W,U); end % 3.2.5. Representation figure hold on box on ax=gca; ax.FontSize = 20; ax.FontName = 'Arial'; h3=plot(E0,U0_th*1e6,'-k','Linewidth',1.5,'color',colorG); h4 = plot(E0, u0_WR*1e6,'^r','Linewidth',1.25,'Color',colorR,... 'MarkerFaceColor',colorR); legend([h3 h4],'Theoretical','Experiments','Location','SE') xlabel('Electric field magnitude (V/m)') ylabel('Slip Velocity (\mum)') title('3 \mum Particles. Voltage Sweep 50 Hz') set(gca,'XScale','log') %% Supporting functions function U = UvsP(p, L,l,h) % Calculo de la maxima velocidad en un canal de seccion rectangular con un flujo de % Poiseuille dada una diferencia de presion % Todas las unidades en el SI % Se para la evaluacion del sumatorio al tener una discrepancia entre % sucesivos sumandos por debajo del 1% % p -- PRESION % L -- LONGITUD DEL CANAL % l -- ALTURA DEL CANAL % h -- ANCHURA DEL CANAL % Por si en algun momento se utiliza como funcion de la posicion en el % canal, se define y = h/2; z=l/2; eta = 1e-3; % Pa*s Viscosidad del agua G = p/L; U2D = G/2/eta*y*(h-y); factor = 4*G*h^2/eta/pi^3; N = 1; n=1; U = U2D; while N>0.01 beta = (2*n-1)*pi/h; Un = -factor*(sinh(beta*z)+sinh(beta*(l-z)))*sin(beta*y)/sinh(beta*l)/(2*n-1)^3; N = abs(Un/U); U = U+Un; n=n+1; end end function Uslip = Uslip(k,Du,zeta0) %alphaE_NoConv % Direct result of Mathematica Computation Uslip= (1/40).*(1+4.*Du).^(-1).*abs(2+2.*k+k.^2+Du.*(2+k).^2).^(-2).*(( ... -2).*abs(zeta0).*(30.*Du.^2.*(1+2.*Du).*real((1+k).*(2+2.*conj( ... k)+conj(k).^2))+abs(2+2.*k+k.^2+Du.*(2+k).^2).^2.*(10.*Du.* ... real(Du.*exp(1).^k.*k.^3.*(2+2.*k+k.^2+Du.*(2+k).^2).^(-1).*(2+2.* ... conj(k)+conj(k).^2+conj(Du).*(2+conj(k)).^2) ... .^(-1).*((-1).*(2+2.*conj(k)+conj(k).^2).*(3.*k.^3.* ... igamma((-6),k)+3.*k.^3.*igamma((-5),k)+2.*k.^3.*igamma((-4),k)+(-6).* ... igamma((-3),k)+(-6).*igamma((-2),k)+(-2).*igamma((-1),k))+2.* ... conj(Du).*(3.*k.^3.*igamma((-6),k)+3.*k.^3.*igamma((-5),k)+2.* ... k.^3.*igamma((-4),k)+12.*igamma((-3),k)+12.*igamma((-2),k)+4.*igamma(( ... -1),k)+conj(k).^2.*(3.*k.^3.*igamma((-6),k)+3.*k.^3.*igamma(( ... -5),k)+2.*k.^3.*igamma((-4),k)+3.*igamma((-3),k)+3.*igamma((-2),k)+ ... igamma((-1),k))+conj(k).*(3.*k.^3.*igamma((-6),k)+3.*k.^3.* ... igamma((-5),k)+2.*k.^3.*igamma((-4),k)+12.*igamma((-3),k)+12.*igamma(( ... -2),k)+4.*igamma((-1),k)))))+(-1).*(1+4.*Du).*(real(Du.*exp(1).^k.* ... k.^(-2).*(2+2.*k+k.^2+Du.*(2+k).^2).^(-1).*(12.*k.^5.*igamma((-3), ... k)+12.*k.^5.*igamma((-2),k)+4.*k.^5.*igamma((-1),k)+exp(1).^((-1).* ... k).*k.^3.*(2+2.*conj(k)+conj(k).^2+conj(Du).*(2+ ... conj(k)).^2).^(-1).*((-2)+(-2).*conj(k)+(-1).*conj( ... k).^2+2.*conj(Du).*(1+conj(k)+conj(k).^2)).*(6+6.* ... exp(1).^k.*k.^5.*igamma((-6),k)+6.*exp(1).^k.*k.^5.*igamma((-5),k)+ ... 4.*exp(1).^k.*k.^5.*igamma((-4),k)+9.*exp(1).^k.*igamma((-1),k)+9.* ... exp(1).^k.*igamma(0,k))+18.*igamma(2,k)+18.*igamma(3,k)+6.*igamma(4,k) ... ))+(-3).*real(Du.*exp(1).^k.*k.*(2+2.*k+k.^2+Du.*(2+k).^2).^(-1).*(( ... -6).*k.^2.*igamma((-3),k)+(-6).*k.^2.*igamma((-2),k)+(-2).*k.^2.* ... igamma((-1),k)+(-1).*exp(1).^((-1).*k).*(2+2.*conj(k)+ ... conj(k).^2+conj(Du).*(2+conj(k)).^2).^(-1).*((-2)+( ... -2).*conj(k)+(-1).*conj(k).^2+2.*conj(Du).*(1+ ... conj(k)+conj(k).^2)).*((-2)+3.*exp(1).^k.*k.^5.*igamma(( ... -6),k)+3.*exp(1).^k.*k.^5.*igamma((-5),k)+2.*exp(1).^k.*k.^5.* ... igamma((-4),k)+(-3).*exp(1).^k.*igamma((-1),k)+(-3).*exp(1).^k.* ... igamma(0,k))+6.*k.^(-3).*igamma(2,k)+6.*k.^(-3).*igamma(3,k)+2.*k.^( ... -3).*igamma(4,k))))))+180.*Du.^2.*abs(1+k).^2.*(Du.*abs(zeta0)+(-2) ... .*(1+2.*Du).*log(cosh((1/4).*zeta0))+(-1).*(1+4.*Du).*sech((1/2).* ... zeta0).*sinh((1/4).*zeta0).^2)+5.*(96.*Du.^3.*log(cosh((1/4).* ... zeta0)).*real((1+k).*(2+2.*conj(k)+conj(k).^2))+abs(2+2.* ... k+k.^2+Du.*(2+k).^2).^2.*((-16).*Du.*log(cosh((1/4).*zeta0)).*real( ... Du.*exp(1).^k.*k.^3.*(2+2.*k+k.^2+Du.*(2+k).^2).^(-1).*(2+2.* ... conj(k)+conj(k).^2+conj(Du).*(2+conj(k)).^2) ... .^(-1).*((-1).*(2+2.*conj(k)+conj(k).^2).*(3.*k.^3.* ... igamma((-6),k)+3.*k.^3.*igamma((-5),k)+2.*k.^3.*igamma((-4),k)+(-6).* ... igamma((-3),k)+(-6).*igamma((-2),k)+(-2).*igamma((-1),k))+2.* ... conj(Du).*(3.*k.^3.*igamma((-6),k)+3.*k.^3.*igamma((-5),k)+2.* ... k.^3.*igamma((-4),k)+12.*igamma((-3),k)+12.*igamma((-2),k)+4.*igamma(( ... -1),k)+conj(k).^2.*(3.*k.^3.*igamma((-6),k)+3.*k.^3.*igamma(( ... -5),k)+2.*k.^3.*igamma((-4),k)+3.*igamma((-3),k)+3.*igamma((-2),k)+ ... igamma((-1),k))+conj(k).*(3.*k.^3.*igamma((-6),k)+3.*k.^3.* ... igamma((-5),k)+2.*k.^3.*igamma((-4),k)+12.*igamma((-3),k)+12.*igamma(( ... -2),k)+4.*igamma((-1),k)))))+9.*(1+4.*Du).*real(Du.*(1+k).*(2+2.*k+ ... k.^2+Du.*(2+k).^2).^(-1).*(2+2.*conj(k)+conj(k).^2+2.* ... conj(Du).*(1+conj(k))).*(2+2.*conj(k)+conj(k) ... .^2+conj(Du).*(2+conj(k)).^2).^(-1)).*tanh((1/2).*abs( ... zeta0))))); end function Uslip_conv = Uslip_conv(k,Du,zeta0,Diff) %alphaE_NoConv % Direct result of Mathematica Computation Uslip_conv = (-1).*(10+40.*Du).^(-1).*(Du.*abs(zeta0)+(-2).*(1+2.*Du).*log( ... cosh((1/4).*zeta0))).*(real(Diff.^(-1).*Du.*k.^(-2).*(2+2.*k+k.^2+ ... Du.*(2+k).^2).^(-1).*(2+2.*conj(k)+conj(k).^2+conj( ... Du).*(2+conj(k)).^2).^(-1).*conj((2+2.*k+k.^2+2.*Du.*(1+ ... k)).*zeta0+8.*Du.*(1+k).*log(cosh((1/4).*zeta0))).*(6.*k.^3+6.* ... exp(1).^k.*k.^8.*igamma((-6),k)+6.*exp(1).^k.*k.^8.*igamma((-5),k)+ ... 4.*exp(1).^k.*k.^8.*igamma((-4),k)+(-12).*exp(1).^k.*k.^5.*igamma(( ... -3),k)+(-12).*exp(1).^k.*k.^5.*igamma((-2),k)+9.*exp(1).^k.*k.^3.* ... igamma((-1),k)+(-4).*exp(1).^k.*k.^5.*igamma((-1),k)+9.*exp(1).^k.* ... k.^3.*igamma(0,k)+(-18).*exp(1).^k.*igamma(2,k)+(-18).*exp(1).^k.* ... igamma(3,k)+(-6).*exp(1).^k.*igamma(4,k)))+3.*real(Diff.^(-1).*Du.* ... k.^(-2).*(2+2.*k+k.^2+Du.*(2+k).^2).^(-1).*(2+2.*conj(k)+ ... conj(k).^2+conj(Du).*(2+conj(k)).^2).^(-1).* ... conj((2+2.*k+k.^2+2.*Du.*(1+k)).*zeta0+8.*Du.*(1+k).*log( ... cosh((1/4).*zeta0))).*((-2).*k.^3+3.*exp(1).^k.*k.^8.*igamma((-6), ... k)+3.*exp(1).^k.*k.^8.*igamma((-5),k)+2.*exp(1).^k.*k.^8.*igamma(( ... -4),k)+(-6).*exp(1).^k.*k.^5.*igamma((-3),k)+(-6).*exp(1).^k.* ... k.^5.*igamma((-2),k)+(-3).*exp(1).^k.*k.^3.*igamma((-1),k)+(-2).* ... exp(1).^k.*k.^5.*igamma((-1),k)+(-3).*exp(1).^k.*k.^3.*igamma(0,k)+ ... 6.*exp(1).^k.*igamma(2,k)+6.*exp(1).^k.*igamma(3,k)+2.*exp(1).^k.* ... igamma(4,k)))); end function [U0] = U0vsY_v4(a,y,y0,L,W,U) % SI units f=@(yp) 1/1536*(18*(1-2*yp).^2-9*(1-2*yp).^4+2*(1-2*yp).^6-12*log(-1+2*yp)); U0 = 32/3*U*W^3/L/a^2*(f(y/W)-f(y0/W)); end