Skip to content
Snippets Groups Projects
Commit 04bcfd6a authored by Sebastian Hofmann's avatar Sebastian Hofmann :point_up_tone3:
Browse files

add vel_vector_1mm2_bin for St calc and vel_uncert

parent 9988c609
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
data/code/3.3_stokes_number/vel_vectors_1mm2_bin.m 0 → 100644
+ 351
− 0
View file @ 04bcfd6a
% ** coded by Sebastian Hofmann **
% ** s.m.hofmann@mailbox.org ** sebastian.hofmann@tuhh.de **
%
% Creation of data for Figure 6 and Figures 19 and 20, velocity profile over reactor height
% 4 graphs due to respective azimuthal velocity profiles in resp. radius r
% Load .mat created by "loadDat_", experimental or simulation, one after
% another!!
% (1) is done for all data sets, (2) combines all data sets to one .mat file
% (3) first plot, (4) second plot, (5) third plot
%
% INFORMATION
% 1. Define dimensions in l.33-34, define azimuthal positions and
% distances in l.53-63 and enter filename in l.127-128
% 2. Used data her 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!
clc
close all
clear all
%% (1) Load .mat file with separated tracks from pie_piece and calculate velocity and std for resp. area
[FileName,FilePath] = uigetfile('*.mat','Select the .mat file',...
'C:\','MultiSelect', 'off');
load(strcat(FilePath,FileName));
% 1st step - Define dimensions
d_2 = 0.036; % ENTER DIAMETER STIRRER HERE, m
% General calculation for all subs, radial distance
r = sqrt((x).^2+(z).^2); % calculation of all x and z positions to "r"
% Calculate center of bins of respective r-axis subdivision
r_spaces = linspace(5,65,60); % zu kompliziert geschrieben, bins sind "doppelt gerechnet...
[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
% center_edges_r = center_edges_r';
% Calculate center of bins of respective y-axis subdivision
y_spaces = linspace(0,228,228); % zu kompliziert geschrieben, bins sind "doppelt gerechnet...
[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 velocities
u_avg_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1); % Create zero vector for u velocity component
v_avg_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1); % Create zero vector for v velocity component
w_avg_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1); % Create zero vector for w velocity component
v_mag_bin = zeros(size(edges_y,2)-1,size(edges_r,2)-1);
%hbar = parfor_progressbar(size(edges_r,2)-1,'load data...');
m = 0;
for k = 1:size(edges_r,2)-1
m = m + 1;
ur = u(r>edges_r(1,m) & r<edges_r(1,m + 1)); % Look for all velocities in entire r bin
vr = v(r>edges_r(1,m) & r<edges_r(1,m + 1));
wr = w(r>edges_r(1,m) & r<edges_r(1,m + 1));
v_magr = v_mag(r>edges_r(1,m) & r<edges_r(1,m + 1));
y_ur = y(r>edges_r(1,m) & r<edges_r(1,m + 1)); % Look for respective y coordinates to prev. determined velocity components
y_vr = y(r>edges_r(1,m) & r<edges_r(1,m + 1));
y_wr = y(r>edges_r(1,m) & r<edges_r(1,m + 1));
y_v_magr = y(r>edges_r(1,m) & r<edges_r(1,m + 1));
u_avg_bin_r = zeros(size(edges_y,2)-1,1); % Define zero vector for upcoming for loop for avg
v_avg_bin_r = zeros(size(edges_y,2)-1,1);
w_avg_bin_r = zeros(size(edges_y,2)-1,1);
v_mag_avg_bin_r = zeros(size(edges_y,2)-1,1);
u_avg_bin_r_std = zeros(size(edges_y,2)-1,1); % Define zero vector for upcoming for loop for std
v_avg_bin_r_std = zeros(size(edges_y,2)-1,1);
w_avg_bin_r_std = zeros(size(edges_y,2)-1,1);
v_mag_avg_bin_r_std = zeros(size(edges_y,2)-1,1);
j = 0;
for i = 1:size(edges_y,2)-1 % Now, start scanning from bottom to top to corr. defined y bins
j = j + 1;
%[a ~] = find(y_ur(:,1) > edges_y(1,j) & y_ur(:,1) < edges_y(1,j + 1)); % Here: with find function, to write resp. row number in a
a = y_ur(:,1) > edges_y(1,j) & y_ur(:,1) < edges_y(1,j + 1);
u_avg_bin_r(j,1) = mean(ur(a,:),'omitnan'); % Write avg velocity in resp. r bin in resp. y bin, goes with j+1 for entire reactor height
v_avg_bin_r(j,1) = mean(vr(a,:),'omitnan');
w_avg_bin_r(j,1) = mean(wr(a,:),'omitnan');
v_mag_avg_bin_r(j,1) = mean(v_magr(a,:),'omitnan');
% for calc of standard deviation
u_avg_bin_r_std(j,1) = std(ur(a,:),'omitnan'); % Write std in resp. r bin in resp. y bin, goes with j+1 for entire reactor height
v_avg_bin_r_std(j,1) = std(vr(a,:),'omitnan');
w_avg_bin_r_std(j,1) = std(wr(a,:),'omitnan');
v_mag_avg_bin_r_std(j,1) = std(v_magr(a,:),'omitnan');
% for calc of uncertainty determine counts in resp. bin
u_avg_bin_r_std_count(j,1) = length(ur(a,:)) - sum(isnan(ur(a,:)));
v_avg_bin_r_std_count(j,1) = length(vr(a,:)) - sum(isnan(vr(a,:)));
w_avg_bin_r_std_count(j,1) = length(wr(a,:)) - sum(isnan(wr(a,:)));
v_mag_avg_bin_r_std_count(j,1) = length(v_magr(a,:)) - sum(isnan(v_magr(a,:)));
% calc of uncertainty with std and count
u_avg_bin_r_delta(j,1) = u_avg_bin_r_std(j,1) ./ sqrt(u_avg_bin_r_std_count(j,1));
v_avg_bin_r_delta(j,1) = v_avg_bin_r_std(j,1) ./ sqrt(v_avg_bin_r_std_count(j,1));
w_avg_bin_r_delta(j,1) = w_avg_bin_r_std(j,1) ./ sqrt(w_avg_bin_r_std_count(j,1));
v_mag_avg_bin_r_delta(j,1) = v_mag_avg_bin_r_std(j,1) ./ sqrt(v_mag_avg_bin_r_std_count(j,1));
end
u_avg_bin(:,m) = u_avg_bin_r; % Write entire created avg column for r bin in prepared zero vector, goes with m+1
v_avg_bin(:,m) = v_avg_bin_r;
w_avg_bin(:,m) = w_avg_bin_r;
v_mag_avg_bin(:,m) = v_mag_avg_bin_r;
u_avg_bin_std(:,m) = u_avg_bin_r_std; % Write entire created std column for r bin in prepared zero vector, goes with m+1
v_avg_bin_std(:,m) = v_avg_bin_r_std;
w_avg_bin_std(:,m) = w_avg_bin_r_std;
v_mag_avg_bin_std(:,m) = v_mag_avg_bin_r_std;
u_avg_bin_std_count(:,m) = u_avg_bin_r_std_count;
v_avg_bin_std_count(:,m) = v_avg_bin_r_std_count;
w_avg_bin_std_count(:,m) = w_avg_bin_r_std_count;
v_mag_avg_bin_std_count(:,m) = v_mag_avg_bin_r_std_count;
u_avg_bin_delta(:,m) = u_avg_bin_r_delta;
v_avg_bin_delta(:,m) = v_avg_bin_r_delta;
w_avg_bin_delta(:,m) = w_avg_bin_r_delta;
v_mag_avg_bin_delta(:,m) = v_mag_avg_bin_r_delta;
%hbar.iterate(1); % update progress by one iteration
end
%close(hbar);
% Save results for dataset
fname = sprintf('vel_vectors_all_1mm2_sim_200LX_1024_740_252rpm_pie_uncert'); % Enter filename here!
%fname = sprintf('vel_vectors_all_1mm2_exp_1024_740_252rpm_pie_uncert'); % Enter filename here!
fpath = sprintf('C:/'); % Enter filepath here!
save(strcat(fpath,'./',fname,'.mat'),...
'u_avg_bin','v_avg_bin','w_avg_bin','v_mag_avg_bin',...
'u_avg_bin_std_count','v_avg_bin_std_count','w_avg_bin_std_count','v_mag_avg_bin_std_count',...
'u_avg_bin_delta','v_avg_bin_delta','w_avg_bin_delta','v_mag_avg_bin_delta',...
'edges_r','edges_y')
fprintf('Done saving!\n')
%% Combine all sets and save again all together (for later Stokes evaluation)
clear all
%n..1-6 252 rpm: sim flowtracer - sim 180 1000 - sim 740 1024 - exp 180 1000 - exp 740 1024
%n..7-12 450 rpm: sim flowtracer - sim 180 1000 - 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;
vel_avg_bin{n} = cat(4,[u_avg_bin],[v_avg_bin],[w_avg_bin],[v_mag_avg_bin]);
vel_avg_bin_std{n} = cat(4,[u_avg_bin_std],[v_avg_bin_std],[w_avg_bin_std],[v_mag_avg_bin_std]);
end
fname = sprintf('graph_vel_vectors_all_1mm2_252rpm_450rpm_pie'); % Enter filename here!
fpath = sprintf('C:/'); % Enter filepath here!
save(strcat(fpath,'./',fname,'.mat'),...
'vel_avg_bin','vel_avg_bin_std',...
'edges_r','edges_y')
%% Combine all sets and save again all together (just for error analysis)
%n..1-2 252 rpm: exp 180 1000 - exp 740 1024
n = 0;
for i=1:2
[FileName,FilePath] = uigetfile('*.mat','Select the .mat file', ...
'C:/','MultiSelect', 'off');
load(strcat(FilePath,FileName));
n = n+1;
vel_avg_bin{n} = cat(4,[u_avg_bin],[v_avg_bin],[w_avg_bin],[v_mag_avg_bin]);
vel_avg_bin_std_count{n} = cat(4,[u_avg_bin_std_count],[v_avg_bin_std_count],[w_avg_bin_std_count],[v_mag_avg_bin_std_count]);
vel_avg_bin_delta{n} = cat(4,[u_avg_bin_delta],[v_avg_bin_delta],[w_avg_bin_delta],[v_mag_avg_bin_delta]);
end
fname = sprintf('graph_vel_vectors_all_1mm2_252rpm_pie_uncert'); % Enter filename here!
fpath = sprintf('C:/'); % Enter filepath here!
save(strcat(fpath,'./',fname,'.mat'),...
'vel_avg_bin','vel_avg_bin_std_count','vel_avg_bin_delta',...
'edges_r','edges_y')
%% Plot 1 - PE part, exp
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(vel_avg_bin{1}(:,:,4)),1),[0 0.4]);
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(vel_avg_bin{1}(:,:,4)),1),[0:0.05:0.4],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'$v_\mathrm{mag,\, PE}$ / $m \, s^{-1}$'};
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 particles}';'\textit{Velocity magnitude}';'252 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_PE_part_exp_252_vmag'),...
% '-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(vel_avg_bin_delta{1}(:,:,4)),1),[0 0.01]);
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(vel_avg_bin_delta{1}(:,:,4)),1),[0:0.002:0.01],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'$\Delta_{\mathrm{PE}}$ / $m \, s^{-1}$'};
cl.Label.Interpreter = 'latex';
cl.FontName='Times';
cl.FontSize=22;
xlabel('Reactor radius \textit{r} / mm','FontSize',17,'interpreter','latex');
title({'\textbf{PE particles}';'\textit{Velocity magnitude uncertainty}';'252 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_PE_part_exp_252_vmag_delta'),...
% '-dpng','-r600');
%% Plot 2 - alg beads, exp
fig = figure(1)
set(0,'DefaultAxesFontName', 'Calibri')
set(0,'DefaultAxesFontSize', 18)
set(gcf,'PaperPositionMode','auto')
set(fig,'units','centimeters','position',[2,2,8.7,27])
imagesc(imgaussfilt(flipud(vel_avg_bin{2}(:,:,4)),1),[0 0.4]);
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(vel_avg_bin{2}(:,:,4)),1),[0:0.05:0.4],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'$v_\mathrm{mag, \, alg}$ / $m \, s^{-1}$'};
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{Alginate beads}';'\textit{Velocity magnitude}';'252 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_alg_beads_exp_252_vmag'),...
% '-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(vel_avg_bin_delta{2}(:,:,4)),1),[0 0.01]);
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(vel_avg_bin_delta{2}(:,:,4)),1),[0:0.002:0.01],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Label.String = {'$\Delta_{\mathrm{alg}}$ / $m \, s^{-1}$'};
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}';'\textit{Velocity magnitude uncertainty}';'252 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_alg_beads_exp_252_vmag_delta'),...
% '-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
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