rotate_gradu.tex

\documentclass[12pt]{article}
\usepackage{amssymb,latexsym,amsmath,bm}

\newif\ifpdf
\ifx\pdfoutput\undefined
\pdffalse % we are not running PDFLaTeX
\else
\pdfoutput=1 % we are running PDFLaTeX
\pdftrue
\fi
\ifpdf
\usepackage[pdftex]{graphicx}
\else
\usepackage{graphicx}
\fi
\ifpdf
\DeclareGraphicsExtensions{.pdf, .jpg, .tif}
\else
\DeclareGraphicsExtensions{.eps, .jpg}
\fi



\textwidth = 6.5 in
\textheight = 9 in
\oddsidemargin = 0.0 in
\evensidemargin = 0.0 in
\topmargin = 0.0 in
\headheight = 0.0 in
\headsep = 0.0 in
\parskip = 0.2 in
\parindent = 0.0 in


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\newcommand{\khat}{\hat{\mathbf k}}
\newcommand{\uv}{\mathbf u}
\newcommand{\up}{\mathbf u'}
\newcommand{\w}{\mathbf w}
\newcommand{\grad}{\nabla}
\newcommand{\curlp}{\gradp \times}
\newcommand{\curl}{\grad \times}
\newcommand{\gradp}{\nabla'}

\newcommand{\commentout}[1]

\DeclareMathOperator{\Span}{span}

\title{Rotation matrix for transforming the strain tensor}

\begin{document}
%\maketitle
\section*{Shear tensor}
There are 9 components of the shear (strain rate) $S_{ij} = \grad {\bm u} =
\partial_j u_i$. We can compute all these components in the
grid-coordinate system ($\bm{{\hat x}}_1,\bm{{\hat x}}_2, \bm{{\hat
x}}_3$). We would like to be able to say what these components are in
a coordinate system whose $x_1$-axis is along $\bm{{\hat r}}$, one of
the 73 directions we compute in our angle averaging procedure.

\section*{New coordinate system}
Let $\bm{{\hat r}}$, $\bm{{\hat t}}_1$
and $\bm{{\hat t}}_2$ (which are determined by the existing code in the grid-coordinate system) form our new coordinate system.
The axis of the new system $\bm{{\hat x}}_1',\bm{{\hat x}}_2',
\bm{{\hat x}}_3'$ are obviously 
given in the transformed coordinate system as $\bm{{\hat x}}_1' =
\left(\begin{array}{c} 1 \\ 0 \\0
\end{array}\right)$, $\bm{{\hat x}}_2' = \left(\begin{array}{c} 0 \\ 1
\\0 \end{array}\right)$ and $\bm{{\hat x}}_3' = \left(\begin{array}{c}
0 \\ 0 \\1 \end{array}\right)$.

The transformation between the two coordinate systems is governed by
the rotation matrix ${\bf A}$ which can be written in terms of the
known components of the new coordinate system. That is, we have 

$$
\bm{{\hat x}}_1' = \left(\begin{array}{c} 1 \\ 0 \\0\end{array}\right) = {\bf A}\left(\begin{array}{c} \bm{{\hat r}}\cdot \bm{{\hat x}}_1 \\ \bm{{\hat r}}\cdot \bm{{\hat x}}_2 \\\bm{{\hat r}}\cdot \bm{{\hat x}}_3\end{array}\right);~~
\bm{{\hat x}}_2' = \left(\begin{array}{c} 0 \\ 1 \\0\end{array}\right) = {\bf A}\left(\begin{array}{c}\bm{{\hat t}}_1\cdot \bm{{\hat x}}_1 \\ \bm{{\hat t}}_1 \cdot\bm{{\hat x}}_2 \\\bm{{\hat t}}_1 \cdot \bm{{\hat x}}_3 \end{array}\right);~~
\bm{{\hat x}}_3' = \left(\begin{array}{c} 0 \\ 0 \\1 \end{array}\right)= {\bf A}\left(\begin{array}{c}\bm{{\hat t}}_2\cdot \bm{{\hat x}}_1 \\ \bm{{\hat t}}_2 \cdot\bm{{\hat x}}_2 \\\bm{{\hat t}}_2 \cdot \bm{{\hat x}}_3 \end{array}\right)
$$

