# Nonlinear finite elements/Homework 6/Solutions/Problem 1/Part 10

## Problem 1: Part 10

Assume that the velocities at nodes 5 and 6 are ${\displaystyle \mathbf {v} _{5}}$  and ${\displaystyle \mathbf {v} _{6}}$ , respectively. Compute the velocity gradient at the blue point.

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 node 1 of the element (global node 5), the ${\displaystyle \mathbf {T} }$  matrix is

${\displaystyle \mathbf {T} _{1}={\begin{bmatrix}1&0&0.0866\\0&1&-0.05\\1&0&-0.0866\\0&1&0.05\end{bmatrix}}~.}$

Therefore, the velocities of the slave nodes are

${\displaystyle {\begin{bmatrix}v_{1-}^{x}\\v_{1-}^{y}\\v_{1+}^{x}\\v_{1+}^{y}\end{bmatrix}}={\begin{bmatrix}1.0866\\1.95\\0.9134\\2.05\end{bmatrix}}~.}$

For node 2 of the element (global node 6), the ${\displaystyle \mathbf {T} }$  matrix is

${\displaystyle \mathbf {T} _{2}={\begin{bmatrix}1&0&0.05\\0&1&-0.0866\\1&0&-0.05\\0&1&0.0866\end{bmatrix}}~.}$

Therefore, the velocities of the slave nodes are

${\displaystyle {\begin{bmatrix}v_{2-}^{x}\\v_{2-}^{y}\\v_{2+}^{x}\\v_{2+}^{y}\end{bmatrix}}={\begin{bmatrix}2.05\\0.9134\\1.95\\1.0866\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}}}

Plugging in the values of the nodal velocities and the shape functions, we get

