A generalized finite element interface method for mesh reduction of composite materials simulations

Abstract

This paper proposes interface and polynomial enrichments using the generalized finite element method (IGFEM) for the material interface in composite materials without matching the finite element mesh to the boundaries of different materials. Applications in structural members such as laminated beams and heterogeneous composites (matrix and inclusions) employing coarse and fine meshes are employed. The results were compared with conventional GFEM and analytical solutions. Verification and simulations proved the efficiency of the suggested framework for solving problems with discontinuous gradients resulting from a material interface. The proposed method allows flexibility in mesh generation for composite materials by letting the interface be embedded in an element without the need to match the mesh to the material interface. This improves the computational efficiency over conventional methods.

Keywords
GFEM; Composite materials; Interface

Graphical Abstract

1 INTRODUCTION

The development of non-conventional formulations required to solve complex problems remains a central issue in computational mechanics. Composite materials, for example, that are formed by two or more stages with different mechanical properties, require modeling with the Finite Element Method (FEM), fits or discretization (mesh) to the contours of the various phases of the material. This can significantly increase modeling's computational cost, making micromechanical analysis and multi-scale techniques prohibitively expensive.

In recent years, various non-conventional methods for solving specific problems in computational mechanics have grown in importance and popularity. An et al. (2013) used the numeric manifold method to model cracked material interfaces. The weak discontinuity through the material interface is described by introducing two types of customized functions. Zhang and Wang (2015) presented an isogeometric enriched quasi-convex meshfree formulation with material interface modeling. Phan and Mukherjee (2009) presented a multi-domain boundary contour method for interfaces and different materials that offers additional dimensionality reduction. Ahmadi et al. (2010) presented a method for studying steady state heat conduction in anisotropic and heterogeneous materials that does not use a local mesh. Soghrati (2014), Soghrati and Ahmadian (2015), Soghrati and Barrera (2016) and Aragón et al. (2020) proposed a hierarchical interface-enriched finite element method (HIFEM) for simulating various interfaces of materials in proximity or contact while capturing the gradient discontinuity. Soghrati et al. (2017, 2018) and Nagarajan and Soghrati (2018) proposed a method called Conforming to Interface Structured Adaptive Mesh Refinement (CISAMR) to simulate interfaces and intersecting boundaries with complex morphologies.

Based on the works of Babuška et al. (1994), Melenk and Babuška (1996) and Babuška and Melenk (1997) another meshfree methods were suggested: Special Finite Elements Methods and the Partition of Unity Finite Element Method. A similar strategy was also suggested by Duarte and Oden (1995), Duarte and Oden (1996) and Oden et al. (1998), under the names Clouds hp and Finite Element Method based on Clouds hp. In the works of Duarte et al. (2000, 2001) and Strouboulis et al. (2000) an approximation named GFEM was proposed that combines enrichment functions from a priori knowledge of the boundary value problem with the FEM shape functions. The method uses interface enrichment and polynomial, singular, and discontinuous enrichment functions to improve the quality of the approximation and make the modeling more flexible. Applications of GFEM are found in Evangelista Jr and Moreira (2020), Wolenski et al. (2020), Gomes et al. (2021) and Campos et al. (2022). The papers Duarte and Kim (2008), Kim et al. (2008), Evangelista Jr et al. (2020), Monteiro et al. (2020) and Novelli et al. (2020) applied a GFEM strategy with global-local enrichment functions built numerically (GFEM^g-l). Another approach was proposed by Belytschko and Black (1999) and Moës et al. (1999) under the name Extended Finite Element Method (XFEM). Those methods add special functions in the domain of interest.

In the case of composite materials, the method can allow the generated mesh to not conform to the physical interface between materials, as it does in FEM, allowing for the modeling and simulation of complex, multi-scale, optimization, stochastic, and other problems at a lower computational cost. This is because the use of specific enrichment functions by this method does not require the generation of successive conforming meshes to achieve the desired approximation quality. Moës et al. (2003) implemented enrichment functions for the material interface problems so that the finite element mesh used does not have to conform to the material surface. These functions are discontinuous in the first derivative, which describes the interface's existing jump behavior. This method relies on the implicit description of the location of the interface plane, which cuts the finite element. This description is provided by the level-set function, which is a distance function mapping and interpolating the interface location. Simone et al. (2006) developed the formulation and implementation of GFEM using discontinuous enrichment functions. Furthermore, the method requires an elaborate map of elements cut by the planes of the interfaces, which can become a disadvantage in problems with a large number of degrees of freedom. The application of the essential boundary conditions in this methos also requires special treatment of the enriched nodes.

Soghrati et al. (2012) developed a method for solving problems with discontinuous gradient fields involving heat transfer problems for analysis in two-dimensional domains. The method was called Interface Enrichment with Generalized Finite Element Method (IGFEM). Enrichment functions that are associated with generalized degrees of freedom are created at the intersection of the phase interface with the edges of the element. These functions were built from the linear combination of Lagrangian functions in the integration element. Among these advantages can be mentioned: the lowest computational cost, easy implementation, and simple handling of the Dirichlet boundary conditions. In another study conducted by Soghrati and Geubelle (2012), the implementation has been extended to a three-dimensional domain. Both works, however, use the IGFEM to solve field problems, specifically heat transfer.

