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

Problem 1: Part 9

Construct a laminar coordinate system at the blue point. (Use equations 9.3.16 and 9.3.17 from the book chapter). Assume that the blue point is at the center of element 5.

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

$\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} = \begin{bmatrix} 0.55000 \\ 0.95263 \\ 0.95263 \\ 0.55000 \end{bmatrix}$

$\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} = \begin{bmatrix} 0.50000 \\ 0.86603 \\ 0.86603 \\ 0.50000 \end{bmatrix}$

$\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} = \begin{bmatrix} 0.60000 \\ 1.0392 \\ 1.0392 \\ 0.60000 \end{bmatrix}$

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

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

The isoparametric map is

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

Plugging in the numbers, we get

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

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

$\begin{align} x_{,\xi} = \frac{\partial x}{\partial \xi} &= 0.20131 + 0.018301 \eta \\ y_{,\xi} = \frac{\partial y}{\partial \xi} &= -0.20131 - 0.018301 \eta~. \end{align}$

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

$\begin{align} x_{,\xi} = \frac{\partial x}{\partial \xi} &= 0.20131 \\ y_{,\xi} = \frac{\partial y}{\partial \xi} &= -0.20131 ~. \end{align}$

Therefore,

$\mathbf{x}_{,\xi} = 0.20131 \mathbf{e}_x - 0.20131 \mathbf{e}_y, ~~ \lVert\mathbf{x}_{,\xi}\rVert_{}= \sqrt{0.20131^2 + 0.20131^2} = 0.28470 ~.$

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

$\widehat{\mathbf{e}}_x = \cfrac{\mathbf{x}_{,\xi}}{\lVert\mathbf{x}_{,\xi}\rVert_{}} = 0.70711 \mathbf{e}_x - 0.70711 \mathbf{e}_y ~.$

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

$\widehat{\mathbf{e}}_y = \mathbf{e}_z\times\widehat{\mathbf{e}_x} = 0.70711 \mathbf{e}_x + 0.70711 \mathbf{e}_y ~.$

A plot of the vectors is shown in Figure 13.

Figure 13. Laminar base vectors.

The Maple code for this calculation is shown below.

> #
> # Shape functions
> #
> N1m := 1/4*(1-xi)*(1-eta):
> N2m := 1/4*(1+xi)*(1-eta):
> N1p := 1/4*(1-xi)*(1+eta):
> N2p := 1/4*(1+xi)*(1+eta):
> N := linalg[matrix](1,4,[N1m,N2m,N1p,N2p]);
> #
> # Compute local laminar basis vectors at the blue point
> #
> # Isoparametric map
> #
> x := x1m*N1m + x2m*N2m + x1p*N1p + x2p*N2p;
> y := y1m*N1m + y2m*N2m + y1p*N1p + y2p*N2p;
> #
> # Derivative of Isoparametric map
> #
> dx_dxi := diff(x, xi);
> dy_dxi := diff(y, xi);
> #
> # Jacobian evluated at blue point
> #
> dx_dxi_cen := subs(eta=0, dx_dxi);
> dy_dxi_cen := subs(eta=0, dy_dxi);
> #
> # Local laminar basis vectors at the blue point
> #
> norm_dx_dxi_cen := simplify(sqrt(dx_dxi_cen^2 + dy_dxi_cen^2));
> ehat_x := vector([simplify(dx_dxi_cen/norm_dx_dxi_cen),
> simplify(dy_dxi_cen/norm_dx_dxi_cen), 0]);
> ehat_z := vector([0, 0, 1]);
> ehat_y := crossprod(ehat_z, ehat_x);