# Nonlinear finite elements/Homework 7/Solutions
## Problem 1: Index Notation
### Solutions
#### Part 1

* Determine whether the following expressions are valid in index notation. If valid, identify the free indices and the dummy indices.
1) $A_{ms}=b_{m}(c_{r}-d_{r})$
Invalid. The free indices are $m,s$ on the LHS and $m,r$ on the RHS.
2) $A_{ms}=b_{m}(c_{s}-d_{s})$
Valid. Both $m$ and $s$ are free indices.
3) $t_{i}=\sigma _{ji}n_{j}$
Valid. The free index is $i$ and the dummy index is $j$.
4) $t_{i}=\sigma _{ji}n_{i}$
Invalid. The free index is $i$ on the LHS and $j$ on the RHS.
5) $x_{i}x_{i}=r^{3}$
Valid. The dummy index is $i$. So the sum is a scalar which is equal to $r^{3}$.
6) $B_{ij}c_{j}=3$
Invalid. There is one free index on the LHS and no free index on the RHS.

#### Part 2

Show the following:
1) $\delta _{ii}=3$
    $\delta _{ii}=\delta _{11}+\delta _{22}+\delta _{33}=1+1+1=3\qquad \implies \qquad \delta _{ii}=3$
2) $e_{ijk}e_{pqk}=\delta _{ip}\delta _{jq}-\delta _{iq}\delta _{jp}$
The LHS is
    $e_{ijk}e_{pqk}=e_{ij1}e_{pq1}+e_{ij2}e_{pq2}+e_{ij3}e_{pq3}$
If $i=j$, then we have
    $LHS=({e_{iik}})~~(e_{pqk})=0~;~~RHS=\delta _{ip}\delta _{iq}-\delta _{iq}\delta _{ip}=0$
We get the same result if $p=q$. The only nonzero LHS and RHS occur when $i\neq j$ and $p\neq q$.
* Case 1: $i=1$, $j=2$, $p=1$, $q=2$.
    $LHS=({e_{121}e_{121}})~~+({e_{122}e_{122}})~~+({e_{123}e_{123}})~~=1~;~~RHS=({\delta _{11}\delta _{11}})~~-({\delta _{12}\delta _{12}})~~=1$
* Case 2: $i=1$, $j=2$, $p=2$, $q=1$ or
$i=2$, $j=1$, $p=1$, $q=2$.
    $LHS=({e_{121}e_{211}})~~+({e_{122}e_{212}})~~+({e_{123}e_{213}})~~=-1~;~~RHS=({\delta _{12}\delta _{21}})~~-({\delta _{11}\delta _{22}})~~=-1$
* Case 3: $i=2$, $j=1$, $p=2$, $q=1$.
    $LHS=({e_{211}e_{211}})~~+({e_{212}e_{212}})~~+({e_{213}e_{213}})~~=1~;~~RHS=({\delta _{22}\delta _{11}})~~-({\delta _{21}\delta _{12}})~~=1$
* Case 4: $i=2$, $j=3$, $p=2$, $q=3$.
    $LHS=({e_{231}e_{231}})~~+({e_{232}e_{232}})~~+({e_{233}e_{233}})~~=1~;~~RHS=({\delta _{22}\delta _{22}})~~-({\delta _{23}\delta _{23}})~~=1$
* Case 5: $i=2$, $j=3$, $p=3$, $q=2$ or
$i=3$, $j=2$, $p=2$, $q=3$.
    $LHS=({e_{231}e_{321}})~~+({e_{232}e_{322}})~~+({e_{233}e_{323}})~~=-1~;~~RHS=({\delta _{23}\delta _{32}})~~-({\delta _{22}\delta _{33}})~~=-1$
* Case 6: $i=3$, $j=2$, $p=3$, $q=2$.
    $LHS=({e_{321}e_{321}})~~+({e_{322}e_{322}})~~+({e_{323}e_{323}})~~=1~;~~RHS=({\delta _{33}\delta _{22}})~~-({\delta _{32}\delta _{23}})~~=1$
