% Representation of the experimental dataset for % "Stationary Electroosmotic Flow Driven by AC Fields Around Charged Dielectric Spheres" close all;clear;clc %% Experimental data for U0. % Field: 80 kV/m. Conductivity: 1.7 mS/m % 167Hz U0_167= [262.00653, 305.29099, 319.26687, 317.43353, 205.50513, 234.51935 ... 420.11423 170.15919 289.14377 265.17237 305.54558 398.25553 343.37439 ... 342.16947 207.12036 265.25451 342.24173]*1e-6; %m M_U0_167 = mean(U0_167); STD_U0_167 = std(U0_167); % 292Hz U0_292= [158, 197.87408, 137.38681, 189.01469, 137.43513, 203.56408,160.60123,... 170.81300,179.08539, 99.00018,205.88668,138.64840,181.88285,182.89,183.84270,... 148.36302,207.63641,147.67829,211.59531,174.50382,274.19922,160.07999,... 259.26205,77.29851]*1e-6; %m M_U0_292 = mean(U0_292); STD_U0_292 = std(U0_292); % 527 Hz U0_527= [99, 63, 147.69812,148.96679,110.49295,91.59467,171.31164,72.51698,... 103.88572,92.79771,139.46999,92.93605,143.7545,99.35827,116.81864,123.29044,... 97.71257,99.53143,135.96428,80.24131,107.88546,151.29214,73.02991,118.57834,... 176.43192,121.15181,118.59056,89.09187,78.09237,118.47169,166.45603]*1e-6; %m M_U0_527 = mean(U0_527); STD_U0_527 = std(U0_527); % 1.7kHz U0_1700= [21.00830,77.88052,102.49395,84.94266,75.74054,105.37700,115.88873,... 81.66619,122.96756,74.97226,75.53630,94.36877,88.88602,92.11911,109.52405,... 145.38195,129.77701,111.07995,51.58619,103.43681,88.74308,104.60295,... 122.57149,99.11064]*1e-6; %m M_U0_1700 = mean(U0_1700); STD_U0_1700 = std(U0_1700); % 3kHz U0_3000= [55.96952,53.61539,36.03259,65.89498,66.91511,60.16586,45.26259,... 89.63774,84.90840,30.61107,88.55048,56.55701]*1e-6; %m M_U0_3000 = mean(U0_3000); STD_U0_3000 = std(U0_3000); f_exp = [167 292 527 1700 3000]; U0 = [M_U0_167,M_U0_292,M_U0_527,M_U0_1700,M_U0_3000]; errU0 = [STD_U0_167,STD_U0_292,STD_U0_527,STD_U0_1700,STD_U0_3000]; %% Theoretical values % Comparison with the experimental set f_th = logspace(2,4,20); % Scaling a = 0.5*3*1e-6; % m Particle Radius e = 1.6e-19; % C eta = 1e-3; % m^2/s epsilon = 80*8.85e-12; zeta0Dim = -75*1e-3; % V Zeta potential of particles E0 = 8e4; % V/m D = 2e-9; % m^2/s phith = 25e-3; % V Eth = phith/a; % V/m zeta0 = zeta0Dim/phith; % Scaled zeta potential % Calculation of different Du Du=zeros(1,3); Ks = 1e-9; % S sigma = 1.7e-3; % S/m Du(1) = Ks/sigma/a; conv = 1e-2/0.140823; % molal/(S/m) N_A = 6.022e23; c0 = sigma*conv*997.05*N_A; kappa = sqrt(2*e*c0/(epsilon*phith)); % 1/m delta = 1/kappa/a; alpha = epsilon*Eth^2*a^2/eta/D; Sigma0 = epsilon*kappa*phith; % C/m^2 Sigma_1 = 4*c0*e/kappa*sinh(zeta0/2); % C/m^2 Sigma = Sigma_1/Sigma0; Bi=delta*Sigma; Du(2) = -Bi*(1+2*alpha); Diff = 1/alpha; v0 = epsilon*a*Eth^2/eta; % Velocity scale omega0 = D/a^2; % rad/s f0 = omega0/2/pi; % Hz omega_th = 2*pi*f_th; omega = omega_th/omega0; k = sqrt(1i*omega); % Least square fitting for obtaining Du omega_exp = 2*pi*f_exp/omega0; xdata = sqrt(1i*omega_exp); ydata = U0/v0; fun = @(x,xdata) (E0/Eth)^2*(Uslip(xdata,x,zeta0)+... Uslip_conv(xdata,x,zeta0,Diff)); x0 = 0.1; x = lsqcurvefit(fun,x0,xdata,ydata); Du(3)=x; U0_th = zeros(length(Du),length(f_th)); for j2=1:length(Du) for j=1:length(f_th) U0_th(j2,j) = (E0/Eth)^2*(Uslip(k(j),Du(j2),zeta0)+... Uslip_conv(k(j),Du(j2),zeta0,Diff)); end Progress(j2,length(Du),'Computing theoretical Slip Velocity',0) end %% Representation figure('Name','Comparison Theory-Experiments') hold on box on axis([1e2 1e4 1e1 1e3]) ax=gca; ax.FontSize = 20; ax.LineWidth = 1; ax.FontName = 'Times New Roman'; set(gca, 'XScale','log', 'YScale','log') xlabel('\fontname{Times New Roman} \fontsize{24} f') ylabel('\fontname{Times New Roman} \fontsize{24} u_{slip}') colorR = [236,100,100]/255; colorG = [195,245,199]/255; errorbar(f_exp,U0*1e6,errU0*1e6,'ok','MarkerSize',8,'MarkerFaceColor','k') plot(f_th,U0_th(1,:)*v0*1e6,'-k','LineWidth',1); plot(f_th,U0_th(2,:)*v0*1e6,'--k','LineWidth',1); plot(f_th,U0_th(3,:)*v0*1e6,':k','LineWidth',1.75); legend('Experimental','Du = 0.392','Du = 0.113','Du = 0.048') %% Supporting functions function Uslip_conv = Uslip_conv(k,Du,zeta0,Diff) %Uslip_conv % Direct result of Mathematica Computation for the computation of the % convection term of the slip velocity under the CPEO theory 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 Uslip = Uslip(k,Du,zeta0) %Uslip_conv % Direct result of Mathematica Computation for the computation of the % Pe=0 term of the slip velocity under the CPEO theory 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 Progress(p,N,name,subindent) % Progress.m % Details the evaluation process sp = blanks(4*subindent); lineLength = length(sp)+length(name)+3+3; if p==1 perc = 0; fprintf('%s%s: %3d%%',sp,name, perc); elseif p==N perc = 100; fprintf(repmat('\b',1,lineLength)) lineLength = fprintf('%s%s: %3d%%',sp,name, perc); fprintf(repmat('\b',1,lineLength)) fprintf('\n') else perc = floor(p/N*100); if perc == 100 else fprintf(repmat('\b',1,lineLength)) fprintf('%s%s: %3d%%',sp,name, perc); end end end