Nonlinear finite elements/Homework 6/Hints

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 ${\displaystyle 5}$ be ${\displaystyle 1}$ for node ${\displaystyle 5}$, and ${\displaystyle 2}$ for node ${\displaystyle 6}$.

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

${\displaystyle \theta ={\cfrac {\pi }{4}}~;~~\theta _{1}={\cfrac {\pi }{3}}~;~~\theta _{2}={\cfrac {\pi }{6}}~.}$

Let ${\displaystyle r_{1}=1}$ and ${\displaystyle r_{2}=1.2}$. Since the blue point is midway between the two, ${\displaystyle r=1.1}$.

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

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

Let the velocities for nodes ${\displaystyle 1}$ and ${\displaystyle 2}$ of the element be

${\displaystyle \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 ${\displaystyle x}$ and ${\displaystyle y}$ coordinates of the master and slave nodes are

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

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

${\displaystyle {\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

${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle \xi }$ are

${\displaystyle {\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 ${\displaystyle \xi }$ and ${\displaystyle \eta }$ at that point are zero. Therefore,

${\displaystyle \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 ${\displaystyle {\widehat {\mathbf {e} }}_{x}}$ is given by

${\displaystyle {\widehat {\mathbf {e} }}_{x}={\cfrac {\mathbf {x} _{,\xi }}{\lVert \mathbf {x} _{,\xi }\rVert _{}}}}$

The laminar basis vector ${\displaystyle {\widehat {\mathbf {e} }}_{y}}$ is given by

${\displaystyle {\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

${\displaystyle {\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

${\displaystyle {\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

${\displaystyle {\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

${\displaystyle {\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

${\displaystyle {\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 (${\displaystyle \xi ,\eta }$). We have to convert them to the (${\displaystyle x,y}$) system in order to compute the velocity gradient. To do that we recall that

${\displaystyle {\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

${\displaystyle {\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

${\displaystyle {\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,

${\displaystyle {\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

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

The global base vectors are

${\displaystyle \mathbf {e} _{x}={\begin{bmatrix}1\\0\end{bmatrix}};~\qquad \mathbf {e} _{y}={\begin{bmatrix}0\\1\end{bmatrix}}~.}$

Therefore, the rotation matrix is

${\displaystyle \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

${\displaystyle \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

${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle x}$-direction maps to the composite ${\displaystyle 3}$-direction and the laminar ${\displaystyle y}$-directions maps to the composite ${\displaystyle 1}$-direction. Hence the constitutive equation can be written as

${\displaystyle {\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,

${\displaystyle {\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 ${\displaystyle \sigma _{yy}=0}$ in the laminar coordinate system. We assume that the rate of ${\displaystyle \sigma _{yy}}$ is also zero. In that case, we get

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

or,

${\displaystyle 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

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