Nonlinear finite elements/Kinematics - polar decomposition

Polar decomposition
The Polar decomposition theorem states that any second order tensor whose determinant is positive can be decomposed  uniquely into a symmetric part and an orthogonal part.

In continuum mechanics, the deformation gradient $$\boldsymbol{F}$$ is such a tensor because $$\det(\mathbf{F}) > 0$$. Therefore we can write

\boldsymbol{F} = \boldsymbol{R}\cdot\boldsymbol{U} = \boldsymbol{V}\cdot\boldsymbol{R} $$ where $$\boldsymbol{R}$$ is an orthogonal tensor ($$\boldsymbol{R}\cdot\boldsymbol{R}^T = \boldsymbol{\mathit{1}}$$) and $$\boldsymbol{U}, \boldsymbol{V}$$ are symmetric tensors ($$\boldsymbol{U} = \boldsymbol{U}^T$$ and $$\boldsymbol{V} = \boldsymbol{V}^T$$) called the  right stretch tensor and the  left stretch tensor, respectively. This decomposition is called the polar decomposition of $$\boldsymbol{F}$$.

Recall that the  right Cauchy-Green deformation tensor is defined as

\boldsymbol{C} = \boldsymbol{F}^T\cdot\boldsymbol{F} $$ Clearly this is a symmetric tensor. From the polar decomposition of $$\boldsymbol{F}$$ we have

\boldsymbol{C} = \boldsymbol{U}^T\cdot\boldsymbol{R}^T\cdot\boldsymbol{R}\cdot\boldsymbol{U} = \boldsymbol{U}\cdot\boldsymbol{U} = \boldsymbol{U}^2 $$ If you know $$\boldsymbol{C}$$ then you can calculate $$\boldsymbol{U}$$ and hence $$\boldsymbol{R}$$ using $$\boldsymbol{R} = \boldsymbol{F}\cdot\boldsymbol{U}^{-1}$$.

How do you find the square root of a tensor?
If you want to find $$\boldsymbol{U}$$ given $$\boldsymbol{C}$$ you will need to take the square root of $$\boldsymbol{C}$$. How does one do that?

We use what is called the  spectral decomposition or  eigenprojection of $$\boldsymbol{C}$$. The spectral decomposition involves expressing $$\boldsymbol{C}$$ in terms of its eigenvalues and eigenvectors. The tensor product of the eigenvectors acts as a basis while the eigenvalues give the magnitude of the projection.

Thus,

\boldsymbol{C} = \sum_{i=1}^3 \lambda_i^2~ \boldsymbol{N}_i\otimes\boldsymbol{N}_i $$ where $$\lambda_i^2$$ are the principal values (eigenvalues) of $$\boldsymbol{C}$$ and $$\boldsymbol{N}_i$$ are the principal directions (eigenvectors) of $$\boldsymbol{C}$$.

Therefore,

\boldsymbol{U}^2 = \sum_{i=1}^3 \lambda_i^2~ \boldsymbol{N}_i\otimes\boldsymbol{N}_i $$ Since the basis does not change, we then have

\boldsymbol{U} = \sum_{i=1}^3 \lambda_i~ \boldsymbol{N}_i\otimes\boldsymbol{N}_i $$ Therefore the $$\lambda_i$$ can be interpreted as principal stretches and the vectors $$\boldsymbol{N}_i$$ are the directions of the principal stretches.

Exercise:
If

\boldsymbol{U} = \sum_{i=1}^3 \lambda_i~ \boldsymbol{N}_i\otimes\boldsymbol{N}_i $$ show that

\boldsymbol{U}^2 = \boldsymbol{U}\cdot\boldsymbol{U} = \sum_{i=1}^3 \lambda_i^2~ \boldsymbol{N}_i\otimes\boldsymbol{N}_i ~. $$

Example of polar decomposition
Let us assume that the motion is given by

