Mathematical Formulation and Numerical Implementation of a Finite Element with Anisotropic Geometry

Email: davide.tumino@unikore.it Abstract: A plane quadrangular element with geometric anisotropy has been developed to perform 2D Finite Elements Analyses in cases where high stress concentrations, varying with very different laws along two orthogonal directions, are present. The element has been implemented into a finite element code. To validate the element behavior, analyses in the adhesive of a single lap joint and in a bimaterial interface have been performed, comparing the stress fields obtained with those get from different methodologies (analytical, experimental and numerical with very fine meshes). It has been found that, using the same number of nodes, the analyses with the developed element returned better results with respect to the ones obtained with standard geometric isotropic elements.


Introduction
The problem of the stress concentration analysis in mechanical parts has been faced using different methodologies. From the first analytical approaches based on the Theory of Elasticity to the experimental and the recent numerical ones of the last decades, the analysis of the stress fields in parts characterized by geometric discontinuity or singularities dues to the matching of different materials has reached important results. In numerical methods based on Finite Elements (FE) it's a common practice to use mesh with high refinements in the areas where the stress concentrations are present (Zienkiewicz and Taylor, 1994). This leads to a good simulation of the problem but also to an increased computational effort, using the elements contained in the element libraries of the commercial FE codes (Reddy, 1993). The latter contain generally elements provided with geometric isotropy (Bathe, 1996), where the shape functions describing the displacement field within the element are polynomials of the same order in the three cartesian directions. Removing the condition of geometric isotropy it is possible to obtain elements where the shape functions describe the internal displacement field with laws of different form in the three cartesian directions. Therefore, in the cases where the quantities to estimate (displacements, strains, stresses, etc…) change with very different laws along the three directions, an element without geometric isotropy, appropriately oriented, is able to lead to reliable results also using a relative reduced number of elements. Then, it would be useful to have in the libraries of the FE programs also this kind of element to study problems with high stress concentrations varying with very different laws along the cartesian directions, with the purpose to reduce the storage resources and the computational time. Other approaches based on Lagrangian multipliers can be found in literature able to solve problems with anisotropic material with inextensible fibers (Wriggers et al., 2016).
In this study a quadrangular, plane, eight-noded finite element without geometric isotropy has been developed. The element has been implemented into a home made FE code, compiled in Fortran. With the aim to verify the efficiency of the element, some FE analyses of problems with high stress concentrations have been done and the results compared with the ones obtained using isotropic element with the same meshes and by other analytical or experimental approaches.

Analytical Formulation of the Anisotropic Element
Experimental, analytical and numerical approaches should be followed when investigation of analytical formulation of the anisotropic element is taken into account. However, experimental methods fail to give information about loading history and the changes introduced during the course of penetration process, while being well expensive. Due to their high reliability, theoretically and analytically methods are a simple and fast method for achieving the desired results (Hedayati and Vahedi, 2017;2018).
Following the typical FE approach in Zienkiewicz and Taylor (1994), displacements within the element are commonly expressed by the matrix relation: where, {s} is the displacement vector of a generic point inside the element, [N] is the shape function matrix, {q} is the displacement vector of the element nodes.
Using the generalized coordinates (Bathe, 1996) the most general expressions of the in-plane displacement of a generic point are the u and v components in the x and y directions, respectively: where, i are coefficients called generalized coordinates.
For a given element, Equation (2) have to be particularized so that the compatibility and completeness criteria, which ensure the convergence of the solution as the element number grows, are satisfied. The condition of geometric isotropy is fulfilled if the x and y terms are chosen symmetrically with respect to the axis from the Pascal's triangle in Fig. 1.

■■
In the case of the plane quadrangular eight-noded isotropic element in Fig. 2, called hereinafter ISO8, choosing from Pascal's triangle the terms circled in Fig. 1 As a consequence, the u and v displacements in (3) have the same analytical shape with respect to the x and y coordinates.
In the case of the plane quadrangular eight-noded anisotropic element in Fig. 3, called hereinafter ANISO8, displacements can be described by Equation (4) The terms chosen from Pascal's triangle are shown in Fig. 4. Now displacements are described by polynomials where the x coordinate appears with linear terms, while y coordinate appears with terms up to cubic.
From the expressions of the displacement functions (2), imposing displacement values at nodes {q}, the matrix of the shape functions can be obtained by: where, [A] is the matrix of the generalized coordinates and [C] is the matrix of the nodal coordinates. By associating the generic element to a square centred with respect to the origin of a local s and t reference system, with sides described by the segments s = 1 and t = 1 (Fig. 5), expressions of the shape functions, using natural s, t coordinates, can be obtained rewriting the element shape functions in generalized coordinates in s, t terms and performing the matrix operation (5). In this way, shape functions of easier usage are obtainable, because they are dependent not from the element geometry but only from the element type. Figure 6 shows, in the natural domain s, t, the shape functions associated, for both elements, to a corner node ((a) for ISO8 and (c) for ANISO8) and to an intermediate node ((b) for ISO8 and (d) for ANISO8). It can be seen how ANISO8 element is able to describe displacement field with different laws along two perpendicular directions. The element stiffnes matrix, [k], normally expressed in the global x, y coordinates, must be given in s and t terms. To compute it, the Gauss quadrature method is usually used (Zienkiewicz and Taylor, 1994), where an analytical integration is substituted with a numerical sum. It can be written: where, G(s,t) is a polynomial whose order depends on the number of degrees of freedom of the element along the s and t directions. Following the numerical integration procedure, G(s,t) is computed in a specified number of points, called Gauss points, inside the natural domain, whose number is equal to ns  nt; the obtained values are multiplied by some coefficients H named weights and then summation is executed. The number of Gauss points necessary to the exact integration of (6) is established by the order of G(s,t). Moreover, while with isotropic element is ns = nt, with anisotropic element ns and nt depend on the order of the coordinates s and t in G(s,t), respectively. With n Gauss points along a fixed direction is possible to achieve the exact integration of a 2n-1 order polynomial in the corresponding variable. For the eightnoded isotropic element 33 = 9 Gauss points are necessary, while for the eight-noded anisotropic element 24 = 8 Gauss points are enough. In Fig. 7 the positions of the Gauss points for the ANISO8 element are shown and in Table 1 the associated coordinates and weights.

Applications
To evaluate the reliability of the element formulation, the developed procedure has been used to estimate the stress field in two problems characterized by high stress concentrations and the results have been compared with the ones obtained by other ways.

Analysis of a Single Lap Joint
The determination of the stress field in a single lap joint has been the object of several studies with different methodologies. Preliminary studies can be found in Hovgaard (1930) and later in Volkersen (1938) and in Goland and Reissner (1944). The attention was focused on the determination of the stress field in the adhesive film, according with the assumption that adhesive failure is the greater reason to the joint breakdown. Numerical studies begin with Wooley and Carver (1971) where the value of the stress concentration factor was estimated for different joint geometries. At the same time, Adams and Peppiatt (1974) appraised the reliability of analytical results performing numerical and experimental analyses. Among the later works is remarkable Andruet et al. (2000), where the adherends mesh has been obtained with beam elements, while the adhesive mesh with quadrangular four-noded elements.
In the present paper, analyses of the joint depicted in Fig. 8 are performed with different methods: Analytical using the approach in Goland and Reissner (1944), numerical with the code ANSYS APDL using eightnoded quadrangular isotropic Element (PLANE82) with a very fine mesh and numerical with the home made FE software using the ISO8 and ANISO8 elements.
The analyses have been performed under these hypotheses:  Plane strain problem  Homogeneous and isotropic materials  Linear-elastic behaviour  In the analytical formulation x is the only strain component unequal to zero in the adherends The elastic characteristics of materials are reported in Table 2.
To refine the mesh in the regions where there is the higher stress concentration, just the part of the joint interested by the overlap has been modelled, replacing the adherends upstream and downstream the overlap with the normal force T and the bending moment M0 in the sections a and b (Fig. 9a). The bending moment M0 = 36.8 Nmm is computed according to the formulation in Goland and Reissner (1944). The mesh created in the Fortran software is shown in Fig. 9b together with the orientation in the x,y plane of the anisotropic elements. With the same FE model, analyses have been performed both with ISO8 element and with ANISO8 element, using for both cases 288 elements. A finer mesh has been employed in ANSYS using 1956 elements, in order to obtain a solution as closer as possible to the exact one. In Table 3       Near x = 0 (Fig. 10a) the curve related to y obtained by ANSYS is slightly different from the one obtained by the analytical approach owing to the simplified hypotheses about x introduced by Goland and Reissner; the values calculated by means of the ANISO8 element stay in the middle between the two curves mentioned above and furnish an error referred to the maximum value smaller than the one obtained by the ISO8 element with respect to ANSYS (3.3% versus 9.9%, respectively)  The shear stress (Fig. 10b) obtained by means of ANISO8 are almost equal to the results obtained by ANSYS (error referred to the minimum value close to 0.5% versus 13.6% related to the analysis using ISO8); this behaviour is appreciable also near the free surface, where isotropic element shows greater errors; according to analytical formulation (Goland and Reissner, 1944) in x = 0 xy should have a minimum, while, by means of numerical analyses, results practically xy = 0, because of the principle of reciprocity, the external surface being free from loads  The error related to the maximum value of x ( Fig.   10c) estimable from the results obtained by the anisotropic element, with respect to ANSYS, is equal to 1.35%, while the one given using the isotropic element is equal to 11.8% Therefore clearly appears that the anisotropic element gives better results than the ones given by the isotropic element ISO8 with the same number of elements and really close to the ones obtained by ANSYS using very fine meshes having six times the total number of elements.

Analysis of a Bimaterial Joint
The use of joints between two or more different materials, in the last decades, has begun a very common practice in several applications. In this study, the case of a joint between two plates made with isotropic and homogeneous materials and with rectilinear interface is considered. The stress field near the interface tip caused by mechanical loads (either static or dynamic) or thermal loads strictly depends on the materials of the joint (Kelly et al., 1992). The more different are the materials, the more the stress concentration increases near the interface tip. This often brings the appearance of cracks or plasticisations. Thus, it's evident how the exact knowledge of the stresses field plays a relevant role in the design of such a joint. The main analytical studies on the topic have been reported in Bogy (1968;1971), in the case of plane problems. Relating to this study, two parameters  and  have been introduced in Dundurs (1969) depending on the elastic characteristics of both materials, denoted with the subscripts 1 and 2 and defined by: and where G e v are respectively the shear modulus and the Poisson's ratio. Dundurs supposes that the stress field pattern near the interface tip depends just on the values assumed by  and  (Kelly et al., 1992): It is possible to have singular stress fields (with logarithmic or hyperbolic pattern) or non-singular.
In the most diffused case of hyperbolic singularity, the order of singularity derives only from the elastic characteristics of the joined materials. Later studies (Yang and Munz, 1997;Cirello and Zuccarello, 2001) have shown that it is possible to completely describe the hyperbolic stress field in a polar coordinate system (Fig. 11), in a restricted area around the origin, through the relation: where,  is the order of singularity, K is the stress intensity factor connected with the elastic characteristics of materials, with the joint geometry and with the load system, fij() are functions derived from the characteristics of the two materials. The original theory formulated by Bogy allows to calculate the order of singularity by solving the implicit relation: Where: Once  has been computed, the functions fij() can be obtained considering that they are the eigenfunctions associated to the eigenvalue . However, analytical formulation permits to obtain just the order of singularity  but not the value of the factor K, for which it is necessary to accomplish also numerical or experimental analyses. It has to be reminded that both numerical and experimental analyses are often difficult to perform. In fact the high stress concentration in a very restricted area, near the singularity point, needs, in numerical analyses, extremely refined meshes and the results can be affected by the type and the shape of the element used. Experimental approach presents difficulties due to the random errors associated to this kind of analysis, not only about the values obtained by the investigation but also about the exact measurement of the polar coordinates (Cirello and Zuccarello, 2001). However, to study the zone regulated by singularity, experimental analysis remain preferable, allowing to consider the eventual effects of hardening and/or plasticisations, the latter difficult to verify and simulate with numerical analyses. The bimaterial joint in Fig. 12 has been studied. The hypotheses of the analysis are:  Plane stress problem  Isotropic and homogeneous materials (aluminium and plexiglass)  Linear-elastic behaviour  Rectilinear interface  The characteristics of materials are reported in Table 4. From Equation (7)

 
Via Equation (10), the theoretical value of the order of singularity, Bogy, is computed: 0.2277.

 
To obtain the order of singularity  via numerical or experimental analyses, Equation (9) has to be rewritten in logarithmic form: To validate the results obtained analytically, these have been compared with the results of the experimental studies performed in Cirello and Zuccarello (2001) about the same couple of materials and the same geometric and loads configuration. These studies have been achieved with the automatic photoelasticity technique using a circular dark field polariscope in white light and a digital camera. Once the images of isochromatic fringes have been digitalized, the retardation  has been computed along a segment oriented in the direction  = 45°. In a bilogarithmic reference system (Fig. 13) the interpolation line of the results is traced; the slope of this line is equal to the order of singularity exp = -0.2273 and is very close to the theoretical value. From Fig. 13, a deviation of the computed points from the interpolation line describing the hyperbolic behaviour predicted by Bogy is noted. This aspect can derive, near the interface tip, from the presence of plastic strains, while, far away of the tip, from the loss of validity of the hyperbolic law.
In the present paper, numerical analyses have been performed using the home made FE software and also ANSYS (with very fine mesh) to validate the results. The model used in the home made software is shown in Fig. 14 considering the symmetry of the problem with respect to the line x = 37.5 mm. Four increasing mesh levels have been used, both with isotropic elements and with anisotropic elements (Table 5).
In the figure the orientation of the anisotropic element on the x,y plane is also depicted: This allows to have more nodes along the direction parallel to the interface.
Doing as in Fig. 13, the r stress along the free edge of aluminium ( = -90°) has been plotted in Fig. 15. To evaluate  the stress values have been selected in the field 0 < log(r) < 0.8, that is 1 mm < r < 6.3 mm, in fact as r0 a scattering of the points is noted (mostly in the analyses using isotropic element), while beyond r = 6.3 mm Bogy's formulation is not valid. In the field r < 6.3 the results obtained by ANSYS, not shown, are in perfect agreement with the theoretical formulation. Figure 15, the greater scattering of points related to the isotropic element can be noticed in particular as r0, while, points computed with the anisotropic element, settle on the interpolation line since the first node (log(r) =-0.61). This different behaviour is more appreciable using coarse meshes. What has been remarked results evident observing Table 5, where the values computed for  with the different refining levels are reported. Since mesh D2 the analyses performed using the anisotropic elements return a value for  equal to the analytical one, while using the isotropic element the mesh D4 is necessary to obtain a fair value. Approximately it can be noticed that to obtain the theoretical singularity order, it has to be used a mesh with twice isotropic elements than anisotropic elements.

Conclusion
A plane quadrangular eight-noded element without geometric isotropy has been developed to perform 2D FE analyses of problems with high stress concentrations characterized by very different laws along two orthogonal directions. The aim is to reduce the total number of elements and the computational time, obtaining the same accuracy in results as with isotropic elements.
Benchmark analyses are performed in some typical stress concentration cases, to compare results given by models meshed with the same number of isotropic and anisotropic elements. Results show that, to reach the same accuracy, a higher number of isotropic elements is in effect needed with respect to anisotropic ones having the same d.o.f. per element.
Moreover, the anisotropic element has also the advantage that needs one less integration point than isotropic element to perform the exact computation of the stiffness matrix; this means that an iteration less of the computing routines is necessary for each element; that is a computational time reduction around 10% with respect to an equivalent analysis using isotropic elements is reachable.
The advantages, in terms of resource saving, employing this kind of element become even more evident in the cases of iterative analyses (optimisation or non-linear problems), where the solution is obtained by several subsequent steps, or in cases of 3D and orthotropic material structures.

Ethics
This article is original and contains unpublished material. The corresponding author confirms that all of