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 =\frac{\pi }{4};{\theta }_1=\frac{\pi }{3};{\theta }_2=\frac{\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

${𝐯}_1=\left[\begin{array}{c}{v}_1^x\\ v_1^y\\ {\omega }_1\end{array}\right]=\left[\begin{array}{c}1\\ 2\\ 1\end{array}\right];{𝐯}_2=\left[\begin{array}{c}{v}_2^x\\ v_2^y\\ {\omega }_2\end{array}\right]=\left[\begin{array}{c}2\\ 1\\ 1\end{array}\right]$

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

$\left[\begin{array}{c}{x}_1\\ y_1\\ x_2\\ y_2\end{array}\right]=\left[\begin{array}{c}r\cos {\theta }_1\\ r\sin {\theta }_1\\ r\cos {\theta }_2\\ r\sin {\theta }_2\end{array}\right]$

$\left[\begin{array}{c}{x}_{1-}\\ y_{1-}\\ x_{2-}\\ y_{2-}\end{array}\right]=\left[\begin{array}{c}{r}_1\cos {\theta }_1\\ r_1\sin {\theta }_1\\ r_1\cos {\theta }_2\\ r_1\sin {\theta }_2\end{array}\right]$

$\left[\begin{array}{c}{x}_{1+}\\ y_{1+}\\ x_{2+}\\ y_{2+}\end{array}\right]=\left[\begin{array}{c}{r}_2\cos {\theta }_1\\ r_2\sin {\theta }_1\\ r_2\cos {\theta }_2\\ r_2\sin {\theta }_2\end{array}\right]$

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

$\begin{array}{cccc}{N}_{1-}(\xi ,\eta )& =\frac{1}{4}(1-\xi )(1-\eta )& N_{2-}(\xi ,\eta )& =\frac{1}{4}(1+\xi )(1-\eta )\\ N_{1+}(\xi ,\eta )& =\frac{1}{4}(1-\xi )(1+\eta )& N_{2+}(\xi ,\eta )& =\frac{1}{4}(1+\xi )(1+\eta ).\end{array}$

The isoparametric map is

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

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

$\begin{array}{c}{x}_{,\xi }=\frac{\partial x}{\partial \xi }\\ y_{,\xi }=\frac{\partial y}{\partial \xi }\end{array}$

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

${𝐱}_{,\xi }=x_{,\xi }{𝐞}_x-y_{,\xi }{𝐞}_y,\Vert {𝐱}_{,\xi }{\Vert }_{}=\sqrt{x_{,\xi }^2+y_{,\xi }^2}.$

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

${\widehat{𝐞}}_x=\frac{{𝐱}_{,\xi }}{\Vert {𝐱}_{,\xi }{\Vert }_{}}$

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

${\widehat{𝐞}}_y={𝐞}_z\times \widehat{{𝐞}_x}$

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

$\left[\begin{array}{c}{v}_{i-}^x\\ v_{i-}^y\\ v_{i+}^x\\ v_{i+}^y\end{array}\right]=\left[\begin{array}{ccc}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{array}\right]\left[\begin{array}{c}{v}_i^x\\ v_i^y\\ {\omega }_i\end{array}\right]$

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

$\left[\begin{array}{c}{v}_{1-}^x\\ v_{1-}^y\\ v_{1+}^x\\ v_{1+}^y\end{array}\right]$

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

$\left[\begin{array}{c}{v}_{2-}^x\\ v_{2-}^y\\ v_{2+}^x\\ v_{2+}^y\end{array}\right]$

The interpolated velocity within the element is given by

$\begin{array}{cc}𝐯(\mathit{\xi },t)& =\frac{D}{Dt}[𝐱(\mathit{\xi },t)]\\ & =\sum _{i-=1}^2\frac{\partial }{\partial t}[{𝐱}_{i-}(t)]N_{i^{-}}(\xi ,\eta )+\sum _{i+=1}^2\frac{\partial }{\partial t}[{𝐱}_{i+}(t)]N_{i^{+}}(\xi ,\eta )\\ & =\sum _{i-=1}^2{𝐯}_{i-}(t)N_{i-}(\xi ,\eta )+\sum _{i+=1}^2{𝐯}_{i+}(t)N_{i+}(\xi ,\eta ).\end{array}$

The velocity gradient is given by

$\mathrm{\nabla }𝐯=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{array}{cc}\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{array}$

In matrix form

$\left[\begin{array}{c}\frac{\partial v_x}{\partial \xi }\\ \frac{\partial v_x}{\partial \eta }\end{array}\right]=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }\end{array}\right]\left[\begin{array}{c}\frac{\partial v_x}{\partial x}\\ \frac{\partial v_x}{\partial y}\end{array}\right]=𝐉\left[\begin{array}{c}\frac{\partial v_x}{\partial x}\\ \frac{\partial v_x}{\partial y}\end{array}\right]$

