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.

Solution

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.

Figure 1. Patch test.

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{aligned}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{aligned}}$

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{aligned}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{aligned}}$

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

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_{i}^{e}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_{i}^{e}=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.

• Stability:

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.

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{aligned}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{aligned}}$

where $s\in [-1,1]$ and $t\in [-1,1]$.

Figure 4. Node numbering for the Plane42 element.

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

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

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{aligned}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{aligned}}$

where $s\in [-1,1]$ and $t\in [-1,1]$.

Figure 5. Node numbering for the Plane82 element.

• 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{aligned}u&=0.001+0.002x+0.003y\\v&=-0.001-0.002x-0.003y\end{aligned}}$

Then, assuming plane strain, the strains in the patch are

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

The stresses are given by

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

If $E=1000$ and $\nu =0.3$, the stresses are

${\begin{aligned}\sigma _{xx}&=9.615384615384616e-01,\qquad \sigma _{yy}=-2.884615384615385e+00,\\\sigma _{zz}&=-5.769230769230769e-01,\qquad \tau _{xy}=3.846153846153846e-01~.\end{aligned}}$

If $E=1000$ and $\nu =0.499999$, the stresses are

${\begin{aligned}\sigma _{xx}&=-1.666651111145334e+05,\qquad \sigma _{yy}=-1.666684444500888e+05,\\\sigma _{zz}&=-1.666664444487555e+05,\qquad \tau _{xy}=3.333335555557037e-01~.\end{aligned}}$

Let us apply the patch test to the patch shown in Figure 6.

Figure 6. A patch of PLANE182 elements showing deformed shape.

The exact solution for the displacements at the nodes is

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

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:

/prep7
! 
! Element type
! 
et,1,182
! 
! Keyopt(1) - 0 = full integration, 1 = uniform reduced integration
! 2 = enhanced strain
! 
!keyopt,1,1,0
!keyopt,1,1,1
keyopt,1,1,2
! 
! Keyopt(3) - 0 = plane stress, 2 = plane strain
! 
keyopt,1,3,2
! 
! Keyopt(6) - 0 = disp (u) formulation, 1 = mixed (u-p) formulation
! 
!keyopt,1,6,0
keyopt,1,6,1
! 
! Material properties
! 
mp,ex,1,1000.0
mp,prxy,1,0.499999
!
! Create nodes
!
n,1,0,0
n,2,2,0
n,3,2,3
n,4,0,2
n,5,0.4,0.4
n,6,1.4,0.6
n,7,1.5,2
n,8,0.3,1.6
!
! Create elements
!
e, 1, 2, 6, 5
e, 2, 3, 7, 6
e, 7, 3, 4, 8
e, 1, 5, 8, 4
e, 5, 6, 7, 8
!
! Apply displacements to the boundary nodes
!
d,1,ux,0.001
d,1,uy,-0.001
d,2,ux,0.005
d,2,uy,-0.005
d,3,ux,0.014
d,3,uy,-0.014
d,4,ux,0.007
d,4,uy,-0.007
fini
!
! Solve (Ku = f)
!
/solu
solve
fini
!
! Write out results
!
/post1
/output,plane182Patch,out
prnsol, u, comp
presol, s, comp 
presol, epel, comp
/output,,
fini

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

Node	$u$ (FEM)	$u$ (Exact)	$v$ (FEM)	$v$ (Exact)
1	0.001	0.001	-0.001	-0.001
2	0.005	0.005	-0.005	-0.005
3	0.014	0.014	-0.014	-0.014
4	0.007	0.007	-0.007	-0.007
5	0.003	0.003	-0.003	-0.003
6	0.0056	0.0056	-0.0056	-0.0056
7	0.01	0.01	-0.01	-0.01
8	0.0064	0.0064	-0.0064	-0.0064

Element	$\varepsilon _{xx}$ (FEM)	$\varepsilon _{xx}$ (Exact)	$\varepsilon _{yy}$ (FEM)	$\varepsilon _{yy}$ (Exact)	$\gamma _{xy}$ (FEM)	$\gamma _{xy}$ (Exact)
1	0.002	0.002	-0.003	-0.003	0.001	0.001
2	0.002	0.002	-0.003	-0.003	0.001	0.001
3	0.002	0.002	-0.003	-0.003	0.001	0.001
4	0.002	0.002	-0.003	-0.003	0.001	0.001
5	0.002	0.002	-0.003	-0.003	0.001	0.001

Element	$\sigma _{xx}$ (FEM)	$\sigma _{xx}$ (Exact) !	$\sigma _{yy}$ (FEM)	$\sigma _{yy}$ (Exact)	$\sigma _{zz}$ (FEM)	$\sigma _{zz}$ (Exact)
	($\times 10^{5}$)	($\times 10^{5}$)	($\times 10^{5}$)	($\times 10^{5}$)	($\times 10^{5}$)	($\times 10^{5}$)
1	-1.6667	-1.666651	-1.6667	-1.666684	-1.6667	-1.666664
2	-1.6667	-1.666651	-1.6667	-1.666684	-1.6667	-1.666664
3	-1.6667	-1.666651	-1.6667	-1.666684	-1.6667	-1.666664
4	-1.6667	-1.666651	-1.6667	-1.666684	-1.6667	-1.666664
5	-1.6667	-1.666651	-1.6667	-1.666684	-1.6667	-1.666664
Element	$\tau _{xy}$ (FEM)	$\tau _{xy}$ (Exact)
1	0.33333	0.33333355
2	0.33333	0.33333355
3	0.33333	0.33333355
4	0.33333	0.33333355
5	0.33333	0.33333355

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.

2-D 4-node element with node numbers.

