COMSOL 6.4 - Stress Analysis of an Elliptic Membrane

In this benchmark, the static stress analysis described in the NAFEMS Test LE1, “Elliptic Membrane”, found on page 5 in Ref. 1 is performed. It is an analysis of a linear elastic plane stress model.

The computed stress level is compared with the values given in the benchmark report.

In addition to the original benchmark, a mesh convergence study is performed.

GEOMETRY

The geometry is an ellipse with an elliptical hole in it. The outer and inner edges are defined by the equations

The thickness (which actually does not influence the analysis) is 0.1 m.

Due to symmetry in load and in geometry, the analysis only includes a quarter of the geometry as shown in Figure 1.

Figure 1: The geometry and load. Only the quarter that is analyzed is shown.

Material

Isotropic with E = 210 GPa and ν = 0.3.

Load

An evenly distributed load of 10 MPa acts along the outward normal of the outer boundary.

Constraints

Symmetry conditions are used along the cuts at X = 0 and Y = 0.

Model Setup

The Solid Mechanics interface with the plane stress assumption is used.

Four meshes are exactly specified in Ref. 1. The “coarse” mesh has 6 quadrilateral or 12 triangular elements. The “fine” mesh has 24 quadrilateral or 48 triangular elements. The triangular elements are created by splitting the quadrilateral elements along a diagonal. The specified meshes are shown in Figure 2 and Figure 3. For the quadrilateral meshes, a reduced integration scheme with hourglass stabilization is also used.

Note that the algorithm used to subdivide the quadrilaterals into triangles tries to place the diagonal so that the two new triangles have the best possible shape. As a result, the triangular mesh used is not identical to the one in the original benchmark, even though the number of elements and even node locations are identical. For very coarse meshes, this can change the stress results appreciably.

For the mesh convergence study, these meshes are uniformly refined using a parameter div. The number of elements along the elliptical boundaries is 3*div and the number of elements along the symmetry cuts is 2*div.

The number of displacement degrees of freedom varies from 48 (div = 1 and quadrilaterals with linear shape order) to 935810 (div = 64 and triangles with cubic shape order).

Figure 2: The meshes as specified in Ref. 1. Left column: “coarse” (div=1). Right column: “fine” (div=2).

Figure 3: A quadrilateral mesh with div=8.

Due to the specification of the benchmark, the modeling differs somewhat from what you would use in practice:

•

The interior boundaries in the model are created for matching the specification of the mesh in the NAFEMS benchmark as close as possible. If you were to solve the problem without these constraints, the modeling would be significantly simplified. Only two ellipses would be needed in the Geometry sequence.

•

The knowledge about where a stress concentration is expected suggests that you should use a mesh such that more elements are present in the region around point D to get optimal accuracy, see Figure 1.

•

Using the possibility to generate a free triangular mesh instead of one where quadrilateral elements are split along the diagonals would also give a mesh with better element quality.

Results and Discussion

The distribution of the direct stress in the Y direction is shown in Figure 4. As can be seen the result has steep gradients toward the point with maximum values.

Figure 4: The distribution of the σy stress component using div=4 and second order quadrilateral elements.

The normal stress σy at the elliptic hole is evaluated at the point D located at X = 2, Y = 0 (see Figure 1). The target value according to Ref. 1 is 92.7 MPa. The value is based on an analytical result. The COMSOL Multiphysics results for the “coarse” and “fine” meshes are given in Table 1.

Table 1: Computed results for the meshes specified in the benchmark.
Study number	Element TYPE	Discretization	Mesh	COMPUTED VALUE	RELATIVE ERROR
1	Quadrilateral	Linear	Coarse	71.0	-23.3%
1	Quadrilateral	Linear	Fine	85.5	-7.7%
2	Quadrilateral	Quadratic	Coarse	89.0	-4.0%
2	Quadrilateral	Quadratic	Fine	92.7	0.0%
3	Quadrilateral	Cubic	Coarse	93.7	+1.1%
3	Quadrilateral	Cubic	Fine	92.5	-0.2%
4	Triangle	Linear	Coarse	33.9	-63.4%
4	Triangle	Linear	Fine	54.1	-41.6%
5	Triangle	Quadratic	Coarse	74.2	-20.0%
5	Triangle	Quadratic	Fine	87.7	-5.4%
6	Triangle	Cubic	Coarse	86.0	-7.2%
6	Triangle	Cubic	Fine	91.8	-0.9%
9	Quadrilateral	Linear, reduced	Coarse	47.7	-48.5%
10	Quadrilateral	Linear, reduced	Fine	62.3	-32.7%
11	Quadrilateral	Quadratic, reduced	Coarse	83.1	-10.3%
12	Quadrilateral	Quadratic, reduced	Fine	91.2	-1.54%
13	Quadrilateral	Cubic, reduced	Coarse	87.6	-5.5%
14	Quadrilateral	Cubic, reduced	Fine	91.0	-1.7%