and

$\left[\begin{array}{c}\frac{\partial v_y}{\partial \xi }\\ \frac{\partial v_y}{\partial \eta }\end{array}\right]=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }\end{array}\right]\left[\begin{array}{c}\frac{\partial v_y}{\partial x}\\ \frac{\partial v_y}{\partial y}\end{array}\right]=𝐉\left[\begin{array}{c}\frac{\partial v_y}{\partial x}\\ \frac{\partial v_y}{\partial y}\end{array}\right].$

Therefore,

$\left[\begin{array}{c}\frac{\partial v_x}{\partial x}\\ \frac{\partial v_x}{\partial y}\end{array}\right]={𝐉}^{-1}\left[\begin{array}{c}\frac{\partial v_x}{\partial \xi }\\ \frac{\partial v_x}{\partial \eta }\end{array}\right]\text{and}\left[\begin{array}{c}\frac{\partial v_y}{\partial x}\\ \frac{\partial v_y}{\partial y}\end{array}\right]={𝐉}^{-1}\left[\begin{array}{c}\frac{\partial v_y}{\partial \xi }\\ \frac{\partial v_y}{\partial \eta }\end{array}\right].$

The rate of deformation is given by

$𝑫=\frac{1}{2}[\mathrm{\nabla }𝐯+(\mathrm{\nabla }𝐯{)}^T].$

The global base vectors are

${𝐞}_x=\left[\begin{array}{c}1\\ 0\end{array}\right];{𝐞}_y=\left[\begin{array}{c}0\\ 1\end{array}\right].$

Therefore, the rotation matrix is

${𝐑}_{\text{lam}}=\left[\begin{array}{cc}{𝐞}_x\bullet \widehat{{𝐞}_x}& {𝐞}_x\bullet \widehat{{𝐞}_y}\\ {𝐞}_y\bullet \widehat{{𝐞}_x}& {𝐞}_y\bullet \widehat{{𝐞}_y}\end{array}\right]$

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

${𝐃}_{\text{lam}}={𝐑}_{\text{lam}}^T𝐃{𝐑}_{\text{lam}}$

The rate constitutive relation of the material is given by

$\frac{D}{Dt}\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{23}\\ {\sigma }_{31}\\ {\sigma }_{12}\end{array}\right]=\left[\begin{array}{cccccc}{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{array}\right]\left[\begin{array}{c}{D}_{11}\\ D_{22}\\ D_{33}\\ D_{23}\\ D_{31}\\ D_{12}\end{array}\right].$

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

$\frac{D}{Dt}\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{33}\\ {\sigma }_{31}\end{array}\right]=\left[\begin{array}{ccc}{C}_{11}& C_{13}& 0\\ C_{13}& C_{33}& 0\\ 0& 0& C_{44}\end{array}\right]\left[\begin{array}{c}{D}_{11}\\ D_{33}\\ D_{31}\end{array}\right].$

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

$\frac{D}{Dt}\left[\begin{array}{c}{\sigma }_{yy}\\ {\sigma }_{xx}\\ {\sigma }_{xy}\end{array}\right]=\left[\begin{array}{ccc}{C}_{11}& C_{13}& 0\\ C_{13}& C_{33}& 0\\ 0& 0& C_{44}\end{array}\right]\left[\begin{array}{c}{D}_{yy}\\ D_{xx}\\ D_{xy}\end{array}\right].$

Rearranging,

$\frac{D}{Dt}\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right]=\left[\begin{array}{ccc}{C}_{33}& C_{13}& 0\\ C_{13}& C_{11}& 0\\ 0& 0& C_{44}\end{array}\right]\left[\begin{array}{c}{D}_{xx}\\ D_{yy}\\ D_{xy}\end{array}\right].$

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=\frac{D}{Dt}({\sigma }_{yy})=C_{13}{D}_{xx}+C_{11}Dyy$

or,

$D_{yy}=-\frac{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

$\frac{D}{Dt}(\mathit{\sigma })={𝐑}_{\text{lam}}\frac{D}{Dt}({\mathit{\sigma }}_{\text{lam}}){𝐑}_{\text{lam}}^T;𝐃={𝐑}_{\text{lam}}{𝐃}_{\text{lam}}{𝐑}_{\text{lam}}^T.$