/prep7
!change this line to et,1,182 for PLANE182
et,1,42
!change this line to keyopt,1,2,0 to include the extra displacement
keyopt,1,2,1
mp,ex,1,1e3
mp,prxy,1,0.3
k,1,0,0
k,2,2,0
k,3,2,3
k,4,0,2
k,5,0.4,0.4
k,6,1.4,0.6
k,7,1.5,2
k,8,0.3,1.6
!create lines and set their number of division to 1 for each line
l,1,2,1
l,2,3,1
l,3,4,1
l,4,1,1
l,5,6,1
l,6,7,1
l,7,8,1
l,8,5,1
l,1,5,1
l,2,6,1
l,3,7,1
l,4,8,1
a,1,2,6,5
a,3,4,8,7
a,1,5,8,4
a,5,6,7,8
a,2,3,7,6
!AMAP is used instead of AMESH in this case to ensure the quad elements
AMAP,1,1,2,6,5
AMAP,2,3,4,8,7
AMAP,3,1,5,8,4
AMAP,4,5,6,7,8
AMAP,5,6,2,3,7
!apply displacement conditions to the surrounding nodes
dk,1,all,0
dk,2,ux,0.004
dk,2,uy,0
dk,3,ux,0.004
dk,3,uy,-0.0018
dk,4,ux,0
dk,4,uy,-0.0012
fini

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

2-D 8-node element with node numbers.

/prep7
!change this line to et,1,82 for PLANE82
et,1,183
!sets the element for using the mixed formulation (u-p).
keyopt,1,6,1
mp,ex,1,1e3
mp,prxy,1,0.4999
k,1,0,0
k,2,2,0
k,3,2,3
k,4,0,2
k,5,0.4,0.4
k,6,1.4,0.6
k,7,1.5,2
k,8,0.3,1.6
!create lines and set their number of division to 1 for each line
l,1,2,1
l,2,3,1
l,3,4,1
l,4,1,1
l,5,6,1
l,6,7,1
l,7,8,1
l,8,5,1
l,1,5,1
l,2,6,1
l,3,7,1
l,4,8,1
a,1,2,6,5
a,3,4,8,7
a,1,5,8,4
a,5,6,7,8
amesh,all
!this area is meshed separately to force a quad mesh
a,2,3,7,6
AMAP,5,6,2,3,7
!apply displacement conditions to the surrounding nodes
d,1,all,0
d,2,ux,0.004
d,2,uy,0
d,9,ux,0.004
d,9,uy,-0.0018
d,10,ux,0
d,10,uy,-0.0012
d,3,ux,0.002
d,3,uy,0
d,20,ux,0.004
d,20,uy,-0.0006*1.5
d,11,ux,0.002
d,11,uy,-0.0006*2.5
d,18,ux,0
d,18,uy,-0.0006
fini

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.

Node	Coordinates		Displacements		Stresses	Test B
$i$	$x_{i}$	$y_{i}$	$u_{i}$	$v_{i}$	$\sigma _{x}$
1	0.0	0.0	0.0000	0.0000	2.0000	pass
2	2.0	0.0	0.40000E-02	0.0000	2.0000	pass
3	1.4	0.6	0.28000E-02	-0.36000E-03	2.0000	pass
4	0.4	0.4	0.80000E-03	-0.24000E-03	2.0000	pass
5	2.0	3.0	0.40000E-02	-0.18000E-02	2.0000	pass
6	0.0	2.0	0.0000	-0.12000E-02	2.0000	pass
7	0.3	1.6	0.60000E-03	-0.96000E-03	2.0000	pass
8	1.5	2.0	0.30000E-02	-0.12000E-02	2.0000	pass

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.

Node	Coordinates		Displacements		Stresses	Test B
$i$	$x_{i}$	$y_{i}$	$u_{i}$	$v_{i}$	$\sigma _{x}$
1	0.0	0.0	0.0000	0.0000	2.2664	pass
2	2.0	0.0	0.40000E-02	0.0000	2.2664	pass
3	1.0	0.0	0.20000E-02	0.0000	2.2664	pass
4	1.4	0.6	0.28000E-02	-0.36000E-03	2.2664	pass
5	1.7	0.3	0.34000E-02	-0.18000E-03	2.2664	pass
6	0.4	0.4	0.80000E-03	-0.24000E-03	2.2664	pass
7	0.9	0.5	0.18000E-02	-0.30000E-03	2.2664	pass
8	0.2	0.2	0.40000E-03	-0.12000E-03	2.2664	pass
9	2.0	3.0	0.40000E-02	-0.18000E-02	2.2664	pass
10	0.0	2.0	0.0000	-0.12000E-02	2.2664	pass
11	1.0	2.5	0.20000E-02	-0.15000E-02	2.2664	pass
12	0.3	1.6	0.60000E-03	-0.96000E-03	2.2664	pass
13	0.15	1.8	0.30000E-03	-0.10800E-02	2.2664	pass
14	1.5	2.0	0.30000E-02	-0.12000E-02	2.2664	pass
15	0.9	1.8	0.18000E-02	-0.10800E-02	2.2664	pass
16	1.75	2.5	0.35000E-02	-0.15000E-02	2.2664	pass
17	0.35	1.0	0.70000E-03	-0.60000E-03	2.2664	pass
18	0.0	1.0	0.0000	-0.60000E-03	2.2664	pass
19	1.45	1.3	0.29000E-02	-0.78000E-03	2.2664	pass
20	2.0	1.5	0.40000E-02	-0.90000E-03	2.2664	pass

Hence, we conclude that these four elements used by ANSYS, PLANE42, PLANE82, PLANE182, and PLANE183 passed the patch test.

Nonlinear finite elements/Homework 3