In order to bridge the gap between exact and discretized geometries of curved material interfaces, Soghrati et al. (2015) present the IGFEM. For the 3D computational design of actively cooled panels, Tan and Geubelle (2017) present a reduced thermal model and a gradient-based shape optimization technique. Malekan and Barros (2018) investigated the behavior of micro-defects and inhomogeneities around a main crack through the global-local strategy of the GFEM. The same strategy was applied by Monteiro et al. (2018), who studied the behavior of a material (inclusion) within a semi-infinite domain of another material (matrix).

The primary goal of this paper is to extend the IGFEM method to vector quantities for the analysis of computational mechanics problems involving composite materials. This paper suggests and develops a method for the analysis of structures in the linear elastic regime by incorporating polynomial and interface enrichment functions. The performance, effectiveness, and convergence of the proposed method are investigated for both mesh refinement and polynomial degree enrichment. The results show that the proposed framework makes the mesh generation process more flexible for composite materials since the method can efficiently embed the interface between two materials and eliminate the need for a conforming mesh. The proposed method also has the advantage of only assigning enrichment fictitious nodes at interfaces, which resolves a common problem with existing methods that involves assigning boundary conditions and mapping boundary continuity through the elements.

2 FORMULATION OF THE PROPOSED METHOD

The proposed method addresses the concept of influence clouds (support) between the elements in the same way as Duarte and Oden (1996) hp-Cloud method, which uses clouds defined by their radius to discretize the problem domain under analysis. The clouds in the simulated method, on the other hand, are formed by a set of finite elements that share a node (Figure 1). For each cloud, the partition of unity (POU) is determined by the set of shape functions of the finite element that comprise the cloud associated with the base node. To improve the quality of the approximation, multiply the POU functions already established by linearly independent functions at each node cloud, also known as local approximation functions or enrichment functions:

${\mathfrak{I}}_j=\left\{L_{j1}\left(\mathit{x}\right),\dots ,L_{ji}\left(\mathit{x}\right),\dots ,L_{jq}\left(\mathit{x}\right)\right\}$ with $\left\{L_{j1}\left(\mathit{x}\right)=1\right\}$ (1)

in which $L_{ji}$ are the enrichment functions defined in each node ${\mathit{x}}_j$ cloud $w_j$.

Figure 1
Construction of the proposed enrichment method's shape functions (a) continuous; (b) high-order discontinuous (Evangelista Jr. et al., 2020).

The product of POU by the functions of equation (1) results in the so-called product function or shape function of the proposed method, shown in equation (2):

${\left\{{\varphi }_{ji}\left(\mathit{x}\right)\right\}}_{i=1}^q=N_j\left(\mathit{x}\right){\left\{L_{ji}\left(\mathit{x}\right)\right\}}_{i=1}^q$ (2)

in which $N_j$ denotes the functions POU, also known as shape functions of the FEM and ${\varphi }_{ji}$ are the shape functions of the GFEM, clearly presented in Figure 1 for a 2D domain. The POU functions defined by four elements sharing a node, ${\mathit{x}}_j$, are shown at the top of Figure 1. The enrichment functions to characterize the local problem are defined in the central section, with continuous enrichment functions (polynomial functions) and discontinuous enrichment functions (Heaviside function) shown in Figures 1a and 1b, respectively. The shape function of the GFEM result of the product of the POU and the enrichment function is shown at the bottom.

The combination of the POU functions and enrichment functions established in Figure (1) gives an approximation shown in equation (3):

$\tilde{\mathit{u}}\left(\mathit{x}\right)=\sum _{j=1}^n{N}_j\left(\mathit{x}\right)\left\{{\mathit{u}}_j+\sum _{i=1}^q{L}_{ji}\left(\mathit{x}\right){\mathit{b}}_{ji}\right\}$ (3)

in which ${\mathit{u}}_j$ are the degrees of freedom of the structure tied to the knot ${\mathit{x}}_j$ cloud $w_j$ and ${\mathit{b}}_{ji}$ are additional degrees of freedom in correspondence with each component of the enriched functions.

The structures with smooth solutions (Figure 1a) are accurately simulated using coarse meshes when polynomial enrichments of different orders are selected for approximation $\tilde{\mathit{u}}\left(\mathit{x}\right)$. The proposed method makes use of interface enrichment. Figure 2 shows the definition of the proposed method with interface enrichment in a domain with two materials.

Figure 2
Approaches in a domain with 2 materials. (a) Conventional FEM interpolation; (b) Interpolation with a non-conforming element; (c) Interpolation using the proposed method with interface enrichment.

In Figure 2a, the two materials domain is discretized by two elements, one for each material, where FEM performs interpolation using Lagrangian shape functions as interpolation functions.

The FEM approximation is shown in equation (4):$\mathit{}$

$u(x)=N_j^{((}1))u_j+N_{(j+1)}^{((}2))u_{(j+1)}+(N_{(j+2)}^{((}1))+N_{(j+2)}^{((}2)))u_{(j+2)}$ (4)