\begin{align} x_1 &= \cfrac{1}{4} \left[4~X_1 + (9 - 3~X_1 - 5~X_2 - X_1~X_2)~t\right] \\ x_2 &= X_2 + (4 + 2~X_1)~t \end{align} $$ The adjacent figure shows how a unit square subjected to this motion evolves over time.

Deformation gradient
The deformation gradient is given by

\boldsymbol{F} = \frac{\partial \mathbf{x}}{\partial \boldsymbol{X}} \quad \implies \quad F_{ij} = \frac{\partial x_i}{\partial X_j} $$ Therefore

\begin{align} F_{11} &= \frac{\partial x_1}{\partial X_1} = \cfrac{1}{4}\left[4 + (- 3 - X_2)~t\right] \\ F_{12} &= \frac{\partial x_1}{\partial X_2} = \cfrac{1}{4}\left[(- 5 - X_1)~t\right] \\ F_{21} &= \frac{\partial x_2}{\partial X_1} = 2~t \\ F_{22} &= \frac{\partial x_2}{\partial X_2} = 1 \end{align} $$ At $$t = 1$$ at the position $$\boldsymbol{X} = (0,0)$$ we have

\mathbf{F} = \begin{bmatrix} \frac{\partial x_1}{\partial X_1} & \frac{\partial x_1}{\partial X_2} \\ \frac{\partial x_2}{\partial X_1} & \frac{\partial x_2}{\partial X_2} \end{bmatrix} = \cfrac{1}{4} \begin{bmatrix} 1 & -5 \\ 8 & 4        \end{bmatrix} $$ You can calculate the deformation gradient at other points in a similar manner.

Right Cauchy-Green deformation tensor
We have

\boldsymbol{C} = \boldsymbol{F}^T\cdot\boldsymbol{F} $$ Therefore,

\mathbf{C} = \mathbf{F}^T~\mathbf{F} = \cfrac{1}{16} \begin{bmatrix} 65 & 27 \\ 27 & 41 \end{bmatrix} $$ To compute $$\boldsymbol{U}$$ we have to find the eigenvalues and eigenvectors of $$\boldsymbol{C}$$. The eigenvalue problem is

(\mathbf{C} - \lambda^2~\mathbf{I})\mathbf{N} = \mathbf{0} $$ where

\mathbf{I} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} $$ To find the eigenvalues we solve the characteristic equation

\det(\mathbf{C} - \lambda^2~\mathbf{I}) = 0 $$ Plugging in the numbers, we get

\det\begin{bmatrix} \cfrac{65}{16} - \lambda^2 & \cfrac{27}{16} \\ \cfrac{27}{16} & \cfrac{41}{16} - \lambda^2 \end{bmatrix} = 0 $$ or

\lambda^4 - \cfrac{53}{8}~\lambda^2 + \cfrac{121}{16} = 0 $$ This equation has two solutions

\begin{align} \lambda_1^2 & = \cfrac{53}{16} + \cfrac{3}{16}~\sqrt{97} = 5.159 \\ \lambda_2^2 & = \cfrac{53}{16} - \cfrac{3}{16}~\sqrt{97} = 1.466 \end{align} $$ Taking the square roots we get the values of the principal stretches

\lambda_1 = 2.2714 \qquad \lambda_2 = 1.2107 $$ To compute the eigenvectors we plug into the eigenvalues into the eigenvalue problem to get

\left\{ \begin{bmatrix} 65 & 27 \\ 27 & 41  \end{bmatrix} - \lambda^2_1 ~ \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \right\}~\begin{bmatrix} N_1^{(1)} \\ N_2^{(1)} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} $$ Because this system of equations is not linearly independent, we need another equation to solve this system of equations for $$N^{(1)}_1$$ and $$N^{(1)}_2$$. This problem is eliminated by using the following equation (which implies that $$\mathbf{N}$$ is a unit vector)

N^{(1)}_2 = \sqrt{1 - (N^{(1)}_1)^2} $$ Solving, we get

