# Nonlinear finite elements/Homework 6/Hints

## Hints for Homework 6: Problem 1: Section 8

Use Maple to reduce your manual labor.

The problem becomes easier to solve if we consider numerical values of the parameters. Let the local nodes numbers of element $5$ be $1$ for node $5$, and $2$ for node $6$.

Let us assume that the beam is divided into six equal sectors. Then,

$\theta = \cfrac{\pi}{4}~;~~ \theta_1 = \cfrac{\pi}{3}~;~~ \theta_2 = \cfrac{\pi}{6} ~.$

Let $r_1 = 1$ and $r_2 = 1.2$. Since the blue point is midway between the two, $r = 1.1$.

Also, let the components of the stiffness matrix of the composite be

$C_{11} = 10;~~ C_{33} = 100;~~ C_{12} = 6;~~ C_{13} = 60; ~~ C_{44} = 30~.$

Let the velocities for nodes $1$ and $2$ of the element be

$\mathbf{v}_1 = \begin{bmatrix} v_1^x \\ v_1^y \\ \omega_1 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \\ 1 \end{bmatrix} ~;~~ \mathbf{v}_2 = \begin{bmatrix} v_2^x \\ v_2^y \\ \omega_2 \end{bmatrix} = \begin{bmatrix} 2 \\ 1 \\ 1 \end{bmatrix}$

The $x$ and $y$ coordinates of the master and slave nodes are