in which, $N_j^{\left(i\right)}$ are the Lagrangian interpolation functions at the node $j$ associated with the element $i$, ${\mathit{u}}_j$ are displacements at each nodal point $j$.

However, if the two elements of Figure 2a are grouped together and processed as a single element in non-compliance, a domain with two materials belonging to the same element will be created, as seen in Figure 2b. The interface region between the materials of this element is unable to reconstruct the discontinuity gradient and the jump in the displacement value does not occur in the interface region. Figure 2c shows how the proposed method's strategy with interface enrichment can be used to retrieve the missing part in the interpolation of the displacement field.

The strategy is to use a domain with two materials but only one element on which the extract Lagrangian shape functions $N_j^{\left(p\right)}$ are called the parent element for each node $j$, i.e., the element which houses the interface region. Consequently, it creates nodes at the intersection of the edge element with the interface, and it increases the number of degrees of freedom, but these nodes do not belong to the original mesh (Soghrati et al., 2012).

With nodes entered in the interface junction with the edge of the element, it creates new elements, called integration elements, as shown in Figure 2c, where each integration element belongs to each material, having Lagrangian shape functions associated with each nodal point of the same, which are responsible for the calculation of numerical integration on each material separately. Finally, it connects each integration element to the parent element in order to mount the element stiffness matrix. Equation (5a) depicts the proposed method's approximation $\mathit{u}\left(\mathit{x}\right)$ with interface enrichment:

$\mathit{u}\left(\mathit{x}\right)=N_j^{\left(p\right)}{\mathit{u}}_j+N_{j+1}^{\left(p\right)}{\mathit{u}}_{j+1}+{\psi }_i{\alpha }_i$

${\psi }_i=N_{j+2}^{\left(1\right)}+N_{j+2}^{\left(2\right)}$ (5,a,b,c)

${\alpha }_i={\mathit{u}}_{j+2}-{\mathit{u}}_{j+2}^I$

in which, $N_j^{\left(p\right)}$ are the Lagrangian shape functions of the parent element in each node $j$, ${\mathit{u}}_j$ are the displacements at each nodal point $j$ of the parent element, ${\psi }_i$ are the interface enrichment functions for each node $i$ of interface, ${\alpha }_i$are the degrees of freedom generalized for each interface node $i$, ${\mathit{u}}_j^I$ are the displacements obtained in the interface region between materials when simulated by FEM with only one element, i.e. without reconstruction gradient discontinuity in the interface region as mentioned previously, and $\mathit{u}\left(\mathit{x}\right)$ is the generic approximation of displacement by proposed method with interface enrichment. Enrichment functions are constructed by combining Lagrangian shape functions that correspond to the node shared by descendant elements intersecting the same interface region. (Soghrati and Geubelle, 2012). In Figure 2c, as two nodes were inserted, two enrichment functions have been added. Equation (6) shows a generic way to describe this approximation $\mathit{u}\left(\mathit{x}\right)$ of the proposed method with interface enrichment:

$\mathit{u}\left(\mathit{x}\right)=\sum _{j=1}^n{N}_j\left(\mathit{x}\right){\mathit{u}}_j+\sum _{i=1}^{n_{en}}{s\psi }_i\left(\mathit{x}\right){\alpha }_i$ (6)

where $s$ is a scaling factor related to the aspect ratio between the integration element and the parent element, and n_en is the total number of nodes enriched with interface functions. The interface enrichment functions are evaluated by splitting the parent element into a minimum number of integration elements to obtain a quadrature need.

In Figure 3, the triangular parent element is intercepted by the interface and divided into three triangular integration elements. The enrichment functions corresponding to each node of the interface of Figure 3 are described by equations (7a,b).

Figure 3
Evaluation of enrichment functions in the proposed method with interface enrichment for triangular elements.

