Nonlinear finite elements/Homework 3/Solutions

Problem: The Patch Test
 Given:

Taylor, R.L., Simo, J.C., Zienkiewicz, O.C., and Chan, A.C.H, 1986, The patch test - a condition for assessing FEM convergence,  International Journal for Numerical Methods in Engineering, 22}, pp. 39-62.

Part 1
What is the patch test and why is it used?


 * What is the patch test?

The standard patch test involves setting up a mesh of elements as shown in Figure 1(a) and checking the predicted solution against a known solution. This mesh is called a patch.

The elements should be arbitrary quadrilaterals rather than rectangles. No body forces should be applied. Material properties should be uniform and linear elastic.

The displacements of nodes $$1$$, $$2$$, $$3$$, and $$4$$ are prescribed. The prescribed value depends on the order of the polynomial that is to be represented exactly by the elements. A similar process is followed for three-dimensional elements.

Suppose we wish to have both constant and linear displacements be represented exactly by the elements. Then the prescribed displacement field is

\text{(2)} \qquad { \begin{align} u_x(x, y) & = c_0 + c_1 x + c_2 y \\ u_y(x, y) & = d_0 + d_1 x + d_2 y \end{align} } $$ where ($$c_0, c_1, c_2$$) and ($$d_0, d_1, d_2$$) are constants set by the user. These constants should be nonzero.

Therefore, the prescribed nodal displacements are

\text{(3)} \qquad \begin{align} u_x(x_i, y_i) & = c_0 + c_1 x_i + c_2 y_i \\ u_y(x_i, y_i) & = d_0 + d_1 x_i + d_2 y_i \end{align} $$ where $$(x_i, y_i)$$ are the coordinates of node $$i$$. The finite element solution for these prescribed displacements should be given by equation (2) throughout the patch.

Hence the displacements at nodes $$5$$, $$6$$, $$7$$, and $$8$$ should be given by equation (3)}.

The element strains should be constant and be given by

\begin{align} \varepsilon_{xx} & = \frac{\partial u_x}{\partial x} = c_1 \\ \varepsilon_{yy} & = \frac{\partial u_y}{\partial y} = d_2 \\ \varepsilon_{xy} & = \frac{1}{2}\left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}\right) = c_2 + d_1 \end{align} $$

The element stresses should also be constant.

These results should be satisfied with a high degree of precision.

The paper talks about an additional test that can be used to assess stability (see Figure 1(b)).


 * Why is it used?

For a finite element analysis, we seek elements (and therefore shape functions) that behave such that as the mesh is refined the approximate solution converges to the exact solution.

The patch test is a necessary condition for convergence. If a finite element solution is to converge with refinement, the finite elements must pass the patch test.

In practice, the patch test checks whether an element passes the completeness test.

Suppose the approximate solution over an element is given by

u_h = \sum_{i=1}^n u^e_i N_i ~. $$ Suppose that we would like the element to represent exactly linear and constant functions. Then, in two dimensions, the shape functions ($$N_i$$) are said to be  complete if

u^e_i = c_0 + c_1 x^e + c_2 y^e \qquad \text{implies that} \qquad u_h(x,y) = c_0 + c_1 x + c_2 y $$ where $$c_0$$, $$c_1$$ and $$c_2$$ are arbitrary constants.

This means that a  complete element can represent  exactly all rigid motions (constant displacements) and all constant strain states.

The patch test is also used to check the correct implementation of a finite element program.

Part 2
What is the consistency requirement?What is the stability condition?


 * Consistency:

The term consistency condition comes from the finite difference literature where it refers to the formal consistency of the difference equations with the correct differential equation.

In the paper, the consistency requirement implies that as the size of the elements tends to zero, the approximate equation $$\mathbf{K} \mathbf{a} = \mathbf{f}$$ will represent the governing differential equation exactly.

By stability, the authors mean that the solution of the finite element system of equations is unique. For linear systems, this means that the stiffness matrix has to be nonsingular.
 * Stability:

Part 3
What do the authors mean by convergence?

Let the finite element approximation be

u_h = \mathbf{N}^T\mathbf{a} $$ where $$\mathbf{N}$$ is the vector of shape functions, and $$\mathbf{a}$$ is the vector of nodal displacements.