$\begin{bmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{bmatrix} = \begin{bmatrix} r\cos\theta_1 \\ r\sin\theta_1 \\ r\cos\theta_2 \\ r\sin\theta_2 \end{bmatrix}$
$\begin{bmatrix} x_{1-} \\ y_{1-} \\ x_{2-} \\ y_{2-} \end{bmatrix} = \begin{bmatrix} r_1\cos\theta_1 \\ r_1\sin\theta_1 \\ r_1\cos\theta_2 \\ r_1\sin\theta_2 \end{bmatrix}$
$\begin{bmatrix} x_{1+} \\ y_{1+} \\ x_{2+} \\ y_{2+} \end{bmatrix} = \begin{bmatrix} r_2\cos\theta_1 \\ r_2\sin\theta_1 \\ r_2\cos\theta_2 \\ r_2\sin\theta_2 \end{bmatrix}$

Since there are two master nodes in an element, the parent element is a four-noded quad with shape functions

\begin{align} N_{1-}(\xi,\eta) & = \cfrac{1}{4}(1-\xi)(1-\eta) & N_{2-}(\xi,\eta) & = \cfrac{1}{4}(1+\xi)(1-\eta) \\ N_{1+}(\xi,\eta) & = \cfrac{1}{4}(1-\xi)(1+\eta) & N_{2+}(\xi,\eta) & = \cfrac{1}{4}(1+\xi)(1+\eta) ~. \end{align}

The isoparametric map is

\begin{align} x(\xi, \eta) & = x_{1-}~N_{1^-}(\xi,\eta) + x_{2-}~N_{2^-}(\xi,\eta) + x_{1+}~N_{1^+}(\xi,\eta) + x_{2+}~N_{2^+}(\xi,\eta)\\ y(\xi, \eta) & = y_{1-}~N_{1^-}(\xi,\eta) + y_{2-}~N_{2^-}(\xi,\eta) + y_{1+}~N_{1^+}(\xi,\eta) + y_{2+}~N_{2^+}(\xi,\eta)~. \end{align}

Therefore, the derivatives with respect to $\xi$ are

\begin{align} x_{,\xi} = \frac{\partial x}{\partial \xi} \\ y_{,\xi} = \frac{\partial y}{\partial \xi} \end{align}

Since the blue point is at the center of the element, the values of $\xi$ and $\eta$ at that point are zero. Therefore,

$\mathbf{x}_{,\xi} = x_{,\xi} \mathbf{e}_x - y_{,\xi} \mathbf{e}_y, ~~ \lVert\mathbf{x}_{,\xi}\rVert_{}= \sqrt{x_{,\xi}^2 + y_{,\xi}^2} ~.$

The local laminar basis vector $\widehat{\mathbf{e}}_x$ is given by

$\widehat{\mathbf{e}}_x = \cfrac{\mathbf{x}_{,\xi}}{\lVert\mathbf{x}_{,\xi}\rVert_{}}$

The laminar basis vector $\widehat{\mathbf{e}}_y$ is given by

$\widehat{\mathbf{e}}_y = \mathbf{e}_z\times\widehat{\mathbf{e}_x}$

To compute the velocity gradient, we have to find the velocities at the slave nodes using the relation

$\begin{bmatrix} v^x_{i-} \\ v^y_{i-} \\ v^x_{i+} \\ v^y_{i+} \end{bmatrix} = \begin{bmatrix} 1 & 0 & - (y_{i-}-y_i) \\ 0 & 1 &(x_{i-}-x_i)\\ 1 & 0 & - (y_{i+}-y_i) \\ 0 & 1 &(x_{i+}-x_i) \end{bmatrix} \begin{bmatrix} v_i^x \\ v_i^y \\ \omega_i \end{bmatrix}$

For master node 1 of the element (global node 5), the velocities of the slave nodes are

$\begin{bmatrix} v^x_{1-} \\ v^y_{1-} \\ v^x_{1+} \\ v^y_{1+} \end{bmatrix}$

For master node 2 of the element (global node 6), the velocities of the slave nodes are

$\begin{bmatrix} v^x_{2-} \\ v^y_{2-} \\ v^x_{2+} \\ v^y_{2+} \end{bmatrix}$

The interpolated velocity within the element is given by

\begin{align} \mathbf{v}(\boldsymbol{\xi}, t) & = \cfrac{D}{Dt}[\mathbf{x}(\boldsymbol{\xi},t)] \\ & = \sum_{i- = 1}^2 \frac{\partial }{\partial t}[\mathbf{x}_{i-}(t)]~N_{i^-}(\xi,\eta) + \sum_{i+ = 1}^2 \frac{\partial }{\partial t}[\mathbf{x}_{i+}(t)]~N_{i^+}(\xi,\eta)\\ & = \sum_{i- = 1}^2 \mathbf{v}_{i-}(t)~N_{i-}(\xi,\eta) + \sum_{i+ = 1}^2 \mathbf{v}_{i+}(t)~N_{i+}(\xi,\eta) ~. \end{align}

The velocity gradient is given by

$\boldsymbol{\nabla} \mathbf{v} = v_{i,j} = \frac{\partial v_i}{\partial x_j} ~.$

The velocities are given in terms of the parent element coordinates ($\xi,\eta$). We have to convert them to the ($x,y$) system in order to compute the velocity gradient. To do that we recall that

\begin{align} \frac{\partial v_x}{\partial \xi} & = \frac{\partial v_x}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial v_x}{\partial y}\frac{\partial y}{\partial \xi} \\ \frac{\partial v_x}{\partial \eta} & = \frac{\partial v_x}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial v_x}{\partial y}\frac{\partial y}{\partial \eta} \\ \frac{\partial v_y}{\partial \xi} & = \frac{\partial v_y}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial v_y}{\partial y}\frac{\partial y}{\partial \xi} \\ \frac{\partial v_y}{\partial \eta} & = \frac{\partial v_y}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial v_y}{\partial y}\frac{\partial y}{\partial \eta} ~. \end{align}

In matrix form

$\begin{bmatrix} \frac{\partial v_x}{\partial \xi} \\ \frac{\partial v_x}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \\ \end{bmatrix} \begin{bmatrix} \frac{\partial v_x}{\partial x} \\ \frac{\partial v_x}{\partial y} \end{bmatrix} = \mathbf{J} \begin{bmatrix} \frac{\partial v_x}{\partial x} \\ \frac{\partial v_x}{\partial y} \end{bmatrix}$

and

$\begin{bmatrix} \frac{\partial v_y}{\partial \xi} \\ \frac{\partial v_y}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \\ \end{bmatrix} \begin{bmatrix} \frac{\partial v_y}{\partial x} \\ \frac{\partial v_y}{\partial y} \end{bmatrix} = \mathbf{J} \begin{bmatrix} \frac{\partial v_y}{\partial x} \\ \frac{\partial v_y}{\partial y} \end{bmatrix}~.$

Therefore,