${\psi }_{1^{\mathrm{\text{'}}}}=N_3^{\left(1\right)}+N_1^{\left(3\right)}$

${\psi }_{2^{\mathrm{\text{'}}}}=N_2^{\left(1\right)}+N_3^{\left(2\right)}+N_2^{\left(3\right)}$ (7a,b)

Figure 4 shows the quadrilateral parent element divided by the interface into two quadrilateral integration elements. The enrichment functions for each node of Figure 4 interface are described by equations (8a,b).

Figure 4
Evaluation of enrichment functions in the proposed method with interface enrichment for quadrilateral elements.

${\psi }_1=N_4^{\left(2\right)}+N_1^{\left(1\right)}$

${\psi }_2=N_3^{\left(2\right)}+N_2^{\left(1\right)}$ (8a,b)

According to Soghrati and Geubelle (2012), the subelements are only created to evaluate the enrichment functions in numerical integration. Because these integration elements have high aspect ratios, the gradient of enrichment functions has high values. This can cause result in the formation of an ill-conditioned stiffness matrix. To avoid this problem, a scale factor $s$ is implemented, which appears in the second part of the equation (6). This factor can be used for any aspect ratio value or when it is less than a specific value found in Soghrati et al. (2012). The proposed method uses polynomial and interface enrichment functions together to get solutions through mesh refinement and/or polynomial degree enrichment.

3 COMPUTER MODELING OF COMPOSITE MATERIALS

This section presents the numerical simulations of structures with linear elastic behavior in a two-material analysis domain used to validate the proposed method with interface enrichment. A comparison with the conventional methods presented in the preceding sections is also provided. The first validation example is a beam subjected to simple bending composed of two different materials and discretized by triangular elements. The second example corresponds to a cantilever sheet subjected to a tensile force constituted by composite material that is discretized by quadrilateral elements.

3.1 Application in laminated beams

In this example, a cantilever beam composed of two distinct materials is analyzed with the following dimensions: L = 10, b = 2, and unitary thickness (dimensions in consistent units). As illustrated in Figure 5, the beam is subjected to simple bending due to a stress applied to its free edge.

Figure 5
Geometry and loading of the two-material-composed beam subjected to simple bending.

As for mechanical properties, the beam has a Reference Modulus of Elasticity E_R = 30E+06, Poisson ratio ν = 0.25. The loading was regarded as an applied stress ${\tau (0,y)=\tau }_0=150$ (consistent units).

The E_R serves as a comparison baseline to which the values of E₁ and E₂ are related. The objective is to compare and contrast the convergence to strain energy values using the conventional method with polynomial enrichment functions and the proposed method with interface and polynomial enrichment functions implemented together to calculate approximations. The analytical solution for strain energy, $U$, is given by:

$U=\frac{{\tau }_0^2{b}^2{L}^3}{6E_{tr}{I}_{tr}}\left[1+3(1+\nu )/L^2\right]$ (9)

in which, $E_{tr}$ is the longitudinal modulus of elasticity of the material to be transformed and $I_{tr}$ is the inertia of the transformed cross section. In this model, both the conventional method and the proposed method with interface enrichment are applied to triangular elements with polynomial enrichment functions.

As shown in Figure 6a for the conventional method and Figure 6b for the proposed method with interface enrichment, these analyses are conducted on five different discretization of the structure with respect to the mesh.

Figure 6
Different mesh configurations for the two-materials-composited beam: (a) the conventional method; (b) the proposed method with interface enrichment.

Table 1 summarizes the number of elements and nodes for each mesh displayed for the conventional method (Figure 6a) and for the proposed method with interface enrichment (Figure 6b), respectively.

Table 1
Number of elements and number of nodes in each mesh for the conventional method and the proposed method with interface enrichment.

Configuration Mesh	Conventional method		Proposed method
	Number of elements	Number of nodes	Number of elements	Number of nodes
I	8	9	4	6
II	32	25	24	20
III	128	81	112	72
IV	512	289	480	272
V	2048	1089	1984	1056

Figure 6a depicts the configuration meshes for the conventional method that are generated to conform to the interface material, i.e., the interface does not cut the element, but is instead aligned with the edges of the elements that form it. In the proposed method with interface enrichment (Figure 6b), the meshes do not conform to the interface and cut all housed elements.

Regarding the displacements on the beam example, different boundary conditions for the conventional and proposed method were used. In the conventional method, all of the cantilever's nodes are used for horizontal constraint, while only the central node is used for vertical constraint. In the proposed method, the beam's central coordinate at the cantilever edge has an interface node rather than a geometric mesh node. Alternatively, it constrains the interface node and the two nearest mesh nodes. The displacements of the nodes located along the cantilever's edge were not enriched.

The relationship between materials for numerical simulation in which the elasticity modulus has values equal to E₁= 10E_R for the material 1, and E₂ = E_R for the material 2. The ratio E₁/E₂ = 10 is the calculated reference value for the relationship between elasticity moduli, and U = 0.2841 is the calculated reference value for strain energy.

3.1.1 Analysis of mesh refinement convergence

An initial analysis for both methods was to obtain the strain energy while taking mesh refinement into account, starting with mesh I (coarser) and progressing to mesh V (finer), as shown in Figure 6a,b for different degrees of polynomial enrichment: P0 to P3. Figure 7 shows corresponding graphics to relative error with respect to number of degrees of freedom (ndof), obtained through mesh refinement for each degree of approximation enriched separately, analyzed for strain energy for the conventional method (Figure 7a) and the proposed method with interface enrichment (Figure 7b), and a ratio E₁/E₂=10. The increased ndof in these graphs is due to an increase in the number of elements for each polynomial degree used.

Figure 7
Relative error ${\epsilon }_r$, for the strain energy in relation to the ndof for different degrees of enrichment functions to the relationship E₁/E₂ = 10.

In Figure 7a it is observed that the solution quality only is strengthened, when using polynomial enrichment on the approximation provided by conventional method. Figure 7a shows that the solution quality is only improved when polynomial enrichment is applied to the conventional method's approximation. The approximation without enrichment (P0) reduces ${\epsilon }_r$ approximately from 94% (Mesh I) to 7% (Mesh V), but does not reach an acceptable convergence. This approach is carried out by the traditional FEM. It can be seen that for every enrichment of polynomial degree P1 to P3, the convergence of the solution is guaranteed, which occurs according to the mesh refinement. The approximation with enrichment P1 has ${\epsilon }_r$ ≈ 8% in mesh I, which decreases to ${\epsilon }_r$ ≈ 1% in mesh V. It is also evidenced that in mesh I with final approximations P2 and P3, the ${\epsilon }_r$ is small on the order of 3% and 2%, respectively, decreasing to around 1% in mesh V. This demonstrates that the polynomial enrichment in the conventional method has good precision even with coarse meshes. Figure 7b shows that, as with the conventional method (Figure 7a), a reduction in ${\epsilon }_r$ is obtained with mesh refinement for the simulation of the proposed method. In the proposed method, the P0 approximation has ${\epsilon }_r$ ≈ 88% in the mesh I and goes to ${\epsilon }_r$ ≈ 7% in the mesh V, whereas the P1 approximation has ${\epsilon }_r$≈ 8% in the mesh I and goes to ${\epsilon }_r$ ≈ 1% in the mesh V. The resulting approximations P2 and P3 enhanced the solution quality in the desired way for mesh I with ${\epsilon }_r$ ≈ 2% and ${\epsilon }_r$ ≈ 1% in mesh V. The proposed method is effective in approximating achieving satisfactory results, as is the approximation of the conventional method for the relationship between the established Young's modulus and the polynomial enrichment degree.

3.1.2 Analysis of the convergence of polynomial degree of enrichment

Another analysis was performed in Figure 6a and b to determine the convergence of polynomial enrichment degrees for each mesh. It was initially used the linear approximation of POU of the FEM without adding enrichment, yielding a final approximation P0, and in other cases the POU approximation with enrichment P1, P2, and P3.

Figure 8 shows the relationship between ${\epsilon }_r$ and ndof for each mesh, with added enrichment functions P0 to P3 to analyze the strain energy, using the conventional method (Figure 8a) and the proposed method with interface enrichment (Figure 8b). The increased ndof in these graphs is due to the increase in the degree of the approximation for each mesh.

Figure 8
Relative error, ${\epsilon }_r$, with respect to the ndof for each mesh by adding enrichment functions at different degrees on each mesh for determining the strain energy.

Figure 8a illustrates how, for all mesh settings, ${\epsilon }_r$ decreased as polynomial degree enrichment was added, resulting in an ${\epsilon }_r$ ≈ 1%. for all meshes. The use of mesh I with only eight triangular elements produced better results in simulations with polynomial enrichment P3 than simulations with more elements and nodes, resulting in a satisfactory gain in computational cost, demonstrating the capability of the conventional method in the magnification of the space of approximation. Figure 8b presents the simulations using the proposed method. By including the polynomial enrichment in each mesh configuration, a reduction in ${\epsilon }_r$ can be identified (${\epsilon }_r$ ≈ 1%). Note that the convergence curves are very similar to the conventional method. Therefore, the proposed method retrieves the convergence rates from the conventional method even without matching the mesh to the material boundaries.

The method with interface enrichment is advantageous because the interface cuts some elements, and with the addition of polynomial enrichments, it produces accurate results regardless of the number of elements and nodes in the chosen mesh. Because it generates the mesh only once, with only the interface coordinates changing, the proposed method with interface enrichment is very useful in problems involving changes in material form.

3.2 Application in composite panel with matrix and inclusions

In this model, a piece set on one end and subjected to a tensile force at the other end is analyzed. The simulated piece consists of a heterogeneous composite material formed by the inclusions with an elastic constitutive constant ${\mathrm{C}}_{\mathrm{I}}^{11}$ and a matrix composed of an elastic constitutive constant ${\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$, both related to the loading direction. Figure 9 depicts the composite material in three different settings in terms of volume.

Figure 9
Geometry and loading of the composite model for different volume relations

The piece has a width a = 1, length b = 1, and a thickness h = 0.036. An elastic constitutive constant reference, $C_{\mathrm{R}}^{11}$ = 30E+06, and a Poisson ratio, ν = 0 are used as mechanical properties for the sheet. It has a uniformly distributed force on the free end (all units are consistent).

The goal in this example is to analyze the displacement responses in order to define an equivalent constitutive constant ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$, which changes the value of the inclusions ${\mathrm{C}}_{\mathrm{I}}^{11}$ in relation to the matrix ${\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$, which remains constant for analysis by the proposed interface enrichment method. Figure 10 depicts four different mesh settings with V₁/V₂ = 0.16.

Figure 10
Mesh settings for V₁/V₂ = 0.16.

FEM was used to simulate the mesh I, with the elements arranged in the interface between the matrix and the inclusions. The proposed method with interface enrichment simulates the meshes II, III, and IV, with the interface of inclusions cutting some elements in the discretization. Along with the interface enrichment functions, polynomial enrichment functions of linear and quadratic degree are employed in mesh IV.

In the analyses carried out by the proposed method with interface enrichment to the mesh (Figure 10), the parent element is divided into two or three integration elements for the calculation of numerical integration. Figure 11 shows the addition of a corner region where the parent element has been divided into three integration elements. Equations 10 a, b, and c show the evaluation of enrichment functions for each interface node.

Figure 11
Enrichment functions in 2D of the proposed method with the creation of integration elements for four-node quadrilateral elements, cut by an interface made by two line segments intercepted.

${\psi }_1=N_3^{\left(1\right)}+N_3^{\left(2\right)}$

${\psi }_2=N_4^{\left(1\right)}+N_2^{\left(2\right)}+N_3^{\left(3\right)}$ (10a,b,c)

${\psi }_3=N_1^{\left(2\right)}+N_4^{\left(3\right)}$

Beyond the nodes of the edge of the element with the interface, nodes are created within the element and this procedure applies when the objective is to reduce errors caused by the geometry relating to various forms of interface between materials. Figure 12 shows the normalized values ${}^{\ast }C{}_{I+II}^{11}$ of the equivalent constitutive constant ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ for various mesh configurations, as well as the types of simulations for different ratios ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$ to the value ${\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}={\mathrm{C}}_{\mathrm{R}}^{11}$. The ${}^{\ast }C{}_{I+II}^{11}$ values are the normalization of each ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ in relation to the ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ of the mesh I FEM case for the relation ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}=1$.

Figure 12
Normalized ${}^{\ast }C{}_{I+II}^{11}$ values for different meshes and types of simulations with different $C_I^{11}/C_{II}^{11}$relations.

Figure 12 depicts a significant proximity of the ${}^{\ast }C{}_{I+II}^{11}$ between mesh I simulated by FEM (reference value) and mesh II, simulated by the proposed method with enrichment of the interface with the same number of elements for all established ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$ relationships. The mesh II obtained an ${\epsilon }_r$ = 0% for the ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}=1$ ratio, a result obtained for all meshes and the largest ${\epsilon }_r$ was 0.06%, achieved in the ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}=2000$ ratio ensuring excellent approximations for the proposed method with interface enrichment in refined mesh (convergence h). The difference ${}^{\ast }C{}_{I+II}^{11}$ for the mesh increases as the number of elements decreases and the value of the ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$ ratio between the matrix and the inclusions increases, as can be note in meshes III and IV. However, the intermediate mesh III of the proposed method with interface enrichment obtains satisfactory results with the highest ${\epsilon }_r$ found of 0.47% for ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}=2000$ ratio.

When the mesh IV is analyzed for different enrichments (P0 to P2), it is clear that polynomial enrichment improves the solution, but not in the desired way for the quality of the result in the proposed method with interface enrichment. The P0 approximation has ${\epsilon }_r$ ≈ 1.7% for the ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}=10$ ratio and reduces to ${\epsilon }_r$ ≈ 1.2% in the P1 and P2 approximations. The P0 solutions for the ratios ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}>10$ have ${\epsilon }_r$ ≈ 2.5% which drops to ${\epsilon }_r$ ≈ 2% in the P1 and P2 solutions.

Figure 13 shows ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ values using the mesh II for various rigidities relations (${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11})$ and volume relations. It can be seen a constant result value for ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ in the case V₁/V₂ = 0 ratio. This is due to the absence of the inclusions, obtaining only the value of the constitutive constant of the matrix.

Figure 13
${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ Value in mesh for different relations of $C_I^{11}/C_{II}^{11}$ and proportions between volumes.

Figures 12 and 13 show how the curves ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ change relative to the reference as the inclusion volumes in relation to the total part volume increase. This occurs because the percentage of inclusions increases or decreases with respect to the matrix that serves as the reference value, separating the values of ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ from the matrix's elastic constitutive constant (${\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$) and approaching to the inclusions' elastic constitutive constant (${\mathrm{C}}_{\mathrm{I}}^{11}$). As a result, the values of ${\mathrm{C}}_{\mathrm{I}+\mathrm{I}\mathrm{I}}^{11}$ for any relationship between the studied intervals can be extracted. This is useful for materials for composite design that specify stiffness as well as determining the relationship between constitutive elastic constants ${\mathrm{C}}_{\mathrm{I}}^{11}/{\mathrm{C}}_{\mathrm{I}\mathrm{I}}^{11}$ and volumes V₁/V₂ for matrix and inclusions using the graphs in Figure 13.

This example highlights the effectiveness of using the proposed method with interface enrichment for optimization problems in composite materials. The proposed technique minimizes the generation of the finite element mesh for the discretization of different boundary problems in situations where the volume of inclusions is consistently changing.

4 CONCLUSIONS

This study focuses on the formulation and implementation of unconventional FEM methods, specifically a proposed method that uses polynomial and interface enrichment functions for applications in two-dimensional composite material structures. To validate the proposed formulations and implementations, two models were numerically simulated. The proposed enrichment strategy for composite materials demonstrated great potential for solving problems with discontinuous gradients caused by the stiffness differences between the materials without using a mesh that fits the interface edges between two materials. The main feature of the proposed method is that the degrees of freedom of enrichment are assigned only to nodes with fictitious interfaces. This formulation change eliminates the problems encountered in some enrichment functions in the conventional method for awarding essential boundary conditions of enriched nodes. Furthermore, the enrichment functions are simply constructed by linear combination of Lagrangian shape functions in integration elements, which reduces the cost and simplifies its application.

The method was validated for composite laminates. The proposed method's convergence in the laminated beams with the interface enrichment, the polynomial degree of enrichment (p-refinement) and mesh size refinement (h-refinement) is satisfactory and comparable to the conventional method. This demonstrates that the proposed method can be used with reduced meshes without compromising solution precision. The proposed method with interface enrichment has the advantage of producing good results for approximation with the use of elements that embed the interface without remeshing.

The proposed method with interface enrichment efficiently captures the values of the equivalent constitutive constants in composite problems and provides good h-refinement. The results demonstrated the utility of these analyses in the development of composite materials by optimizing the equivalent stiffness of the composite through the simultaneous optimization of the rigidity and volume of the inclusions.

Finally, it is concluded that using the method improves the quality of the final approximation, with the rational use of enrichment functions improving the computational efficiency. The application of the proposed interface functions significantly reduces the problems in the modeling of problems with discontinuous gradient fields, particularly in problems involving more than one material, such as composite materials.

Acknowledgment

The authors are grateful to CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), FAPDF (Fundação de Apoio à Pesquisa do Distrito federal), UFPI (Universidade Federal do Piauí), UnB (Universidade de Brasília), and IFPI (Instituto Federal do Piauí).

• https://doi.org/10.1590/1679-78257273

References

• Ahmadi, I., Sheikhy, N., Aghdam, M.M., Nourazar, S.S., (2010). A new local meshless method for steady-state heat conduction in heterogeneous materials, Engineering Analysis with Boundary Elements 34 (12): 1105-1112.
• An, X.M., Zhao, Z.Y., Zhang, H.H., He, L., (2013). Modeling bimaterial interface cracks using the numerical manifold method, Engineering Analysis with Boundary Elements 37 (2): 464-474.
• Aragón, A.M., Liang, B., Ahmadian, H., Soghrati, S., (2020). On the stability and interpolating properties of the hierarchical interface-enriched finite element method, Computer Methods in Applied Mechanics and Engineering 362: 112671.
• Babuška, I., Caloz, G., Osborn, J.E., (1994). Special finite element method for a class second order elliptic problems with rough coefficients, SIAM Journal on Numerical Analysis 31 (4): 945 - 981.
• Babuška, I., Melenk, J. M., (1997). The partition of unity method, International Journal for Numerical Methos in Engineering 40: 727-758.
• Belytschko, T., Black, T., (1999). Elastic crack growth in finite elements with minimal remeshing, International Journal for Numerical Methods in Engineering 45 (5): 601-620.
• Campos, B. C., Barros, F. B., Penna, S. S. (2022). Numerical integration in G/XFEM analysis of 2-D fracture mechanics problems for physically nonlinear material and cohesive crack propagation,Engineering Computations, 39 (3): 1134-1160.
• Duarte, C. A., Babuška, I., Oden, J. T., (2000). Generalized finite element methods for three-dimensional Structural mechanics problems, Computers & Structures 77 (2): 215-232.
• Duarte, C. A., Hamzeh, O. N., Liszka, T. J., Tworzydlo, W. W., (2001). A generalized finite element method for the simulation of three dimensional dynamic crack propagation, Comput. Methods Appl. Mech. Engrg., 190: 2227-2262.
• Duarte, C.A., Kim, D. J. (2008). Analysis and applications of a generalized finite element method with global-local enrichment functions, Comput. Methods Appl. Mech. Engrg., 197: 487-504.
• Duarte, C. A., Oden, J. T., (1995). Hp clouds - a meshless method to solve boundary value problem, Technical Report 9505, TICAM, University of Texas at Austin.
• Duarte, C. A., Oden, J. T., (1996). An h-p adaptive method using clouds, Computer Methods in Applied Mechanics and Engineering 139 (1-4): 237-262.
• Evangelista Jr, F., Alves, G. S., Moreira, J. F. A and Paiva, G.O. F., (2020). A global-local strategy with the generalized finite element framework for continuum damage models, Comput. Methods Appl. Mech. Engrg., 363: 112888.
• Evangelista Jr, F., Moreira, J. F. A., (2020). A novel continuum damage model to simulate quasi-brittle failure in mode I and mixed-mode conditions using a continuous or a continuous-discontinuous strategy, Theoretical and Applied Fracture Mechanics, 1: 102745.
• Gomes, L. L., Barros, F. B., Penna, S. S. and Pitangueira, R. L. d. S. (2021). Geometrically nonlinear analysis by the generalized finite element method,Engineering Computations, 38 (1): 266-288.
• Kim, D. J., Duarte, C. A., Pereira, J. P. (2008). Analysis of interacting cracks using the generalized finite element method with global local enrichment functions, J. Appl. Mech., 75 (5): 763-813.
• Malekan, M., Barros, F.B., (2018). Numerical Analysis of a main crack interactions with micro-defects/inhomogeneities using two-scale generalized/extended finite element method, Comput. Mech., 62 (4): 783-801.
• Melenk, J. M., Babuška, I., (1996). The partition of unity finite element method: Basic theory and applications, Computer Methods in Applied Mechanics and Engineering 139 (1-4): 289-314.
• Moës, N., Dolbow, J., Belytschko, T., (1999). A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineeriing 46: 131-150.
• Moës, N., Cloiree, M., Cartraud, P., Remacle, J. F., (2003). A computational approach to handle complex microstructure geometries, Computer methods in applied mechanics and engineering 192: 3163-3177.
• Monteiro, H. A. S., Fonseca, G. M., Novelli, L., Pitangueira, R. L. S., Barros, F. B., (2018). Análisis global-local de medios lineales utilizando el Método de Elementos Finitos Generalizados, Mecánica Computacional, v. XXXVI. Number 32 : p. 1539-1368.
• Monteiro, H. A. S., Novelli, L., Fonseca, G. M., Pitangueira, R. L. S., Barros, F. B. (2020). A new approach for physically nonlinear analysis of continuum damage mechanics problems using the generalized/extended finite element method with global-local enrichment, Engineering Analysis with Boundary Elements, 113: 277-295.
• Nagarajan, A., Soghrati, S., (2018). Conforming to interface structured adaptive mesh refinement: 3D algorithm and implementation, Computational Mechanics 62 (5): 1213-1238.
• Novelli, L., Oliveira, T. S., Monteiro, H. A. S., Fonseca, G. M., Pitangueira, R. L. S., Barros, F. B. (2020). Stable Generalized/eXtended Finite Element Method with global-local enrichment for material nonlinear analysis, Comput. Methods Appl. Mech. Engrg., 372: 113429.
• Oden, J. T., Duarte, C.A., Zienkiewicz, O. C., (1998). A new cloud-basedhp finite element method, Computer methods in applied mechanics and engineering 153: 117-126.
• Phan, A.-V, Mukherjee, S., (2009). The multi-domain boundary contour method for interface and dissimilar material problems, Engineering Analysis with Boundary Elements 33(5): 668-677.
• Simone, A., Duarte, C. A., Giessen, E. V., (2006). A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries, International Journal for Numerical Methods in Engineering 67 (8): 1122-1145.
• Soghrati, S., (2014). Hierarchical interface-enriched finite element method: An automated technique for mesh-independent simulations, Journal of Computational Physics 275: 41-52.
• Soghrati, S., Ahmadian, H., (2015). 3D hierarchical interface-enriched finite element method: Implementation and applications, Journal of Computational Physics 299: 45-55.
• Soghrati, S., Aragón, A. M., Duarte, C.A., Geubelle, P. H., (2012). An interface-enriched generalized FEM for problems with discontinuous gradient fields, International Journal for numerical methods in engineering 89: 991-1008.
• Soghrati, S., Barrera, J.L., (2016). On the application of higher‐order elements in the hierarchical interface‐enriched finite element method, International Journal for Numerical Methods in Engineering 105 (6): 403-415.
• Soghrati, S., Chen, Y., Mai, W., (2018). A non-iterative local remeshing approach for simulating moving boundary transient diffusion problems, Finite Elements in Analysis and Design 140: 23-37.
• Soghrati, S., Duarte, C.A., Geubelle, P. H., (2015). An adaptive interface-enriched generralized FEM for the treatment of problems with curved interfaces. International Journal for Numerical Methods in Engineering. 102 (6): 1352-1370.
• Soghrati, S., Geubelle, P. H., (2012). An 3D interface-enriched generalized finite element method for weakly discontinuous problems with complex internal geometries, Comput. Methods Appl. Mech. Engrg (217-220): 46-57.
• Soghrati, S., Nagarajan, A., Liang, B., (2017). Conforming to interface structured adaptive mesh refinement: New technique for the automated modeling of materials with complex microstructures, Finite Elements in Analysis and Design 125: 24-40.
• Strouboulis, T., Babuška, I., Copps, K., (2000). The design and analysis of the Generalized Finite Element Method, Computer Methods in Applied Mechanics and Engineering 181 (1-3): 43-69.
• Tan, M.H.Y., Geubelle, P.H., (2017). 3D dimensionally reduced modeling and gradient-based optimization of microchannel cooling networks, Computer Methods in Applied Mechanics and Engineering 323: 230-249.
• Wolenski, A. R. V., Monteiro, A. B., Penna, S. S., Pitangueira, R. L. S., Barros, F. B., (2020). Damage propagation using novel G/XFEM strategies: computational aspects and numerical investigations, Journal of the Brazilian Society of Mechanical Sciences and Engineering, 42: 1-14.
• Zhang, H., Wang, D., (2015). An isogeometric enriched quasi-convex meshfree formulation with application to material interface modeling, Engineering Analysis with Boundary Elements 60: 37-50.

Edited by

Editor: Marco L. Bittencourt

Publication Dates

• Publication in this collection
  07 Nov 2022
• Date of issue
  2022

History

• Received
  05 Sept 2022
• Reviewed
  07 Oct 2022
• Accepted
  10 Oct 2022