As can be expected, the coarse mesh is not able to capture the stress concentration unless elements with high order are used. Generally the quadrilaterals perform better than the corresponding triangles.

The mesh that is denoted as “fine” is more similar to what you would use in an analysis of a larger structure in a case where you are not specifically interested in a high resolution of the stress concentration. Still, with quadratic shape order elements, the accuracy is good enough for most engineering purposes. With the current mesh, the triangular elements will have a small angle at the stress evaluation point, hence the less accurate result.

Using elements with linear shape functions for structural analysis is commonly avoided in the finite element community.

The results of the mesh convergence study are shown in Figure 5. The element size h is defined as 0.417[m]/div, which is the length of an edge in the element where the stress is measured.

The target value in Ref. 1, 92.7 MPa, is given with only three digits. This is not accurate enough for the convergence study here. Instead, the error is measured relative to the value 92.65817 MPa, toward which σy converges.

The convergence behavior is as expected since it is faster for elements with a higher shape function order. It can also be seen that quadrilaterals are somewhat more accurate than triangles for quadratic and cubic elements. The quadrilateral elements with a reduced ingratiation rule perform similar to the triangles for first order shape functions. For higher order shape functions, they are almost on par with their fully integrated counterparts.

The other two in-plane stress components σx and τxy should both be zero at point D since the boundary is free. In Figure 6 and Figure 7 similar convergence graphs are shown for these stress components.

Figure 5: Error with respect to the stress target value as a function of the element size h.

Figure 6: Error in the stress σx. The values are normalized with the target for σy.

Figure 7: Error in the stress τxy. The values are normalized with the target for σy.

Since elements with different shape function orders are used, a comparison based only on element size may not be fair when efficiency is considered. The number of degrees of freedom in the model varies a lot for the same element size, and so does the solution time. In Figure 8, the error is shown as a function of the number of degrees of freedom. Also when compared this way, the elements with cubic shape functions have the best performance. This is usually true as long as the solutions are smooth, but it may not be true, for example, when solving nonlinear problems.

However, for coarser meshes, which are more representative or real life mesh resolution, the second order quadrilateral elements are quite efficient when compared versus the number of degrees of freedom.

Note that somewhat fewer degrees-of-freedom are used for when a reduced integration rule is used. This follows from the lower order of the auxiliary shape functions used for the out-of-plane displacement derivative.

Figure 8: Error with respect to the stress target value as a function of the number of degrees of freedom.

Reference

1. G.A.O. Davies, R.T. Fenner, and R.W. Lewis, Background to Benchmarks, NAFEMS, Glasgow, 1993.

Structural_Mechanics_Module/Verification_Examples/elliptic_membrane

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
div	1	1	Mesh refinement factor
sy_ref	92.65817[MPa]	9.2658E7 Pa	Target stress

Geometry 1

Ellipse 1 (e1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 90.

In the text field, type 3.25.

In the text field, type 2.75.

Create an extra mesh control ellipse.

Ellipse 2 (e2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.417.

In the text field, type 1.583.

In the text field, type 90.

Ellipse 3 (e3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the objects and only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the object only.

Click

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Locate the section. In the text field, type 1.783 and to 2.3.

Locate the section. In the text field, type 1.165 and to 0.812.

Line Segment 2 (ls2)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Locate the section. In the text field, type 2.833 and to 1.348.

Locate the section. In the text field, type 1.783 and to 0.453.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	210[GPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.3	1	Young’s modulus and Poisson’s ratio
Density	rho	0	kg/m³	Basic