* Case 7: $i=3$, $j=1$, $p=3$, $q=1$.
    $LHS=({e_{311}e_{311}})~~+({e_{312}e_{312}})~~+({e_{313}e_{313}})~~=1~;~~RHS=({\delta _{33}\delta _{11}})~~-({\delta _{31}\delta _{13}})~~=1$
* Case 8: $i=3$, $j=1$, $p=1$, $q=3$ or
$i=1$, $j=3$, $p=3$, $q=1$.
    $LHS=({e_{311}e_{131}})~~+({e_{312}e_{132}})~~+({e_{313}e_{133}})~~=-1~;~~RHS=({\delta _{31}\delta _{13}})~~-({\delta _{33}\delta _{11}})~~=-1$
* Case 9: $i=1$, $j=3$, $p=1$, $q=3$.
    $LHS=({e_{131}e_{131}})~~+({e_{132}e_{132}})~~+({e_{133}e_{133}})~~=1~;~~RHS=({\delta _{11}\delta _{33}})~~-({\delta _{13}\delta _{31}})~~=1$
Hence the $~e-\delta$ relation is satisfied for all cases.
3) $\delta _{ij}e_{ijk}=0$
    ${{\delta _{ij}e_{ijk}=e_{iik}=0}.}$
4) $e_{qrs}d_{q}d_{s}=0$
    ${\begin{aligned}e_{qrs}d_{q}d_{s}=&\sum _{q=1}^{3}\sum _{s=1}^{3}e_{qrs}d_{q}d_{s}\\=&({e_{1r1}d_{1}d_{1}})~~~+e_{1r2}d_{1}d_{2}+e_{1r3}d_{1}d_{3}+\\&e_{2r1}d_{2}d_{1}+({e_{2r2}d_{2}d_{2}})~~~+e_{2r3}d_{2}d_{3}+\\&e_{3r1}d_{3}d_{1}+e_{3r2}d_{3}d_{2}+({e_{3r3}d_{3}d_{3}})\end{aligned}}$
For $r=1$,
    $e_{qrs}d_{q}d_{s}=({e_{112}})~~~d_{1}d_{2}+({e_{113}})~~~d_{1}d_{3}+({e_{211}})~~~d_{2}d_{1}+({e_{213}})~~~d_{2}d_{3}+({e_{311}})~~~d_{3}d_{1}+({e_{312}})~~~d_{3}d_{2}=0$
For $r=2$,
    $e_{qrs}d_{q}d_{s}=({e_{122}})~~~d_{1}d_{2}+({e_{123}})~~~d_{1}d_{3}+({e_{221}})~~~d_{2}d_{1}+({e_{223}})~~~d_{2}d_{3}+({e_{321}})~~~d_{3}d_{1}+({e_{322}})~~~d_{3}d_{2}=0$
For $r=3$,
    $e_{qrs}d_{q}d_{s}=({e_{132}})~~~d_{1}d_{2}+({e_{133}})~~~d_{1}d_{3}+({e_{231}})~~~d_{2}d_{1}+({e_{233}})~~~d_{2}d_{3}+({e_{331}})~~~d_{3}d_{1}+({e_{332}})~~~d_{3}d_{2}=0$
Hence shown.

#### Part 3

The elasticity tensor is given by
    ${\boldsymbol {\mathsf {C}}}=\lambda ~{\boldsymbol {\mathit {1}}}\otimes {\boldsymbol {\mathit {1}}}+2\mu ~{\boldsymbol {\mathsf {I}}}$
where $\lambda ,\mu$ are Lame constants, ${\boldsymbol {\mathit {1}}}$ is the second order identity tensor, and ${\boldsymbol {\mathsf {I}}}$ is the fourth-order symmetric identity tensor. The two identity tensors are defined as
    ${\begin{aligned}{\boldsymbol {\mathit {1}}}&=\delta _{ij}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\\{\boldsymbol {\mathsf {I}}}&={\frac {1}{2}}[\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}]~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\otimes \mathbf {e} _{k}\otimes \mathbf {e} _{l}\end{aligned}}$
