Interpolation functions of a hybrid-Trefftz perforated super-element featuring nodes on the hole boundary

Finite Elements in Analysis and Design 91 (2014) 40–47

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Interpolation functions of a hybrid-Trefftz perforated super-element featuring nodes on the hole boundary C. Hennuyer a,b,n, N. Leconte b, B. Langrand a, E. Markiewicz b a b

ONERA – The French Aerospace Lab, F-59045 Lille, France LAMIH, UMR CNRS 8201, UVHC, Le Mont Houy, F-59313 Valenciennes Cedex 9, France

art ic l e i nf o

a b s t r a c t

Article history: Received 15 January 2014 Received in revised form 29 April 2014 Accepted 12 July 2014

The topic of the paper is related to the macroscopic modelling of riveted assemblies for full-scale aircraft structures subjected to crash/impact loadings. “Super-elements” have been developed to model the rivet and its rupture. However, stress concentrations in riveted assemblies can also initiate cracks from hole boundaries that propagate to the next ones and can lead to the catastrophic loss of the airplane by dismantlement of structural parts. An hybrid-Trefftz perforated super-element has been previously developed by the authors. However, the hole boundary of this super-element is a traction-free analytical boundary, preventing any connection with rivet elements. As a necessary first step to make the connection with the rivet possible, it is proposed to formulate additional nodes on the hole boundary of the perforated super-element. However, increasing the number of nodes implies increasing the superelement interpolation functions series order. Their ability to accurately describe the strain and stress fields is evaluated. It is found that when the series order is increased, the mechanical fields are still well described, but are a little less accurate than lower order interpolation functions. This evolution is discussed and seems to be linked to the ill-conditioning of the system. The accuracy is however still sufficient for the formulation of a 16-node super-element featuring eight nodes on the hole boundary. & 2014 Elsevier B.V. All rights reserved.

Keywords: Riveted assemblies Hybrid-Trefftz FE Interpolation functions

1. Introduction Finite element analysis (FEA) of airframe under high velocity impacts hardly succeeds in representing the failure of the structure when it occurs in riveted joint areas. Computational and experimental results were compared for different types of impact loadings [1,2]. The analysis shows that the macroscopic plastic strains are not sufficiently localised within the shell finite elements (without holes), to which rivets macro-finite elements are connected [1,3], so as to initiate and propagate the failure along rivet lines. In fact, the structural embrittlement, caused by holes in the sheets due to the riveting process, is not introduced in the standard shell FE formulation that is used for structural computations [4]. Modelling the geometrical defects (holes) with a really fine mesh is not suitable for full-scale aircraft crashworthiness, as an aircraft can feature more than hundreds of thousands of riveted assemblies [5]. Therefore, there is a rising need for alternative FE methods that would facilitate structural modelling and diminish

n Corresponding author at: ONERA – The French Aerospace Lab, F-59045 Lille, France. Tel.: þ33 3 20 49 69 00. E-mail address: [email protected] (C. Hennuyer).

http://dx.doi.org/10.1016/j.finel.2014.07.008 0168-874X/& 2014 Elsevier B.V. All rights reserved.

the computation cost by enabling one to build “super-elements” that model geometrical defects. Special finite elements have been developed to take into account the influence of geometrical defects (hole, crack) or inclusions (rigid, elastic). Tong et al. [6] were the first to formulate a “super-element” for the analysis of elastic stress intensity factors at a crack-tip. In the sense of Tong et al., a “super-element” is defined as a finite element that is able to represent complex mechanical fields with few degrees of freedom, contrary to conventional linear or quadratic finite elements. Following Tong et al., Piltner [7] formulated a super-element featuring an elliptical hole. One should note that Piltner's formulation is somehow general, because it can either model an elliptical hole or an internal crack. Leconte et al. [8] proposed a hybrid-Trefftz superelement with an internal circular hole based on the Piltner's formulation. Qin et al. [9] proposed an extension of Piltner's element in order to consider arbitrarily oriented elliptical holes. As an alternative to Piltner's formulation, Wang and Qin [10] developed special elliptical-hole elements based on the fundamental solutions [11]. Zhang et al. [12–15] developed n-sided polygonal elements featuring inclusions (rigid or elastic) or elliptical holes for the micromechanical analysis of heterogeneous materials. Their works are based on the Voronoï cell finite element method (VCFEM) proposed by Ghosh et al. [16,17] where the

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47

