Find Q-criterion and Lambda2 from velocity values

Hello everyone, hope you are doing great.

I have a matlab code which reads a file that has x, y and z velocity values, as well as corresponding time for each value this is what it looks like:

clear all; clc

[file_list, path_n] = uigetfile('.txt', 'Multiselect', 'on');
filesSorted = natsortfiles(file_list);

if iscell(filesSorted) == 0;
    filesSorted = (filesSorted);
end


for i = 1:length(filesSorted);
    filename = filesSorted{i};
    data = load([path_n filename]);
    
    % Define x and y from the uploaded files
    x =  data (:,1);
    y =  data (:,2);
    
    % Define Time for each value of the uploaded files
    Time(i) = data (1,3);
    Time_tp = Time';
    
    
end

Is there a way to get Lambda2 "λ" and Q-criterion values from the velocity data?

Thank you so much, stay safe

Answer 1

[Martino Pinto]

Hi! Here's a couple of functions that will do it:

function lambda2 = computeL2Criterion(u,v,w,grid_h)
    % Computes the lambda2 criterion from the cartesian velocity 
    % field u,v,w and the grid spacing grid_h

    [du_dx, du_dy, du_dz] = gradient(u,grid_h);
    [dv_dx, dv_dy, dv_dz] = gradient(v,grid_h);
    [dw_dx, dw_dy, dw_dz] = gradient(w,grid_h);
    
    lambda2 = zeros(size(u));
    for i = 1:size(u, 1)
        for j = 1:size(u, 2)
            for k = 1:size(u, 3)
                % Calculate the symmetric and anti-symmetric parts of the velocity gradient tensor
                S = [du_dx(i,j,k), (du_dy(i,j,k) + dv_dx(i,j,k)) / 2, (du_dz(i,j,k) + dw_dx(i,j,k)) / 2;
                     (dv_dx(i,j,k) + du_dy(i,j,k)) / 2, dv_dy(i,j,k), (dv_dz(i,j,k) + dw_dy(i,j,k)) / 2;
                     (dw_dx(i,j,k) + du_dz(i,j,k)) / 2, (dw_dy(i,j,k) + dv_dz(i,j,k)) / 2, dw_dz(i,j,k)];
                
                Omega = [0, (du_dy(i,j,k) - dv_dx(i,j,k)) / 2, (du_dz(i,j,k) - dw_dx(i,j,k)) / 2;
                         (dv_dx(i,j,k) - du_dy(i,j,k)) / 2, 0, (dv_dz(i,j,k) - dw_dy(i,j,k)) / 2;
                         (dw_dx(i,j,k) - du_dz(i,j,k)) / 2, (dw_dy(i,j,k) - dv_dz(i,j,k)) / 2, 0];
                
                M = S*S + Omega*Omega;   % Combine deformation and rotation
                
                eigenvalues = eig(M);
                eigenvalues = sort(eigenvalues, 'descend');
                lambda2(i,j,k) = eigenvalues(2);
            end
        end
    end

end



function Q = computeQCriterion(U, V, W, grid_h)
    % Computes the lambda2 criterion from the cartesian velocity 
    % field u,v,w and the grid spacing grid_h

    [dxU, dyU, dzU] = gradient(U, grid_h);
    [dxV, dyV, dzV] = gradient(V, grid_h);
    [dxW, dyW, dzW] = gradient(W, grid_h);

    S_xx = dxU;
    S_yy = dyV;
    S_zz = dzW;
    S_xy = 0.5 * (dyU + dxV);
    S_xz = 0.5 * (dzU + dxW);
    S_yz = 0.5 * (dzV + dyW);

    Omega_xy = 0.5 * (dyU - dxV);
    Omega_xz = 0.5 * (dzU - dxW);
    Omega_yz = 0.5 * (dzV - dyW);

    Q = 0.5 * ((Omega_xy.^2 + Omega_xz.^2 + Omega_yz.^2) - ...
               (S_xx.^2 + S_yy.^2 + S_zz.^2 + 2*(S_xy.^2 + S_xz.^2 + S_yz.^2)));
end

[William Thielicke]

How would computeQCriterion look like if the data is only 2-dimensional (U and V)?