By convergence, the authors mean that the approximate solution $$u_h \rightarrow u$$ as $$h \rightarrow 0$$, where $$u$$ is the exact solution, and $$h$$ is the element size.

Part 4
List the 2-D structural elements available in ANSYS and describe what shape functions and degrees of freedom each uses.

Listed below are some example of the 2-D structural elements used by ANSYS. The descriptions and the shape functions of each of these elements are from the  ANSYS Element Manual.
 * PLANE42: 2-D 4-Node Structural Solid
 * The PLANE42 element can be used either as a plane element (plane stress or plane strain) or as an axisymmetric element. The element is defined by four nodes (see Figure 4) with two translational degrees of freedom at each node. The shape functions for this element are

\begin{align} N_1 & = \cfrac{1}{4}(1-s)(1-t),& N_2 & = \cfrac{1}{4}(1+s)(1-t) \\ N_3 & = \cfrac{1}{4}(1+s)(1+t),&N_4 & = \cfrac{1}{4}(1-s)(1+t) \end{align} $$
 * where $$s \in [-1,1]$$ and $$t \in [-1,1]$$.


 * However, if incompatible elements are used (KEYOPT(2)=0), then the approximate solution is of the form

\begin{align} u & = u_1 N_1 + u_2 N_2 + u_3 N_3 + u_4 N_4 + u_a (1 - s^2) + u_b (1 - t^2) \\ v & = v_1 N_1 + v_2 N_2 + v_3 N_3 + v_4 N_4 + v_a (1 - s^2) + v_b (1 - t^2) \\ \end{align} $$
 * where nodes $$a$$ and $$b$$ are incompatible nodes. The forces at the incompatible nodes are zero. Hence the element stiffness matrix can be statically condensed and $$a$$ and $$b$$ do not have to correspond to any physical nodes.


 *  PLANE82: 2-D 8-Node Structural Solid


 * The PLANE82 element is defined by eight nodes (see Figure 5) with two translational degrees of freedom at each node. The element may be used as a plane elementor as an axisymmetric element. The shape functions of this element are

\begin{align} N_1 & = -\cfrac{1}{4}(1-s)(1-t)(1+s+t),& N_2 & = -\cfrac{1}{4}(1+s)(1-t)(1-s+t) \\ N_3 & = -\cfrac{1}{4}(1+s)(1+t)(1-s-t),& N_4 & = -\cfrac{1}{4}(1-s)(1+t)(1+s-t) \\ N_5 & = \cfrac{1}{2}(1-s^2)(1-t),& N_6 & = \cfrac{1}{2}(1+s)(1-t^2) \\ N_7 & = \cfrac{1}{2}(1-s^2)(1+t),& N_8 & = \cfrac{1}{2}(1-s)(1-t^2) \end{align} $$
 * where $$s \in [-1,1]$$ and $$t \in [-1,1]$$.


 *  PLANE182: 2-D 4-Node Structural Solid


 * PLANE182 is a four-noded quadrilateral element and has the same displacement shape functions as PLANE42. In addition, the element allows for pressure degrees of freedom if  KEYOPT(6) = 1 is chosen (mixed u-p formulation). The hydrostatic pressure has an interpolation function one order lower than the displacement interpolation. If  KEYOPT(1) = 0 is chosen, a $$\bar{B}$$ method with full integration is used. One pressure degree of freedom isused per element with this option. If  KEYOPT(1) = 1 is chosen, a $$\bar{B}$$ method with selective reduced integration (and hourglass control) is used in the computations. If  KEYOPT(1) = 2 is chosen, an  enhanced strain method isused and the pressure is interpolated with linear shape functions. There are three pressure degrees of freedom in the element.


 *  PLANE183 - 2-D 8-Node Structural Solid


 * PLANE183 is a two-dimensional eight-node element similar to PLANE82. This element can also include a mixed u-p formulationwhen  KEYOPT(6) is set to 1.

Part 5
Check whether each of these elements passes the patch test. Show details of your calculations.

 Patch test for PLANE182 elements

In this test we assume that the displacements on the patch have the form