$\begin{bmatrix} \frac{\partial v_x}{\partial x} \\ \frac{\partial v_x}{\partial y} \end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \frac{\partial v_x}{\partial \xi} \\ \frac{\partial v_x}{\partial \eta} \end{bmatrix} \qquad \text{and} \qquad \begin{bmatrix} \frac{\partial v_y}{\partial x} \\ \frac{\partial v_y}{\partial y} \end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \frac{\partial v_y}{\partial \xi} \\ \frac{\partial v_y}{\partial \eta} \end{bmatrix}~.$

The rate of deformation is given by

$\boldsymbol{D} = \frac{1}{2}[\boldsymbol{\nabla} \mathbf{v} + (\boldsymbol{\nabla} \mathbf{v})^T] ~.$

The global base vectors are

$\mathbf{e}_x = \begin{bmatrix} 1 \\ 0 \end{bmatrix}; ~\qquad \mathbf{e}_y = \begin{bmatrix} 0 \\ 1 \end{bmatrix}~.$

Therefore, the rotation matrix is

$\mathbf{R}_{\text{lam}} = \begin{bmatrix} \mathbf{e}_x\bullet\widehat{\mathbf{e}_x} & \mathbf{e}_x\bullet\widehat{\mathbf{e}_y} \\ \mathbf{e}_y\bullet\widehat{\mathbf{e}_x} & \mathbf{e}_y\bullet\widehat{\mathbf{e}_y} \end{bmatrix}$

Therefore, the components of the rate of deformation tensor with respect to the laminar coordinate system are

$\mathbf{D}_{\text{lam}} = \mathbf{R}^T_{\text{lam}}\mathbf{D}\mathbf{R}_{\text{lam}}$

The rate constitutive relation of the material is given by

$\cfrac{D}{Dt} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{31} \\ \sigma_{12} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & (C_{11}-C_{12})/2 \end{bmatrix} \begin{bmatrix} D_{11} \\ D_{22} \\ D_{33} \\ D_{23} \\ D_{31} \\ D_{12} \end{bmatrix}~.$

Since the problem is a 2-D one, the reduced constitutive equation is

$\cfrac{D}{Dt} \begin{bmatrix} \sigma_{11} \\ \sigma_{33} \\ \sigma_{31} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{13} & 0 \\ C_{13} & C_{33} & 0\\ 0& 0& C_{44} \end{bmatrix} \begin{bmatrix} D_{11} \\ D_{33} \\ D_{31} \end{bmatrix}~.$

The laminar $x$-direction maps to the composite $3$-direction and the laminar $y$-directions maps to the composite $1$-direction. Hence the constitutive equation can be written as

$\cfrac{D}{Dt} \begin{bmatrix} \sigma_{yy} \\ \sigma_{xx} \\ \sigma_{xy} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{13} & 0 \\ C_{13} & C_{33} & 0\\ 0& 0& C_{44} \end{bmatrix} \begin{bmatrix} D_{yy} \\ D_{xx} \\ D_{xy} \end{bmatrix}~.$

Rearranging,

$\cfrac{D}{Dt} \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{bmatrix} = \begin{bmatrix} C_{33} & C_{13} & 0\\ C_{13} & C_{11} & 0 \\ 0& 0& C_{44} \end{bmatrix} \begin{bmatrix} D_{xx} \\ D_{yy} \\ D_{xy} \end{bmatrix}~.$

The plane stress condition requires that $\sigma_{yy} = 0$ in the laminar coordinate system. We assume that the rate of $\sigma_{yy}$ is also zero. In that case, we get

$0 = \cfrac{D}{Dt}(\sigma_{yy}) = C_{13}~D_{xx} + C_{11}~D{yy}$

or,

$D_{yy} = - \cfrac{C_{13}}{C_{11}}~D_{xx} ~.$

To get the global stress rate and rate of deformation, we have to rotate the components to the global basis using

$\cfrac{D}{Dt}(\boldsymbol{\sigma}) = \mathbf{R}_{\text{lam}} \cfrac{D}{Dt}(\boldsymbol{\sigma}_{\text{lam}}) \mathbf{R}^T_{\text{lam}} ~; \qquad \mathbf{D} = \mathbf{R}_{\text{lam}} \mathbf{D}_{\text{lam}} \mathbf{R}^T_{\text{lam}} ~.$