heterogeneous domain is modelled by a network of Voronï cells depending on the location, shape and size of the second phase. However, according to Zhang et al. the technique of transformation strain used by Ghosh et al. to take into account the effect of the second phase on the mechanical response of the material is not enough predictive [18]. In Refs. [6–15], the interpolation functions used to take into account the influence of such discontinuities (crack, hole, or inclusions) are based on Kolosov–Muskhelishvili's formalism [19]. It allows obtaining an analytical expression for displacements and stresses according to two complex potentials. In Refs. [10,11], the displacement and stress fields are expressed as a linear combination of fundamental solutions at a number of source points thanks to this formalism. In Refs. [6– 8,12–15], the two complex functions are expanded in Laurent series by taking into account the boundary conditions applied to the hole or the inclusion boundary (traction-free analytical boundary). In the case of a super-element featuring an inclusion, a displacement continuity condition is introduced at the interface between the matrix and the inclusion. The paper focuses on a particular type of finite elements: the hybrid-Trefftz displacement finite elements [20]. These elements are based on (i) a hybrid-Trefftz displacement variational principle, which enables the compatibility with conventional finite elements via an inter-element boundary displacement field, and (ii) interpolation functions satisfying a priori the governing equations of the problem, which enables to reduce the integrations to the element boundaries only. Most of the above-mentioned elements featuring geometrical discontinuities are hybrid-Trefftz displacement finite elements because, with this approach, complex mechanical fields near the defect can be represented without the need to mesh finely the area. The hybrid-Trefftz approach is not only efficient for the stress concentration analysis near a defect, but also for other various domains. The more recent works on hybrid-Trefftz finite elements concern the finite element analysis of Helmholtz problems [21], the analysis of axisymmetric potential problems [22], the elastostatic analysis of biphasic media [23], or even the analysis of thin plates [24]. Some of the authors of this paper developed a hybrid-Trefftz super-element featuring a hole to improve the riveted joints modelling [8]. The aim of the research presented here is to enable interactions between the rivet macro-element and the perforated super-elements. In its current version, the hole boundary is a traction-free analytical boundary. Consequently, loadings cannot be exchanged between the rivet macro-model and the perforated plane element. In order to establish a permanent dialogue between these models, it is proposed to formulate nodes on the hole boundary in order to allow the application of interaction forces between the rivet nodes and the ones of the perforation. The increase in the number of nodes of the perforated superelement implies an increase in the number of series terms of the interpolation functions (i.e., an increase in the order of the Laurent series). Moreover, the literature survey shows that the evolution of the interpolation functions accuracy depending on the number of series terms remains difficult to predict [25–27], and it is necessary to study it. Therefore, the capability of interpolation functions to calculate the displacement and stress fields inside the perforated super-element when it features additional nodes placed on the hole boundary, in particular, is studied in the paper. Given that a higher order has to be used in the Laurent series, close attention is paid to evaluate the accuracy of the interpolation functions. The paper is organised as follows. First, the formulation of the perforated super-element together with the expression of the associated interpolation functions is briefly reminded. Then, the numerical study of the interpolation functions of the superelement is performed for different orders of the Laurent series.

41

Finally, the main features of the interpolation functions are discussed.

