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}} ~. $$