Commit 4c3efd40 authored by Sebastian Hofmann's avatar Sebastian Hofmann ☝🏽
Browse files

updated vel_vector_1mm2_bin

parent 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
% This calculation serves two purposes:
% (1) Firstly: Conversion of Lagrangian velocity data from 4D-PTV to Eulerian data
% in 1 mm x 1 mm bin (via azimuthal data) to calculate St number for characteristic 2D contour plot
% (2) Secondly: Creation of data for Figures 23 and 24 to show 2D contour plot of velocity
% magnitude and corresponding uncertainty throughout the reactor (2x 4 graphs, 252 and 450 rpm)
%
% Load .mat created by "pie_piece_from_loadDat_exp_sim", 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
% (1) is done for all data sets, (2) combines all data sets to one .mat file for later St calculation
% (3) combines just relevant "error" plot data into one (see Figures 23 & 24),
% (4) plot for Figure 23 and (5) plot for Figure 24
%
% 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. Enter filename in l.154-155
% 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
......@@ -23,21 +27,20 @@ clc
close all
clear all
%% (1) Load .mat file with separated tracks from pie_piece and calculate velocity and std for resp. area
%% (1) Load .mat file with separated tracks from pie_piece and calculate velocity, std and uncertainty for resp. bin
[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"
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
r_spaces = linspace(5,65,60); % reactor radius total 130mm; half is 65 mm, but inner 5 mm
% skipped due to shaft obscuration-> 60 mm,
% division into 60 bins to gain 1 mm grid size
[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);
......@@ -46,11 +49,10 @@ for i = 1:size(edges_r,2)-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
y_spaces = linspace(0,228,228); % reactor height with 228 mm ->
% division into 228 bins to gain 1 mm grid height
[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);
......@@ -110,13 +112,13 @@ for k = 1:size(edges_r,2)-1
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,:)));
u_avg_bin_r_std_count(j,1) = length(ur(a,:)) - sum(isnan(ur(a,:))); % Write number of counts in resp. r bin in resp. y bin, for entire reactor height
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));
% calc of uncertainty with std and sqrt(count)
u_avg_bin_r_delta(j,1) = u_avg_bin_r_std(j,1) ./ sqrt(u_avg_bin_r_std_count(j,1)); % Write uncertainty in resp. r bin in resp. y bin, for entire reactor height
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));
......@@ -128,17 +130,17 @@ for k = 1:size(edges_r,2)-1
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
u_avg_bin_std(:,m) = u_avg_bin_r_std; % Write entire created std column for r bins, 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;
u_avg_bin_std_count(:,m) = u_avg_bin_r_std_count; % Write entire created "count" column for r bins, goes with m+1
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;
u_avg_bin_delta(:,m) = u_avg_bin_r_delta; % Write entire created uncertainty column for r bins, goes with m+1
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;
......@@ -149,9 +151,9 @@ 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!
%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_450rpm_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',...
......@@ -162,7 +164,7 @@ save(strcat(fpath,'./',fname,'.mat'),...
fprintf('Done saving!\n')
%% Combine all sets and save again all together (for later Stokes evaluation)
%% (2) Combine all single saved sets and save them again together (for later Stokes calculation)
clear all
......@@ -188,13 +190,16 @@ save(strcat(fpath,'./',fname,'.mat'),...
'edges_r','edges_y')
%% Combine all sets and save again all together (just for error analysis)
%% (3) From here on: Only for figures for appendix for error analysis - combine all sets and save again all together
%n..1-2 252 rpm: exp 180 1000 - exp 740 1024
%n..3-4 450 rpm: exp 180 1000 - exp 740 1024
n = 0;
for i=1:2
for i=1:4
[FileName,FilePath] = uigetfile('*.mat','Select the .mat file', ...
'C:/','MultiSelect', 'off');
load(strcat(FilePath,FileName));
......@@ -213,7 +218,9 @@ save(strcat(fpath,'./',fname,'.mat'),...
'edges_r','edges_y')
%% Plot 1 - PE part, exp
%% (4) Figure 23: Experimentally determined velocity magnitudes v_mag and corresponding uncertainty \Delta for pie_piece for
% PE particles and alginate beads for 252 rpm
% Plot 1.1 - PE parts, exp, 252 rpm
fig = figure(1)
set(0,'DefaultAxesFontName', 'Calibri')
......@@ -270,7 +277,7 @@ contour(imgaussfilt(flipud(vel_avg_bin_delta{1}(:,:,4)),1),[0:0.002:0.01],'ShowT
% '-dpng','-r600');
%% Plot 2 - alg beads, exp
% Plot 1.2 - alg beads, exp, 252 rpm
fig = figure(1)
set(0,'DefaultAxesFontName', 'Calibri')
......@@ -327,6 +334,123 @@ contour(imgaussfilt(flipud(vel_avg_bin_delta{2}(:,:,4)),1),[0:0.002:0.01],'ShowT
% '-dpng','-r600');
%% (5) Figure 24: Experimentally determined velocity magnitudes v_mag and corresponding uncertainty \Delta for pie_piece for
% PE particles and alginate beads for 450 rpm
% Plot 2.1 - PE part, exp, 450 rpm
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{3}(:,:,4)),1),[0 0.8]);
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{3}(:,:,4)),1),[0:0.1:0.8],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Ticks = [0,0.4,0.8]
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}';'450 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_PE_part_exp_450_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{3}(:,:,4)),1),[0 0.04]);
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{3}(:,:,4)),1),[0:0.005:0.04],'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}';'450 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_PE_part_exp_450_vmag_delta'),...
% '-dpng','-r600');
% Plot 2.2 - alg beads, exp, 450 rpm
fig = figure(3)
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{4}(:,:,4)),1),[0 0.8]);
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{4}(:,:,4)),1),[0:0.1:0.8],'ShowText','on','LineStyle','-','Color','r')
cl = colorbar('Southoutside');
cl.Ticks = [0,0.4,0.8]
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}';'450 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_alg_beads_exp_450_vmag'),...
% '-dpng','-r600');
fig = figure(4)
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{4}(:,:,4)),1),[0 0.04]);
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{4}(:,:,4)),1),[0:0.005:0.04],'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}';'450 rpm, exp'},'FontSize',17,'interpreter','latex')
hold off
% print(strcat('C:/','./','uncert_alg_beads_exp_450_vmag_delta'),...
% '-dpng','-r600');
%SH%%#%%%%%%%%%%%%%%%%#%%%%%%%####
......
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