2. Common rules for the formulation of hybrid-Trefftz finite elements featuring hole 2.1. Variational principle The perforated super-element [8] is based on a hybrid-Trefftz displacement variational principle [28]. It derives from the total potential energy principle [29]. A potential is added to ensure the displacement continuity between the super-element and the neighbouring conventional finite elements [28]. By applying the divergence theorem [30], and, by supposing that the interpolation functions a priori satisfy the equilibrium equations (Trefftz method) [7], the interior domain integrations are reduced to boundary integrations (1) Z Z Z 1 Π HT ðui ; T i ; u~ i Þ ¼ T i ui dS  Tb i ui dS þ ðu~ i ui ÞT i dS ð1Þ 2 S St S u [ Si where u is the vector of displacement field, T is the traction vector b is the prescribed traction vector, and, u~ is the vector of (T ¼ σ :n), T independently assumed boundary displacement. The fields u and T are related by the governing equations of elasticity (strain compatibility equations, linear elastic material behaviour) [8]. The boundaries Su and St represent the prescribed displacement boundary and the prescribed traction boundary, respectively. The boundary Si is the inter-element boundary, i.e., the boundary where the displacement continuity is ensured between the hybrid-Trefftz super-element and the neighbouring conventional finite elements. The boundary S is the disjoint union of the three above-mentioned boundaries S ¼ Su [ St [ Si , S being the external boundary of the super-element, and the perforation being a traction-free analytical boundary of the super-element. The fields u and T are expressed in terms of special interpolation functions and generalised degrees of freedom (see Section 2.2), and the inter-element displacement field u~ is approximated by means of conventional shape functions (linear or quadratic) and nodal degrees of freedom. The stiffness matrix and the load vector of the super-element are finally obtained by introducing these expressions in Eq. (1) and by taking the stationary conditions of the variational principle [8]. Note that each hole of a FE problem is discretised by a single super-element, which can be surrounded by conventional finite elements. 2.2. Interpolation functions The interpolation functions of this super-element comes from Kolosov–Muskhelishvili's formalism [19] to take account of the influence of the hole within the element. The displacements and stresses for a load-free perforated plate featuring a circular hole of radius a are expressed in the form of a Laurent series which depend on ðαj ; βj Þ parameters (2)–(6). These fields are known as a homogeneous solution of the considered problem uhx ¼

1 M j ∑ α ½ðkR þR  j Þ cos ðjθÞ  jðRj  Rj  2 Þ cos ðj  2Þθ 2μ j ¼  N j þ β j ½  ðkR þ R  j Þ sin ðjθÞ þ jðRj  Rj  2 Þ sin ðj 2Þθ j

uhy ¼

ð2Þ

1 M j ∑ α ½ðkR þR  j Þ sin ðjθÞ þ jðRj  Rj  2 Þ sin ðj 2Þθ 2μ j ¼  N j þ β j ½ðkR þ R  j Þ cos ðjθÞ þ jðRj  Rj  2 Þ cos ðj 2Þθ j

ð3Þ

42

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47

dofs are placed (in particular along the internal or external boundaries).

1 M ∑ αj j½2Rj  1 cos ðj  1Þθ R  j  1 cos ðj þ 1Þθ a j ¼ N

σ hxx ¼

 ½ðj  1ÞRj  1 ðj  2ÞRj  3  cos ðj 3Þθ þ βj j½  2Rj  1 sin ðj  1Þθ þ R  j  1 sin ðjþ 1Þθ þ ½ðj  1ÞRj  1 ðj  2ÞRj  3  sin ðj 3Þθ

σ hyy ¼

3. Study of the interpolation functions ð4Þ

M

1 ∑ α j½2Rj  1 cos ðj  1Þθ þ R  j  1 cos ðj þ 1Þθ a j ¼ N j

þ ½ðj  1ÞRj  1 ðj  2ÞRj  3  cos ðj 3Þθ þ βj j½  2Rj  1 sin ðj  1Þθ  R  j  1 sin ðjþ 1Þθ  ½ðj  1ÞRj  1 ðj  2ÞRj  3  sin ðj 3Þθ

σ hxy ¼

ð5Þ

1 M ∑ α j½ R  j  1 sin ðj þ 1Þθ a j ¼ N j

3.1. Methodology

þ ½ðj  1ÞRj  1 ðj  2ÞRj  3  sin ðj 3Þθ þ βj j½  R  j  1 cos ðjþ 1Þθ þ ½ðj  1ÞR

j1

j3

ðj  2ÞR

 cos ðj 3Þθ

ð6Þ

with ðR; θÞ being the polar coordinates (the radial coordinate is dimensionless R ¼r/a), μ is the shear modulus of the material, k is Muskhelishvili's constant (k ¼ ð3  νÞ=ð1 þ νÞ in plane stress, k ¼ 3 4ν in plane strain, where ν is Poisson's ratio), and N and M are the truncation orders of the series. Note that the superscript h denotes homogeneous solution. The fields uh and σ h can be written more concisely by introducing the matrices N and Q , which are the interpolation functions associated to displacements and stresses, respectively uh ¼ Nc

and

σ h ¼ Qc

The formulation of a hybrid-Trefftz perforated super-element has been reminded. It has been shown that increasing the number of nodes means increasing the order of truncation of the series appearing in the interpolation functions of the super-element. With the aim of formulating a hybrid-Trefftz perforated superelement featuring nodes on its hole boundary (Fig. 1(b)), it is proposed to study the interpolation functions associated with this super-element. The presented study assesses, more generally, these interpolation functions according to the series order.

ð7Þ

where c ¼ ðαj ; β j Þ represents the vector of the generalised degrees of freedom (dofs) of the super-element. For the sake of conciseness, “the perforated plate analytical solution calculated thanks to Kolosov–Muskhelishvili's formalism” will be simply denoted “K–M solution” in the forthcoming sections of the paper.

The methodology presented in this section has already been used to study the interpolation functions of an existing 8-node super-element in which series are truncated to N ¼ M ¼ 4 [8]. The purpose is to evaluate the capability of the interpolation functions, truncated to N ¼ M ¼ 8 in particular, to calculate the displacement and stress fields (inside the element) from input values of displacements extracted at selected points (in particular placed in the same way as the nodes of the studied element Fig. 1(b)). The objective is to study the influence of the order of the Laurent series, and the number and the position of the considered data points on the accuracy of the reinterpolated fields. The method consists in three steps once analytical reference solutions are set up for different elementary load cases (i.e., uniaxial or biaxial tension, simple or pure shear). Firstly, displacement data are collected on a reference solution. Secondly, a linear system (Ac ¼ b) can be written assuming that the displacements computed by K–M functions (Ac) should be equal to those collected on the reference solution (b). In the method, the linear system Ac ¼ b is built as follows:

 c is a row vector composed of the parameters ðαj ; βj Þ, which are the unknowns of the system,

2.3. Series truncation

 b is a row vector built from displacement data collected on the

The upper and lower bounds of the series appearing in Eqs. (2)–(6) are fixed in accordance with the number of nodes of the element thanks to relation (8) [26]

 A is a matrix constituted by the displacement equations, Eqs.

nu Znσ Z nq  r

ð8Þ

where nu and nσ represent the sufficient number of parameters to define displacements and stresses, respectively, nq is the number of physical dofs, and r is the number of rigid body modes. In our case, the displacement and stress fields in Eq. (7) being both expressed with the vector c, the relation between nu and nσ can be obtained such that nu ¼ nσ . The problem being bidimensional, nq is twice the number of nodes and the number of rigid body modes is r ¼3. Thanks to relation (8), the minimum number of required parameters is set relatively to the number of nodes of the formulated element. In the literature [7,26,28,31], it is recommended to take nσ approximatively equal to nq  r. This minimum number does not always guarantee a stiffness matrix with full rank. Full rank may be achieved by suitably increasing the number of parameters. However, too many parameters would make the element become too stiff [26,28,31]. The optimal value of parameters for a given type of element should be found by numerical experimentation [32,33]. For example, following these theoretical requirements, N and M are set to 4 for the formulation of a 8-node perforated element (Fig. 1(a)), and to 8 for a 16-node perforated element (Fig. 1(b)). Note that these rules do not vary wherever the

reference solution, (2) and (3).

The parameters ðαj ; β j Þ are finally identified after solving this system (c ¼ A  1 b). Consequently the K–M solution is completely defined: all fields are calculated by introducing the identified parameters into the previous analytical expressions (2)–(6). Finally, the computed solution is compared with the reference solution, by analysing the relative error between reference and computed displacement and stress components ((ux, uy) and (σxx, σyy, σxy), respectively). 3.2. Parameters of the numerical study A perforated plate featuring a centred circular hole of radius a¼ 0.2 whose material properties are E¼ 74,000 MPa and ν ¼0.3 is considered throughout the analysis. The value of Muskhelishvili's constant k and the shear modulus μ, required for the calculation of K–M solution, are k¼ 2.0769 (plane stress state) and μ ¼ 28; 462 MPa. A previous study showed that the distance between the perforation center and the data points (i.e. the parameter L appearing in Fig. 2) has no influence on the K–M solution [8]. The value of parameter L can be arbitrary fixed. In the following, L¼ 1 is considered. Finally, the studied parameters are:

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47

43

Fig. 1. (a) A 8-node perforated super-element and (b) a 16-node perforated super-element featuring nodes on the hole boundary.

Fig. 2. Configuration of data points collected for the calculation of the K–M solution. The data points are represented by the symbol ♦: (a) 8 points on the square, (b) 16 points on the square, (c) 8 points on the square and 8 on the perforation.

 The load case: uniaxial tension, biaxial tension, pure shear,   

simple shear. The truncation order of the series in the K–M solution (2)–(6): N and M, with N ¼M. The number of data points (where displacements are collected on the reference solution). The location of the data points (Fig. 2).

We have chosen to study two truncation orders of the interpolation functions: N ¼ M ¼ 4 and N ¼ M ¼ 8. The order 4 means that the K–M solution is driven by 16 parameters ðαj ; β j Þ, while 32 parameters ðαj ; β j Þ are available when the order is equal to 8. Because Eq. (8) set up a relationship between the number of parameters and the number of element nodes, it is at least necessary to collect data in eight points for N ¼ M ¼ 4, and 16

points for N ¼ M ¼ 8. It is proposed to analyse the three configurations presented in Fig. 2. They feature various data positions and numbers: either all the data are collected far from the perforation (Fig. 2(a) and (b)) or the data are collected on the perforation and away from the perforation (on a square) as shown in Fig. 2(c). In the following, each configuration will be represented by a triplet fns =np =no g where ns is the number of points collected on the square, np is the number of points collected on the perforation, and no is the interpolation functions' order. For example, if the interpolation functions are truncated at order 4 and if the data points are set up as shown in Fig. 2(c), then the configuration is denoted by f8=8=4g. It can be noticed that in the case of an overdetermined system (i.e., when data are collected on 16 points and the order is 4) a resolution via the least square method is performed.

44

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47

Fig. 3. Comparison between the analytical reference solution of Kirsch and the K–M solution for the configuration f8=8=8g: (a) σ VM distribution computed from the K–M solution, (b) σ VM along the y-axis, (c) ux along the x-axis, and (d) σyy along the hole boundary.

3.3. Results The accuracy of the interpolation functions is studied thanks to a first reference test-case [8]. This test-case corresponds to the analytical solution of an elastic infinite plate featuring a circular hole subjected to a far-field uniaxial tensile load σ 1 which was established by Kirsch [34]. The displacements and stresses expressions are given in Eqs. (9). The equivalent von Mises stress field obtained for this test-case is presented in Fig. 3(a):    2    σ1 a a4 4a2 ð1  νÞr þ ð1 þ νÞ þ r  3 cos ð2θÞ þ cos ð2θÞ ur ðr; θÞ ¼ 2E r r r 

   σ1 a4 2a2 ð1 þ νÞ 1 þ 4 r þ ð1  νÞ sin ð2θÞ uθ ðr; θÞ ¼  2E r r     σ a2 σ 3a4 4a2 σ rr ðr; θÞ ¼ 1 1  2 þ 1 1 þ 4  2 cos ð2θÞ 2 2 r r r     2 4 σ a σ 3a σ θθ ðr; θÞ ¼ 1 1 þ 2  1 1 þ 4 cos ð2θÞ 2 2 r r   σ 3a4 2a2 σ rθ ðr; θÞ ¼  1 1  4 þ 2 sin ð2θÞ 2 r r

ð9Þ

The K–M solution built from the computed coefficients ðαj ; β j Þ is compared to the reference solution inside the whole perforated domain fðx; yÞ A ½  1; 12 g using the maximal percent error. The results are given in particular along three relevant lines: around the perforation, and along fx ¼ 0g and fy ¼ 0g. The maximal percent errors, inside the perforated domain ½  1; 12 , are reported in Table 1 for each displacement and stress component, and for all the studied configurations. The first line of Table 1 corresponds to

Table 1 Maximal percent error (%) observed between the K–M solution and the Kirsch solution for different configurations (order of magnitude). Configuration f8=0=4g f8=8=4g f16=0=4g f16=0=8g f8=8=8g

ux

σxx

uy  14

10 10  14 10  14 10  12 10  8

 13

10 10  13 10  13 10  12 10  7

 13

10 10  14 10  14 10  11 10  7

σyy  13

10 10  13 10  13 10  10 10  5

σxy 10  11 10  12 10  14 10  9 10  4

the interpolation functions of the 8-node super-element, series truncated at order 4 [8]. Firstly, one can notice that, whatever the studied configuration for series truncated at order 4, the accuracy does not decrease (Table 1). Increasing the number of data points (with some of them on the hole boundary) has no influence on the accuracy of interpolation functions when truncated at order 4 (see configurations fns =np =4g). Secondly, the analysis of results obtained at the order 8 shows a small decrease in the accuracy compared with the reference results. For example, the percent error between the reference and K–M solutions is about 10  9 % for the configuration f16=0=8g and the σxy component of stresses instead of 10  11 % as previously (Table 1). When the data are collected on the hole boundary, the maximal percent error between the reference and the K–M solution built with the coefficients associated to order 8 is more important (for example the percent error is about 10  4 % for the σxy component of stresses in Table 1). This configuration is the less accurate among all those studied. It seems to be counter intuitive to decrease the accuracy while the series order increases, but,

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47

increasing the order means, in fact, increasing the number of parameters. Most of them should be equal to zero (for simple load cases) but they are in fact close to zero (and not exactly). A little bias may be consequently introduced for each contribution added to the solution given by Eqs. (2)–(6), which generates discrepancies between the exact solution and the K–M solution. Nevertheless, the results at order 8 with data points on the hole boundary (configuration f8=8=8g) remain very close to the reference solution (Fig. 3); the maximal percent error being to the worst for the component σxy and is about 10  4 %. It has been verified that changing the position of the data points around the perforation has no influence on the accuracy of the K–M solution. Indeed, another analysis has been done with data points on the perforation rotated with an angle of π =8, for example, and results obtained at order 8 for the configuration f8=8=8g are similar to those obtained previously (Table 1). The reference solution for the other load cases (i.e., biaxial tension, simple or pure shear) are given by the K–M solution. Indeed, one can observe that only few parameters ðαj ; βj Þ have non-zero values for elementary loads. For example, the previous uniaxial tensile loading is defined by the only parameters α  1 and α1 (α  1 ¼ 10 and α1 ¼ 5 in Eq. (10) when σ 1 ¼ 100 MPaÞ:      1 k 1 1 þ R cos θ þ  3 cos 3θ uhx ðR; θÞ ¼ α  1 2μ R R R    1 2 cos θ ð10Þ þ α1 ðk  1ÞR þ 2μ R The parameters ðαj ; βj Þ are analysed using FE results, given by [8], for the other load cases. It is shown that, using the same method as in Section 3.1, the biaxial tension, the simple and pure shear are defined by parameters α1, (β  1, β1) and β  1, respectively. The analytical expression of x-displacement, for example, for the three above-mentioned load cases is thus given by Eqs. (11)–(13)    1 2 uhx ðR; θÞ ¼ α1 ðk 1ÞR þ cos θ ð11Þ 2μ R uhx ðR; θÞ ¼

uhx ðR; θÞ ¼



    k 1 1 þ R sin θ þ  3 sin 3θ R R R

1 β 2μ  1 1 þ β 1 ½  ðk þ 1ÞR sin θ 2μ 1 β 2μ  1

    k 1 1 þ R sin θ þ  3 sin 3θ R R R

4. Discussion The analysis presented above shows clearly the sensitivity of the accuracy of the interpolation functions depending on the number of parameters ðαj ; β j Þ. This problem has been encountered in the literature by Rezaiee et al. [33] who developed a 4-node hybrid stress membrane element with drilling degrees of freedom. The interpolation functions of this FE are expressed in the form of a series which depends on Airy stress functions and unknown coefficients, and, the highest accuracy is provided by the minimum number of coefficients satisfying relation (8). This subject had also been discussed by Dhanasekar et al. [25], and Piltner [26,27]. The two proposed solutions have been evaluated and compared at one or two positions near the hole for the calculation of hoop stress. Table 2 Maximal percent error (%) observed between the K–M solution and the new reference solutions for different configurations (order of magnitude). Test-case

Configuration

ux

σxx

uy

σyy

 14

 13

 14

 13

σxy

Biaxial tension

f8=0=4g f8=8=4g f16=0=4g f16=0=8g f8=8=8g

10 10  14 10  14 10  12 10  7

10 10  14 10  14 10  12 10  8

10 10  14 10  14 10  12 10  7

10 10  14 10  14 10  12 10  7

10  14 10  12 10  12 10  10 10  5

Simple shear

f8=0=4g f8=8=4g f16=0=4g f16=0=8g f8=8=8g

10  14 10  13 10  14 10  12 10  7

10  13 10  13 10  13 10  11 10  6

10  12 10  12 10  12 10  9 10  5

10  12 10  11 10  12 10  9 10  5

10  13 10  13 10  13 10  10 10  6

Pure shear

f8=0=4g f8=8=4g f16=0=4g f16=0=8g f8=8=8g

10  14 10  13 10  14 10  11 10  6

10  14 10  13 10  14 10  11 10  6

10  12 10  12 10  14 10  9 10  4

10  12 10  12 10  14 10  9 10  4

10  13 10  13 10  14 10  10 10  5

ð13Þ

Each reference solution in Fig. 4 is set-up with the parameters' value equal to:  α1 ¼ 17:64 for the biaxial tension,  ðβ  1 ; β1 Þ ¼ ð9:15;  4:20Þ for the simple shear,  β  1 ¼ 19:02 for the pure shear.

The methodology, already presented for the uniaxial tension test-case and performed for the three configurations in Fig. 2 and the two orders 4 and 8, is undertaken for these new test-cases. Results are presented in Table 2. The conclusions are similar compared to the previous uniaxial tensile test-case: whatever the test-case considered, the maximal percent error increases when the order is equal to 8 and particularly when data points are located on the hole boundary. The maximal percent error is observed for components σxx and σyy for the pure shear test-case. It should be noted that for every configuration all the significant parameters have values which match the reference values previously given for the test-cases. All the other parameters of the K– M solution have non-zero values but are very close to zero.

ð12Þ



45

Fig. 4. Equivalent von Mises stress σ VM (MPa) distribution obtained from analytical reference solutions: (a) biaxial tension, (b) simple shear, and (c) pure shear.

46

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47

Table 3 Information about the linear system Als c ¼ bls . Configuration

Size of Als

Size of c

Size of bls

condðAls Þ

f8=0=4g f8=8=4g f16=0=4g f16=0=8g f8=8=8g

16  16 16  16 16  16 32  32 32  32

16  1 16  1 16  1 32  1 32  1

16  1 16  1 16  1 32  1 32  1

1:5025e7 4:7238e6 8:7729e6 2:0136e14 3:9873e14

Piltner assures that there is no loss of accuracy with the increasing number of parameters ðαj ; βj Þ [26]. Besides, Piltner has presented more accurate results with a 24-node element than a 8-node element (cf. Tables 2–5 in Ref. [26]). But the conclusion is based on the results observed at two positions only. The study presented in the present paper endeavours to analyse each obtained solution in the whole domain ½ 1; 12 , and not only at few points. It should be noted that close to the hole, all configurations give good results too and in particular for the value of the hoop stress (the maximal percent error in Tables 1 and 2 are never observed at the hole boundary). The parameters ðαj ; β j Þ of the K–M solution are identified by solving the linear system: Ac ¼ b as presented in Section 3.1. The study of the condition number of the linear system can thus be relevant because it measures the sensitivity of the solution of a linear system to errors introduced in the matrix, A þ ΔA, or, in the second member vector, b þ Δb. A little bias introduced in the matrix, or in the second member vector could imply a large error in the identified parameters of the vector c. Given that a least square method is applied, the linear system built is Als c ¼ bls where Als ¼ A  1 A and bls ¼ A  1 b. Note that the matrix A is the same whatever the load case because it is completely defined by the coordinates of the data points and the order of the considered Laurent series, only vector b differs in each test-case. The value of the condition number of the matrix Als is given in Table 3 according to each configuration. It is shown that the identification of the parameters ðαj ; βj Þ is more sensitive to numerical inaccuracies at order 8 than at order 4. It can be also noticed that overdetermined systems may diminish a little the condition number of the matrix Als . Numerical inaccuracies of the interpolation functions, at order 8 in particular, can be explained by an ill-conditioned linear system to be solved. Numerical methods such as preconditioning methods can be used to improve the condition number of the linear system.

5. Conclusions The present paper deals with the formulation of a hybridTrefftz perforated super-element featuring nodes on the hole boundary (rivet/perforated plate interaction), and focuses especially on the interpolation functions of this new super-element. The formulation of the current version of the hybrid-Trefftz perforated super-element which features eight nodes (i.e., the series appearing in interpolation functions are truncated at order 4) is reminded. The current line of research consists in formulating eight additional nodes located on the hole boundary of the superelement. But, the increase in the number of nodes implies an increase in the series order appearing in the interpolation functions of this special FE. Therefore, it is proposed to evaluate the interpolation functions of a 16-node perforated super-element featuring eight nodes on its perforation, whose the series are now truncated at order N ¼ M ¼ 8. The numerical study presented in the paper shows that the interpolation functions provide a little

less accurate results when the series terms are computed for order 8 instead of the ones computed for order 4 (interpolation functions of the 8-node perforated super-element). The worst results are obtained for the configuration where data points are collected on the hole boundary (i.e., the configuration f8=8=8g), which is the most interesting configuration for our problem. Nevertheless, even if the solution computed from K–M functions remains still very close to the reference solution for all considered load cases, the paper highlights the fact that the obtained results are counterintuitive. Indeed, the accuracy of the computed solution decreases when increasing the series order, due to the increased sensitivity. This phenomenon of sensitivity has already been subject of discussion in the literature, and has also been observed in recent works for Airy stress functions, but it has never been established for K–M functions. The ill-conditioning of the system is suggested as being the origin of the sensitivity of these interpolation functions. Finally, these inaccuracies at order 8 are not significant, and they are not a priori an obstacle for the forthcoming developments. Therefore, the theoretical formulation of a 16-node perforated super-element featuring eight nodes on its perforation (i.e., the variational principle and the stiffness matrix) will be presented in a future paper. The next step will consist in modelling the interaction between the fastener model and the perforated super-element. In order to do so, the influence of forces applied on the hole boundary will be taken into account in the super-element formulation.

Acknowledgements The present research work has been supported by the International Campus on Safety and Intermodality in Transportation, the Nord-Pas-de-Calais Region, the European Community, the Regional Delegation for Research and Technology, the Ministry of Higher Education and Research, and the National Center for Scientific Research. The authors gratefully acknowledge the support of these institutions.

References [1] B. Langrand, E. Deletombe, E. Markiewicz, P. Drazétic, Riveted joint modelling for numerical analysis of airframe crashworthiness, Finite Elem. Anal. Des. 38 (2001) 21–44. [2] B. Langrand, A.S. Bayart, Y. Chauveau, E. Deletombe, Assessment of multiphysics FE methods for bird impact modelling—application to a metallic riveted airframe, Int. J. Crashworthiness 7 (4) (2002) 415–428. [3] A. Combescure, F. Delcroix, L. Caplain, S. Espanol, P. Eliot, A finite element to simulate the failure of weld points on impact, Int. J. Impact Eng. 28 (2003) 783–802. [4] B. Langrand, E. Deletombe, A model of structural embrittlement: from plastic location to failure—theoretical basis and experimental characterisation on perforated plates, in: Proceedings of the 12th International Symposium on Plasticity (PLASTICITY 2006), Halifax, Nova Scotia, Canada, 2006. [5] R. Ortiz, J.L. Charles, J.F. Sobry, Structural loading of a complete aircraft under realistic crash conditions: generation of a load database for passenger safety and innovative design, in: Proceedings of the 24th Congress of International Council of the Aeronautical Sciences, Yokohama, Japan, 2004, 〈http://www. icas.org/ICAS_ARCHIVE/ICAS2004/ABSTRACTS/063.HTM〉 (accessed 17.12.13). [6] P. Tong, T.H.H. Pian, S.J. Lasry, A hybrid-element approach to crack problems in plane elasticity, Int. J. Numer. Methods Eng. 7 (1973) 297–308. [7] R. Piltner, Special finite elements with holes and internal cracks, Int. J. Numer. Methods Eng. 21 (1985) 1471–1485. [8] N. Leconte, B. Langrand, E. Markiewicz, On some features of a plate hybridTrefftz displacement element containing a hole, Finite Elem. Anal. Des. 46 (2010) 819–828. [9] Q.H. Qin, X.Q. He, Special elliptic hole elements of Trefftz FEM in stress concentration analysis, J. Mech. MEMS 1 (2009) 335–348. [10] H. Wang, Q.H. Qin, A new special element for stress concentration analysis of a plate with elliptical holes, Acta Mech. 223 (6) (2012) 1323–1340. [11] H. Wang, Q.H. Qin, Fundamental-solution-based hybrid FEM for plane elasticity with special elements, Comput. Mech. 48 (5) (2011) 515–528.

C. Hennuyer et al. / Finite Elements in Analysis and Design 91 (2014) 40–47 [12] J. Zhang, N. Katsube, A hybrid finite element method for heterogeneous materials with randomly dispersed rigid inclusions, Int. J. Numer. Methods Eng. 38 (1995) 1635–1653. [13] J. Zhang, N. Katsube, A hybrid finite element method for heterogeneous materials with randomly dispersed elastic inclusions, Finite Elem. Anal. Des. 19 (1995) 45–55. [14] J. Zhang, N. Katsube, A polygonal element approach to random heterogeneous media with rigid ellipses or elliptical voids, Comput. Methods Appl. Mech. Eng. 148 (1997) 225–234. [15] J. Zhang, P. Dong, A hybrid polygonal element method for fracture mechanics analysis of resistance spot welds containing porosity, Eng. Fract. Mech. 59 (1998) 815–825. [16] S. Ghosh, S.N. Mukhopadhyay, A material based finite element analysis of heterogeneous media involving Dirichlet tessellations, Comput. Methods Appl. Mech. Eng. 104 (1993) 211–247. [17] S. Ghosh, R.L. Mallett, Voronoi cell finite elements, Comput. Struct. 50 (1994) 33–46. [18] J. Zhang, N. Katsube, Problems related to application of transformation strains in a finite element analysis, Int. J. Numer. Methods Eng. 37 (1994) 3185–3193. [19] N.I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, Noordhoff, Groningen, Holland, 1953. [20] J.A. Teixeira de Freitas, Hybrid finite element formulations for elastodynamic analysis in the frequency domain, Int. J. Solids Struct. 36 (1999) 1883–1923. [21] K.Y. Sze, Y.K. Cheung, A hybrid-Trefftz finite element model for Helmholtz problem, Commun. Numer. Methods Eng. 24 (2008) 2047–2060. [22] K. Wang, L. Zhang, P. Li, A four-node hybrid-Trefftz annular element for analysis of axisymmetric potential problems, Finite Elem. Anal. Des. 60 (2012) 49–56.

47

[23] I.D. Moldovan, T.D. Cao, J.A. Teixeira de Freitas, Hybrid-Trefftz finite elements for biphasic elastostatic, Finite Elem. Anal. Des. 66 (2013) 68–82. [24] M. Rezaiee-Pajand, M. Karkon, Two higher order hybrid-Trefftz elements for thin plate bending analysis, Finite Elem. Anal. Des. 85 (2014) 73–86. [25] M. Dhanasekar, J. Han, Q.H. Qin, A hybrid-Trefftz element containing an elliptic hole, Finite Elem. Anal. Des. 42 (2006) 1314–1323. [26] R. Piltner, Some remarks on finite elements with an elliptic hole, Finite Elem. Anal. Des. 44 (2008) 767–772. [27] R. Piltner, Comment on: “A hybrid-Trefftz element containing an elliptic hole”, by Manicka Dhanasekar, Jianjun Han and Qinghua Qin, Finite Elem. Anal. Des. 42 (2006) 1314–1323, Finite Elem. Anal. Des. 43 (2007) 975. [28] P. Tong, New displacement hybrid finite element models for solid continua, Int. J. Numer. Methods Eng. 2 (1970) 73–83. [29] N. Leconte, B. Langrand, E. Markiewicz, Hybrid-displacement FE formulations including a hole, Struct. Eng. Mech. 31 (2009) 439–451. [30] C. Felippa, Three-dimensional linear elastostatics, in: Advanced Finite Element Method, University of Colorado Courses, Spring, USA 2013, 〈http://www. colorado.edu/engineering/CAS/courses.d/AFEM.d/〉 (accessed 17.12.13). [31] T.H.H. Pian, D. Chen, On the suppression of zero energy deformation modes, Int. J. Numer. Methods Eng. 19 (1983) 1741–1752. [32] Q.H. Qin, The Trefftz Finite and Boundary Element Method, WIT Press, Southampton, 2000. [33] M. Rezaiee-Pajand, M. Karkon, An effective membrane element based on analytical solution, Eur. J. Mech.—A/Solids 39 (2013) 268–279. [34] G. Kirsch, Die theorie der elastizitat und die bedurfnisse der festigkeitslehre, Z VDI 42 (1898) 797–807.