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

Problem 1: Part 10

Assume that the velocities at nodes 5 and 6 are $\mathbf {v} _{5}$ and $\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

${\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 $\mathbf {T}$ matrix is

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

${\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 $\mathbf {T}$ matrix is

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

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

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

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

${\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 isoparametric map is

${\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 $\xi$ and $\eta$ are

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

$\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 ($\xi ,\eta$) = (0,0), we can compute $\mathbf {J}$ at this point

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

The inverse is

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

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

${\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 ($\xi ,\eta$) = (0,0), we get

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

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

${\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]]);

Nonlinear finite elements/Homework 6/Solutions/Problem 1