Commit 391be32d authored by Sebastian Hofmann's avatar Sebastian Hofmann ☝🏽
Browse files

added metric_anle_sim.m for figure 14

parent 6004e4e5
% ** coded by Christian Weiland **
% ** christian.weiland@tuhh.de **
%
% This script "metric_angle_sim.m" creates data for figure 14.
% It calculates the metric d(u,v) and angle phi(u,v) for each particle
% according to equation 30 and equation 29, respectively. These values
% are averaged in selected bins over time.
%
% ADDITIONAL INFORMATION
% 1. The input files must conform to the standard output structure of the
% particle data by the option "Moving Probes" of M-Star CFD, manipulated by
% the .bash script "skriptSortieren.bash".
clear variables
close all
tic
%Radius of the reactor
R = 0.065;
%Height of the reactor
H = 0.25;
disp('Loading data...')
dataidGS = importdata('probeGS.txt');
dataidKL = importdata('probeKL.txt');
disp('... ready with loading data!')
GS = dataidGS.data;
KL = dataidKL.data;
%Particle-ID starts at 0, should be 1
GS(:,1) = GS(:,1) + 1;
KL(:,1) = KL(:,1) + 1;
%Sort
GS = sortrows(GS);
KL = sortrows(KL);
%Program puts (0,0,0) at outer point, should be at the stirrer shaft
% -> linear translation in x- und z direction
GS(:,3) = GS(:,3) - R;
GS(:,5) = GS(:,5) - R;
KL(:,3) = KL(:,3) - R;
KL(:,5) = KL(:,5) - R;
%Calculation of radial position of particles
rGS = sqrt(GS(:,3).^2 + GS(:,5).^2);
rKL = sqrt(KL(:,3).^2 + KL(:,5).^2);
%Radius is 14th entry
GS = [GS rGS];
KL = [KL rKL];
clear data*
%How many particles? Which are the time steps?
schritte = find(GS(:,1)==1, 1, 'last');
numPart = max(GS(:,1));
zeit = GS(1:schritte,2)-GS(1,2);
%Velocity particle
uGS = [GS(:,7) GS(:,8) GS(:,9)];
uKL = [KL(:,7) KL(:,8) KL(:,9)];
%Velocity fluid
vGS = [GS(:,11) GS(:,12) GS(:,13)];
vKL = [KL(:,11) KL(:,12) KL(:,13)];
disp('Calculate the angles...')
%Calculation of angle between particle and fluid velocity (not used for paper)
winkelGS = acos( sum(uGS.*vGS,2) ./ ( sqrt(sum(uGS.*uGS,2)) .* sqrt(sum(vGS.*vGS,2)) ));
winkelKL = acos( sum(uKL.*vKL,2) ./ ( sqrt(sum(uKL.*uKL,2)) .* sqrt(sum(vKL.*vKL,2)) ));
disp('Calculate the velocity metric')
% Modified sign function
vorzeichenGS = sign(sum(uGS.*vGS,2));
vorzeichenGS(vorzeichenGS==0) = 1;
vorzeichenKL = sign(sum(uKL.*vKL,2));
vorzeichenKL(vorzeichenKL==0) = 1;
% Calculation of metric (eq 30) for each particle at each recorded time step
differenzGS = vorzeichenGS.*sqrt(sum((uGS - vGS).*(uGS - vGS),2));
differenzKL = vorzeichenKL.*sqrt(sum((uKL - vKL).*(uKL - vKL),2));
radius = 0:0.001:R;
hoehe = 0:0.001:H;
phiGS = zeros(length(hoehe), length(radius));
phiKL = zeros(length(hoehe), length(radius));
metrikGS = zeros(length(hoehe), length(radius));
metrikKL = zeros(length(hoehe), length(radius));
disp('... now for loop...')
for r = 1:length(radius)-1
for h = 1:length(hoehe)-1
% Assign positions to metric and angles. Average over time.
reihe = find( GS(:,14) >= radius(r) & GS(:,14) < radius(r+1) & GS(:,4) >= hoehe(h) & GS(:,4) < hoehe(h+1) );
phiGS(h, r) = mean(winkelGS(reihe), 'omitnan');
metrikGS(h,r) = mean(differenzGS(reihe), 'omitnan');
reihe = find( KL(:,14) >= radius(r) & KL(:,14) < radius(r+1) & KL(:,4) >= hoehe(h) & KL(:,4) < hoehe(h+1) );
phiKL(h, r) = mean(winkelKL(reihe), 'omitnan');
metrikKL(h,r) = mean(differenzKL(reihe), 'omitnan');
end
disp(['Ready with: radius = ', num2str(radius(r)), ' m'])
end
dauer = toc;
disp(['... ready! The script took ', num2str(dauer), ' s.'])
%CW%%#%%%%%%%%%%%%%%%%#%%%%%%%####
%%%%%%%%%%%%#%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