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

Problem 1: Part 10

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 node 1 of the element (global node 5), the $𝐓$ matrix is

${𝐓}_1=\left[\begin{array}{ccc}1& 0& 0.0866\\ 0& 1& -0.05\\ 1& 0& -0.0866\\ 0& 1& 0.05\end{array}\right].$

Therefore, 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]=\left[\begin{array}{c}1.0866\\ 1.95\\ 0.9134\\ 2.05\end{array}\right].$

For node 2 of the element (global node 6), the $𝐓$ matrix is

${𝐓}_2=\left[\begin{array}{ccc}1& 0& 0.05\\ 0& 1& -0.0866\\ 1& 0& -0.05\\ 0& 1& 0.0866\end{array}\right].$

Therefore, 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]=\left[\begin{array}{c}2.05\\ 0.9134\\ 1.95\\ 1.0866\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}$

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

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

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

The derivatives with respect to $\xi$ and $\eta$ are

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

Therefore,

$𝐉=\left[\begin{array}{cc}0.20131+0.018301\eta & -0.20131-0.018301\eta \\ 0.0683+0.018301\xi & 0.0683-0.018301\xi \end{array}\right].$

Since we are only interested in the point ($\xi ,\eta$) = (0,0), we can compute $𝐉$ at this point

$𝐉=\left[\begin{array}{cc}0.20131& -0.20131\\ 0.0683& 0.0683\end{array}\right].$

The inverse is

${𝐉}^{-1}=\left[\begin{array}{cc}2.4837& 7.3205\\ -2.4837& 7.3205\end{array}\right].$

The derivatives of $𝐯$ with respect to $\xi$ and $\eta$ are

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

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

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

Therefore,

$\left[\begin{array}{c}\frac{\partial v_x}{\partial x}\\ \frac{\partial v_x}{\partial y}\end{array}\right]=\left[\begin{array}{c}0.74184\\ -1.7418\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}-0.74184\\ 1.7418\end{array}\right].$

Therefore, the velocity gradient at the blue point is

$\mathrm{\nabla }𝐯=\left[\begin{array}{cc}0.74184& -1.7418\\ -0.74184& 1.7418\end{array}\right]$

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]]);