The stress-strain relation is
    ${\boldsymbol {\sigma }}={\boldsymbol {\mathsf {C}}}:{\boldsymbol {\varepsilon }}$
Show that the stress-strain relation can be written in index notation as
    $\sigma _{ij}=2\mu \varepsilon _{ij}+\lambda \varepsilon _{kk}\delta _{ij}~.$
Write the stress-strain relations in expanded form.
    ${\begin{aligned}{\boldsymbol {\sigma }}&={\boldsymbol {\mathsf {C}}}:{\boldsymbol {\varepsilon }}\\&=\left(\lambda ~{\boldsymbol {\mathit {1}}}\otimes {\boldsymbol {\mathit {1}}}+2\mu ~{\boldsymbol {\mathsf {I}}}\right):{\boldsymbol {\varepsilon }}\\&=\left(\lambda ~{\boldsymbol {\mathit {1}}}\otimes {\boldsymbol {\mathit {1}}}\right):{\boldsymbol {\varepsilon }}+2\mu ~{\boldsymbol {\mathsf {I}}}:{\boldsymbol {\varepsilon }}\end{aligned}}$
Now, in dyadic notation
    ${\boldsymbol {\mathit {1}}}\otimes {\boldsymbol {\mathit {1}}}=(\delta _{ij}~\mathbf {e} _{i}\otimes \mathbf {e} _{j})\otimes (\delta _{kl}~\mathbf {e} _{k}\otimes \mathbf {e} _{l})=\delta _{ij}\delta _{kl}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\otimes \mathbf {e} _{k}\otimes \mathbf {e} _{l}~{\text{and}}~{\boldsymbol {\varepsilon }}=\varepsilon _{kl}~\mathbf {e} _{k}\otimes \mathbf {e} _{l}$
Therefore
    ${\begin{aligned}({\boldsymbol {\mathit {1}}}\otimes {\boldsymbol {\mathit {1}}}):{\boldsymbol {\varepsilon }}&=\delta _{ij}\delta _{kl}\varepsilon _{kl}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\\&=\delta _{ij}\varepsilon _{kk}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\end{aligned}}$
Similarly,
    ${\begin{aligned}{\boldsymbol {\mathsf {I}}}:{\boldsymbol {\varepsilon }}&=\left({\frac {1}{2}}[\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}]\varepsilon _{kl}\right)~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\\&={\frac {1}{2}}\left(\delta _{ik}\varepsilon _{kj}+\delta _{il}\varepsilon _{jl}\right)~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\\&={\frac {1}{2}}\left(\varepsilon _{ij}+\varepsilon _{ji}\right)~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\\&=\varepsilon _{ij}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\qquad {\text{(symmetry)}}\end{aligned}}$
The stress-strain law becomes
    ${\begin{aligned}{\boldsymbol {\sigma }}&=\lambda ~\delta _{ij}\varepsilon _{kk}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}+2\mu ~\varepsilon _{ij}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}\\&=\left(\lambda \delta _{ij}\varepsilon _{kk}+2\mu \varepsilon _{ij}\right)\mathbf {e} _{i}\otimes \mathbf {e} _{j}\end{aligned}}$
Expanding the left hand side, we get
    $\sigma _{ij}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}=\left(\lambda \delta _{ij}\varepsilon _{kk}+2\mu \varepsilon _{ij}\right)\mathbf {e} _{i}\otimes \mathbf {e} _{j}$
Therefore,
    $\sigma _{ij}=\lambda \delta _{ij}\varepsilon _{kk}+2\mu \varepsilon _{ij}$

## Problem 2: Rotating Vectors and Tensors

Let ($\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3}$) be an orthonormal basis. Let ${\boldsymbol {A}}$ be a second order tensor and $\mathbf {u}$ be a vector with components
    ${\begin{aligned}{\boldsymbol {A}}&=5~\mathbf {e} _{1}\otimes \mathbf {e} _{1}-4~\mathbf {e} _{2}\otimes \mathbf {e} _{1}+2~\mathbf {e} _{3}\otimes \mathbf {e} _{3}\\\mathbf {u} &=-2~\mathbf {e} _{1}+3~\mathbf {e} _{3}\end{aligned}}$

