% Model rheometer and predict behavior for different liquids % plot experimental acceleration ratios % use plot to find viscosity that matches experimental data %% User instructions for rheometer % 1. run accelprocessor_foruser for your accelerometer data % 2. enter the mass of the sphere below mass = 29.67*10^-3; %kg % 3. enter the radius of the sphere below radius = 0.015; %m % 4. enter the spring constant of the spring k = 1050.76; % Nm % 5. enter the density of the fluid density = 1261; % kg/m^3 % 6. enter your guess for the viscosity of the fluid eta = 100; % cP %% processing code % DEFINE CONSTANTS IN FORCE EQUATION (FROM DOLFO PAPER) % convert viscosity from cP to Pa*s eta_guess = eta*10^-3; % create matrix for range of omega values omega = linspace(100,800,1000); % viscous penetration depth (m) delta = sqrt(2*eta_guess./(density.*omega)); % calculate damping coefficient, Nm/s c = -6*pi*eta_guess*radius.*(1+radius./delta); % natural frequency, rad/s omega_n = sqrt(k/mass); % calculate viscous damping factor zeta = c./(2*mass*omega_n); % ratio of frequency to natural frequency r = omega/omega_n; % set up matrix to store calculated acceleration ratios Accel_ratio = []; % calculate acceleration ratios at every frequency for i = 1:length(omega) Accel_ratio(i) = 1/sqrt((1-r(i)^2)^2+(2*zeta(i)*r(i))^2); end % set up matrix of frequencies at which experimental data was collected omega_exp = 15:5:75; omega_exp = 2*pi*omega_exp; %% calculate best approximation of viscosity % comment this out if you want to brute force the approximation % find maximum acceleration ratio in exp. data max_A = max(ratio_mat); % find maximum acceleration ratio in estimated calulations A_approx = max(Accel_ratio); % find % difference between maximum accel ratios percent_dif_w = abs((max_A-A_approx))/max_A; % define the increment that the approximate viscosity should increase by % to be 10% of its estimated value increment = 0.1*eta_guess; % this can be adjust by changing the fraction % continue to increase estimated viscosity until max amplitude give % % difference of less than 1 while percent_dif_w > 0.01 eta_guess = eta_guess + increment; % increase estimate % recalculate constants delta = sqrt(2*eta_guess./(density.*omega)); c = -6*pi*eta_guess*radius.*(1+radius./delta); zeta = c./(2*mass*omega_n); for i = 1:length(omega) Accel_ratio(i) = 1/sqrt((1-r(i)^2)^2+(2*zeta(i)*r(i))^2); end [A_approx, omegan_approx] = max(Accel_ratio); % recalculate max accel ratio percent_dif_w = abs((max_A-A_approx))/max_A; % recalculate % dif end %% Plot experimental and calculated acceleration ratios figure hold on % plot experimental and calculated accel ratios vs frequency scatter(omega_exp, ratio_mat) plot(omega, Accel_ratio) xlabel('Omega (rad/s)','fontsize', 16) ylabel('Acceleration Ratio','fontsize', 16) title('Expected vs Calculated Magnification','fontsize',18) xlim([100 800]) %ylim([0 10]) legend('Experimental response','Predicted response') % print calculated viscosity visc = ['Viscosity = ', num2str(eta_guess*10^3)]; disp(visc)