From these relations we can immediately write down 
${\bf A}^{-1}$.  Since the vectors are ortho-normal, 
${\bf A A}^{T} = {\bf I}$, so ${\bf A} = {\bf A}^{-T}$ and thus
\begin{equation}
{\bf A} =  \left(\begin{array}{ccc}  
\bm{{\hat r}}\cdot \bm{{\hat x}}_1 & \bm{{\hat r}}\cdot \bm{{\hat x}}_2 &\bm{{\hat r}}\cdot \bm{{\hat x}}_3 \\
\bm{{\hat t}}_1\cdot \bm{{\hat x}}_1 & \bm{{\hat t}}_1 \cdot\bm{{\hat x}}_2 & \bm{{\hat t}}_1 \cdot \bm{{\hat x}}_3 \\
\bm{{\hat t}}_2\cdot \bm{{\hat x}}_1 & \bm{{\hat t}}_2 \cdot\bm{{\hat x}}_2 & \bm{{\hat t}}_2 \cdot \bm{{\hat x}}_3 
\end{array} 
\right)
\end{equation}


So the strain tensor in the new coordinate system is  given by 
${\bf S}' = {\bf A S A}^{T}$, or component-wise
$S_{ij}' = A_{ik} S_{kl} A_{jl}$.

We cycle through the 73 different $\bm{{\hat r}}$ in this way to find the one 
which give the largest off-diagonal
contribution to $S_{1j}'$.  That is, we find the 
direction with the largest value of $S_{12}'^2 + S_{13}'^2$.  

We then construct the {\em tilde} coordinate system which consists
of $\bm{{\hat r}}$ and new tangent vectors $( \bm{{\tilde t}}_1, 
\bm{{\tilde t}}_2 )$ given by a rotation of $( \bm{{\hat t}}_1, \bm{{\hat t}}_2 )$
about the $\bm{{\hat r}}$ axis.  The rotation is chosen so that the 
tangential derivative of $u_1'$ is entirely in the $\bm{{\tilde t}}_1$
direction.  The strain tensor in this coordinate system
is given by ${\bf \tilde S}$, and the velocity vector by 
$( \tilde u_1, \tilde u_2, \tilde u_3)$. 

The tilde-tangent vectors are given by 

\begin{equation}
\bm{{\tilde t}}_1 = \frac{S_{12}'  \bm{{\hat t}}_1 +  S_{13}'  \bm{{\hat t}}_2 }
                         {\sqrt{ S_{12}'^2 +  S_{13}'^2 } } 
\qquad
\bm{{\tilde t}}_2 = \bm{{\hat r}} \times \bm{{\tilde t}}_1
\end{equation}
so that
\begin{equation}
\tilde S_{12} = 
\begin{pmatrix}  S_{11}' \\
                 S_{12}' \\
                 S_{13}' 
\end{pmatrix}
\cdot \bm{{\tilde t}}_1 = 
 \sqrt{ S_{12}'^2 + S_{13}'^2 }
\end{equation}
and 
\begin{equation}
\tilde S_{13} = 
\begin{pmatrix}  S_{11}' \\
                 S_{12}' \\
                 S_{13}' 
\end{pmatrix}
\cdot \bm{{\tilde t}}_2 =  0
\end{equation}
THE ABOVE TWO EQUATIONS (3-4) ARE WRONG. SEE KURIEN SECTION.

We then compute the mixed structure function $\langle \Delta \tilde u_1 (r) 
\Delta \tilde u_2 (r)\rangle$ and 
$\langle \Delta \tilde u_1 (r) \Delta \tilde u_3 (r)\rangle$ for that direction of $\bm{{\hat r}}$.

\section*{Redundancies}

In the above procedure, we looked for the direction $ \bm{{\hat r}}$
which gave the maximum value of $S_{12}'^2 + S_{13}'^2$.  
Is this sufficient, or should we also consider the other
components $S_{21}'^2 + S_{23}'^2$ and 
$S_{31}'^2 + S_{32}'^2$?  

In isoave.F90, we construct the set of vectors $\bm{{\hat r}}$ 
by taking all directions expressible as a vector with integer
coefficients of norm less than $\sqrt{11}$. Call this set
of directions $V$.   If we consider all the components mentioned
above, this amounts to maximizing $S_{12}'^2 + S_{13}'^2$ over
all directions in $V$ as well as all directions given by
all of the $\bm{{\hat t}}_1$ and $\bm{{\hat t}}_2$ vectors.

In isoave.F90, from each $\bm{{\hat r}}$ we construct its associated
$\bm{{\hat t}}_1$ by permutation of indices of $\bm{{\hat r}}$.
This means that $\bm{{\hat t}}_1$ is in a direction given by integer coefficients
with norm less than $\sqrt{11}$, and so $\bm{{\hat t}}_1 \in V$.  
Thus it is redundant to maximize
both $S_{12}'^2 + S_{13}'^2$ and $S_{21}'^2 + S_{23}'^2$.

