Nonlinear finite elements/Homework 6/Hints

< Nonlinear finite elements < Homework 6

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

This article is issued from Wikiversity - version of the Tuesday, January 08, 2008. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.