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{aligned}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{aligned}}$

The isoparametric map is

${\begin{aligned}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{aligned}}$

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

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

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_{i-}^{x}\\v_{i-}^{y}\\v_{i+}^{x}\\v_{i+}^{y}\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_{1-}^{x}\\v_{1-}^{y}\\v_{1+}^{x}\\v_{1+}^{y}\end{bmatrix}}$

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

${\begin{bmatrix}v_{2-}^{x}\\v_{2-}^{y}\\v_{2+}^{x}\\v_{2+}^{y}\end{bmatrix}}$

The interpolated velocity within the element is given by

${\begin{aligned}\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{aligned}}$

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{aligned}{\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{aligned}}$

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} _{\text{lam}}^{T}\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} _{\text{lam}}^{T}~;\qquad \mathbf {D} =\mathbf {R} _{\text{lam}}\mathbf {D} _{\text{lam}}\mathbf {R} _{\text{lam}}^{T}~.$

Nonlinear finite elements/Homework 6