\mathbf{N}_1 = \begin{bmatrix} N_1^{(1)} \\ N_2^{(1)} \end{bmatrix} = \begin{bmatrix} 0.8385 \\ 0.5449 \end{bmatrix} $$ We can do the same thing for the other eigenvector $$\mathbf{N}_2$$ to get

\mathbf{N}_2 = \begin{bmatrix} N_1^{(2)} \\ N_2^{(2)} \end{bmatrix} = \begin{bmatrix} -0.5449 \\0.8385 \end{bmatrix} $$ Therefore,

\boldsymbol{N}_1\otimes\boldsymbol{N}_1 = \mathbf{N}_1~\mathbf{N}_1^T = \begin{bmatrix} 0.8385 \\ 0.5449 \end{bmatrix} \begin{bmatrix} 0.8385 & 0.5449 \end{bmatrix} = \begin{bmatrix} 0.7031 & 0.4569 \\ 0.4569 & 0.2969 \end{bmatrix} $$ and

\boldsymbol{N}_2\otimes\boldsymbol{N}_2 = \mathbf{N}_2~\mathbf{N}_2^T = \begin{bmatrix} -0.5449 \\0.8385 \end{bmatrix} \begin{bmatrix} -0.5449 & 0.8385 \end{bmatrix} = \begin{bmatrix} 0.2969 & -0.4569 \\ -0.4569 & 0.7031 \end{bmatrix} $$ Therefore,

\boldsymbol{C} = \lambda_1^2~\boldsymbol{N}_1\otimes\boldsymbol{N}_1 + \lambda_2^2~\boldsymbol{N}_2\otimes\boldsymbol{N}_2 \quad \implies \quad \mathbf{C} = 5.159~ \begin{bmatrix} 0.7031 & 0.4569 \\ 0.4569 & 0.2969 \end{bmatrix} + 1.466~      \begin{bmatrix} 0.2969 & -0.4569 \\ -0.4569 & 0.7031 \end{bmatrix} $$ We usually don't see any problem to calculate $$\boldsymbol{C}$$ at this point and go straight to the right stretch tensor.

Right stretch
The right stretch tensor $$\boldsymbol{U}$$ is given by

\boldsymbol{U} = \lambda_1~\boldsymbol{N}_1\otimes\boldsymbol{N}_1 + \lambda_2~\boldsymbol{N}_2\otimes\boldsymbol{N}_2 \quad \implies \quad \mathbf{U} = 2.2714~ \begin{bmatrix} 0.7031 & 0.4569 \\ 0.4569 & 0.2969 \end{bmatrix} + 1.2107~      \begin{bmatrix} 0.2969 & -0.4569 \\ -0.4569 & 0.7031 \end{bmatrix} $$ or

\mathbf{U} = \begin{bmatrix} 1.9565 & 0.4846 \\ 0.4846 & 1.5256 \end{bmatrix} $$ We can invert this matrix to get

\mathbf{U}^{-1} = \begin{bmatrix} 0.5548 & -0.1762 \\ -0.1762 & 0.7114 \end{bmatrix} $$

Rotation
We can now find the rotation matrix by using th relation

\boldsymbol{R} = \boldsymbol{F}\cdot\boldsymbol{U}^{-1} $$ In matrix form,

\mathbf{R} = \cfrac{1}{4} \begin{bmatrix} 1 & -5 \\ 8 & 4 \end{bmatrix} \begin{bmatrix} 0.5548 & -0.1762 \\ -0.1762 & 0.7114 \end{bmatrix} = \begin{bmatrix} 0.3590 & -0.9334 \\ 0.9334 & 0.3590 \end{bmatrix} $$

You can check whether this matrix is orthogonal by seeing whether $$\mathbf{R}~\mathbf{R}^T = \mathbf{R}^T~\mathbf{R} = \mathbf{I}$$.

You thus get the polar decomposition of $$\boldsymbol{F}$$. In an actual calculation you have to be careful about floating point errors. Otherwise you might not get a matrix that is orthogonal.