metric_angle_sim.m 4.31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
% ** 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