Finite Element Analysis

• Finite Element Analysis
  • What is FE Method?
  • Fundamentals
  • Direct Finite Element Method (Matrix Stiffness Method)
• \end{bmatrix}
  • Truss Analysis (2D Direct FEM)
  • Strong Form
  • Weak Form
  • Approximations
  • FEA derivaiton
  • 2D Heat Conductance (Steady State)
  • 2D Approximations

References: Belytschko, Hartmann, Matrix Partitioner/Symmetry Maker

Assignments: 1234 5

What is FE Method?

FE is a method of:

• restriction: it is based on restricting the structure into shape functions
• substitution: replaces original load with work-equivalent load
• projection: FE solution is the projection of solution onto strain energy exact solution. In a vertical projection the length of a shadow is always less than the length of the original vector (see Bessel's inequality [175]); this implies that the strain energy of the FE solution is always less than the strain energy of the exact solution. An engineer would say that the FE solution overestimates the stiffness of the structure.
• energy method: based on principle of virtual work
• method of appeoximate influence functions. Displacement and stress at point x is scalar product of the applied load and corresponding influence (green's) function $u(x)=\int_0^lG_0(y,x)p(y)dy$. FE replaces exact Green's functions with approimate Green's functions $G_o^h$ , such that $u(x)-u_h(x)=\int_0^l[G_0(y,x)-G_o^h(y,x)]p(y)dy$

Fundamentals

The finite element method (FEM) consists of the following five steps:

1. Preprocessing: subdividing the problem domain into finite elements
2. Element formulation: development of equations for elements
3. Assembly: obtaining the equations of the entire system from the equations of individual elements.
4. Solving the equations.
5. Postprocessing: determining quantities of interest, such as stresses and strains, and obtaining visualizations of the response

Matrix Partitioning: how to simplify a matrix so that it is easier to solve by breaking it down into partitions. With the function $Kd=f+r$, the d vector, or the vector of nodal degrees of freedom or displacement, is broken down into two types. $d^E$ identifies essential or fixed boundary conditions and $d^F$ identifies free boundary conditions. Steps are simply as follows:

1. Place columns in K matrix in the corresponding order to the order of degrees of freedom, given that the essential (fixed) degrees of freedom will be moved to the top
2. Reorder the rows to make it symmetric again, similarly to as before

Direct Finite Element Method (Matrix Stiffness Method)

Nodal internal forces are denoted as $F_1^e$ and $F_2^e$, while nodal displacements are denoted by $u_1^e$ and $u_2^e$. For an elastic bar, assuming linear analysis, elastic material, and constant stiffness EA, the following are used to relate $F$ and $u$:

1. Equilibrium of the element $F_1^e=F_2^e$
2. Hooke's Law $\sigma_x=E\varepsilon$ ($F_2^e=\sigma_xA$)
3. Strain-displacement relationship $\varepsilon=\frac{\Delta{L}}{L}=\frac{u_2^e-u_1^e}{L}$

Thus, $F_1^e=\frac{AE}{L}u_1^e-\frac{AE}{L}u_2^e$ and $F_2^e=\frac{AE}{L}(-u_1^e)+\frac{AE}{L}u_2^e$. This is simply:

$$\begin{bmatrix} F_1^e \\ F_2^e \end{bmatrix} = \frac{AE}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} u_1^e \\ u_2^e \end{bmatrix}$$

The first $F^e$ matrix is the element internal force matrix, the other the element stiffness matrix $K^e$, and the final $d^e$ matrix the element degree of freedom matrix. Compact, $F^e=K^ed^e$. This is normally denoted with a tilda on the bottom to show that they are matrices, not scalars.

Considering the analysis of several connected bar elements, such that a mesh with $n^e$ element and nn nodes such that $nn=n^e+1$, the governing principle is equilibrium. Enforcing equilibrium at each node can provide a series of local equations. These can then be mapped from local to global with the transformation matrix $L^e$, such that $d^-=L^td$. $L^e$ is gather matrix, which uses 1s and 0s to select the local displacements from the global set. Since $L^e$ is different for each element and is a series of 1 and 0 and $(L^e)^T$ transforms $F^e$ to be the same size of the boundary conditions, we have, we have $K=\sum_{e=1}^nL_e^TKL^e$.

The connectivity matrix shows which nodes are connected per element. See here for example:

(derivation)

Truss Analysis (2D Direct FEM)

Trusses can be analysed element by element. Assuming only axial deformation cause stress, we get the following system, with $\phi$ being the counterclockwise angle of rotation of node 1 to node 2 in radians:

$$\begin{bmatrix}F_{1x}\\F_{1y}\\F_{2x}\\F_{2y}\\\end{bmatrix}= \frac{AE}{L}\begin{bmatrix} \cos^2\theta&\cos\theta\sin\theta&-\cos^2\theta&-\cos\theta\sin\theta\\ \cos\theta\sin\theta&\sin^2\theta&-\cos\theta\sin\theta&-\sin^2\theta\\ -\cos^2\theta&-\cos\theta\sin\theta&\cos^2\theta&\cos\theta\sin\theta\\ -\cos\theta\sin\theta&-\sin^2\theta&\cos\theta\sin\theta&\sin^2\theta\\ \end{bmatrix} \begin{bmatrix}u_{1x}\\u_{1y}\\u_{2x}\\u_{2y}\\\end{bmatrix}$$

This is alternatively the following:

$$K^e=T^T\frac{AE}{L}\begin{bmatrix} 1&-1\\-1&1\\\end{bmatrix} T$$

$$T=\begin{bmatrix} \cos&\sin&0&0\\ 0&0&\cos&\sin\\\end{bmatrix}$$

To assemble truss elements, reactions and loads are taken positive in the positive x-dir and y-dir. The following steps shall be applied (See example):

1. Number nodes and elements globally
2. Define conncetivity of the mesh (connectivity matrix)
3. Build equilibrium equations at each node
4. Repeat for all nodes

It is notable that the strain $\varepsilon=\frac{\Delta{u}}{L}=B^ed^e$. This refers to displacement at an angle. In technicality, $\Delta{u}=u_{2x}^e\prime-u_{1x}^e\prime=B^ed^e$. This are defined below:

$$\begin{bmatrix} u_{1x}\prime\\ u_{2x}\prime\\\end{bmatrix}= \begin{bmatrix} \cos\phi{u_{1x}^e}+\sin\phi{u_{1y}^e}\\ \cos\phi{u_{2x}^e}+\sin\phi{u_{2y}^e}\\\end{bmatrix}$$

$$B^e=\frac{\begin{bmatrix}-1 & 1\\\end{bmatrix}T^e}{L}$$

$$T^e=\begin{bmatrix} \cos\phi & \sin\phi &0&0\\ 0&0&\cos\phi&\sin\phi\\\end{bmatrix}$$

See full derivation

Strong Form

FEA follows the following progression:

Within this progression, functions are part of a continuity.

• $C^2$ function space is a set of functions w/ continuous 2nd derivatives
• $C^1$ function space is a set of functions w/ continuous 1st derivatives
• $C^0$ function space is a set of functions that are continuous
• $C^{-1}$ function space is a set of piecewise continuous functions

The strong form is the governing differential equation of concern, with boundary conditions and initial conditions. In contrast, the weak form is an integral form of these equations, which is needed to formulate the finite element method. Approximation functions are combined with the weak form to obtain the discrete finite element equations. Finite element method requires:

1. Strong form
2. Weak form
3. Approximation functions

In the following, we aspire to find the stress distribution $\sigma{(x)}$ in the bar. This stress is caused by deformation of the body, or displacements u(x). The bar is subject to:

• Body force b(x), which is force along the length of the bar (units N/m). This could be gravity of thermal heating
• Traction force $\bar{t}$, which is the force applied at the end of the bar (units $\frac{N}{m^2}$).

Conditions to satisfy are as follows:

1. Must be in equilibrium
2. Must satisfy Hooke's law -relation of stress and strain
3. Must have compatible displacement field
4. Must satisfy strain-displacement equation

Steps considering equilibrium of internal force p(x) and external force b(x):

1. $-p(x)+b(x+\frac{\Delta{x}}{2})\Delta{x}+p(x+\Delta{x})=0$
2. $\frac{p(x+\Delta{x})-p(x)}{\Delta{x}}+b(x+\frac{\Delta{x}}{2})=0$
3. $\frac{dp(x)}{dx}+b(x)=0$
4. We know $p(x)=A(x)\sigma{(x)}$, as well as that $\varepsilon{(x)}=\frac{du}{dx}$ and $\sigma{(x)}=E(x)\varepsilon{(x)}$. Plugging in, the form $\frac{d}{dx}(AE\frac{du}{dx})+b=0,0<x<l$ is determined
   • This assumes linearity in the definition of the strain (3.3) and the stress–strain law (3.4). Compatibility is satisfied by requiring the displacement to be continuous. More will be said later about the degree of smoothness, or continuity, which is required.
5. To solve, the boundary conditions at the ends of the bar $\sigma{(\textrm{end})}=E\frac{du}{dx}=-\bar{t}$ and $u(\textrm{other end})=\bar{u}$ are used
   • $\bar{t}$ has the same units as stress (force/area), but its sign is positive when it acts in the positive x-direction regardless of which face it is acting on, whereas the stress is positive in tension and negative in compression, so that on a negative face a positive stress corresponds to a negative traction. Note that either the load or the displacement can be specified at a boundary point, but not both.

The end solution, with differential equation and boundary conditions, is the strong form.

In heat flow, there is no such thing as heat. It is heat difference. The goal of the following is to determine the heat temp distribution. We are given:

• A(x): area normal to heat flow direction
• s(x): heat generated per unit thickness of wall (l). This is the heat source.
  • rate of heat generation per unit length = energy/time*s= $\frac{J}{m\cdot{s}}$
• Heat flux q(x): rate of heat flow across a surface. Unit: $\frac{J}{m^2\cdot{s}}=\frac{W}{m^2}$

Energy in=energy out. We have the following procedure

1. $s(x+\Delta{x}/2)\Delta{x}+q(x)A(x)-q(x+\Delta{x})A(x+\Delta{x})=0$
2. $\frac{q(x+\Delta{x})A(x+\Delta{x})-q(x)A(x)}{\Delta{x}}=s(x+\Delta{x}/2)$
3. $\frac{d(qA)}{dx}=s$
4. We know that $q=\frac{-kdT}{dx}$, with T as the temperature and h asthe thermal conductivity. $\therefore\frac{d}{dx}(Ah\frac{dT}{dx})+s=0$, $0\leq{x}<1$
5. Boundary conditions= temperature $T-\bar{T}$ at end, and flux $q=-\bar{q}$ at other end. Flux $\bar{q}$ positive if heat flows out of bar

1D Elastic Body looks something like this:

$$AEu_{xx}+b=0\tag{Elastic Strong Form}$$

By equilibrium, $-\frac{d(Aq)}{dx}+S=0$. Since by Fourier's Law $q=-k\frac{dT}{dx}$ (slope negative because it goes from hot to cold), we have

$$AkT_{xx}+S=0\tag{Heat Transfer Strong Form}$$

See derivation here

Weak Form

The weak form is defined using weight functions for convenience. They also transform the result into only having first derivatives, as this would make the trial solutions easier to determine.

To find the weak form, follow the procedure below (example.

1. Define Equation 1 as your strong form, Equation 2 as an essential boundary condition and Equation 3 as a traction boundary condition (related to stress at the end of the bar from applied force)
2. Multiply EQN 1 by w(x) and integrate over the domain. Let w(x) be any function of sufficient smoothness (meaning integrals exist and are finite)
3. Apply integration by parts
4. Using traction boundary condition, we can evaluate the boundary term @ x=L
5. Limit the scope of w(x). Wherever there is no displacement, let w(x)=0. For instance, w(0)=0 is a possibility
6. Find the weak form is; remember it must belong to a function space with boundary conditions satisfied. Note if a function is in $C^0$, it is likewise in $H\prime$

Alternative procedure:

1. Multiply governing equation and traction boundary condition by an arbitrary function w(x). w(x) is a weight function. $\forall{w}$ denotes that it is arbitrary. Weight function is enforcer; whatever it multiplies is enforced to be zero by its arbitrariness.
   • $w[\frac{d}{dx}(AE\frac{du}{dx})+b]dx=0,\forall{w}$
   • Note: not applied to displacement condition, as it is easy to figure out
2. Integrate over domain (likely $0\to{l}$). Note: just multiply traction coefficient by area
   • $wAE\frac{du}{dx}+wA\bar{t}|_{x=0},\forall{w}$
   • $\int_0^Lw[\frac{d}{dx}(AE\frac{du}{dx})+b]dx=0,\forall{w}$
3. Convenient for w(l)=0. Impose restriction on set of weight functions
4. For convenience, transform $\int_0^Lw[\frac{d}{dx}(AE\frac{du}{dx})+b]dx=0,\forall{w}$ to only contain first derivatives using integration by parts
5. Substitute boundary condition $E\frac{du}{dx}=\sigma$

From the derivation, we state that the trial solution that satisfies the below for all smooth w(x) with w(l)=0 is the solution. Thus, the problem is to find u(x) among the smooth functions that satisfy $u(l)=\bar{u}$ such that the following is satisfied:

$$\int_0^l{\frac{dw}{dx}AE\frac{du}{dx}dx}=(wA\bar{t})_{x=0}+\int_0^l{wbdx},\forall{w}\epsilon{w(l)=0}\tag{Elastic Weak Form}$$

In other words, the trial solutions u(x) that satisfy this are the strong form. A trial solution that is smooth and satisfies the essential boundary conditions (displacement) are called admissible, and weight functions that are smoth and vanish on essential boundaries are admissible

Approximations

For FEM to converge, must be continuous and complete

• Continuous: trial solutions and weight functions are sufficiently smooth. Trial solutions are the set of admissible solutions u(x) which satisfy certain conditions
• Completeness: the series must approximate the smooth function with enough accuracy. This is sufficient if as the element sizes approach zero, the trial solutions and weight functions and their derivatives up to and including the highest order derivative appearing in the weak form are capable of assuming constant values. For elasticity, this requires that the displacement field and its derivative can take constant values so that the finite elements can represent rigid body motion and constant strain states exactly.

Approximate functions are of the form $u(x)=\alpha_1+\alpha_2x_1+\alpha_3x_1^2\cdots$. In matrix form, we have $\begin{bmatrix}u_1^e\\u_2^e\end{bmatrix}=\begin{bmatrix}1&x_1^e\\1&x_2^e\end{bmatrix}\begin{bmatrix}\alpha_1^e\\\alpha_2^e\end{bmatrix}$, or $u^e=M^e\alpha^e$. In other words, $u_1^e(x)=\begin{bmatrix}1&x\end{bmatrix}\begin{bmatrix}\alpha_1^e\\\alpha_2^e\end{bmatrix}$, or $u_1^e=p(x)\alpha^e$. From here, the shape function $N^e(x)=p(x)(M^e)^{-1}$ is derived, with $u^e(x)=N^e(x)d^e$

$N_1^e(x_1^e)=1$, while $N_1^e(x_2^e)=0$. Generally, $N_I^e(x_J^e)=\delta_{IJ}$, with $\delta_{IJ}=\begin{Bmatrix}1&\textrm{ if I=J}\\0&\textrm{ if I not J}\end{Bmatrix}$

From the weak form, derivatives of trial solutions and weight functions were required. For the linear approximation, $\frac{du^e}{dx}=\frac{d(N^ed^e)}{dx}=\begin{bmatrix}\frac{dN_1^e}{dx}&\frac{dN_2^e}{dx}\end{bmatrix}\begin{bmatrix}u_1^e\\u_2^e\end{bmatrix}=B^ed^e$. Since we know the slopes of the shape functions, we can say:

$$B^e=\begin{bmatrix}\frac{dN_1^e}{dx}&\frac{dN_2^e}{dx}\end{bmatrix}=\frac{1}{l_e}\begin{bmatrix}-1&+1\end{bmatrix}$$

The above values of N are for inear approximation. The derivation is similar for quadratic approximations, except that the expressions with alpha are to the second order. The functions of N are determined in such a way that they are zero at every node except the ones where they should equal one.

Weight functions can be approximated similarly (Galerkin FEM), such that $w^e(x)=N^e(x)w^e$ and $\frac{dw^e}{dx}=B^ew^e$

Globally, we have:

$$u^h=\sum_{e=1}^{n_{el}}N^ed^e=(\sum_{e=1}^{n_{el}}N^eL^e)d\tag{Global Trial Solutions}$$

$$w^h=\sum_{e=1}^{n_{el}}N^ew^e=(\sum_{e=1}^{n_{el}}N^eL^e)w\tag{Global Weight functions}$$

$$N=\sum_{e=1}^{n_{el}}N^eL^e\tag{Global Shape functions}$$

FEA derivaiton

To solve weak form, it is to be evaluated for the finite element trial solutions and weight functions. Then, by invoking the arbitrariness of the weight functions, a set of linear equations will be produced. This is demonstrated here

Steps to Galerkin:

1. Approximate $w(r)=N^e(r)w^e$ and dependent variable $C(r)=N^e(r)d$.
   • Remember, $\frac{dw}{dr}=\frac{dN^e(r)}{dr}w^e=B^ew^e=w^{e^T}B^{e^T}$, $\frac{dC}{dr}=B^ed$
2. Sub these values into the weak form
3. Set up a form of summing the system for all elements, while integrating each integral over $x_2$ and $x_1$ of their respective element
4. Sub in equations for N and $B^e$, remembering that $x_2-x_1=h^e$
5. Solve

For L3 quadratic elements, we instead have:

$$\frac{du}{dx}=\begin{bmatrix}N_{1x}^e&N_{2x}^e&N_{3x}^e\end{bmatrix}\begin{bmatrix}u_{1}^e\\u_2^e\\u_3^e\end{bmatrix}$$

For convenience, parent coordinates $\zeta$ are used. See below:

From this, we see:

$$N_1^e(\zeta)=\frac{(\zeta-0)(\zeta-1)}{(-1-0)(-1-1)}=\frac{\zeta{(\zeta-1)}}{2}$$

$$N_2^e(\zeta)=(\zeta+1)(1-\zeta)$$

$$N_3^e(\zeta)=\frac{\zeta{(\zeta+1)}}{2}$$

We know $\frac{dN_I^e}{dx}=\frac{dN_I^e}{d\zeta}\frac{d\zeta}{dx}=\frac{\frac{dN_I^e}{d\zeta}}{h^e}=\frac{N_1^e}{h^e}$. Furthermore, the acclivity of N steepens. $N_{1x}^e=\frac{\zeta-0.5}{h^e}$, $N_{2x}^e=\frac{-2\zeta}{h^e}$, and $N_{3x}^e=\frac{\zeta+0.5}{h^e}$. Altogether,

$$K^e=\int_{-1}^1{B^{e^T}(\zeta)AEB^e(\zeta)\frac{dx}{d\zeta}d\zeta}=\int_{-1}^1{B^{e^T}(\zeta)AEB^e(\zeta)h^ed\zeta}$$

Eventually, it is seen that $K^e=\frac{AE}{12h^e}\begin{bmatrix}14&-16&2\\-16&32&-16\\2&-16&14\end{bmatrix}$. This matrix has the following properties:

• sum of rows andcolumns equals zero
• symmetric

In the quadratic realm, the body force vector is similarly generated. This vector starts as $f_{\Omega}^e=\int_{x_1^e}^{x_2^e}{N^{e^T}b(x)dx}$. b(x) must be approximated by $b^h$ as follows:

$$b^h(\zeta)=\sum_1^{\zeta}{N_I^e(\zeta)b_I^e=N^eb^e},b^e=\begin{bmatrix}b_1^e\\b_2^e\\b_3^e\end{bmatrix}$$

From there, the terms are subbed together to form $f_{\Omega}^e=\int_{-1}^1{N^{e^T}(\zeta)N^e(\zeta)b^e\frac{dx}{d\zeta}d\zeta}=f_{\Omega}^e=\int_{-1}^1{N^{e^T}(\zeta)N^e(\zeta)b^eh^ed\zeta}$. This generally leads to $f_{\Omega}^e=M^eb^e$, with the integral leading to $M^e=\frac{h^e}{15}\begin{bmatrix}4&2&-1\\2&16&2\\-1&2&4\end{bmatrix}$

There is always error in finite element analysis. The convergence is how much the error goes to zero.

The L2 error, or the normalized error, is related to the displacement. See below:

$$e_{L2}=\sqrt{\int_0^L{[u^{ex}(x)-u^h]^2dx}}$$

The energy error is related to the strain energy. See below:

$$e_{L2}=\sqrt{\int_0^L{[\varepsilon^{ex}(x)-\varepsilon^h]^2dx}}=\sqrt{\int_0^L{[\frac{du}{dx}^{ex}-\frac{du}{dx}^h]^2dx}}$$

$e^{L2}\alpha{h^{p+1}}$. In other words, L2 error varies with the measure of the element size, with p as the polynomial order of approximation. From the below, we can see that if we reduce the linear element size by 2, the error reduces by $2^{p+1}=2^2=4$:

The energy error is 1 order of magnitude lower than L2 error. Decreases by error value to the power of p. This implies that error in $\frac{d_u^h}{dx}$, such as strain, stress, and flux, decreases much more slowly than that of $u^h(x)$. See below:

2D Heat Conductance (Steady State)

See notes on 2D heat conductance.

In steady state, the temp in domain doesn't change over time. This is a fundamental assumption. See below:

Recalling Fourier's Law, $q=-D\nabla{T}$. In matrix form, this is equal to $\begin{bmatrix}q_x\\q_y\end{bmatrix}=-\begin{bmatrix}k_{xx}&k_{xy}\\k_{yx}&k_{yy}\end{bmatrix}\begin{bmatrix}\frac{\partial{T}}{\partial{x}}\\\frac{\partial{T}}{\partial{y}}\end{bmatrix}$

The domain $\Omega$ has boundaries $\Gamma_T$ (temperature) and $\Gamma_q$ (heat flux). See below:

The goal is to find $T(x,y)$ such that $\nabla{(k\nabla{T})}+S=0$, x and y $\epsilon\Omega$. Assuming an isotropic material, $\vec{q}=-k\nabla{T}$. The boundary conditions are $T(x,y)=\bar{T}(x,y),(x,y)\epsilon\Gamma_T$ (essential BC) and ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at position 12: \vec{q}\dot\̲v̲e̲c̲{n}=\bar{q},(x,… (traction BC)

Note that $\Gamma_T\cup\Gamma_q=\Gamma$ and $\Gamma_T\cap\Gamma_q=0$; there is no overlapping.

$\bar{q}$ is magnitude in normal direction. $\vec{q}$ is in any arbitrary direction, and $\vec{n}$ transforms it. $\vec{q}=q\vec{i}+q\vec{j}$. We have the following:

To obtain the weak form of the 2D heat conductance strong form:

1. Multiply strong form by weight function w(x,y) and integrate over domain: $\int_{\Omega}w\nabla{(k\nabla{T})d\Omega}+\int_{\Omega}wsd\Omega=0\forall{w}$
2. Use green's Formula to modify the first term to move derivates over to w
3. Move to matrix form
4. Divide integral over $\Gamma$ to 2 integrals over $\Gamma_T$ and $\Gamma_q$

Strong form and weak form are presented below:

2D Approximations

See notes on approximations

For triangular elements,

$$d^e=\begin{bmatrix}1&x_1^e&y_1^e\\1&x_2^e&y_2^e\\1&x_3^e&y_3^e\end{bmatrix}\begin{bmatrix}\alpha_0\\\alpha_1\\\alpha_2\end{bmatrix}$$

$$N_1^e(x,y)=\frac{1}{2A^2}(x_2^ey_3^e-y_2^ex_3^e-y_{23}x+x_{32}y)$$

$$N_2^e(x,y)=\frac{1}{2A^2}(x_3^ey_1^e-y_3^ex_1^e-y_{31}x+x_{13}y)$$

$$N_3^e(x,y)=\frac{1}{2A^2}(x_1^ey_2^e-y_1^ex_2^e-y_{12}x+x_{21}y)$$

$$A^e=\frac{1}{2}(x_{32}y_{12}-x_{12}y_{32})$$

Note: $x_{IJ}=x_I-x_J$. Given that $T^h(x,y)=N^ed^e=\sum_{I=1}^3N_I^e(x,y)T_I^e$, the gradient of $T^h$ is as follows:

Properties of the shape functions are that the sum of $N_I^e(x,y)=1$. Also, the sum of derivatives of shape function with respet to x or y = 0. Therefore, the sum of rows in $B^e$ are zero

For higher order, there are several formulations. Quadrilateral elements are 4-node rectangular elements, with edges parallel to the axes. See below

The shape functions of these are products of the linear shape functions along the edges.

$$N_1(x,y)=N_1(x)N_1(y)=\frac{1}{A^e}(x-x_2)(y-y_2)$$

$$N_2(x,y)=N_2(x)N_2(y)=\frac{1}{A^e}(x_1-x)(y-y_3)$$

$$N_3(x,y)=N_2(x)N_2(y)=\frac{1}{A^e}(x-x_1)(y-y_1)$$

$$N_4(x,y)=N_1(x)N_2(y)=\frac{1}{A^e}(x-x_2)(y_1-y)$$

The curve of $T_h$ along the quadrilateral looks like the following due to the bilinear term xy (not quadratic):

The accuracy if said Q4 element is between a T3 and T6, but still converges to a linear approximation accuracy. See below:

The case of isoparametric elements are more tricky. For instance, a 4-node quadrilateral with arbitrary edge/side orientation. This must be analysed with a parent coordinate system, using $\eta$ and $\zeta$. See below:

Evidently, $N_1^{L2}(\zeta)=\frac{1}{2}(1-\zeta)$ and $N_1^{L2}(\eta)=\frac{1}{2}(1-\eta)$ etc. Altogether, we have:

$$N_1^e=\frac{1}{4}(\zeta-1)(\eta-1)$$

Using isoparametric mapping, we have the following:

$$GN^{Q4}=\begin{bmatrix}N_{1,\zeta}&N_{2,\zeta}&N_{3,\zeta}&N_{4,\zeta}\\N_{1,\eta}&N_{2,\eta}&N_{3,\eta}&N_{4,\eta}\end{bmatrix}$$

$$J^e=GN^{Q4}\begin{bmatrix}x_1^e&y_1^e\\x_2^e&y_2^e\\x_3^e&y_3^e\\x_4^e&y_4^e\end{bmatrix}$$

$$d^e=\begin{bmatrix}T_1^e\\T_2^e\\T_3^e\\T_4^e\end{bmatrix}$$

$$\nabla{T^h}={J^{e}}^{-1}GNd^e=B^ed^e$$

As we derived a Q4 shape function by multiplying L2 shape functions, we can do similarly for a Q9 by multiplying L3 shape functions and Q16 by multiplying Q4 shape functions (see here). Such a representation is displayed below:

One can imagine that having fewer triangles (as in a T6) is more accurate than a T3. This is correct because it allows for curved surfaces. Such is also seen in Q9. The Q9 has better accuracy than T6, though, but convergence is stoll only 2nd order in energy.

See notes for great Galerkin example

We also have:

$$f_{\Gamma}=\int_{\Gamma_q}{N^e}^T\bar{q}d\Gamma\tag{Boundary Flux}$$

• $\Gamma$ is normally defined along edge using local coordinate system
• If $s(x,y)=N^es$ and $s=\begin{bmatrix}s_1^e\\s_2^e\\s_3^e\end{bmatrix}$, then ParseError: KaTeX parse error: Expected '}', got '&' at position 45: …egin{bmatrix}{2&̲1&1\\1&2&1\\1&1…. Also, $f_{\Omega}^e=\frac{A_s}{3}\begin{bmatrix}1\\1\\1\end{bmatrix}$