### Solution
#### Part 1

Write out ${\boldsymbol {A}}$ and $\mathbf {u}$ in matrix notation.
    ${\mathbf {A} ={\begin{bmatrix}5&0&0\\-4&0&0\\0&0&2\end{bmatrix}};\qquad \mathbf {u} ={\begin{bmatrix}-2\\0\\3\end{bmatrix}}}$

#### Part 2

Find the components of the vector $\mathbf {v} ={\boldsymbol {A}}\bullet \mathbf {u}$ in the basis ($\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3}$).
The components of $\mathbf {v}$ are
    $v_{i}=A_{ij}u_{j}$
Therefore
    $v_{1}=(5)(-2)=-10;~~v_{2}=(-4)(-2)=8;~~v_{3}=(2)(3)=6$
    ${\mathbf {v} =-10~\mathbf {e} _{1}+8\mathbf {e} _{2}+6\mathbf {e} _{3}}$

#### Part 3

Find the components of the vector $\mathbf {w} =\mathbf {v} \times \mathbf {u}$ in the basis ($\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3}$).
The cross product is given by
    $\mathbf {w} =\mathbf {v} \times \mathbf {u} ={\begin{vmatrix}\mathbf {e} _{1}&\mathbf {e} _{2}&\mathbf {e} _{3}\\v_{1}&v_{2}&v_{3}\\u_{1}&u_{2}&u_{3}\end{vmatrix}}={\begin{vmatrix}\mathbf {e} _{1}&\mathbf {e} _{2}&\mathbf {e} _{3}\\-10&8&6\\-2&0&3\end{vmatrix}}=(8)(3)\mathbf {e} _{1}-[(-10)(3)-(6)(-2)]\mathbf {e} _{2}-(8)(-2)\mathbf {e} _{3}$
Therefore,
    ${\mathbf {w} =24~\mathbf {e} _{1}+18\mathbf {e} _{2}+16\mathbf {e} _{3}}$

#### Part 4

Find the components of the tensor ${\boldsymbol {C}}=\mathbf {v} \otimes \mathbf {w}$ in the orthonormal basis.
The tensor product is given by
    $\mathbf {v} \otimes \mathbf {w} =v_{i}w_{j}~\mathbf {e} _{i}\otimes \mathbf {e} _{j}$
Hence, in matrix notation
    $\mathbf {C} ={\begin{bmatrix}v_{1}\\v_{2}\\v_{3}\end{bmatrix}}{\begin{bmatrix}w_{1}&w_{2}&w_{3}\end{bmatrix}}={\begin{bmatrix}-10\\8\\6\end{bmatrix}}{\begin{bmatrix}24&18&16\end{bmatrix}}$
    ${\mathbf {C} ={\begin{bmatrix}-240&-180&-160\\192&144&128\\144&108&96\end{bmatrix}}}$

#### Part 5