\begin{align} u &= 0.001 + 0.002 x + 0.003 y \\ v &= -0.001 - 0.002 x - 0.003 y \end{align} $$ Then, assuming plane strain, the strains in the patch are

\begin{align} \varepsilon_{xx} &= \frac{\partial u}{\partial x} = 0.002 \\ \varepsilon_{yy} &= \frac{\partial v}{\partial y} = -0.003 \\ \gamma_{xy} &= \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} = 0.003 - 0.002 = 0.001 \\ \end{align} $$ The stresses are given by

\begin{align} \sigma_{xx} &= \cfrac{E}{(1+\nu)(1-2\nu)} \left[(1-\nu)\varepsilon_{xx} + \nu\varepsilon_{yy}\right] \\ \sigma_{yy} &= \cfrac{E}{(1+\nu)(1-2\nu)} \left[(1-\nu)\varepsilon_{yy} + \nu\varepsilon_{xx}\right] \\ \sigma_{zz} &= \cfrac{E}{(1+\nu)(1-2\nu)} \left[\nu(\varepsilon_{xx}+\varepsilon_{yy})\right] \\ \tau_{xy} &= \cfrac{E}{2(1+\nu)}\gamma_{xy} \end{align} $$ If $$E = 1000$$ and $$\nu = 0.3$$, the stresses are

\begin{align} \sigma_{xx} & = 9.615384615384616e-01, \qquad \sigma_{yy} = -2.884615384615385e+00, \\ \sigma_{zz} & = -5.769230769230769e-01, \qquad \tau_{xy} = 3.846153846153846e-01 ~. \end{align} $$ If $$E = 1000$$ and $$\nu = 0.499999$$, the stresses are

\begin{align} \sigma_{xx} & = -1.666651111145334e+05, \qquad \sigma_{yy} = -1.666684444500888e+05, \\ \sigma_{zz} & = -1.666664444487555e+05, \qquad \tau_{xy} = 3.333335555557037e-01 ~. \end{align} $$ Let us apply the patch test to the patch shown in Figure 6. The exact solution for the displacements at the nodes is

\begin{align} u_1 & = 0.001, \qquad u_2 = 0.005, \qquad u_3 = 0.014, \qquad u_4 = 0.007 \\ u_5 &= 0.003, \qquad u_6 = 0.0056, \qquad u_7 = 0.01, \qquad u_8 = 0.0064 \\ v_1 &= -0.001, \qquad v_2 = -0.005, \qquad v_3 = -0.014, \qquad v_4 = -0.007 \\ v_5 &= -0.003, \qquad v_6 = -0.0056, \qquad v_7 = -0.01, \qquad v_8 = -0.0064 \end{align} $$ We apply the exact displacements at nodes 1, 2, 3, and 4 and compute the remaining displacements and the stresses and strains using ANSYS. The ANSYS script we have used is shown below:

The results from ANSYS are compared with the exact values in the tables below.

These results clearly show that the element is passing the patch test. Similar results are obtained by changing the  keyopt values in the input file.

Patch test for PLANE42 The patch test (Test B) is applied to the model shown in Figure 5a in the paper. The same solution given in the handout is considered for this analysis.

u=0.002x\quad\mbox{and}\quad v=-0.0006y \qquad \text{(10)} $$ Listed below are the inputs for the 4-node element. To analyze the incompressible materials, $$\nu$$ is set to 0.4999 and the element type is changed to PLANE182. The  keyopt command should be changed to reflect the options available of PLANE182.

Patch test for PLANE82

For an 8-node element, the similar model is created. However, the displacements must also be applied to the midnodes. Only the inputs for PLANE183 are shown here. To run the analysis using PLANE82, the  ET line should be changed to  et,1,82, and the  keyopt line should be changed to reflect the options available for PLANE82. The input for the model using PLANE183 is shown below

 Results The displacements obtained from ANSYS are compared to (10). In order pass the patch test, the displacement at each node must match (10), and stresses have to be constant everywhere. Similarly for PLANE182 and PLANE183, the results shown in the table below use $$\nu=0.4999$$ to simulate the incompressible material. The verification of the results is also similar to the previous cases. Hence, we conclude that these four elements used by ANSYS, PLANE42, PLANE82, PLANE182, and PLANE183 passed the patch test.