Wrong implementation equation for strain rate in postprocessing

[mgvisakh]

The equation for strain rate given in PIVLab is ∂u/∂x - ∂v/∂y and is referred to in the paper "Stamhuis, E. J. (2006), Aquatic Ecology, 40(4), pp. 463–479.". However, this equation seems inconsistent with the classical definition in the literature (kindly refer to section 3.4 in Fluid Mechanics by Cohen and Kundu (edition 5) or section 2.4 of Viscous Fluid Flow by Papanastasiou or even Wikipedia https://en.wikipedia.org/wiki/Strain-rate_tensor). Strain rate is a tensor quantity, and cannot be described by a single scalar, to begin with.
Strain rate tensor is defined as S = 1/2 (∇u + ∇u^T). It is a symmetric tensor, the off-diagonal elements are the same for 2D case, and is equal to 1/2(∂v/∂x + ∂u/∂y) which is the shear strain.
The trace (sum of the diagonal elements) is the bulk /volumetric strain rate ∇.u = ∂u/∂x + ∂v/∂y for 2D case and is equal to the divergence of the velocity field, which is also the first invariant of the strain rate tensor. The diagonal components of the strain rate tensor are the linear (or extensional) strain rates i.e, for x-direction ∂u/∂x and for y-direction, ∂v/∂y separately. Their sum for 2D flows is the volumetric strain rate for compressible flows and is 0 for 2D incompressible flows. For 3D incompressible flows, this sum is -∂w/∂z which is the extensional rate in the out-of-plane direction.
To summarize, the post-processing variables could be renamed as components of the strain rate tensor and could be as follows:
Shear strain rate = (∂v/∂x + ∂u/∂y)
Linear (extensional) strain rate in x-direction: ∂u/∂x; for y-direction: ∂v/∂y .
Divergence of 2D field = ∇.u = ∂u/∂x + ∂v/∂y
Linear (extensional) strain rate in out-of-plane direction (z) (only for incompressible flows): = -∇.u = -(∂u/∂x + ∂v/∂y)
One can look up these definitions in the above-mentioned references or in section 7.5 of the book on Particle Image Velocimetry by Raffel et. al (3rd edition)

[Shrediquette]

You are suggesting to use the following variables:

• Shear strain rate
• Linear (extensional) strain rate in x-direction
• Linear (extensional) strain rate in y-direction
• Divergence of 2D field
• Linear (extensional) strain rate in out-of-plane direction (z)

However, I only have space for three of these (otherwise I need to do some more major changes).

1: AFAIK divergence should be ok in PIVlab (it uses Matlabs built-in function), and I would keep it as it is.
2: Then I would modify the currently named "simple strain rate" to the shear strain rate.
3: Which of the remaining parameters is most useful to users?

How would you code this in Matlab? The current implementation is as follows:

function out=strain(x,y,u,v)
hx = x(1,:);
hy = y(:,1);
[px, junk] = gradient(u, hx, hy);
[junk, qy] = gradient(v, hx, hy); %#ok<*ASGLU>
out = px-qy;

function out=shear(x,y,u,v)
hx = x(1,:);
hy = y(:,1);
[junk, py] = gradient(u, hx, hy);
[qx, junk] = gradient(v, hx, hy);
out= qx+py;

Thanks for your help and all the best,
William

[Shrediquette]

Changed the calculation in ab2497d
following the proposal of @CSfrom

[Shrediquette]

`function out=shear(x,y,u,v)
%original code until 3.10
%{
hx = x(1,:);
hy = y(:,1);
[junk, py] = gradient(u, hx, hy);
[qx, junk] = gradient(v, hx, hy);
out= qx+py;
%}

%New code, 'Shear rate (magnitude of the rate-of-strain tensor)'
%based on discussion with https://github.com/CSFrom
%x is a 2D matrix containing the x-coordinates of every node
%y is a 2D matrix containing the y-coordinates of every node
%u is a 2D matrix containing the velocities in x direction at every node
%v is a 2D matrix containing the velocities in y direction at every node

dx=x(1,:);
dy=y(:,1);

[dux_dx, dux_dy] = gradient(u, dx, dy);
[duy_dx, duy_dy] = gradient(v, dx, dy);

% where the components of the rate-of-strain tensor components
Dxx = dux_dx;
Dyy = duy_dy;

Dxy = 0.5 * (dux_dy + duy_dx); % Only Dxy since it is symmetric by definition (also the sum of diagonals, Dxy + Dyx = dux_dy + duy_dx).

%Compute the magnitude of the rate-of-strain tensor, i.e.,
out = sqrt(2 * (Dxx.^2 + Dyy.^2 + 2 * Dxy.^2));
`

[alexlib]

`function out=shear(x,y,u,v) %original code until 3.10 %{ hx = x(1,:); hy = y(:,1); [junk, py] = gradient(u, hx, hy); [qx, junk] = gradient(v, hx, hy); out= qx+py; %}

%New code, 'Shear rate (magnitude of the rate-of-strain tensor)' %based on discussion with https://github.com/CSFrom %x is a 2D matrix containing the x-coordinates of every node %y is a 2D matrix containing the y-coordinates of every node %u is a 2D matrix containing the velocities in x direction at every node %v is a 2D matrix containing the velocities in y direction at every node

dx=x(1,:); dy=y(:,1);

[dux_dx, dux_dy] = gradient(u, dx, dy); [duy_dx, duy_dy] = gradient(v, dx, dy);

% where the components of the rate-of-strain tensor components Dxx = dux_dx; Dyy = duy_dy;

Dxy = 0.5 * (dux_dy + duy_dx); % Only Dxy since it is symmetric by definition (also the sum of diagonals, Dxy + Dyx = dux_dy + duy_dx).

%Compute the magnitude of the rate-of-strain tensor, i.e., out = sqrt(2 * (Dxx.^2 + Dyy.^2 + 2 * Dxy.^2)); `

I'm confused by the definitions: shear rate is not the magnitude of the rate-of-strain tensor as far as I know.

[mgvisakh]

As I understand, in PIVLab 3.10, simple strain rate is Dxy and shear rate is used to mention the magnitude of the strain rate tensor. But I think its more appropriate to call Dxy as "shear strain rate" and the latter as just magnitude of strain rate tensor. The square of the latter is related to the dissipation rate through viscosity.
In any case, I think it will be better to show the user the equations used, while selecting it in the postprocessing dropdown menu.