However, $\bm{{\hat t}}_2 = \bm{{\hat r}} \times \bm{{\hat t}}_1$.
Most of these vectors are also in $V$, but there are 32 of them 
which are not in $V$.  Thus it is not redundant to maximize over
both $S_{12}'^2 + S_{13}'^2$ and $S_{31}'^2 + S_{32}'^2$.
However, to compute the needed structure functions associated with
a direction which had a maximum value of $S_{31}'^2 + S_{32}'^2$, we
would need to compute increments in the $\bm{{\hat t}}_2$ direction.
This direction is not in $V$, and our code isoave.F90 can only 
compute increments in directions in $V$.  Computing increments in
other directions would be very expensive, requiring much more communication
and interpolation to non-grid point locations.

In summary, we are thus capable of computing the mean shear in
the 73 directions in $V$, as well as a few extra directions given by the
$\bm{{\hat t}}_2$ vectors.   But we are only capable of computing
structure functions for increments in directions in $V$, and thus
we limit our search to maximizing $S_{31}'^2 + S_{32}'^2$ over
all directions in $V$.  


\section*{Kurien's approach}

The above definition of $\bm{{\tilde t}}_1$ is confusing to me because
it is seems to be a definition which mixes two coordinate systems. On
the one hand, $S'_{1j}$ is the $j$-component (in the coordinate system
where $\hat{\bm{ r}} \rightarrow 1$-component, $\hat{\bm{t}}_1
\rightarrow 2$-component and $\hat{\bm{ t}}_2 \rightarrow
3$-component) of the gradient of the velocity along $\hat{\bm{
r}}$. In general $S'_{1j}$ is a three component vector expressed in
the ($\hat{\bm{ r}}$, $\hat{\bm{t}_1}$, $\hat{\bm{t}_2}$) coordinate
system. On the other hand $\hat{\bm{t}}_1$ and $\hat{\bm{t}}_2$ are
expressed (in the code) in the grid-coordinate system.

We must first find what the $vector$ $S'_{1j}$ looks like in the
grid-coordinate system. So it seems to me, that we should first
express the vector $S'_{1j}$ in the grid-coordinate system. Denote
this representation by (superscript $g$ for $grid$)

\begin{equation}S^g_{1j}  =  A_{kj} S'_{1k}\end{equation}

NOTE that $S^g_{1j} \neq S_{1j}$(!), the latter is the $j$-component
of the shear in the velocity along the grid-coordinate axis
$\bm{x}_1$.  

Now, we can proceed with defining the tilde transverse coordinates as before but now using the $S^g_{1j}$ components.
 
The tilde-tangent vectors are given by 
\begin{equation}
\bm{{\tilde t}}_1 = \frac{(S^g_{1}\cdot \hat{\bm{t}}_1) \bm{{\hat t}}_1 +  (S^g_{1}\cdot\hat{\bm{t}}_2) \bm{{\hat t}}_2 }
                         {\sqrt{ (S^g_{1}\cdot \hat{\bm{t}}_1)^2 +  (S^g_{1}\cdot\hat{\bm{t}}_2)^2 }}
\qquad
\bm{{\tilde t}}_2 = \bm{{\hat r}} \times \bm{{\tilde t}}_1
\end{equation}

The following is the relationship between $S_1^g$ and $S_1'$:
\begin{eqnarray}
S^g_{1}\cdot \hat{\bm {r}} = S'_{11}\\
S^g_{1}\cdot \hat{\bm {t}}_1 = S'_{12}\\
S^g_{1}\cdot \hat{\bm {t}}_2 = S'_{13}
\end{eqnarray}

Therefore in the tilde coordinate system we have 
\begin{equation}
\tilde S_{12} =  S^g_{1} \cdot \bm{{\tilde t}}_1 = 
 \sqrt{ S_{12}'^2 + S_{13}'^2 } 
\end{equation}
and 
\begin{equation}
\tilde S_{13} = S^g_{1}\cdot \bm{{\tilde t}}_2 =  0
\end{equation}



We then compute the mixed structure function $\langle \Delta \tilde u_1 (r) 
\Delta \tilde u_2 (r)\rangle$ and 
$\langle \Delta \tilde u_1 (r) \Delta \tilde u_3 (r)\rangle$ for that direction of $\bm{{\hat r}}$.

\section*{Implementation in the code}
Since we are only looking for maximization of $S_{12}$ in the
maximal shear direction, all we need to compute is the following:

