Commit 60d9a3f5 authored by Sebastian Hofmann's avatar Sebastian Hofmann ☝🏽
Browse files

add acc_vecor_1mm2_bin_force_calc

parent eada120b
% ** coded by Sebastian Hofmann **
% ** s.m.hofmann@mailbox.org ** sebastian.hofmann@tuhh.de **
%
% Creation of data for Figure 9, force distribution throughout reactor as 2D contour plot in 1 mm x 1 mm grid
% Load .mat created by "pie_piece_from_loadDat_exp_sim.m", experimental or simulation, one after
% another!!
% (1) is done for all data sets, (2) combines all data sets to one .mat file
% (3) calculate forces from acc and mass, (4) 2D contour plot of force distribution
%
% INFORMATION
% 1. Define dimensions in l.29-30 and file name in l.117-118
% 2. Used data here are as follows:
% 1-6: 252 rpm: sim flowtracer - sim 180 1000 - sim 180 1024 - sim 740 1024 - exp 180 1000 - exp 740 1024
% 7-12: 450 rpm: sim flowtracer - sim 180 1000 - sim 180 1024 - sim 740 1024 - exp 180 1000 - exp 740 1024
% meaning 180 or 740 µm in particle size and 1000 or 1024 kg/m^-3 particle density
% 3. Data set of "sim 180 1024" is not used for the publication!
close all
clear all
%% (1) Load .mat file with separated tracks and pie_piece (front view)
[FileName,FilePath] = uigetfile('*.mat','Select the .mat file', ...
'C:\','MultiSelect', 'off');
load(strcat(FilePath,FileName));
% 1st step - Define particle densities and sizes
rho_p = [1000; 1024]; % part densities in kg/m^3 - PE and alginate
d_p = [180; 732]; % part diameter in µm - PE and alginate
% General calculation for all subs, radial distance
r = sqrt((x).^2+(z).^2); % calculation of "r" position in reactor from x and y
% Calculate center of bins of respective r-axis subdivision - definition of grid size
r_spaces = linspace(5,65,60); % reactor radius with 130 mm divided in 2 = 65 mm ->
% leave inner 5 mm due to shaft obscuration (1mm x 1mm grid)
[N, edges_r]=histcounts(r_spaces,60); % edges_r used later for calculation
% j=0;
% center_edges_r = zeros(size(edges_r,2)-1,1);
% for i = 1:size(edges_r,2)-1
% j = j + 1;
% center_edges_r(j,1) = (edges_r(1,j) + edges_r(1,j + 1))./2;
% end
% Calculate center of bins of respective y-axis subdivision - definition of grid size
y_spaces = linspace(0,228,228); % reactor height with 228 mm ->
% division into 228 bin to gain 1mm x 1mm grid
[N, edges_y]=histcounts(y_spaces,228); % edges_y used later for calculation
% j=0;
% center_edges_y = zeros(size(edges_y,2)-1,1);
% for i = 1:size(edges_y,2)-1
% j = j + 1;
% center_edges_y(j,1) = (edges_y(1,j) + edges_y(1,j + 1))./2;
% end
% Calculate entire grid y-r with accelerations
ax_avg_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1); % Create zero vector for ax acc component
ay_avg_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1); % Create zero vector for ay acc component
az_avg_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1); % Create zero vector for az acc component
a_mag_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1);
m = 0; % counter for radius bin
for k = 1:size(edges_r,2)-1 % starts at bottom left bin and goes to top right bin of overall graph
m = m + 1;
axr = u(r>edges_r(1,m) & r<edges_r(1,m + 1)); % Look for all accelerations in entire r bin
ayr = v(r>edges_r(1,m) & r<edges_r(1,m + 1));
azr = w(r>edges_r(1,m) & r<edges_r(1,m + 1));
a_magr = v_mag(r>edges_r(1,m) & r<edges_r(1,m + 1));
y_axr = y(r>edges_r(1,m) & r<edges_r(1,m + 1)); % Look for respective y coordinates to prev. determined acceleration components
y_ayr = y(r>edges_r(1,m) & r<edges_r(1,m + 1));
y_azr = y(r>edges_r(1,m) & r<edges_r(1,m + 1));
y_a_magr = y(r>edges_r(1,m) & r<edges_r(1,m + 1));
ax_avg_bin_r = zeros(size(edges_y,2)-1,1); % Define zero vector for upcoming for loop for avg
ay_avg_bin_r = zeros(size(edges_y,2)-1,1);
az_avg_bin_r = zeros(size(edges_y,2)-1,1);
a_mag_avg_bin_r = zeros(size(edges_y,2)-1,1);
ax_avg_bin_r_std = zeros(size(edges_y,2)-1,1); % Define zero vector for upcoming for loop for std
ay_avg_bin_r_std = zeros(size(edges_y,2)-1,1);
az_avg_bin_r_std = zeros(size(edges_y,2)-1,1);
a_mag_avg_bin_r_std = zeros(size(edges_y,2)-1,1);
j = 0; % counter for axial/y bin
for i = 1:size(edges_y,2)-1 % Start scanning from bottom to top to corr. defined y bins
j = j + 1;
[a b] = find(y_axr(:,1) > edges_y(1,j) & y_axr(:,1) < edges_y(1,j + 1)); % Here: with find function, to write resp. row number in a
ax_avg_bin_r(j,1) = mean(axr(a,:),'omitnan'); % Write avg velocity in resp. r bin in resp. y bin, goes with j+1 for entire reactor height
ay_avg_bin_r(j,1) = mean(ayr(a,:),'omitnan');
az_avg_bin_r(j,1) = mean(azr(a,:),'omitnan');
a_mag_avg_bin_r(j,1) = mean(a_magr(a,:),'omitnan');
ax_avg_bin_r_std(j,1) = std(axr(a,:),'omitnan'); % Write std in resp. r bin in resp. y bin, not used in this study
ay_avg_bin_r_std(j,1) = std(ayr(a,:),'omitnan');
az_avg_bin_r_std(j,1) = std(azr(a,:),'omitnan');
a_mag_avg_bin_r_std(j,1) = std(a_magr(a,:),'omitnan');
end
ax_avg_bin(:,m) = ax_avg_bin_r; % Write entire created avg column for r bin in prepared zero vector, goes with m+1
ay_avg_bin(:,m) = ay_avg_bin_r;
az_avg_bin(:,m) = az_avg_bin_r;
a_mag_avg_bin(:,m) = a_mag_avg_bin_r;
ax_avg_bin_std(:,m) = ax_avg_bin_r_std; % Write entire created std column for r bin in prepared zero vector, goes with m+1
ay_avg_bin_std(:,m) = ay_avg_bin_r_std;
az_avg_bin_std(:,m) = az_avg_bin_r_std;
a_mag_avg_bin_std(:,m) = a_mag_avg_bin_r_std;
end
toc
% Save results for dataset
%fname = sprintf('acc_vectors_all_1mm2_sim_200LX_1024_740_450rpm_pie'); % Enter filename here!
fname = sprintf('acc_vectors_all_1mm2_exp_1024_740_450rpm_pie'); % Enter filename here!
fpath = sprintf('C:/'); % Enter filepath here!
save(strcat(fpath,'./',fname,'.mat'),...
'ax_avg_bin','ay_avg_bin','az_avg_bin','a_mag_avg_bin',...
'ax_avg_bin_std','ay_avg_bin_std','az_avg_bin_std','a_mag_avg_bin_std',...
'edges_r','edges_y')
fprintf('Done saving!\n')
%% (2) Reload data from single mat file and saves them with resp. index for plotting
clear all
%n..1-6 252 rpm: sim flowtracer - sim 180 1000 - sim 180 1024 - sim 740 1024 - exp 180 1000 - exp 740 1024
%n..7-12 450 rpm: sim flowtracer - sim 180 1000 - sim 180 1024 - sim 740 1024 - exp 180 1000 - exp 740 1024
n = 0;
for i=1:12
[FileName,FilePath] = uigetfile('*.mat','Select the .mat file', ...
'C:/','MultiSelect', 'off');
load(strcat(FilePath,FileName));
n = n+1;
acc_avg_bin{n} = cat(4,[ax_avg_bin],[ay_avg_bin],[az_avg_bin],[a_mag_avg_bin]);
acc_avg_bin_std{n} = cat(4,[ax_avg_bin_std],[ay_avg_bin_std],[az_avg_bin_std],[a_mag_avg_bin_std]);
end
fname = sprintf('graph_acc_vectors_all_1mm2_252rpm_450rpm_pie'); % Enter filename here!
fpath = sprintf('C:/'); % Enter filepath here!
save(strcat(fpath,'./',fname,'.mat'),...
'acc_avg_bin','acc_avg_bin_std',...
'edges_r','edges_y')
%% (3) Calculate forces from acc_avg of PE and alginate particle for 252 and 450 rpm
F_avg_PE_exp_252 = acc_avg_bin{5}(:,:,4) * rho_p(1,1) * 4/3 * pi * ((d_p(1,1)./2)*10^(-6))^3 * 10^(9); % nNewton
F_avg_alg_exp_252 = acc_avg_bin{6}(:,:,4) * rho_p(2,1) * 4/3 * pi * ((d_p(2,1)./2)*10^(-6))^3 * 10^(9); % nNewton
F_avg_PE_exp_450 = acc_avg_bin{11}(:,:,4) * rho_p(1,1) * 4/3 * pi * ((d_p(1,1)./2)*10^(-6))^3 * 10^(9); % nNewton
F_avg_alg_exp_450 = acc_avg_bin{12}(:,:,4) * rho_p(2,1) * 4/3 * pi * ((d_p(2,1)./2)*10^(-6))^3 * 10^(9); % nNewton
%% (4) Figure 9: Acting forces on PE particles and alginate beads for 252 and 450 rpm calculated
% by means of Lagrangian acceleration magnitude data
fig = figure(1)
set(0,'DefaultAxesFontName', 'Calibri')
set(0,'DefaultAxesFontSize', 18)
set(gcf,'PaperPositionMode','auto')
set(fig,'units','centimeters','position',[2,2,10,27])
imagesc(imgaussfilt(flipud(F_avg_PE_exp_252),1),[0 1.0]);
hold on
fig_set = gca
set(gca,'XTick',[-5:10:65])
set(gca,'XTickLabel',[0:10:65])
yticks([29:50:228])
yticklabels(fliplr([0:50:228]))
set(gca,'ydir','reverse')
contour(imgaussfilt(flipud(F_avg_PE_exp_252),1),[0:0.2:1.0],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'\textit{F} / nN'};
cl.Label.Interpreter = 'latex';
cl.FontName='Times';
cl.FontSize=22;
ylabel("Reactor height \textit{h} / mm",'FontSize',25,'interpreter','latex');
xlabel('Reactor radius \textit{r} / mm','FontSize',17,'interpreter','latex');
title('\textbf{PE particle}','252 rpm, exp','FontSize',17,'interpreter','latex')
hold off
%print(strcat('C:/','./','F_plot_1_pie'),...
% '-dpng','-r600');
fig = figure(2)
set(0,'DefaultAxesFontName', 'Calibri')
set(0,'DefaultAxesFontSize', 18)
set(gcf,'PaperPositionMode','auto')
set(fig,'units','centimeters','position',[15,2,8.7,27])
imagesc(imgaussfilt(flipud(F_avg_alg_exp_252),1),[0 100]);
hold on
fig_set = gca
set(gca,'XTick',[-5:10:65])
set(gca,'XTickLabel',[0:10:65])
yticks([29:50:228])
yticklabels(fliplr([0:50:228]))
set(gca,'ydir','reverse')
contour(imgaussfilt(flipud(F_avg_alg_exp_252),1),[0:20:100],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'\textit{F} / nN'};
cl.Label.Interpreter = 'latex';
cl.FontName='Times';
cl.FontSize=22;
%ylabel("Reactor height / mm",'FontSize',25,'interpreter','latex');
xlabel('Reactor radius \textit{r} / mm','FontSize',17,'interpreter','latex');
title('\textbf{Alginate beads}','252 rpm, exp','FontSize',17,'interpreter','latex')
hold off
%print(strcat('C:/','./','F_plot_2_pie'),...
% '-dpng','-r600');
fig = figure(3)
set(0,'DefaultAxesFontName', 'Calibri')
set(0,'DefaultAxesFontSize', 18)
set(gcf,'PaperPositionMode','auto')
set(fig,'units','centimeters','position',[20,2,8.7,27])
imagesc(imgaussfilt(flipud(F_avg_PE_exp_450),1),[0 2]);
hold on
fig_set = gca
set(gca,'XTick',[-5:10:65])
set(gca,'XTickLabel',[0:10:65])
yticks([29:50:228]) %an example
yticklabels(fliplr([0:50:228])) %an example
set(gca,'ydir','reverse')
contour(imgaussfilt(flipud(F_avg_PE_exp_450),1),[0:0.4:2],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'\textit{F} / nN'};
cl.Label.Interpreter = 'latex';
cl.FontName='Times';
cl.FontSize=22;
%ylabel("Reactor height / mm",'FontSize',25,'interpreter','latex');
xlabel('Reactor radius \textit{r} / mm','FontSize',17,'interpreter','latex');
title('\textbf{PE particles}','450 rpm, exp','FontSize',17,'interpreter','latex')
hold off
%print(strcat('C:/','./','F_plot_3_pie'),...
% '-dpng','-r600');
fig = figure(4)
set(0,'DefaultAxesFontName', 'Calibri')
set(0,'DefaultAxesFontSize', 18)
set(gcf,'PaperPositionMode','auto')
set(fig,'units','centimeters','position',[30,2,8.7,27])
imagesc(imgaussfilt(flipud(F_avg_alg_exp_450),1),[0 200]);
hold on
fig_set = gca
set(gca,'XTick',[-5:10:65])
set(gca,'XTickLabel',[0:10:65])
yticks([29:50:228]) %an example
yticklabels(fliplr([0:50:228])) %an example
set(gca,'ydir','reverse')
contour(imgaussfilt(flipud(F_avg_alg_exp_450),1),[0:40:200],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'\textit{F} / nN'};
cl.Label.Interpreter = 'latex';
cl.FontName='Times';
cl.FontSize=22;
%ylabel("Reactor height / mm",'FontSize',25,'interpreter','latex');
xlabel('Reactor radius \textit{r} / mm','FontSize',17,'interpreter','latex');
title('\textbf{Alginate beads}','450 rpm, exp','FontSize',17,'interpreter','latex')
hold off
%print(strcat('C:/','./','F_plot_4_pie'),...
% '-dpng','-r600');
%SH%%#%%%%%%%%%%%%%%%%#%%%%%%%####
%%%%%%%%%%%%#%v~~~~~~\%%%#%%%%%%%%
%%%%%%%#%%%%v' ~~~~\%%%%%%%
%%%%#%%%%%%v'dHHb a%%%#%%%%%%
%%%%%%%#%%v'dHHHA :%%%%%%#%%%%
%%%%%#%%%v' VHHHHaadHHb:%#%%%%%%%%
%%%%%%%#v' `VHHHHHHHHb:%%%%%#%%%
%%%#%%%v' `VHHHHHHH:%%%#%%%%%
%%%%%%%' dHHHHHHH:%%#%%%%%%
%%%%#%% dHHHHHHHH:%%%%%%%%%
%%%%%%% dHHHHHHHHH:%%#%%%%%%
%%#%%%% VHHHHHHHHH:%%%%%#%%%
%%%%%%# b HHHHHHHHV:%%%#%%#%%
%%%%%%% Hb HHHHHHHV'%%%%%%%%%%
%%%%#%% HH dHHHHHHV'%%%#%%%%%%%
%%%#%%% VHbdHHHHHHV'#%%%%%%%%%%%
%%%%%#% VHHHHHHHV'%%%%%#%%#%%%%
%%%#%%%% VHHHHHHH:%%#%%#%%#%%#%
%%#%%#%%#%%#%%#%%#%%#%%#%%#%%#%%#%
%%#%%#%%#%%#%%#%%#%%#%%#%%#%%#%IMS
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment