A technique to develop mesh-distortion immune finite elements

Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A technique to develop mesh-distortion immune finite elements S. Rajendran * School of Mechanical and Aerospace Engineering, Nanyang Technological University, Nanyang Avenue, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 22 September 2007 Received in revised form 13 November 2009 Accepted 20 November 2009 Available online 3 December 2009 Keywords: Finite element Mesh distortion Element distortion 8-Node quadrilateral Quadratic quadrilateral Geometric distortion

a b s t r a c t Three sufficient conditions are proposed for an iso-parametric finite element formulation to be immune to mesh-distortion effects, even for extremely distorted element geometries for which the Jacobian geometric mapping goes negative! Illustrative element formulations are presented to test the sufficiency of the proposed conditions. Numerical results confirm that the illustrative elements are indeed capable of giving mesh-distortion immune solution for extreme mesh distortions. This has been possible partly because the stiffness integral is free of the Jacobian term in the new formulation. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Iso-parametric elements have been time-tested and generally give reliable performance for regular meshes. However, in the presence of severe mesh distortions, the performance of these elements may deteriorate drastically. Mesh distortions are not completely avoidable in practice. Simpler geometries like squares and rectangles can easily be meshed with regular (undistorted) elements thus avoiding distorted meshes altogether. However, meshes for complex geometrical shapes cannot be completely free from mesh distortions. Mesh distortions also naturally arise in the simulation of large deformation process like forging and extrusion. Even if the initial mesh is carefully chosen to have very small or little geometric distortions, severe distortions may eventually set in as the material deforms considerably. Mesh distortions involving negative Jacobian is often possible in such situations. The classical isoparametric elements are not designed to handle such situations. To continue simulation, one may have to resort to remeshing. This, however, involves additional cost as well as attendant difficulties such as how to map the internal forces to the new mesh. In this sense, formulations free from the effects of Jacobian are of much interest, but are currently scarce in the literature. Although, the problem of mesh-distortion sensitivity has been known to researchers for three decades, a permanent remedy to the problem is yet to appear. The problem was first recognised in the seventies [1–7]. Lee and Bathe [8] attributed the mesh-distortion sensitivity to the loss of polynomial completeness of shape * Tel.: +65 6790 5865; fax: +65 6791 1859. E-mail address: [email protected] 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.11.017

functions under distorted element geometries. When the mesh is distorted, the shape functions lose the ability to reproduce exactly the various Cartesian monomial terms present in the displacement interpolation. This observation has motivated researchers to develop newer elements that have quadratic polynomial completeness in physical (Cartesian) space [9–13]. The elements based on quadrilateral area coordinate shape functions developed by Soh et al. [9] and Chen et al. [10,12] as well as the incompatible element proposed by Huang and Li [11] attempt to ensure quadratic completeness in physical space. As a result, these elements emerge to be highly insensitive to mesh distortions. The elements give excellent results for several benchmark problems. Nevertheless, for element shapes involving curved geometries, the area coordinate elements generally do not yield distortion-insensitive solution. Rajendran and his co-researchers [13–16] developed a family of quadratic elements (viz., US-QUAD8 [13,14], US-HEXA20 [15] and TRIA6 [16]) called unsymmetric elements which use two different and BR matrices appearing in types of shape functions for R the BL e the stiffness integral, K ¼ V e BTL DBR dV . Although these elements give rise to unsymmetric stiffness matrices, the performance of these elements is highly tolerant to mesh distortions. Unsymmetric elements are capable of reproducing exact solution for quadratic displacement problems even under severe mesh distortions involving quadratic element geometries. The elements pass not only the constant strain patch-test but also the linear strain patch-test in the strong form [17]. Recently, Prathap et al. [18] studied the unsymmetric formulation in the light of best-fit paradigm. They observed that this formulation gives excellent results only because the use of parametric

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

shape functions for the test functions helps to satisfy the continuity conditions across element edges while the use of metric shape functions for the trial functions ensures that the actual variation of stress in the metric space is retrieved correctly in a best-fit manner even when elements are distorted. In another paper [19] on unsymmetric formulation for a 1D bar element, Prathap et al. made an interesting observation that the Jacobian no longer appears in the stiffness integral. The absence of Jacobian in the stiffness integral is a more than welcome feature because it may help to prevent the mesh-distortion effects from entering the element stiffness via iso-parametric mapping of element geometry. This observation forms a motivation for this paper. The present paper focuses on the root cause of mesh-distortion sensitivity and proposes three sufficient conditions, which when satisfied, ensure mesh-distortion immunity of the element. There may exist other sets of conditions leading to mesh-distortion immune elements. Nevertheless, the sufficient conditions proposed in this paper are perhaps a first attempt towards characterising the requirements for a mesh-distortion immune formulation. The concepts and developments presented in this paper centre around polynomial interpolation functions. Nevertheless, the concepts can be extended equally well to non-polynomial interpolation fields such as the trigonometric displacement field around the crack-tip in fracture mechanics problems. For the classical 8-node and 9-node quadratic quadrilateral elements, the interpolation errors and convergence aspects have been discussed by Arnold et al. [20,21] for affine and bilinear element geometries. For a recent 8-node element derived from 9-node element [22], interpolation error estimates have been obtained for bilinear geometries [23]. Such an analysis for the element discussed in this paper would be of much interest, but seems intractable as the mesh distortions handled in this paper are extremely severe involving even negative Jacobian. For the purpose of this paper, the term mesh-distortion tolerance will be interpreted as the ability of an element to produce generally accurate results even under distorted meshes. However, the term mesh-distortion immunity will be interpreted as the ability of an element to produce results that are completely free from mesh-distortion effects. More formal definitions of mesh-distortion immunity will be attempted in Section 2. Sufficient conditions will be presented in Section 3 and illustrative element formulations in Section 4. Numerical verification of sufficiency of the conditions will be provided in Section 5 and major conclusions from the present work will be summarised in Section 6. 2. Definitions of mesh-distortion immunity Mesh-distortion immunity of an element, in simple terms, is its ability to give consistent performance with no sensitivity to mesh distortions. A higher order element usually exhibits mesh-distortion immunity for lower order displacement problems, but not in general for the higher. In this paper, however, we shall brand an element mesh-distortion immune only when it gives mesh-distortion independent performance, not only for the lower order displacement problems but also for the higher order ones (as high as the order of the element itself), and hence the following definition: Definition 0 (Mesh-distortion immunity of an nth order element). An nth order element is said to have mesh-distortion immunity, if it can reproduce, for all possible distorted meshes, the exact solution to an nth order displacement problem. For iso-parametric elements, this definition is rather ambiguous with respect to the term ‘‘order”. Traditionally, the order of an element refers to the order of parametric polynomial (i.e., a polynomial

1045

written in local/parametric coordinates) displacement field assumed. In contrast, the exact solution for many practical problems is often a Cartesian polynomial. It is thus difficult to relate the order of the element to the order of exact solution. Furthermore, it is also unclear whether the phrase ‘all possible distorted meshes’ includes element geometries with non-positive Jacobian or it just refers to only geometries with positive Jacobian. Hence, for iso-parametric elements, a more specific definition may be proposed as follows: Definition 1 (Jacobian-limited mesh-distortion immunity). An isoparametric element is said to have Jacobian-limited mesh-distortion immunity, if it can exactly reproduce, for all admissible element geometries, a Cartesian polynomial displacement field similar in form to the parametric polynomial displacement interpolation of the element. (The phrase ‘admissible element geometries’ refers to element geometries limited to the ones with positive Jacobian at Gauss points.) Consider, for illustration, the classical iso-parametric quadratic quadrilateral (QUAD8) element with parametric polynomial displacement interpolation of the form

u ¼ a0 þ a1 n þ a2 g þ a3 n2 þ a4 ng þ a5 g2 þ a6 n2 g þ a7 ng2 :

ð1Þ

For this element to exhibit Jacobian-limited mesh-distortion immunity, the element should be able to reproduce, for all admissible element geometries, a Cartesian polynomial displacement field of the form

u ¼ b0 þ b1 x þ b2 y þ b3 x2 þ b4 xy þ b5 y2 þ b6 x2 y þ b7 xy2

ð2Þ

which is similar in form to Eq. (1) with the coordinates n and g replaced by x and y, respectively. One can also think of degenerate versions of Definition 0. In reality, a given element may not have mesh-distortion immunity with respect all the monomial terms, but only with respect to a limited number of them such as linear terms {1, x, y} or up to bilinear terms {1, x, y, xy} or up to quadratic terms {1, x, y, x2, xy, y2}. One may then possibly define partial mesh-distortion immunity with respect to linear field, bilinear field and quadratic field, and so on. QUAD8 element does not indeed possess mesh-distortion immunity with respect to all the terms in Eq. (2). The number of monomial terms that the element can actually reproduce depends on the element geometry, i.e., whether the element is used in linear or bilinear or quadratic element geometry (see e.g., Lee and Bathe [8]). In such cases, one may define partial mesh-distortion immunity with respect to linear/bilinear/quadratic field under linear/bilinear/ quadratic element geometry. In this paper, we shall deal with elements that have mesh-distortion immunity with respect all the terms in Eq. (2), not only with respect to admissible element geometries, but also element geometries with non-positive Jacobian. Hence, the following definition: Definition 2 (Jacobian-unlimited mesh-distortion immunity). An iso-parametric element is said to have Jacobian-unlimited meshdistortion immunity, if it can exactly reproduce, for all possible element geometries, a Cartesian polynomial displacement field similar in form to the parametric polynomial displacement interpolation of the element. (The phrase ‘all possible element geometries’ refers to all element geometries possible in an isoparametric mapping, without limiting to those with positive Jacobian at Gauss points. The cases of coincident nodes are, however, excluded.) Jacobian-unlimited mesh-distortion immunity is the ultimate property one can dream of an element. Elements exhibiting Jacobian-unlimited mesh-distortion immunity are scarce in the literature. The element formulation proposed in this paper, however,

1046

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

exhibits this property. This has been possible partly because, as we shall see later, the stiffness integral of the new formulation is free from the Jacobian term.

3. Sufficient conditions for Jacobian-unlimited mesh-distortion immunity The expression for the element stiffness matrix may be written as:

KðeÞ ¼

Z Ve

e

BTL DBR dV ¼

Z

1

1

Z

1

1

BTL DBR jJjt dn dg;

literature satisfy condition C1 but not C2 in general. On the other hand, the unsymmetric (Petrov–Galerkin) elements [13–16] use two sets of shape functions to obtain BL and BR matrices, in which case, the two sets of shape functions can be so chosen as to satisfy the conditions C1 and C2 individually. This makes the unsymmetric formulation particularly attractive for the purposes of developing mesh-distortion immune elements. Formulations satisfying only conditions C1 and C2, but not C3 possess Jacobian-limited mesh-distortion immunity. Those satisfying all the three conditions exhibit Jacobian-unlimited mesh-distortion immunity, which forms the focus of this paper.

ð3Þ

where BL and BR are the strain–displacement matrices corresponding to test and trial functions, respectively, D is the material constitutive matrix, J is the Jacobian matrix, t is the thickness of the element, and Ve is the volume of the element. Two different strain–displacement matrices, BL and BR, are shown in Eq. (3) so that Eq. (3) can cater to either unsymmetric (Petrov–Galerkin) formulation or symmetric (Galerkin) formulation. For unsymmetric formulation, BL and BR are obtained using two different sets of shape functions. For symmetric formulation, the same set is used, i.e., BR  BL. An important observation from Eq. (3) is that BL, BR and J are the only three quantities which can breed mesh-distortion effects; i.e., these are the only source terms that can inject mesh-distortion sensitivity into a formulation. This observation immediately suggests that the stiffness matrix K can be rendered mesh-distortion immune, provided BL, BR and J can be rendered free of mesh-distortion effects. In order to render BL and BR mesh-distortion immune, the shape functions (used to construct the matrices BL and BR, respectively) should be carefully chosen so that they are immune to mesh-distortion effects. The unsymmetric formulation [13] gives some guidelines for the same.

3.2. MDI-IEC property A displacement interpolation is said to have mesh-distortion independent inter-element compatibility property if continuity of displacements at the element interfaces is maintained for all possible mesh distortions. The displacement interpolation of iso-parametric elements is well-known to possess this MDI-IEC property. This is because their shape functions take a value of unity at that their respective nodal positions and zero at all edges not containing the respective nodes. Some shape functions have mere Kronecker-delta property, i.e., they take a value of unity at their respective nodal positions and zero at other nodal positions, but not necessarily all along the edges not containing the respective nodes. Such shape functions are unsuitable to ensure MDI-IEC property. The metric shape functions [13– 16] (alias Cartesian shape functions) form an example of this category. In summary, to ensure MDI-IEC property, parametric shape functions form a good choice and hence will be used in this paper; Metric shape functions are a bad choice for this purpose. 3.3. MDI-C property

3.1. Sufficient conditions Based on hind sight [13–19], the sufficient conditions for the stiffness matrix to be mesh-distortion immune even for meshes with non-positive Jacobian at Gauss points are stated as follows: C1: The shape functions used in BL matrix must have mesh-distortion independent inter-element compatibility (MDI-IEC) property, i.e. they must ensure inter-element compatibility whatever be the mesh distortion. The MDI-IEC property will be elaborated in Section 3.2. C2: The shape functions used in BR matrix must satisfy mesh-distortion independent completeness property (MDI-C) up to all the terms in the Cartesian displacement polynomial (Eq. (2)) that we want the element to reproduce exactly even under mesh distortion. The MDI-C property will be elaborated in Section 3.3. C3: The Jacobian should not occur, explicitly or implicitly, in the denominator of the stiffness integrand. This is to avoid a rational integrand function and thereby ensure smoothness of variation of integrand function over the element domain under severe mesh distortions. The implication of this condition is rather subtle, and will be elaborated in Sections 3.4 and 5. In the classical iso-parametric elements, BR  BL = B. Hence, for such elements to be mesh-distortion immune, the shape functions used in B-matrix must satisfy both conditions, C1 and C2, simultaneously. This is rather too stringent a requirement to satisfy because the iso-parametric shape functions available in the

A displacement interpolation is said to have mesh-distortion independent completeness if the interpolation is capable of reproducing every term of the Cartesian displacement polynomial that we want to reproduce (Eq. (2)) for all possible element geometries. For a 8-node quadrilateral element for instance, reproducibility conditions of the individual terms are given by [24] 8 X

Ni ¼ 1;

8 X

i¼1 8 X

i¼1

Ni x2i ¼ x2 ;

i¼1

8 X

Ni yi ¼ y;

ð4aÞ

i¼1 8 X

i¼1 8 X

Ni xi ¼ x;

Ni xi yi ¼ xy;

i¼1

Ni x2i yi ¼ x2 y;

8 X

Ni y2i ¼ y2 ;

ð4bÞ

i¼1 8 X

N i xi y2i ¼ xy2 :

ð4cÞ

i¼1

If the shape functions Ni satisfy the above conditions, then all the terms {1, x, y, x2, xy, y2, x2y, xy2} are reproducible for all possible mesh distortions. For iso-parametric shape functions, there is no guarantee that all the above conditions are satisfied for all distorted element geometries. Hence, not all the terms {1, x, y, x2, xy, y2, x2y, xy2} are reproducible for all possible mesh distortions. The terms that are reproducible depends on the element geometry (viz., linear or bilinear or quadratic) and have been listed by Lee and Bathe [8]. Metric shape functions, however, satisfy all the eight conditions above irrespective of the type and extent of mesh distortion, and

1047

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

hence inherently satisfy the MDI-C property. This is because these conditions themselves are solved to get the metric shape functions, as will be seen in Eq. (11) later. In summary, to satisfy MDI-C property, metric shape functions form a good choice and hence will be used in this paper; Parametric shape functions are to be avoided for the purpose.

4.1. Eight-node quadratic quadrilateral element The three conditions C1, C2 and C3 are satisfied as follows: 4.1.1. Satisfaction of condition C1 In order to satisfy condition C1, BL matrix is constructed using the well-known iso-parametric shape functions:

3.4. Sufficiency of the conditions The mesh distortions considered in this paper are extremely severe involving quadratic element geometries with even negative Jacobian values for some meshes. For such severe mesh distortions, a rigorous mathematical proof of the sufficiency of conditions based on approximation theory seems formidable. Nevertheless, the sufficiency of the conditions can be shown to directly emerge from the well-known weak form – strong form equivalence (e.g. see [25,26]): As long as the test and trial functions are sufficiently smooth (satisfying the required compatibility and completeness requirements), a solution to the weak form is also a solution to the strong form, i.e., the exact solution. This is guaranteed by the equivalence theorem [25,26]. Thus, if the compatibility and completeness conditions are satisfied for all possible distorted element geometries, then a solution to the weak form will also be a solution to the strong form for all possible distorted meshes. Put in other words, we are assured of exact solution for all possible distorted element meshes provided the first two sufficiency conditions are satisfied. (In the classical proof of equivalence theorem, the trial space is assumed to be infinite dimensional. In the finite element method, however, the trial space is finite dimensional. Nevertheless, the theorem still holds for a finite dimensional trial space provided the trial space is sufficiently complete, i.e., if the solution vector lies in the trial space. In classical finite element analyses, sufficient completeness is generally achieved usually in the limit of successive mesh refinements. In the present paper, this is ensured by including in the trial basis all the monomial terms of the intended exact solution. For instance, if we intend to capture an exact solution of the form u = b0 + b1x + b2y + b3x2 + b4xy + b5y2 + b6 x2y + b7xy2, then a sufficiently complete trial basis is {1, x, y, x2, xy, y2, x2y, xy2}.) The third condition (and hence the necessity to prove its sufficiency) arises as we use Gaussian quadrature for the evaluation of the stiffness matrix. If it is ever possible to evaluate the integration analytically, this condition can be dispensed with. Gaussian quadrature schemes are guaranteed to give exact numerical integration (for a sufficiently high order of integration) for polynomial functions. For highly distorted iso-parametric elements, the integrand may become a singular function at the Gaussian sampling points. This is due to the appearance of Jacobian in the denominator of the integrand as will be discussed in Section 5. Singular functions are difficult to evaluate exactly by Gaussian quadrature. Inaccurate numerical integration of the stiffness matrix leads to erroneous displacement results. Thus, the third condition has been included to ensure that the stiffness integrand remains a polynomial function so that exact numerical integration is always possible even for extremely severe mesh distortions.

4. Illustrative formulations The rest of the paper is dedicated to illustrate numerically that the three conditions derived in Section 3.1. are sufficient for the design of elements that have mesh-distortion immunity even when Jacobian goes negative. For illustration, we shall consider 8-noded as well as 9-noded quadratic quadrilateral plane stress element.

Ni ¼

1 ð1 þ nni Þð1 þ ggi Þðnni þ ggi  1Þ; 4

ð5Þ

Nj ¼

1 ð1  n2 Þð1 þ ggj Þ; 2

ð6Þ

Nk ¼

1 ð1 þ nnk Þð1  g2 Þ; 2

ð7Þ

where i, j and k refer to corner nodes, mid-side nodes along n = 0 edges, and mid-side nodes along g = 0 edges, respectively. The shape function matrix is defined as

 N¼

N1

0

N2

0

N3

0

   N8

0

0

N1

0

N2

0

N3



0

N8

 ð8Þ

:

BL matrix is given by

2

N1;x

6 BL ¼ 4 0 N1;y

0

N2;x

0

N3;x

0



N8;x

N1;y

0

N 2;y

0

N3;y



0

N1;x

N2;y

N2;x

N3;Y

N3;x

   N8;Y

0

3

7 N8;y 5

ð9Þ

N8;x

Ni,x and Ni,y (i = 1, 2, 3, . . ., 8) are the Cartesian derivatives of the parametric shape functions given by the coordinate transformation:



N i;x Ni;y



¼ ½J1



N i;n

 

Ni;g

Adj½J jJj



N i;n N i;g

 ;

ð10Þ

where Ni,n and Ni,g (i = 1, 2, 3, . . ., 8) are the derivatives of shape functions shown in Eqs. (5)–(7), and jJj is the determinant of Jacobian matrix of coordinate transformation, which will be replaced by unity in Eq. (10) for reasons described in Section 4.1.3. 4.1.2. Satisfaction of condition C2 In order to satisfy condition C2, BR matrix is constructed using the metric (Cartesian) shape functions, Mi, (i = 1, 2, 3, . . . , 8), which are obtained by solving the equations [24]:

2

1

1

1

1

1

1

1

6 x1 x2 x3 x4 x5 x6 x7 6 6 6 y1 y2 y3 y4 y5 y6 y7 6 6 x2 x22 x23 x24 x25 x26 x27 6 1 6 6 x1 y1 x2 y2 x3 y3 x4 y4 x5 y5 x6 y6 x7 y7 6 6 y2 y22 y23 y24 y25 y26 y27 6 1 6 2 4 x1 y1 x22 y2 x23 y3 x24 y4 x25 y5 x26 y6 x27 y7 x1 y21 x2 y22 x3 y23 x4 y24 x5 y25 x6 y26 x7 y27

38 9 8 9 M1 > > 1 > > > > > > > > > > > > > > x > > x8 7 M2 > 7> > > > > > > > > > > > 7> > > > > > y > y8 7> M3 > > > > > > > > 7> > > > > = < 2 7< 2 = x x8 7 M 4 ¼ : 7 x8 y8 7> xy > > > > > M5 > > > > > 7> > > > > > > >M > > y28 7 > > > > 6> > y2 > 7> > > > > > 7> > > > > x2 y > 2 > > > 5 x8 y8 > > > > > M7 > > > ; : 2> ; : 2 x8 y8 M8 xy ð11Þ 1

Note that the equations in the system (11) are nothing but the expanded versions of the completeness conditions shown in Eqs. (4a)–(4c). In other words, the metric shape functions Mi, (i = 1, 2, 3, . . . , 8), are obtained by solving the very completeness conditions themselves. Thus, it is no surprise that metric shape functions are able to satisfy the mesh-distortion independent completeness property with respect to all the eight terms: {1, x, y, x2, xy, y2, x2y, xy2}. For the present element, the x- and y-coordinates appearing in Eq. (11) are measured with respect to a local Cartesian coordinate system whose origin is located at the geometric centre of the eleP P ment (xo, yo) where xo ¼ 18 8i¼1 xi and yo ¼ 18 8i¼1 yi . The iso-parametric mapping of element geometry is carried out as usual

1048

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

P8

using equations, x ¼ i¼1 N i xi and y ¼ i¼1 N i yi , where Ni are the parametric shape functions defined in (5)–(7). The shape function matrix is defined as

M¼



2

P8

M1

0

M2

0

M3

0

0

M1

0

M2

0

M3

   M8 0



0 M8

6 BL ¼ 4 0 N1;y

M 1;x 6 BR ¼ 4 0 M 1;y

ð12Þ

:

N 2;x

0

N3;x

0



N8;x

N1;y

0

N2;y

0

N3;y



0

N1;x

N2;y

N2;x

N3;y

N3;x



N8;y

0

M2;x

0

M3;x

0

M 1;y

0

M 2;y

0

M 3;y



   M 8;x

M1;x

M 2;y

M2;x

M 3;y

M3;x

   M8;y

0

M 8;x

where Mi,x and Mi,y (i = 1, 2, 3, . . . , 8) are the derivatives of Cartesian shape functions, which can be computed at the quadrature points directly by solving the derivative-versions of Eq. (11):

38 9 8 9 1 1 1 1 1 1 1 1 > M1;x > > 0 > > > > > > > > > 7 6 x1 > > > > 1 > > x2 x3 x4 x5 x6 x7 x8 7 > M2;x > 6 > > > > > > > > > > 7> 6 > > > > > > > > 6 y1 0 y2 y3 y4 y5 y6 y7 y8 7> M 3;x > > > > > > > 7> 6 > > > = = < < 6 x2 2 2 2 2 2 2 2 7> x x x x x x x 2x M 4;x 6 1 2 3 4 5 6 7 8 7 ¼ 7 6 > 6 x1 y1 x2 y2 x3 y3 x4 y4 x5 y5 x6 y6 x7 y7 x8 y8 7> y > M5;x > > > > > > > > 7> 6 > > > > > > 6 y2 2 2 2 2 2 2 2 7> > > > > 0 > y2 y3 y4 y5 y6 y7 y8 7> M6;x > > > 6 1 > > > > > > > > 7 6 2 > > > > > > > 4 x1 y1 x22 y2 x23 y3 x24 y4 x25 y5 x26 y6 x27 y7 x28 y8 5> 2xy M > > > > 7;x > > > > ; ; : : 2 2 2 2 2 2 2 2 2 y x1 y1 x2 y2 x3 y3 x4 y4 x5 y5 x6 y6 x7 y7 x8 y8 M8;x 2

ð14Þ 38 9 8 9 M1;y > > 0 > > > > > > > > > > > M2;y > > > > 0 > > x8 7 7> > > > > > > > 7> > > > > > > > > > M3;y > 1 > y8 7> > > > > > > > 7> > > > = < 0 > = < x28 7 7 M4;y ¼ : 7 > > > M5;y > x > x8 y8 7> > > > > > 7> > > > > > >M > > > > y28 7 > > 6;y > > 2y > 7> > > > > > > > > 7 > > > 2 > > > > x x28 y8 5> M > > > > 7;y > > > > ; ; : : 2 x8 y8 M8;y 2xy

3

N8;y 7 5; N8;x

ð17Þ Ni;x

)

Ni;y

3

0 7 M8;y 5;

0

where the Ni;x and N i;y may be computed as follows:

(

ð13Þ

2

0



BR matrix is given by

2

N1;x

¼ Adj½J



Ni;n Ni;g

 ð18Þ

:

However, comparing Eqs. (18) and (10), one can realise that computation of Ni;x and N i;y does not demand any additional effort. One can still use Eq. (10) as usual, but reset the Jacobian to unity, i.e., jJj = 1 in the computer code. Thus, computer implementation of condition C3 is quite simple in an existing computer code. All one has to do is just follow the usual procedure as for computing the iso-parametric BL matrix in the classical implementation (outlined in Section 3.1), but reset the Jacobian to unity in all the calculations. This simple technique mimics the action of cancelling out the jJj terms from the numerator and denominator in Eq. (16). 4.2. Nine-node quadratic quadrilateral element For the second illustration, the development of a 9-node quadratic quadrilateral plane-stress element, called PM-QUAD9*, is considered. For this element, sufficient condition C1 is satisfied by using parametric shape functions for constructing BL matrix. Lagrangian shape functions (e.g. see [28]) are used for the purpose:

ð15Þ

N1 ¼ L1 ðnÞL1 ðgÞ N2 ¼ L3 ðnÞL1 ðgÞ N3 ¼ L3 ðnÞL3 ðgÞ N4 ¼ L1 ðnÞL3 ðgÞ; ð19aÞ N5 ¼ L2 ðnÞL1 ðgÞ N6 ¼ L3 ðnÞL2 ðgÞ N7 ¼ L2 ðnÞL3 ðgÞ N8 ¼ L1 ðnÞL2 ðgÞ; ð19bÞ nð1  nÞ nð1 þ nÞ L1 ðnÞ ¼  L2 ðnÞ ¼ ð1 þ nÞð1  nÞ L3 ðnÞ ¼ ; 2 2 ð19cÞ gð1  gÞ gð1 þ gÞ L1 ðgÞ ¼  L2 ðgÞ ¼ ð1 þ gÞð1  gÞ L3 ðgÞ ¼ ; 2 2 ð19dÞ

Note that in Eqs. (11), (14) as well as (15), the interpolation matrix on the left side is the same and hence for computing the shape functions Mi or the derivatives Mi,x and Mi,y, the interpolation matrix needs to be inverted only once per element.

where the subscripts 1–4 refer to corner nodes and the rest the midside nodes. Sufficient condition C2 is satisfied by using metric shape functions to construct BR matrix exactly as in Section 4.1.2 but with an additional term x2y2 in the basis, viz., {1, x, y, x2, xy, y2, x2y, xy2, x2y2}. Metric shape functions are obtained by solving the system:

1

1

1

1

1

1

1

6 x1 x2 x3 x4 x5 x6 x7 6 6 6 y1 y2 y3 y4 y5 y6 y7 6 6 x2 x22 x23 x24 x25 x26 x27 6 1 6 6 x1 y1 x2 y2 x3 y3 x4 y4 x5 y5 x6 y6 x7 y7 6 6 y2 y22 y23 y24 y25 y26 y27 6 1 6 2 4 x1 y1 x22 y2 x23 y3 x24 y4 x25 y5 x26 y6 x27 y7 x1 y21 x2 y22 x3 y23 x4 y24 x5 y25 x6 y26 x7 y27

1

2 4.1.3. Satisfaction of condition C3 Firstly, we note that computation of BL involves Jacobian jJj (see Eqs. (9) and (10) whereas computation of BR does not (see Eqs. (13)–(15)). Secondly, we note that, by virtue of Eqs. (9) and (10), BL can be rewritten as BL ¼ BL =jJj where the matrix BL is constructed similar to that in Eq. (9), but with the derivatives        Ni;x N i;n N i;x ¼ Adj½J instead of the derivatives defined  N i;y N i;g N i;y by Eq. (10). Hence, Eq. (3) can be written as

ð16Þ Notice that the Jacobian jJj gets cancelled in Eq. (16) and hence the resulting stiffness integral is free of the Jacobian term. As a consequence, the stiffness can be computed without any fear of Jacobian going negative. This interesting property enables this element work even for extremely severe distortions. The expression for BL can be obtained as follows:

1 6 x1 6 6 6 y1 6 2 6 x 6 1 6 6 x1 y 1 6 2 6 y1 6 6 x2 y 6 1 1 6 4 x1 y21 x21 y21 8 > > > > > > > > > > > > > > > < ¼ > > > > > > > > > > > > > > > :

1 1 x2 x3 y2 y3 x22 x23 x2 y2 x3 y3 y22 y23 2 x2 y2 x23 y3 x2 y22 x3 y23 x22 y22 x23 y23 9 1 > > > x > > > > > y > > > > 2 > x > > = xy : > > > y2 > > > 2 > > x y> > > > > xy2 > > ; x2 y2

1 x4 y4 x24 x4 y 4 y24 x24 y4 x4 y24 x24 y24

1 x5 y5 x25 x5 y5 y25 x25 y5 x5 y25 x25 y25

1 x6 y6 x26 x6 y6 y26 x26 y6 x6 y26 x26 y26

1 x7 y7 x27 x7 y 7 y27 x27 y7 x7 y27 x27 y27

1 x8 y8 x28 x8 y 8 y28 x28 y8 x8 y28 x28 y28

38 9 1 > M1 > > > > > > M2 > x9 7 > 7> > > > 7> > > > y9 7> M3 > > > > > 7 > > > > M4 > x29 7 > = < 7> 7 x9 y9 7 M 5 > 7> > > > M6 > y29 7 > > > 7> > > > > 7 2 > > x9 y9 7> M > 7 > > > > 7 > 2 5> > > x9 y9 > M > 8 > > ; : 2 2 x9 y9 M9

ð20Þ

1049

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

y (10,2)

(5,2)

(0,2)

2000

1 (0,1)

(5, 1)

2

(10,1)

B

A (5,0)

(0,0)

(10,0)

2000

x

Fig. 1. A tip-moment cantilever beam modelled with a distorted mesh. Mesh distortion is introduced by moving node A from its original position (5,1). The material is assumed to be linear elastic with material properties E = 1500 and v = 0.25.

The condition C3 is also satisfied as in Section 4.1.3. The numerical results for this element are discussed in Section 5.7.

5. Numerical examples Numerical results for 8-node element will be presented first. The results for 9-node element will be presented in Section 5.7. The following three types of 8-node elements are considered for comparison purposes:

Table 1 Vertical tip displacement at node B for tip-moment cantilever beam (Fig. 1). Coordinates of node A xA

yA

0.1 1 2 3 4 5 6 7 8 9 9.9

0.02 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 1.98

Negative Jacobian at Gauss points?

PP-QUAD8a

MM-QUAD8

PM-QUAD8* (present)

Yes Yes Yes No No No No No Yes Yes Yes

9.011 8.312 15.110 24.857 47.953 100.000 48.227 25.443 18.093 14.571 12.198

97.599 85.740 125.449 107.822 103.401 100.000 100.203 101.247 104.582 106.138 112.534

100.000 100.000 100.000 100.000 100.000 100.000 100.000 100.000 100.000 100.000 100.000

KðeÞ ¼

Z

1

1

Z

1

ðBL ÞT DBL jJjt dn dg;

FðeÞ ¼

Z

NT tt dL;

ð21Þ

L

1

where t is the distributed force vector which is integrated over the loaded edge (L) of the element as usual. Substituting BL ¼ BL =jJj, we can write

Exact solution = 100 a

(i) PP-QUAD8: This is just the same as the classical iso-parametric 8-node quadrilateral element sometimes acronymed in the literature as Q8 or sometimes QUAD8. This is a symmetric element with BL matrix being used both on the left and right side of the D matrix in the stiffness integral. Recall that parametric shape functions are used to construct BL matrix. The prefix ‘PP’ here signifies the fact that parametric shape functions are used on the left as well as the right side of the D matrix. The element stiffness matrix and force vector are computed as

Alias QUAD8.

3000 2000 1000 0 -1000 -2000

(a) PM-QUAD8* (present)

-3000 2000 1000 0 -1000

(b) PP-QUAD8 (alias QUAD8)

-2000 4000 2000 0 -2000 -4000 -6000

(c) MM-QUAD8 Fig. 2. rxx-stress contours for the tip-moment cantilever beam (Fig. 1).

1050

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

20

6

17 5

19

16 7

14

4

8

5

1 4

3

2

13

2

15

3

18

11 12 10

1

9

Fig. 3. A severely distorted 5-element mesh. The nodal coordinates of elements are listed in Table 2.

ðeÞ

K

¼

Z

1

1

Z

1



1

T BL DBL t dn dg; jJj

FðeÞ ¼

Z

NT tt dL;

ð21aÞ

L

which shows the Jacobian effectively appears in the denominator for this element. (ii) MM-QUAD8: This is another possible symmetric element but is unpopular because the shape functions loose compatibility under distorted geometries. Here, BR matrix is used both on the left and right side of the D matrix in the stiffness integral. Recall that metric shape functions are used to construct BR matrix. The prefix ‘MM’ signifies that metric shape functions are used on the left as well as the right side of the D matrix. The element stiffness matrix and force vector are computed as

KðeÞ ¼

Z

1

1

Z

1

ðBR ÞT DBR jJj t dn dg FðeÞ ¼

Z

MT tt dL;

ð22Þ

L

1

where the shape function matrix M is used to evaluate the element force vector. The BR matrix does not involve Jacobian. Hence, the Jacobian appears in the numerator for this element. (iii) PM-QUAD8*: This is the present element, an improved version of US-QUAD8 [13] alias PM-QUAD8 [14]. This is an unsymmetric element (as BR – BL). The prefix ‘PM’ signifies that parametric and metric shape functions are used to obtain BL and BR matrices, respectively. The *signifies that the element satisfies all the three sufficiency conditions, C1, C2 and C3. The element stiffness matrix and force vector are computed as

KðeÞ ¼

Z

1

1

Z

1

1



T BL DBR t dn dg;

FðeÞ ¼

Z

NT tt dL:

ð23Þ

L

Note that the nodal force vector is just the same as for PPQUAD8 element in (i). The Jacobian is absent for this element. For all computations, a 3  3 numerical integration is used to evaluate the element stiffness unless stated otherwise. 5.1. A cantilever beam with a tip-moment load A cantilever beam subjected to a tip moment involves a quadratic displacement field and hence is a good problem to test the mesh-distortion sensitivity of a quadratic element [8,14]. We consider a mesh (Fig. 1) with extreme mesh distortion for which the Jacobian goes negative in some regions of element domain. The distortion introduced in the mesh is controlled by moving node A from its original position of (5,1). The vertical tip deflection at point B is computed for several positions of node A and the results are listed in Table 1. The results demonstrate that the exact solution is captured by PM-QUAD8* irrespective of the mesh distortion introduced. The results obtained with PP-QUAD8 and MM-QUAD8 are also listed for comparison. PP-QUAD8 and MM-QUAD8 produce inferior results.

The contour-plots of rxx-stress are shown in Fig. 2. For these plots, node A was located at (2,0.4). The rxx-stress for this problem has a linear variation across the depth of the beam and hence the contour lines must appear straight and horizontal. PM-QUAD8* reproduces exact stress contours while that of PP-QUAD8 and MM-QUAD8 give very inaccurate ones. The lack of continuity and smoothness of contours at the element interface in the case of PP-QUAD8 and MM-QUAD8 is a measure of their mesh-distortion sensitivity. The problem is also solved using an extremely distorted 5-element mesh as shown in Fig. 3. The nodal coordinates and element connectivity of the mesh are listed in Table 2. The vertical tip displacements at node 18 given by PM-QUAD8*, PP-QUAD8 and MMQUAD8 are 100, 15.534 and 105.576, respectively, against the exact value of 100. PM-QUAD8 element is able to reproduce exact solution. The stress contour-plots are shown in Fig. 4. Here again, PMQUAD8* element is able to capture the exact rxx-stress distribution despite extremely high mesh distortions. 5.2. A cantilever beam with a tip-shear load The same problem as in Section 5.1 is considered here, but with a tip shear loading (Fig. 5). This is a cubic displacement problem. The shear force of 300 units is distributed parabolically, which is

Table 2 Nodal coordinates and element connectivity for the mesh shown in Fig. 3. Node No.

x-coord.

y-coord.

Nodal coordinates 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Element No.

0.00 3.00 1.75 3.20 7.00 0.00 1.65 0.00 10.00 7.00 6.75 7.00 4.00 6.85 7.50 4.00 10.00 10.00 7.00 2.60

0.00 0.70 0.80 1.30 1.00 2.00 1.10 1.00 0.00 0.00 0.70 0.50 0.60 1.20 0.95 1.40 2.00 1.00 1.76 2.00

Element connectivity 1 2 3 4 5

Nodes 1, 1, 9, 2, 4,

3, 2, 5, 4, 7, 6, 8 10, 9, 12, 11, 13, 2, 3 18, 17, 19, 14, 15, 11, 12 13, 11, 15, 14, 16, 4, 5 16, 14, 19, 17, 20, 6, 7

1051

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

3000 2000 1000 0 -1000 -2000 -3000

(a) PM-QUAD8* (present)

4000 2000 0 -2000 -4000

(b) PP-QUAD8 (alias QUAD8) 3000 2000 1000 0 -1000 -2000 -3000

(c) MM-QUAD8 Fig. 4. rxx-stress contours for the tip-moment cantilever beam with the severely distorted 5-element mesh (Fig. 3).

y (10,2)

(5,2)

(0,2) 1 (0,1)

(10,1)

2

(5, 1)

B

300

A (0,0)

(5,0)

(10,0)

x

Fig. 5. A tip shear cantilever beam modelled with a distorted mesh.

simulated by applying equivalent nodal forces of 30, 240 and 30 units at nodes (10, 0), (10, 1) and (10, 2), respectively. The reaction shear force at the fixed end is represented by 30, 240 and 30 units of load at nodes (0, 0), (0, 1) and (0, 2), respectively. The vertical tip deflections computed at point B are listed in Table 3. Since this is a cubic displacement problem, neither PM-QUAD8* nor the other two elements can be expected to reproduce exact solution because they are all quadratic elements. Nevertheless, we notice that the results given by PM-QUAD8* element are less sensitive to distortion as compared to the other two elements. For the severely distorted mesh (Fig. 3) solved with tip-shear loading, the tip displacements obtained with PM-QUAD8*, PPQUAD8 and MM-QUAD8 elements are 14.580, 19.424 and 95.118, respectively, as against the overkill solution of 103.495. At the face of it, the PM-QUAD8* element seems to give very bad solution for this mesh. However, the bad solution is rather due to insufficient order of numerical integration (3  3) than due to the formulation itself. Very high order monomial terms (in n,g coordinates) appear in the stiffness integrand due to severe quadratic (curved-edge) mesh distortions in the mesh. With a 4  4 Gaussian quadrature, the tip displacements given by PM-QUAD8*, PPQUAD8 and MM-QUAD8 are 106.565, 8.205 and 89.611, respectively. (Note that, for a regular mesh without mesh-distortions, a

2  2 quadrature is sufficient to give exact integration for any of the three elements considered, viz., PM-QUAD8*, PP-QUAD8 or MM-QUAD8.)

Table 3 Vertical tip displacement at node B for the tip-shear cantilever beam (Fig. 5). Coordinates of node A xA

yA

0.1 1 2 3 4 5 6 7 8 9 9.9

0.02 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 1.98

Negative Jacobian at Gauss points?

PP-QUAD8a (alias QUAD8)

MMQUAD8

PMQUAD8* (present)

Yes Yes Yes No No No No No Yes Yes Yes

13.994 11.500 20.305 31.322 54.040 100.543 54.094 31.913 24.549 20.693 18.268

91.398 75.258 135.189 110.132 105.585 100.543 99.754 100.021 102.004 106.132 113.343

102.157 97.333 95.243 95.641 98.073 100.543 102.912 104.741 104.790 102.837 97.827

Overkill solution = 103.495b a b

Alias QUAD8. Obtained using ANSYS Plane82 element with a 100  20 mesh.

1052

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Table 4 Convergence of tip displacement at node B with increasing order of Gaussian quadrature for 2-element mesh (Fig. 1) with nodal coordinates A (2, 0.4).

a

Problem

Element type

Order of Gaussian quadrature 22

33

44

55

66

Tip-moment cantilever beam

PP-QUAD8a MM-QUAD8 PM-QUAD8* (present)

75.271 237.400 106.933

15.110 125.449 100.00

12.923 114.519 100.00

10.855 114.031 100.00

15.043 114.018 100.00

Tip shear cantilever beam

PP-QUAD8a MM-QUAD8 PM-QUAD8* (present)

63.293 265.785 94.741

20.305 135.189 95.243

17.439 123.813 93.400

14.828 123.275 93.400

20.614 123.266 93.400

Alias QUAD8.

Table 5 Convergence of tip displacement at node B with increasing order of Gaussian quadrature for 2-element mesh (Fig. 1) with nodal coordinates A (3, 0.6).

a

Problem

Element type

Order of Gaussian quadrature 22

33

44

55

66

Tip-moment cantilever beam

PP-QUAD8a MM-QUAD8 PM-QUAD8* (present)

74.201 134.176 103.188

24.857 107.822 100.000

24.177 107.402 100.000

24.119 107.404 100.000

24.113 107.404 100.000

Tip shear cantilever beam

PP-QUAD8a MM-QUAD8 PM-QUAD8* (present)

69.062 137.691 95.438

31.322 110.132 95.641

30.695 110.246 95.092

30.647 110.243 95.092

30.642 110.243 95.092

Alias QUAD8.

Table 6 Convergence of tip displacement at node 18 with increasing order of Gaussian quadrature for 5-element mesh (Fig. 3). Problem

a

Element type

Order of Gaussian quadrature 22

33

44

55

66

Tip-moment cantilever beam

a

PP-QUAD8 MM-QUAD8 PM-QUAD8* (present)

30.896 108.825 103.400

15.534 105.576 100.000

6.555 102.484 100.000

14.315 102.309 100.000

8.195 102.307 100.000

Tip shear cantilever beam

PP-QUAD8a MM-QUAD8 PM-QUAD8* (present)

16.823 114.900 93.163

19.424 95.118 14.580

8.205 89.611 106.565

58.129 89.272 106.565

10.599 89.258 106.565

Alias QUAD8.

5.3. Order of Gaussian quadrature for PM-QUAD8* The purpose of this section is to show that PM-QUAD8* element gives convergent results with increasing order of numerical integration. The convergence of nodal displacements with the order of Gaussian quadrature is studied. Table 4 lists the tip displacements computed for the 2-element mesh (Fig. 1) with a severe distortion, A (2, 0.4). The converged results are shown in boldface. It is seen that PM-QUAD8* gives convergent results with increasing order of integration while PPQUAD8 exhibits divergence. MM-QUAD8 also seems to give convergent results although the convergence is slower. Apparently, PP-QUAD8 gives divergent results because the Jacobian is zero in some part of element domain and this seems to affect the stiffness integral adversely. In order to verify this conjecture, the computations are repeated for a milder distortion, A (3, 0.6), for which the Jacobian remains positive in the element domains. Table 5 shows the results for this case, which suggests that PP-QUAD8 produces convergent results now. Table 6 shows similar results for the severely distorted 5-element mesh. Here again, we notice that PM-QUAD8* converges faster, MM-QUAD8 converges slowly and PP-QUAD8 gives divergent results. This exercise clearly demonstrates the reliability of PM-QUAD8* element with respect to numerical integration. As the order of inte-

gration is increased, convergence to exact integration is guaranteed. From Tables 4–6, the minimum order required for exact integration is 3  3 and 4  4, respectively, for the tip-moment and tip-shear problems. In order to confirm that the integration beyond this minimum order always is indeed exact, the results of the last row of Table 4 are listed in Table 7 with 16 decimal places displayed. It is seen that for an order 4  4 and above, the computed displacement is the same up to about 10 significant places, which confirms that for orders 4  4 or more, exact integration results. Under severe mesh-distortions, it is possible that the stiffness integrand is a many-to-single valued mapping (i.e., many sets of (n,g) leading to the same integrand value) in some parts of the element domain. An intriguing question then arise: How is it that exact integration of stiffness has been possible for PM-QUAD8*? Why

Table 7 Exact integration for PM-QUAD8*: results of the last row of Table 4 are reproduced with 16 decimal places. Order

Tip displacement

22 33 44 55 66

94.74093149950346 95.24263100009034 93.39955413644687 93.39955413656045 93.39955413624790

1053

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

shown to be a nonzero function. Note that, in the stiffness integral of PP-QUAD8, the determinant of Jacobian matrix appears in the denominator of the integrand (Eq. (21a)). For extreme distortions, the determinant may become zero at some points in the element domain or its neighbourhood, so the integrand function becomes singular at those points. This causes steep gradients in the vicinity of singular points. Steep gradients render exact integration of stiffness integral difficult with a finite order of Gaussian quadrature. Thus, the convergence of numerical integration becomes very slow and sometimes unstable in the case of PP-QUAD8 whenever the mesh distortion is severe, which is evidenced in Tables 4 and 6. (iv) In the stiffness integral of MM-QUAD8, one can see that the Jacobian appears in the numerator (Eq. (22)). Thus, MMQUAD8 does not experience any singularity of stiffness integral as in PP-QUAD8 element. However, the convergence of results is slower compared to PM-QUAD8 element because the presence of Jacobian in the numerator increases the polynomial order of the integrand as compared to PMQUAD8.

is it that exact integration was not possible for PP-QUAD8 element in spite of using higher order integration? Investigations were done on these aspects and the following are some important conclusions. (i) Integration of many-to-single valued functions caused by highly warped element domains does not seem to be a problem in Gaussian quadrature. In order to illustrate this point, a distorted mesh is considered as shown in Fig. 5. This mesh is obtained by taking the coordinates of node A as (0.1, 0.02). Notice that there is overlap of element domains in the region 4-3-2-4 of Fig. 6a. Nevertheless, the areas of the individual elements as computed by Gaussian quadrature are 3.4666. . . and 16.5333. . ., respectively, the sum of which correctly adds up to 20 which is the correct area of the element, i.e., length (10 units) times depth (2 units) of the cantilever beam. It appears that Gaussian quadrature computes the area of the overlapping region 4-3-2-4 (Fig. 6b) as a negative area with respect to element No. 1 and positive area with respect to element No. 2. As a result, the overlapping area does not enter into the area calculation. Thus, Gaussian quadrature over extremely distorted element geometries does not pose any problem. (ii) As long as the integrand of the stiffness integral is a polynomial function (which is the case in PM-QUAD8*), the integrand is guaranteed to be a non-singular function in the element domain (1 6 n 6 1 and 1 6 g 6 1). Hence, the integrand function has no steep gradients or discontinuities, and therefore, exact integration is possible with a finite order of Gaussian quadrature. This explains why exact integration of stiffness matrix has been possible for PM-QUAD8*

The mesh distortions shown in Fig. 6 are, of course, too hypothetical. Although severe mesh distortions are likely in practical applications such as forging and extrusion, it is hard to imagine an application where the mesh distortions are so severe as in Fig. 6. Nevertheless, the mesh in Fig. 6 is useful to demonstrate that the PM-QUAD8* element has Jacobian-unlimited mesh-distortion immunity (i.e., the element’s performance is insensitive to mesh distortion whether or not the Jacobian is positive). 5.4. A special post-processing issue associated with PM-QUAD8*

element even for extremely distorted meshes. (iii) Whenever the integrand function has a rational form (i.e., in the form P(n,g)/Q(n,g)), the non-singularity of integrand cannot be guaranteed, unless the denominator Q(n,g) can be

a

6

7

Although PM-QUAD8* element can give accurate nodal displacement results even for extreme mesh distortions with negative

5

13

12

1 8

2

1

3

4

b

11

2

6

7

9

10

5

1 8

2 1 4

c

3

5

13

11

2

4

3

12

9

10

Fig. 6. Element geometry under severe curved-edge mesh distortion. (a) Mesh, (b) Element No. 1, and (c) Element No. 2.

1054

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Jacobian, difficulties may arise during the post-processing stage, i.e., while attempting to contour-plot the displacements and stresses over the element domains. The extreme distortions in the mesh may cause overlaps (or folds) of element domains due to nonuniqueness of element mapping. Hence, classical post-processing algorithms may experience difficulty in plotting the displacement/stress contours because they often use iso-parametric mapping of element geometry in the process. Fig. 7 shows an example of such a troublesome stress plot. Note that the stress contours are absent in the lower part of element No. 1. This does not mean that stresses are not being calculated for that part of the element. It has been verified that the difficulty is not with computing the stresses, but rather with plotting them. This is because the mapped element geometry has overlaps and the classical contour-plotting algorithm cannot handle the overlaps. However, there is no such problem experienced in element No. 2 as the there is no overlap of element geometry. Thus, for contour-plotting of displacements and stresses, classical plotting algorithms may not work whenever the mapping of element geometry is non-unique.

Better techniques need to be explored. (Such techniques may purely rely on the nodal data rather than element connectivity.) Nevertheless, if the interest is only computing the nodal quantities (viz., nodal displacements/strains/stresses) and not contourplotting over the element domain, this special problem associated with post processing is not an issue. This issue will again be addressed in Section 5.6. 5.5. Nodal rotation patch-test A square panel of size 2  2 and unit thickness is loaded on the right edge as shown in Fig. 8. Only three points on the left edge are restrained as in Fig. 8. Tensile, moment and shear loadings are distributed forces as rx = 1000, rx = 1500y and sxy = 750(1  y2), respectively, so as to give a net tensile load of 2000, a moment of 1000 and a shear load of 1000 as marked in Fig. 8. The equivalent lumped forces applied at the nodes on the right edge are shown in Fig. 8. Reaction forces applied at the left edge are also shown. The above three loadings induce constant strain, linear strain and qua-

3000 2000 1000 0 -1000 -2000

Fig. 7. rxx stress contours for the tip-moment cantilever beam under severe curved-edge mesh distortion. Classical contour plotting algorithms may fail to plot the contours correctly.

A

2000

1000

1000

Equivalent lumped forces applied at the nodes on the right edge Load type Tension Moment Shear

Load type Tension Moment Shear

Direction of force applied

Node number (Figure 9) 8 2000/3 -250 350

9 1000/6 -500 12.5

Reaction forces applied at the nodes on the left edge Direction of force Node number (Figure 9) applied 1 16 15 14 Fx -1000/6 -2000/3 -1000/3 -2000/3 Fx -250 -500 0 250 Fx -500 -1000 0 1000 Fy -12.5 -350 -275 -350

13 -1000/6 500 500 -12.5

Fx Fx Fy

5 1000/6 250 12.5

6 2000/3 500 350

7 1000/3 0 275

Fig. 8. Square panel with tensile, moment and shear loading for the nodal rotation patch test. The material is assumed to be linear elastic with material properties E = 1500 and v = 0.25.

1055

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

13

12

11

14

15

19

20

18

21

16

1

10

17

2

9

13

8

14

7

15

6

3

4

5

12

19

2

12

10

14

15

21

9

13

8

14

7

15

6

16

3

4

5

1

11

10

14 17 15

21 18

2

4

5

12

11

10

9

17

8

18

21

20

7

6

19

2

3

4

5

9

13

8

14

7

15

6

16

5

1

12

11

10

20

9

8

20

17

19

21

7

19

16

1

3

Mesh d (θ = 175 )

Mesh c (θ = 135 ) 12

6

o

o

13

7

20

19

2

17

17

16

1

21

Mesh b (θ = 90o)

11

18

9

8

20

Mesh a (θ = 0o) 13

10

18

16

1

11

3

4

Mesh e (θ = 225 ) o

18

2

6

3

4

5

Mesh f (θ = 270 ) o

Fig. 9. Nodal rotation patch-test to assess mesh-distortion immunity. Distorted meshes are obtained by rotating nodes 17, 18, 19 and 20 of Mesh a about node 21 by an angle h in the anti-clockwise direction. See Fig. 8 for forces applied. Note that Meshes c, d and e have extreme distortions with overlapping element geometries.

dratic strain conditions in the panel, respectively. Six meshes, as shown in Fig. 9, are considered. Distorted meshes are obtained by rotating nodes 17–20 (of Mesh a of Fig. 9) about node 21 by an angle h in the anti-clockwise direction. The crossing of element boundaries in Meshes c–e show the presence of extreme mesh distortions with overlaps in the element geometry.

Mesh-distortions to the extent seen in Meshes c–e may not be encountered in practical applications. However, severe mesh-distortions leading to negative Jacobian is often possible in hot metal forming simulations. The maximum deformation that can be simulated in large deformation processes like forging and extrusion is limited by the occurrence of negative Jacobian. Once zero or neg-

1056

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

ative Jacobian appears, the simulation by a finite element software usually aborts and remeshing may be necessary to continue simulation. Since the Jacobian dependence is avoided in the proposed elements, the maximum deformation that can be simulated can possibly be extended well beyond that of classical elements. Thus, although we do not deal with plastic deformation in this paper, Meshes c–e may give us an indication as to how the proposed element might behave in such situations. The computed vertical displacements at the mid-point (i.e., point A) of the right edge are listed in Tables 8–10. For all computations in this section, the stiffness matrix was integrated with a 4  4 Gaussian quadrature. Tables 8 and 9 show that the present element (PM-QUAD8*) gives exact solution for constant and linear strain patch-tests, respectively, and thus passes the constant and linear strain patch-tests in the strong form. It is interesting to note that even for extreme mesh-distortions with overlapping element geometries (Meshes c–e), the present element able to reproduce the exact solution. For the quadratic strain patch-test (Table 10), the present element shows some sensitivity to mesh distortion, although much less compared to the other two elements. Quadratic strain patch-test is, however, passed by the present element in the weak form as detailed in Section 5.6. 5.6. Convergence test with extremely distorted meshes The convergence of displacements with mesh refinement is studied for the shear loading problem of Fig. 8. Three sequences of meshes with 1  1, 2  2 and 4  4 super-elements are considered (Fig. 10). In order to simulate the shear loading, the equivalent lumped forces applied in the y-direction at the nodes on the right edge of the three meshes are {12.5, 350, 275, 350, 12.5}, {1.5625, 106.25, 96.875, 231.25, 128.125, 231.25, 96.875, 106.25, 1.5625} and {0.195313, 28.9063, 27.7344, 75.7813, 47.2656, 107.031, 58.9844, 122.656, 62.8906, 122.656, 58.9844, 107.031, 47.2656, 75.7813, 27.7344, 28.9063, 0.195313}, respectively. Negative of these values are applied as reactions on the left edge. The reaction forces applied in the x-direction are {500, 1000, 0, 1000, 500}, {250, 750, 250, 250, 0 250, 250, 750, 250} and {125, 437.5, 187.5, 312.5, 125, 187.5, 62.5, 62.5, 0, 62.5, 62.5, 187.5, 125, 312.5, 187.5, 437.5, 125}, respectively. Three sets

Table 8 Patch test results for tensile loading (Fig. 8). Mesh

PP-QUAD8a

MM-QUAD8

PM-QUAD8* (present)

a b c d e f

1.0 1.0 1.0 1.0 1.0 1.0

1.0 1.0 1.0 1.0 1.0 1.0

1.0 1.0 1.0 1.0 1.0 1.0

Exact solution = 1.0 a

Alias QUAD8.

Table 9 Patch test results for moment loading (Fig. 8). PP-QUAD8a

MM-QUAD8

PM-QUAD8* (present)

a b c d e f

1.5000 1.1010 0.9206 1.2447 0.9206 1.1010

1.5000 1.7912 1.6705 1.6051 1.6705 1.7912

1.5000 1.5000 1.5000 1.5000 1.5000 1.5000

Alias QUAD8.

Mesh

PP-QUAD8a

MM-QUAD8

PM-QUAD8* (present)

a b c d e f

3.3069 2.7898 2.5739 4.3555 2.5739 2.7898

3.3069 3.7555 3.4742 3.4712 3.4742 3.7555

3.3069 3.4453 3.3201 3.3391 3.3201 3.4453

Exact solution = 3.3125 a

Alias QUAD8.

of super-element meshes are constructed using Mesh a Mesh b, and Mesh c (of Fig. 8) as super-element. For all computations in this section, the stiffness matrix was integrated with a 4  4 Gaussian quadrature. The computed vertical displacements at point A are shown in Tables 11–13. The energy norm error in strain is computed as

eerror

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u1 R T u ðeexact  eFEM Þ Dðeexact  eFEM Þ dV ¼ t2 V R 1 ðeexact ÞT Deexact dV 2 V

ð24Þ

and is plotted in Fig. 11. PM-QUAD8*not only gives a convergent solution with an order of convergence close to the theoretical value 2, but also more accurate solution compared to the other two elements. The high convergence order of 3.01 given by PP-QUAD8 for mesh sequence B is not reliable; Due to the presence of Jacobian in the denominator, the element in fact gives different orders of convergence with different orders of numerical integration of stiffness matrix. The rxx-stress contours computed for 4  4 super-element meshes are shown in Figs. 12–14. It is seen the stress contours for the distorted meshes given by the present element (Fig. 13) are in close agreement with that of the undistorted mesh (Fig. 12). This suggests that, despite extremely severe mesh distortions, the present element is able to converge to correct stress results. In contrast, the stress contours given the by the classical element (PP-QUAD8) (Fig. 14) are highly discontinuous. Although the displacement results of PP-QUAD8 element are not that bad (Tables 11–13), the stress results are indeed very unreliable. In summary, the convergence test suggests that the folds and overlaps in the element geometry caused by extreme mesh distortions do not hinder the convergence of either displacements or stresses given by the present element. The results also suggest that the present element is capable of capturing quadratic strain variations as the mesh is refined. In other words, the element passes quadratic strain patch-test in the weak form. Furthermore, the smooth stress contours for 4  4 mesh reveal that the classical contour-plotting algorithm (which was used for plotting the stress contours) is indeed able to handle contour-plotting of stresses despite severe mesh-distortions of the element. The special post-processing problem discussed in Section 5.4 does not seem to exist for this fine mesh. 5.7. Numerical results for 9-node quadrilateral elements

Mesh

Exact solution = 1.5 a

Table 10 Patch test results for shear loading (Fig. 8).

Three types of 9-node element, viz., PP-QUAD9, MM-QUAD9 and PM-QUAD9*, are formulated following a procedure similar to that of PP-QUAD8, MM-QUAD8 and PM-QUAD8* introduced at the beginning of Section 5. For all computations in this section, a 5  5 Gaussian quadrature is used to evaluate the element stiffness matrix. The rotational patch-test of Section 5.5 is repeated using the 9node elements. The same meshes as in Fig. 9 are used here with the

1057

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Fig. 10. Sequences of meshes for convergence test. Sequences A, B and C are constructed with Mesh a, Mesh b and Mesh c as the super-element.

Table 11 Convergence results with Mesh a as super-element. Super-element mesh

11 22 44

Table 13 Convergence results with Mesh c as super-element.

Vertical displacement at the mid-point of loaded edge PP-QUAD8a

MM-QUAD8

PM-QUAD8* (present)

3.3069 3.3117 3.3124

3.3069 3.3117 3.3124

3.3069 3.3117 3.3124

Super-element mesh

11 22 44

Exact solution = 3.3125 a

Table 12 Convergence results with Mesh b as super-element.

11 22 44

Vertical displacement at the mid-point of loaded edge PP-QUAD8a

MM-QUAD8

PM-QUAD8* (present)

2.7898 2.9948 3.1795

3.7555 3.7683 3.7718

3.4453 3.3132 3.3125

Exact solution = 3.3125 a

Alias QUAD8.

MM-QUAD8

PM-QUAD8* (present)

2.5737 2.8572 3.0732

3.4742 3.4479 3.5420

3.3201 3.2983 3.3105

Exact solution = 3.3125 a

Alias QUAD8.

Super-element mesh

Vertical displacement at the mid-point of loaded edge PP-QUAD8a

Alias QUAD8.

ninth node of the element located (for simplicity) at the centroid of the four corner nodes of the element. The computed results are shown in Tables 14–16. Tables 14 and 15 show that PM-QUAD9* gives exact solution for constant and linear strain patch-tests, respectively, and thus passes the constant and linear strain patch-tests in the strong form. Even for meshes with overlaps and folds (Meshes c–e), the present element reproduces exact solution. For the quadratic strain patch-test (Table 16), the present element shows much less mesh-distortion sensitivity compared to the other two elements.

1058

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063 -1 -1.2

log10(Energy norm error)

-1.4

PP-QUAD8 PM -QUAD8 * (p res ent)

-1.6 -1.8

R = 1.94 -2 -2.2 -2.4 -2.6 -0.35

-0.15

0.05

0.25

log10(h)

(a) Mesh sequence A 0 -0.2

-1

-1.5

R = 1.07 log10(Energy norm error)

log10(Energy norm error)

-0.5

R = 3.01 R = 2.08

-2

-0.7

-1.2

-1.7

PP-QUAD8

R = 1.88 PP-QUAD8

PM -QUAD8 * (p res ent) -2.5 -0.35

-0.15

0.05

0.25

log10(h)

(b) Mesh sequence B

PM -QUAD8 * (p res ent) -2.2 -0.35

-0.15

0.05

0.25

log10(h)

(c) Mesh sequence C

Fig. 11. Convergence of energy norm strain error for shear loading of square panel (modelled with 8-node elements).

The rxx-stress contours computed for the 4  4 super-element meshes are shown in Figs. 16–18. Here, again, the stress contours given by the present element for the distorted meshes (Fig. 17) are in close agreement with that of the undistorted mesh (Fig. 16). However, the stress contours given the by the classical PP-QUAD9 (Fig. 18) are highly discontinuous. 5.8. Random internal nodes test

Fig. 12. rxx-stress contours given by PP-QUAD8 (alias QUAD8) as well as PMQUAD8* (present) for regular 4  4 mesh.

As a further verification of Jacobian-unlimited mesh-distortion immunity of PM-QUAD8* and PM-QUAD9* elements, all the internal nodes of the mesh shown in Fig. 3 were assigned distinct (noncoincident) random positions within the rectangular domain of cantilever beam using a random number generator, and the tip displacements were computed for the tip-moment problem of Fig. 1 using PM-QUAD8* and PM-QUAD*. The computed tip displacements have been observed to be invariably 100, the exact solution, irrespective of the random positioning of internal nodes. 5.9. Thick cylinder under internal pressure – an example of nonpolynomial displacement field

The results of convergence test (similar to the one in Section 5.6) are shown in Tables 17–19, and plotted in Fig. 15. Here again, the present element (PM-QUAD9*) gives a convergent solution with less error (and with a convergence order close to the theoretical value of 2) as compared the classical PP-QUAD9 element.

The numerical examples considered so far involve polynomial displacement field. In order to investigate how the elements perform under non-polynomial fields, a thick cylinder with internal pressure is considered. Exploiting symmetry of the problem, one

1059

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Fig. 13. rxx-stress contours for distorted meshes given by PM-QUAD8* (present).

Fig. 14. rxx-stress contours for distorted meshes given by PP-QUAD8 (alias QUAD8).

Table 14 Patch test results for tensile loading (Fig. 8).

Table 16 Patch test results for shear loading (Fig. 8).

Mesh

PP-QUAD9a

MM-QUAD9

PM-QUAD9* (present)

Mesh

PP-QUAD9a

MM-QUAD9

PM-QUAD9* (present)

a b c d e f

1.0 1.0 1.0 1.0 1.0 1.0

1.0 1.0 1.0 1.0 1.0 1.0

1.0 1.0 1.0 1.0 1.0 1.0

a b c d e f

3.2979 1.9088 0.4993 2.5498 0.4993 1.9088

3.2979 3.9559 3.8294 3.3657 3.8294 3.9559

3.2979 3.2710 3.3044 3.3677 3.3044 3.2710

Exact solution = 1.0 a

Exact solution = 3.3125 a

Alias QUAD9.

Table 15 Patch test results for moment loading (Fig. 8).

Table 17 Convergence results for 9-node elements with Mesh a as super-element.

Mesh

PP-QUAD9a

MM-QUAD9

PM-QUAD9* (present)

a b c d e f

1.5000 1.2606 1.0372 1.2472 1.0372 1.2606

1.5000 1.7935 1.6398 1.5285 1.6398 1.7935

1.5000 1.5000 1.5000 1.5000 1.5000 1.5000

Exact solution = 1.5 a

Alias QUAD9.

Super-element mesh

11 22 44

Vertical displacement at the mid-point of loaded edge PP-QUAD9a

MM-QUAD9

PM-QUAD9* (present)

3.2979 3.3106 3.3123

3.2979 3.3106 3.3123

3.2979 3.3106 3.3123

Exact solution = 3.3125 a

Alias QUAD9.

Alias QUAD9.

fourth of the geometry is modelled. Two sequences of meshes with regular and distorted elements are considered as in Fig. 19. The ‘1  1 mesh’ of Fig. 19 is used as the super-element for construct-

ing the 2  2 and 4  4 meshes. The mesh-distortion in the sequence B is introduced by rotating the nodes surrounding the central node of the super-element by an angle a in the counterclockwise direction. The thickness of the cylinder is taken as unity

1060

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Table 18 Convergence results for 9-node elements with Mesh b as super-element. Super-element mesh

11 22 44

Vertical displacement at the mid-point of loaded edge PP-QUAD9a

MM-QUAD9

PM-QUAD9* (present)

3.0033 3.1768 3.2661

3.9559 3.9519 3.9598

3.3044 3.3057 3.3116

Exact solution = 3.3125 a

Alias QUAD9.

Table 19 Convergence results for 9-node elements with Mesh c as super-element. Super-element mesh

11 22 44

Vertical displacement at the mid-point of loaded edge PP-QUAD9a

MM-QUAD9

PM-QUAD9* (present)

0.4993 2.9218 3.1531

3.8294 3.9208 3.9660

3.2710 3.3001 3.3110

Exact solution = 3.3125 a

# 2 b ur ¼ 1  m þ 2 ð1 þ mÞ ; uh ¼ 0; 2 r Eðb  a2 Þ ! ! 2 2 a2 p b a2 p b rr ¼ 2 2 1  2 ; rh ¼ 2 2 1 þ 2 ; r r b a b a a2 pr

Alias QUAD9.

and the problem is solved in plane-stress condition. The exact solution to the problem is available in Ref. [29]:

"

ð25Þ

rrh ¼ 0:

For the present computations, the numerical values of the parameters are chosen as p = 1000, a = 1, b = 3. The material is assumed to be linear elastic with material properties E = 1000 and m = 0.25. Fig. 20 shows the convergence plots for this problem. As seen from Fig. 20a, for regular mesh (i.e., mesh sequence A), both elements exhibit almost the same order of convergence, close to the theoretical value of 2 expected of a quadratic element. However, it is seen that the solution error with the present element is larger than with the classical element (viz., PP-QUAD8). The poorer performance of the present element is believed to be due to incompatible (metric) shape functions used in the BR-matrix. The ill-effects of incompatibility seem to show up more seriously under non-polynomial displacement fields than under polynomial fields. The convergence plots for the distorted meshes are shown in Fig. 20b–d, corresponding to nodal rotation angle a = 30°, 45° and 60°, respectively. The convergence order for the present element is higher than that for PP-QUAD8 element for all the angles. How-

-1 -1.2

log10(Energy norm error)

-1.4 -1.6 -1.8 -2

R = 1.93

-2.2 PP-QUAD8 -2.4 -2.6 -0.35

PM -QUAD8 * (p res ent)

-0.15

0.05

0.25

log10(h)

(a) Mesh sequence A 1.2

0

0.7

R = 1.05

log10(Energy norm error)

log10(Energy norm error)

-0.5

-1

-1.5

R = 2.11

-2

-2.5 -0.35

0.2 -0.3

R = 0.85

-0.8 -1.3 -1.8

R = 1.97

PP-QUAD8 PM -QUAD8 * (p res ent)

-0.15

0.05 log10(h)

(b) Mesh sequence B

0.25

ð26Þ

PP-QUAD8

-2.3 -2.8 -0.35

PM -QUAD8 * (p res ent)

-0.15

0.05 log10(h)

(c) Mesh sequence C

Fig. 15. Convergence of energy norm strain error for shear loading of square panel (modelled with 9-node elements).

0.25

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

1061

6. Concluding remarks Three sufficient conditions have been proposed for a displacement-based finite element to be immune to mesh distortion effects even for non-positive Jacobian element geometries: 1. The shape functions used in BL matrix must have mesh-distortion independent inter-element compatibility property. 2. The shape functions used in BR matrix must satisfy mesh-distortion independent completeness property up to all the terms in the Cartesian displacement polynomial that we expect the element to reproduce exactly even under mesh distortion. 3. The Jacobian should not occur in the denominator of the stiffness integrand.

Fig. 16. rxx-stress contours given by PP-QUAD9 (alias QUAD9) as well as PMQUAD9* (present) for regular 4  4 mesh.

ever, the convergence order gradually decreases as the rotation angle increases. Thus quadratic order of convergence is not maintained for this problem. This example seems to suggest that the three conditions C1, C2 and C3 are not sufficient to guarantee quadratic convergence under distorted meshes with non-polynomial displacement fields.

The third condition is needed to avoid the possible instability of Gaussian numerical quadrature of element stiffness under extremely distorted element geometries. PM-QUAD8* and PM-QUAD9* elements have been presented as an illustrative examples satisfying all the three conditions above. Numerical test results show that PM-QUAD8* and PM-QUAD9* reproduce exact nodal displacements for quadratic displacement problems even for extreme mesh distortions for which the Jacobian goes negative. The element can reproduce exact nodal solution for quadratic displacement problem even for random positioning of

Fig. 17. rxx-stress contours for distorted meshes given by PM-QUAD9* (present).

Fig. 18. rxx-stress contours for distorted meshes given by PP-QUAD9 (alias QUAD9).

1062

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

Fig. 19. Super-element mesh sequence for thick cylinder problem.

0

0

-0.5

PP-QUAD8

-0.2

PM-QUAD8* (present)

-0.4

PP-QUAD8 PM-QUAD8* (present)

-0.6 -1 -0.8

-1.5

R=0.95

-1

R=2.1

-1.2 -2

-2.5 -0.35

-1.4

R=1.9 -0.15

0.05

0.25

-1.6

R=1.85

-1.8 -0.35

-0.15

log10(h)

0.25

(b) Distorted mesh with distortion angle, α = 30o

(a) Regular mesh 0

0 -0.2

0.05

log10(h)

PP-QUAD8

-0.2

PM-QUAD8* (present)

PP-QUAD8 PM-QUAD8* (present)

-0.4 -0.4 -0.6 -0.8

-0.6

R=0.94

R=0.98 -0.8

-1 -1.2 -1.4 -0.35

R=1.42 -0.15

-1

0.05

log10(h)

(c) Distorted mesh with distortion angle, α = 45o

0.25

-1.2 -0.35

R=1.21

-0.15

0.05

log10(h)

(d) Distorted mesh with distortion angle, α = 60o

Fig. 20. Convergence of energy norm strain error for thick cylinder problem.

0.25

S. Rajendran / Comput. Methods Appl. Mech. Engrg. 199 (2010) 1044–1063

interior nodes, which is the ultimate test for Jacobian-unlimited mesh-distortion immunity. As regards the numerical integration of the PM-QUAD8* element, exact integration is guaranteed with a sufficiently high order of quadrature. This is because the stiffness integral does not involve the Jacobian, and the integrand turns out to be a polynomial whatever be the extent of mesh distortion. The MM-QUAD8 element requires a higher order quadrature than PM-QUAD* as the stiffness integrand involves the Jacobian term in the numerator. The PP-QUAD8 element requires even higher order quadrature as the Jacobian appears in the denominator of the stiffness integral. For extremely high distortions involving negative Jacobian, the numerical integration for PP-QUAD8 becomes unstable leading to divergence of computed displacements with respect to increasing order of integration. For some choice of nodal positions, the possibility of interpolation matrix (Eqs. (11), (14), (15)) becoming singular cannot be ruled out in the case of metric shape functions. The singularity will lead to interpolation failure which is a kiss of death to the formulation. One way to handle this problem is to sense the impending interpolation failure and bypass it with a nodal perturbation technique or a coordinate transformation as suggested in reference [27]. Investigation into other techniques to handle interpolation failure is a potential area for future research. The technique to develop mesh-distortion immune elements presented in this paper is in its early stage. Further investigations may be necessary to establish its theoretical foundation. Nevertheless, the technique seems to open up exciting possibilities in designing elements with Jacobian-unlimited mesh-distortion immunity. Acknowledgements The research work presented in this paper was supported by the research Grant No. RG 29/08 (M52050104) of Academic Research Fund (AcRF) Tier 1, Ministry of Education, Singapore. References [1] I. Fried, Accuracy of complex finite elements, AIAA J. 10 (1972) 347–349. [2] R.D. Henshell, D. Walters, G.B. Warburton, A new family of curved elements for vibration and stability, J. Sound Vib. 20 (1972) 381–397. [3] I. Fried, Possible loss of accuracy in curved (iso-parametric) finite elements – comments on a paper by Henshell, J. Sound Vib. 23 (1972) 507–510. [4] R.D. Henshell, D. Walters, G.B. Warburton, On possible loss of accuracy in curved finite elements, J. Sound Vib. 23 (1972) 510–513. [5] J.A. Stricklin, W.S. Ho, E.Q. Richardson, W.E. Haisler, On iso-parametric vs. linear strain triangular elements, Int. J. Numer. Methods Engrg. 11 (1977) 1041–1043.

1063

[6] J. Bäcklund, On iso-parametric elements, Int. J. Numer. Methods Engrg. 12 (1978) 731–732. [7] L.N. Gifford, More on distorted iso-parametric elements, Int. J. Numer. Methods Engrg. 14 (1979) 290–291. [8] N.S. Lee, K.J. Bathe, Effects of element distortions on the performance of isoparametric elements, Int. J. Numer. Methods Engrg. 36 (1993) 3553–3576. [9] A.K. Soh, Y. Long, S. Cen, Development of eight-node quadrilateral membrane elements using the area coordinates method, Comput. Mech. 25 (2000) 376– 384. [10] X.M. Chen, S. Cen, Y.Q. Long, Z.H. Yao, Membrane elements insensitive to distortion using the quadrilateral area coordinate method, Comput. Struct. 82 (2004) 35–54. [11] Y.Q. Huang, Q.S. Li, Four-node incompatible plane and axisymmetric elements with quadratic completeness in the physical space, Int. J. Numer. Methods Engrg. 61 (2004) 1603–1624. [12] S. Cen, X.M. Chen, X.R. Fu, Quadrilateral membrane element family formulated by the quadrilateral area coordinate method, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4337–4353. [13] S. Rajendran, K.M. Liew, A novel unsymmetric 8-node plane element immune to mesh distortion under a quadratic displacement field, Int. J. Numer. Methods Engrg. 58 (2003) 1713–1748. [14] S. Rajendran, S. Subramanian, Distortion sensitivity of 8-node plane elasticity elements based on parametric, metric, parametric-metric and metricparametric formulations, Struct. Engrg. Mech. 17 (2004) 767–788. [15] E.T. Ooi, S. Rajendran, J.H. Yeo, A 20-node hexahedral element with enhanced distortion tolerance, Int. J. Numer. Methods Engrg. 60 (2004) 2501–2530. [16] K.M. Liew, S. Rajendran, W. Wang, A quadratic plane triangular element immune to quadratic mesh distortions under quadratic displacement fields, Comput. Methods Appl. Mech. Engrg. 195 (2006) 1207–1223. [17] S. Rajendran, E.T. Ooi, J.H. Yeo, Mesh-distortion immunity assessment of QUAD8 elements by strong-form patch tests, Commun. Numer. Methods Engrg. 23 (2007) 157–168. [18] G. Prathap, S. Mukherjee, Management-by-stress model of finite element computation, Research Report CM 0405, CSIR Centre for Mathematical Modelling and Computer Simulation, Bangalore, 2004. [19] G. Prathap, V. Senthilkumar, S. Manju, Mesh distortion immunity of finite elements and the best-fit paradigm, Sadhana 31 (2006) 505–514. [20] D.N. Arnold, D. Boffi, R.S. Falk, L. Gastaldi, Finite element approximation on quadrilateral meshes, Commun. Numer. Methods Engrg. 17 (2001) 805–812. [21] D.N. Arnold, D. Boffi, R.S. Falk, Approximation by quadrilateral finite elements, Math. Comput. 71 (2001) 909–922. [22] F. Kikuchi, M. Okabe, H. Fujio, Modification of the 8-node serendipity element, Comput. Methods Appl. Mech. Engrg. 179 (1999) 91–109. [23] J. Zhang, F. Kikuchi, Interpolation error estimates of a modified 8-node serendipity finite element, Numer. Math. 85 (2000) 503–524. [24] S. Rajendran, K.M. Liew, Completeness requirements of shape functions for higher order elements, Struct. Engrg. Mech. 10 (2000) 93–110. [25] T.J.R. Hughes, The Finite Element Method, vol. 78, Prentice-Hall, New Jersey, 1987. p. 4. [26] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, New York, 2000. pp. 148–150. [27] E.T. Ooi, S. Rajendran, J.H. Yeo, Remedies to rotational frame dependence and interpolation failures of US-QUAD8 element, Commun. Numer. Methods Engrg. 24 (2008) 1203–1217. [28] T.R. Chandruptla, A.D. Belegundu, Introduction to Finite Elements in Engineering, Prentice Hall, New Jersey, 2002. pp. 220–222. [29] R.J. Roark, W.C. Young, Formulas for Stress and Strain, McGraw-Hill, New York, 1975.