For a given direction $\hat{\bm{r}}$
\begin{equation}
{S}_{12}' = \nabla (\hat {\bm{r}} \cdot \bm{u}) \cdot \hat{\bm{t}}_1
= \hat {\bm{r}}^T  S \cdot \hat{\bm{t}}_1 
= A_{1k} S_{kj} \cdot \hat{\bm{ t}}_1 
\end{equation}
and 
\begin{equation}
{S}_{13}' = \nabla (\hat {\bm{r}} \cdot \bm{u}) \cdot \hat{\bm{t}}_2
= \hat {\bm{r}}^T  S \cdot \hat{\bm{t}}_2 
= A_{1k} S_{kj} \cdot \hat{\bm{ t}}_2 
\end{equation}

Cycle through the 73 directions to find the one for which ${S}_{12}'^2
+ {S}_{13}'^2$ is maximum. Then for this shear tensor in this maximal
strain directions, denoted by
${\bf S}_{max}'$, find the orientation of the tranverse directions (see
definition in previous section) which orient all the strain in the
$(1,2)$ component of the strain,

\begin{equation}
\bm{{\tilde t}}_1 = \frac{S_{max}'_{12} \bm{{\hat t}}_1 +  S_{max}'_{13} \bm{{\hat t}}_2 }
                         {\sqrt{S_{max}'_{12}^2 +  S_{max}'_{13}^2 }}
\qquad
\bm{{\tilde t}}_2 = \bm{{\hat r}} \times \bm{{\tilde t}}_1
\end{equation}

such that (this is the check for correctness):
\begin{equation}
{\tilde S}_{12} = \nabla (\hat{\bm{r}} \cdot \bm{u}) \cdot \tilde{\bm{t}}_1
= \sqrt{S_{max}'_{12}^2 +  S_{max}'_{13}^2 }
\end{equation}
and
\begin{equation}
{\tilde S}_{13} = \nabla (\hat{\bm{r}} \cdot \bm{u}) \cdot \tilde{\bm{t}}_2
= 0
\end{equation}

Return the maximum direction index, $\bm{{\tilde t}}_1$ and $\bm{{\tilde t}}_2$.

\commentout{
Comment by Mark:
  I think there is a mistake in your formula for $\tilde S_{12}$.

\begin{equation}
{\tilde S}_{12} = \nabla (\hat r \cdot u) \cdot {\tilde t}_1
= \hat r^T  S \cdot {\tilde t}_1 
= A_{1k} S_{kj} \cdot {\tilde t}_1 
\end{equation}
Now take
\[
   {\tilde t}_1  = \alpha {\hat t}_1 + \beta {\hat t}_1
\]
\[
   {\tilde t}_2  = \beta {\hat t}_1 - \alpha {\hat t}_1
\]
so that
\begin{equation}
A_{1k} S_{kj} \cdot {\tilde t}_1  = 
\alpha A_{1k} S_{kj} \cdot {\hat t}_1  + 
\beta A_{1k} S_{kj} \cdot {\hat t}_2
\end{equation}
\begin{equation}
 = \alpha S'_{12}  + \beta S'_{13} 
\end{equation}
Similarly,
\[
{\tilde S}_{13} = A_{1k} S_{kj} \cdot {\tilde t}_2  = 
\beta S'_{12}  - \alpha S'_{13} 
\]
Thus, to summerize, 
\[
{\tilde S}_{12} = \alpha S'_{12}  + \beta S'_{13} 
\]
\[
{\tilde S}_{13} = \beta S'_{12}  - \alpha S'_{13} 
\]
and thus we want to take
\[
  \alpha = S'_{12}  \qquad \beta = S'_{13}
\]


Response by Susan:
   I think there is a mistake in going from (9) to (10). I agree completely with (9). But according to (9), shouldn't 
\begin{equation}S^g_{1j} = A_{1k} S_{kj}\quad ?\end{equation}
Eq. (11) which comes from your equation (9), is not the same as
(10). But I checked and (11) IS the same as (5) (since ${\bf A}$ is
orthogonal). That is,
\begin{eqnarray}S^g_{1j} &=& A_{kj} S'_{1k} \quad \mbox{(Kurien definition Eq (5))} \\
&=& A_{kj} A_{1p} S_{pq} A_{kq} \\
&=& \delta_{jq} A_{1p} S_{pq} \quad \mbox{(orthogonality $\Rightarrow A_{kj}A_{kq} == \delta_{jq}$)} \\ 
&=& A_{1p} S_{pj} \quad \mbox{(Taylor definition from Eq (9) (==Eq (11)))} 
\end{eqnarray}
}
\end{document}