{\displaystyle {\begin{aligned}v_{x}&=0.27165(1-\xi )(1-\eta )+0.51250(1+\xi )(1-\eta )\\&+0.22835(1-\xi )(1+\eta )+0.48750(1+\xi )(1+\eta )\\v_{y}&=0.48750(1-\xi )(1-\eta )+0.22835(1+\xi )(1-\eta )\\&+0.51250(1-\xi )(1+\eta )+0.27165(1+\xi )(1+\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 isoparametric map is

{\displaystyle {\begin{aligned}x&=0.12500(1-\xi )(1-\eta )+0.21651(1+\xi )(1-\eta )\\&+0.15000(1-\xi )(1+\eta )+0.25981(1+\xi )(1+\eta )\\y&=0.21651(1-\xi )(1-\eta )+0.12500(1+\xi )(1-\eta )\\&+0.25981(1-\xi )(1+\eta )+0.15000(1+\xi )(1+\eta )~.\end{aligned}}}

The derivatives with respect to ${\displaystyle \xi }$  and ${\displaystyle \eta }$  are

{\displaystyle {\begin{aligned}{\frac {\partial x}{\partial \xi }}&=0.20131+0.018301\eta &{\frac {\partial y}{\partial \xi }}&=-0.20131-0.018301\eta \\{\frac {\partial x}{\partial \eta }}&=0.0683+0.018301\xi &{\frac {\partial y}{\partial \eta }}&=0.0683-0.018301\xi \end{aligned}}}

Therefore,

${\displaystyle \mathbf {J} ={\begin{bmatrix}0.20131+0.018301\eta &-0.20131-0.018301\eta \\0.0683+0.018301\xi &0.0683-0.018301\xi \end{bmatrix}}~.}$

Since we are only interested in the point (${\displaystyle \xi ,\eta }$ ) = (0,0), we can compute ${\displaystyle \mathbf {J} }$  at this point

${\displaystyle \mathbf {J} ={\begin{bmatrix}0.20131&-0.20131\\0.0683&0.0683\end{bmatrix}}~.}$

The inverse is

${\displaystyle \mathbf {J} ^{-1}={\begin{bmatrix}2.4837&7.3205\\-2.4837&7.3205\end{bmatrix}}~.}$

The derivatives of ${\displaystyle \mathbf {v} }$  with respect to ${\displaystyle \xi }$  and ${\displaystyle \eta }$  are

{\displaystyle {\begin{aligned}{\frac {\partial v_{x}}{\partial \xi }}&=0.5+0.018301\eta &{\frac {\partial v_{x}}{\partial \eta }}&=-0.0683+0.018301\xi \\{\frac {\partial v_{y}}{\partial \xi }}&=-0.5+0.018301\eta &{\frac {\partial v_{y}}{\partial \eta }}&=0.0683+0.018301\xi \end{aligned}}}

At (${\displaystyle \xi ,\eta }$ ) = (0,0), we get

{\displaystyle {\begin{aligned}{\frac {\partial v_{x}}{\partial \xi }}&=0.5&{\frac {\partial v_{x}}{\partial \eta }}&=-0.0683\\{\frac {\partial v_{y}}{\partial \xi }}&=-0.5&{\frac {\partial v_{y}}{\partial \eta }}&=0.0683\end{aligned}}}

Therefore,

${\displaystyle {\begin{bmatrix}{\frac {\partial v_{x}}{\partial x}}\\{\frac {\partial v_{x}}{\partial y}}\end{bmatrix}}={\begin{bmatrix}0.74184\\-1.7418\end{bmatrix}};\qquad {\begin{bmatrix}{\frac {\partial v_{y}}{\partial x}}\\{\frac {\partial v_{y}}{\partial y}}\end{bmatrix}}={\begin{bmatrix}-0.74184\\1.7418\end{bmatrix}}~.}$

Therefore, the velocity gradient at the blue point is

${\displaystyle {\boldsymbol {\nabla }}\mathbf {v} ={\begin{bmatrix}0.74184&-1.7418\\-0.74184&1.7418\end{bmatrix}}}$

The Maple code for doing the above calculations is shown below.

> # > # Compute velocity gradient > # > # Map from master to slave nodes > # > T1 := linalg[matrix](4,3,[1,0,y1-y1m, > 0,1,x1m-x1, > 1,0,y1-y1p, > 0,1,x1p-x1]); > T2 := linalg[matrix](4,3,[1,0,y2-y2m, > 0,1,x2m-x2, > 1,0,y2-y2p, > 0,1,x2p-x2]); > v1slave := evalm(T1&*v1master); > v2slave := evalm(T2&*v2master); > # > # Velocity interpolation within element > # > vx := v1slave[1,1]*N[1,1] + v2slave[1,1]*N[1,2] + > v1slave[3,1]*N[1,3] + v2slave[3,1]*N[1,4]; > vy := v1slave[2,1]*N[1,1] + v2slave[2,1]*N[1,2] + > v1slave[4,1]*N[1,3] + v2slave[4,1]*N[1,4]; > # > # Derivatives of velocity within element (at the blue point) > # > dvx_dxi := subs(xi=0,eta=0,diff(vx, xi)); > dvy_dxi := subs(xi=0,eta=0,diff(vy, xi)); > dvx_deta := subs(xi=0,eta=0,diff(vx, eta)); > dvy_deta := subs(xi=0,eta=0,diff(vy, eta)); > # > # Jacobian of the isoparametric map and its inverse > # (at the blue point) > # > dx_dxi := subs(xi=0,eta=0,diff(x, xi)); > dy_dxi := subs(xi=0,eta=0,diff(y, xi)); > dx_deta := subs(xi=0,eta=0,diff(x, eta)); > dy_deta := subs(xi=0,eta=0,diff(y, eta)); > J := linalg[matrix](2,2,[dx_dxi,dy_dxi,dx_deta,dy_deta]); > Jinv := inverse(J); > # > # Compute velocity gradient (at the blue point) > # > dVx_dxi := linalg[matrix](2,1,[dvx_dxi,dvx_deta]); > dVy_dxi := linalg[matrix](2,1,[dvy_dxi,dvy_deta]); > dVxdx := evalm(Jinv&*dVx_dxi); > dVydx := evalm(Jinv&*dVy_dxi); > GradV := linalg[matrix](2,2,[dVxdx[1,1],dVxdx[2,1], >dVydx[1,1],dVydx[2,1]]);