Rotate the basis clockwise by 30 degrees around the $\mathbf {e} _{3}$ direction. Find the components of $\mathbf {u}$, $\mathbf {v}$, $\mathbf {w}$, and ${\boldsymbol {C}}$ in the rotated basis.
The vector transformation rule is
    $v_{i}^{'}=l_{ij}v_{i}$
where $l_{ij}$ are the direction cosines.
In this case, the direction cosines are
    ${\begin{aligned}l_{11}&=\mathbf {e} _{1}^{'}\bullet \mathbf {e} _{1}=\cos(30)=0.866&l_{12}&=\mathbf {e} _{1}^{'}\bullet \mathbf {e} _{2}=\cos(60)=0.5&l_{13}&=\mathbf {e} _{1}^{'}\bullet \mathbf {e} _{3}=\cos(90)=0\\l_{21}&=\mathbf {e} _{2}^{'}\bullet \mathbf {e} _{1}=\cos(120)=-0.5&l_{22}&=\mathbf {e} _{2}^{'}\bullet \mathbf {e} _{2}=\cos(30)=0.866&l_{23}&=\mathbf {e} _{2}^{'}\bullet \mathbf {e} _{3}=\cos(90)=0\\l_{31}&=\mathbf {e} _{3}^{'}\bullet \mathbf {e} _{1}=\cos(90)=0&l_{32}&=\mathbf {e} _{3}^{'}\bullet \mathbf {e} _{2}=\cos(90)=0&l_{33}&=\mathbf {e} _{3}^{'}\bullet \mathbf {e} _{3}=\cos(0)=1\end{aligned}}$
Therefore,
    $\mathbf {u} ^{'}=\mathbf {L} \mathbf {u} ={\begin{bmatrix}0.866&0.5&0\\-0.5&0.866&0\\0&0&1\end{bmatrix}}{\begin{bmatrix}-2\\0\\3\end{bmatrix}}={\begin{bmatrix}-1.732\\1\\3\end{bmatrix}}$
    ${\mathbf {u} ^{'}=-1.732\mathbf {e} _{1}+1\mathbf {e} _{2}+3\mathbf {e} _{3}}$
Similarly,
    $\mathbf {v} ^{'}=\mathbf {L} \mathbf {v} ={\begin{bmatrix}0.866&0.5&0\\-0.5&0.866&0\\0&0&1\end{bmatrix}}{\begin{bmatrix}-10\\8\\6\end{bmatrix}}={\begin{bmatrix}-4.66\\11.93\\6\end{bmatrix}}$
    ${\mathbf {v} ^{'}=-4.66\mathbf {e} _{1}+11.93\mathbf {e} _{2}+6\mathbf {e} _{3}}$
Also,
    $\mathbf {w} ^{'}=\mathbf {L} \mathbf {w} ={\begin{bmatrix}0.866&0.5&0\\-0.5&0.866&0\\0&0&1\end{bmatrix}}{\begin{bmatrix}24\\18\\16\end{bmatrix}}={\begin{bmatrix}29.78\\3.59\\16\end{bmatrix}}$
    ${\mathbf {v} ^{'}=29.78\mathbf {e} _{1}+3.59\mathbf {e} _{2}+16\mathbf {e} _{3}}$
From the handout from Slaughter's book, the tensor transformation rule is
    $T_{ij}^{'}=l_{ip}l_{jq}T_{pq}$
where $l_{ij}$ are the direction cosines.
In matrix form,
    $\mathbf {T} ^{'}=\mathbf {L} \mathbf {T} \mathbf {L} ^{T}$
Therefore the components of ${\boldsymbol {C}}$ in the rotated basis are give by
    $\mathbf {C} ^{'}={\begin{bmatrix}0.866&0.5&0\\-0.5&0.866&0\\0&0&1\end{bmatrix}}{\begin{bmatrix}-240&-180&-160\\192&144&128\\144&108&96\end{bmatrix}}{\begin{bmatrix}0.866&-0.5&0\\0.5&0.866&0\\0&0&1\end{bmatrix}}={\begin{bmatrix}-138.8&-16.7&-74.6\\355.3&42.8&190.9\\178.7&21.5&96\end{bmatrix}}$
    ${\mathbf {C} ^{'}={\begin{bmatrix}-138.8&-16.7&-74.6\\355.3&42.8&190.9\\178.7&21.5&96\end{bmatrix}}}$

## Problem 3: More Beams
### Part A

Consider a beam of length $L$ = 100 in., cross-section 1 in. $\times$ 1 in., and subjected to a uniformly distributed transverse load $q_{0}$ lbf/in. Model one half of the beam using symmetry considerations.

#### Part 1

Hinged-Hinged Beam
The boundary conditions are
    $w_{0}(0)=u_{0}(L/2)=\varphi _{x}(L/2)=0~.$
Compute a plot similar to that shown in Figure 4.3.4 for this case using Beam188 elements. What do you observe?
The result is shown in Table 1.
Table 1. Deflections of a hinged-hinged beam :
| Load | $U_{y}$ at $x=L/2$ |
| 1 | -0.520746 |
| 2 | -1.04086 |
| 3 | -1.55922 |
| 4 | -2.07510 |
| 5 | -2.58768 |
| 6 | -3.09622 |
| 7 | -3.60000 |
| 8 | -4.09835 |
| 9 | -4.59065 |
| 10 | -5.07636 |

#### Part 2

Clamped-Clamped Beam
The boundary conditions are
    $u_{0}(0)=w_{0}(0)=\varphi _{x}(0)=u_{0}(L/2)=\varphi _{x=L/2}=0.$
Compute a plot for this case using Beam188 elements. Comment on your plot.
The result is shown in Table 2.
Table 2. Deflections of a clamped-clamped beam :
| Load | $U_{y}$ at $x=L/2$ |
| 1 | -0.103456 |
| 2 | -0.202476 |
| 3 | -0.294220 |
| 4 | -0.377753 |
| 5 | -0.453387 |
| 6 | -0.521968 |
| 7 | -0.584455 |
| 8 | -0.644185 |
| 9 | -0.696440 |
| 10 | -0.745243 |
Listed below is the ANSYS input code for Problem 3A.1 and 3A.2.

### Part B
#### Part 1

Simulate the unrolling of a cantilever beam from Section 4.1.1 of Ibrahimbegovic (1995) and compare your results with the results shown in the paper.
The result is shown in Table 3.
Table 3. Cantilever free-end displacement components :
| Load | $U_{x}$ | $U_{y}$ | Rotation |
| $M=10\pi$ | -0.040666 | 6.3205 | -3.1287 |
| $M=20\pi$ | 9.9578 | 0.12729 | -6.2577 |
The deformation plots are shown in Figure 4 and 5.
* Figure 4. Deformed shape under $M=10\pi$ for Problem 3B.1.
* Figure 5. Deformed shape under $M=20\pi$ for Problem 3B.1.
The ANSYS input code for this problem is listed below.

#### Part 2

Simulate the clamped-hinged deep circular arch from Example 7.3 of Simo and Vu Quoc (1986) and compare you results with the results shown in the paper.
The inputs are:$R=100$, $\varphi =145^{o}$, $EI=10^{6}$, and $EA=100EI$. We assume a square cross section. Then
    ${\begin{aligned}I&={\cfrac {1}{12}}h^{4}~;&A&=h^{2}\\{\cfrac {EA}{EI}}&={\cfrac {A}{I}}~;&{\cfrac {12h^{2}}{h^{4}}}={\cfrac {12}{h^{2}}}=100\end{aligned}}$
Therefore, $h=0.34641$.
The deformed shape (unconverged) for a load of 905 is shown in Figure 6.
* Figure 6. Deformed shape of arch.
The load-displacement curve (up to the last converged solution) is shown in Figure 7.
* Figure 7. Load-displacement plot for circular arch.
The buckling load is 900.925 compared to 905.28 in Simo and Vu Quoc.
The ANSYS file use for the calculations is shown below.
Here is another version of solution to this problem.
* Figure 8. Force-displacement diagram for Problem 3B.2.
* Figure 9. Shape deformation at the last load step for Problem 3B.2.
The ANSYS input code for Problem 3B.2 is listed below.

#### Part 3

Simulate the buckling of a hinged right-angle frame under both fixed and follower loads from Example 7.4 of Simo and Vu Quoc (1986) and compare your results with those shown in the paper.
The force-displacement diagram is shown in Figure 10. The deformation is illustrated in Figure 11 and 12.
* Figure 10. Force-displacement diagram for Problem 3B.3.
* Figure 11. Deformation (fixed load) at the last load step for Problem 3B.3.
* Figure 12. Force-displacement diagram (follower load) for Problem 3B.3.
The ANSYS input code for Problem 3B.3 is listed below.
Warning: The arc length method no longer converges with Ansys 13. Try using the stabilization option instead of arclen, on: