Topology Optimization – unconventional approaches using the Generalized Finite Element Method and the Stable Generalized Finite Element Method

Abstract

The Structural Optimization process has an increasing importance in industry and academic fields, assisting in the development of designs at the initial stages of a project. Nowadays, the structural optimization methodology can be conducted by Topology Optimization Method (TOM), which is an efficiently combination of the Finite Element Method (FEM) with an optimization algorithm, in order to find the optimized material distribution inside a given design domain subjected to a set of constraints. Application of the FEM in TOM suffers from a series of instability problems, being one of them the checkerboard pattern. This paper investigates the impact of the Generalized Finite Element Method (GFEM) and Stable Generalized Finite Element Method (SGFEM) in the implementation of the TOM algorithm. This work shows that these unconventional FEM formulations are able to solve most of the checkerboard pattern problem when combined with an enriched mesh designed specifically to each example evaluated. Significant improvement in results of the topology optimization is achieved when compared to the conventional formulation of TOM.

Keywords
Generalized Finite Element Method; Stable Generalized Finite Element Method; Topology Optimization; Checkerboard Pattern

Graphical Abstract

1 INTRODUCTION

Basically, the Topology Optimization Method (TOM) is a mathematical methodology that seeks an optimized material distribution in a design domain using the conventional Finite Element Method (FEM) (Bathe, 1996). This process is systematically performed by the identification and removal of unnecessary members and/or connections presented in the initial design domain, thus providing a creative and often non-intuitive solution in the early development phase of a project.

Regardless of its excellent results and broad spectrum of applications, it is possible that due to numerical problems linked with the FEM formulation, the quality of results can be undermined, and leading to an incorrect or undesired output. These classical numerical problems are commonly related to mesh dependency (a domain discretization problem), and checkerboard problems (a numerical instability problem) (Bendsøe and Sigmund, 2003).

The Generalized Finite Element Method (GFEM), the Stable Generalized Finite Element Method (SGFEM), and other unconventional FEM formulations have been developed with the aim of overcoming classical limitations presented by the FEM conventional form. Specifically, the GFEM is a numerical alternative obtained from the composition of positive characteristics of FEM and the Meshless Methods, as well as the SGFEM. Basically, in the GEFM a mesh is employed to define support clouds 𝜔 within, in which the classic FEM shape functions are used as Partition of Unit (PoU). These shape functions can be enriched, within each cloud, by multiplying the PoU by functions of interest (polynomial or not) (Duarte et al., 2000). In the SGFEM structure, enrichment is done in a similar way, with the addition of interpolation functions to each enriched PoU. It is important to emphasize that the result of enrichment is nodal, (Duarte et al., 2000).

Nodal enrichment, in summary, can be defined as the amplification of the approximate bases involved, without the need to introduce new nodal points in the domain. Therefore, this strategy differs from the classical p-adaptive FEM process, which, for a certain class of finite elements, requires the addition of extra nodal points in the finite elements to increase the interpolation degree of the shape functions, (Duarte and Oden, 1995).

In recent years, the number of works regarding Topology Optimization and non-conventional formulations of the FEM has increased. It has been shown that the use of different versions of Smoothed Finite Element Method (Node-Based and Cell Based) may bring interesting results (He et al., 2015; Shobeiri, 2016), in terms of efficiency, although not tackling the checkerboard problem or the use of sensitivity filters. Similarly, the use of the TOM combined with Interface-enriched Generalized Finite Element Method (IGFEM) and level-set functions provide smooth results with well-defined borders and narrow structures. Although this method found improved results compared to classical TOM, it still relies on filtering methods during image processing (van den Boom et al., 2021).

This work aims to contribute with studies of the impacts of the GFEM and SGFEM formulations in solving problems associated specifically to the checkerboard pattern, based on nodal enrichment technique in TOM. In other words, the TOM-GFEM and TOM-SGFEM algorithms are analyzed to better understand the best enrichment conditions to face the checkerboard pattern problem. It aims to do so without relying on sensitivity filters of any kind, as seen on the aforementioned works (He et al., 2015; Shobeiri, 2016; van den Boom et al., 2021). In addition, this work can consolidate the application of GFEM and SGFEM in another unconventional finite element formulations - Hybrid-Mixed Stress Formulation (HMSF) and its combination with GFEM and SGFEM - HMSF-GFEM and HMSF-SGFEM - for the analysis of the stress constraint problem in TOM, for example. The application of HMSF in TOM can be observed in Navarro and Góis (2017).

The paper outline is as follows. In Section 2, the main TOM features are presented. Section 3 describes the fundamental relations of GFEM and SGFEM applied to plane isotropic problems. Since both formulations share much in common, except for one step during enrichment, they were explained together. Section 4 shows the numerical implementation of TOM combined with both methods (GFEM and SGFEM). In Section 5 we highlight the developed numerical simulation results. Finally, in Section 6, the main conclusions of this work are summarized.

2 TOPOLOGY OPTIMIZATION

The main goal of Structural Optimization is the search for optimal material layouts (design variables) to sustain given loads, while satisfying the imposed constraints, (Christensen and Klarbring, 2009). The determination of these optimal values is done according to several criteria, often resulting in complex and original structures. Nowadays, a very robust method for this purpose is the TOM, which allows generating new boundaries (holes) in the design of the structure domain, by removing and adding material in each point of that domain in order to extremize (minimize or maximize) a specified objective function, satisfying given constraints of the optimization problem. The TOM requires no prior knowledge of the layout of the structure, which is obtained systematically from any initial material distribution, allowing to obtain high structural performance increase (Bendsøe and Sigmund, 2003).

The numerical implementation of the material distribution in structural design was first introduced in Bendsøe and Kikuchi (1988), in which the TOM is employed to minimize the mean compliance objective function. The discrete form of optimization problem is given by:

$\underset{\rho }{Minimize}$ $C\left(\rho \right)=U^{T}KU$

Subjected to $KU=F$

${\displaystyle \sum _{i=1}^{ne}{\rho }_i{V}_i\le V}$ (1)

$0\le {\rho }_{\min }\le \rho \le 1$

where the superscript ^T denotes the transpose vector, C(𝛒) is the mean compliance, K is the global stiffness matrix, U and F are the global displacement and load vectors, respectively. ne is the total number of elements in the design domain, V_i is the volume of individual elements, and V is a limitation on the material allowed to be distributed (volume constraint) in the design domain. Finally, the design variables 𝛒 are limited between 0 and 1, however, to avoid numerical problems a very small value (𝛒_min), instead of zero, is adopted to the lower limit of the design variables.

The material model adopted is the SIMP (Solid Isotropic Material with Penalization) (Bendsøe, 1989; Rozvany et al., 1992; Mlejnek, 1992), which defines the material property E (elasticity tensor), at each point of the domain, as follows:

$E={\rho }^p{E}_0$ (2)

where E₀ is the property of the basic material that will be distributed in the domain, ρ is a pseudo-density (design variable) that describes the amount of material at each point in the domain, and p is a penalty factor that must be adjusted (Bendsøe, 1989) to make the topology optimization process be able to achieve binary solutions (interpreted as black and white during image processing).

The solution of the compliance minimization problem, given by Eq. (1), can lead to certain numerical instabilities, generated by two main issues: mesh dependence, and checkerboard patterns, (Sigmund and Petersson, 1998). The mesh dependence refers to the problem of obtaining different solutions for different discretization, respectively. The checkerboard pattern is an instability associated with the application of the classical FEM in the TOM, in which the density distribution alternates from solid to void in a pattern similar to a checkerboard, as shown in Figure 1.

Figure 1
Arising of checkerboard problem.

For further details on the Topology Optimization formulation, as well as, the numerical problems of the method, the reader may refer to the book by Bendsøe and Sigmund (2003). Here, a brief explanation about the checkerboard pattern is given, since the focus of this work is to reduce the presence of this numerical problem in the results obtained by TOM.

The checkerboard pattern is the result of a numerical instability related to the finite element approximation of the original continuous problem, and it is undesirable since it corresponds to a material distribution with artificially high stiffness. Diaz and Sigmund (1995) discuss the reasons for the formation of such patterns in the topology optimization and show that its origin is due to bad numerical modeling that overestimates the stiffness of checkerboards. They demonstrated that the behavior of bilinear four-node elements arranged in a checkerboard fashion is artificially high stiff, which makes the Topology Optimization (TO) algorithm to priorate the formation of checkerboard patterns, when this element is employed, for numerical solution of the mean compliance optimization problem.

Most common methods applied to prevent the checkerboards in the TO design employ filtering techniques from image processing, which are based on heuristic schemes that modify the design sensitivities or design densities used in each iteration of the TO algorithm (Bendsøe and Sigmund, 2003). Although very efficient, filtering techniques introduce changes to values of the design variables (or sensitivities), after (or before) of the optimization step in the TO algorithm, which modify the exact optimized solution. The use of nine-node elements is also suggested to avoid the checkerboard problem, since meshes with these elements are more flexible than the ones with four-node elements. However, this involves high computational cost, because nine-node elements have many more displacement degrees of freedom per material design (variable of optimization problem). Moreover, higher-order finite elements (eight or nine-nodes) does not work well for SIMP approach using large penalty factors (Diaz and Sigmund, 1995) and, thus, it is rarely used in TO design. Thus, this present work contributes by exploring the use of unconventional finite element formulations, called GFEM and SGFEM, as another alternative to solve the problem of the formation of checkerboard in the TO results, as can be seen ahead.

3 GENERALIZED FINITE ELEMENT METHOD (GFEM) AND STABLE GENERALIZED FINITE ELEMENT METHOD (SGFEM) IN 2D ANALYSIS

For both methods (GFEM and SGFEM), consider an isotropic solid of a volume V (the volume of material in a TOM analysis) and regular boundary 𝛤 defined by complementary regions - 𝛤_t and 𝛤_u, where forces and displacements are imposed, respectively. The equilibrium of this solid is expressed by the well-known Principle of Virtual Work, as follows:

${\displaystyle \underset{V}{\int }\delta {\epsilon }^T\sigma }dV={\displaystyle \underset{V}{\int }\delta u^T\overline{b}}dV+{\displaystyle \underset{{\Gamma }_t}{\int }\delta u^T}\overline{p}d{\Gamma }_t,{\forall }^{}\delta u$ (3)

where u is the displacement vector, 𝝈 and 𝜺 are vector columns that hold the components of stress and strain tensors, respectively. 𝛿u is the admissible virtual displacements field, 𝛿𝜺 is the virtual deformations field compatible with 𝛿u, b̅ represents the vector of volumetric forces and p̅ is the vector of surface force prescribed in 𝛤_t.

For the complete characterization of the proposed elastic problem, the quantities involved in this problem (u, 𝝈, 𝜺, 𝛿u and 𝛿𝜺) must still, at each point of V, obey the compatibility equations, constitutive law, and essential boundary conditions in 𝛤_u. In this way, disregarding the volumetric forces, Eq. (3) takes the following form:

${\displaystyle \underset{V}{\int }\delta u^T}{L}^{T}DLudV={\displaystyle \underset{{\Gamma }_t}{\int }\delta u^T}\overline{p}d{\Gamma }_t,\forall \delta u$ (4)

where L is the partial derivative operator and D is the material stiffness matrix. Another presentation for the problem expressed in Eq. (4) can be written as follows:

Determine $u\in E\left(V\right)$ such that $B\left(u,\delta u\right)=F\left(\delta u\right),\forall \delta u\in E\left(V\right)$ (5)

where B(u, 𝛿u) is the bilinear form (energy inner product), E(V) is the energy space with energy norm defined by $|\left|\mathbf{u}\right|{|}_{E\left(V\right)}=\sqrt{B(\mathbf{u},\delta \mathbf{u})}$ and F(𝛿u) is a linear limited form in E(V). The strain energy is also defined as the value obtained from $\frac{1}{2}B(\mathbf{u},\delta \mathbf{u})$. In Eq. (1), it is observed that the mean compliance C(ρ) is the double of the strain energy $\frac{1}{2}B(\mathbf{u},\delta \mathbf{u})$.

Now, considering that the approach adopted to generate an approximation (via GFEM, for example) for the bilinear form presented in Eq. (5) is given by Galerkin, then the u space is equal to the 𝛿u space. In addition, the use of the GFEM for the approximation structure gives a finite dimension to the space E(V).

In this way, the following approximation is defined for u (analogously extendable to 𝛿u):

$u=NU$ (6)

Where N is the matrix that holds the approximation functions of the displacements and U is the vector that collects the weights of this approximation (Duarte et al., 2000). It is worth remembering that the approximation functions that are in N must be kinematically admissible to the finite dimension space that belongs to space E(V).

The following linear equations system is generated by replacing the approximation defined in Eq. (6) into Eq. (5), already specified for finite dimension:

$KU=F$ (7)

with:

$K={\displaystyle \sum _{e=1}^N{\displaystyle \underset{V_e}{\int }{\left(B\right)}^{T}DBdV}}$ and $F={\displaystyle \sum _{e=1}^N{\displaystyle \underset{{\Gamma }_{t_e}}{\int }{\left(N\right)}^T\overline{p}d{\Gamma }_{t_e}}}$ (8)

Where K is the stiffness matrix, V_e and t_e are the volume and the boundary force of the finite element, respectively, B is the matrix that holds the differential operators applied to form functions and is given by the product LN. It is important to note that if a problem has a concentrated force, the load vector F must be completed.

For the application of the GFEM and SGFEM in plane analysis (2D), it is necessary to define supports or clouds formed by four-node quadrilateral isoparametric finite elements. It is important to note that this finite element will be applied in both kinds of analysis. Each cloud 𝝎_j, see Figure 2, is formed by finite elements that share a common node j.

Figure 2
Cloud example for a two-dimensional discretization (quadrilateral elements).

For the finite element clouds consisting of four-node quadrilaterals, the classic bilinear Lagrangian functions are used as PoU for the interpolation of the displacements (parametrized domain).

Now, in order to provide an enrichment of the approximation attached to the nodes j of the finite element mesh, polynomial functions h_jk, k=1, …, I(j), where I(j) represents the number of these functions added to the index node j, are adopted in this work.

Thus, we have a new shape functions family for the displacements field given by:

${\Theta }_n^m=\left\{{\left\{N_j\right\}}_{j=1}^n\cup {\left\{N_j{h}_{jk}\right\}}_{j=1}^n;j=\mathrm{1,}\cdots ,n;k=\mathrm{1,}\cdots ,I\left(j\right)\right\}$ (9)

and used to construct the following approximation:

$u={\displaystyle \sum _{j=1}^n{N}_j}\left\{U_j+{\displaystyle \sum _{i=1}^k{h}_{ji}{e}_{ij}}\right\}$ (10)

where m is the maximum degree of the resulting approximation, n is the total number of nodes in the discretized domain, N_j is the original shape functions (PoU) for the displacement field, U_j are the degrees of freedom in displacement associated with the originals shape functions and e_ji are the new nodal parameters corresponding to each enrichment condition.

Thus, for example, for quadrilateral four-node finite element meshes, the original approach bases for the displacement field can be written as follows:

$N_{V_e}=\left[\begin{array}{cccc}{\phi }_1{\Delta }_1& {\phi }_2{\Delta }_2& {\phi }_3{\Delta }_3& {\phi }_4{\Delta }_4\end{array}\right]$ (11)

where 𝜙_j, j = 1, …, 4 represents the classical bilinear Lagrangian functions attached to the node. V_e represents the finite element domain of the discretization.

In Eq. (11), Δ_j is the polynomial enrichment matrix of the node j in the element domain. This matrix is given by:

${\Delta }_j=\left[\begin{array}{cccccc}{I}_2& h_{j1}{I}_2& \cdots & h_{jk}{I}_2& \cdots & h_{jI\left(j\right)}{I}_2\end{array}\right]$ (12)

where I₂ is the second-order identity matrix.

It is worth mentioning that the enrichment for both GFEM and SGFEM approximation base is performed in the global domain (and not in the parameterized domain) and that the enriching functions are defined in such a way to be null in the enriched nodes to which they are attached. It is evident that if the functions h_jk are null, it is preserved with the I₂ matrix, the conventional structure of the FEM. Figure 3 illustrates the nodal enrichment technique applied to the FEM shape function.

Figure 3
Construction of GFEM shape function (far left), and of SGFEM shape function (middle and far right). N_j is the FEM shape function; h_jk (GFEM) and h̃_jk (SGFEM) are enrichment functions; Θ^m_n (GFEM) and Θ̃^m_n (SGFEM) are the modified PoU – shape functions adapted from (Gupta et al., 2013).

Up to this point, GFEM and SGFEM take advantage of the same enrichment techniques, however SGFEM counts with the additional resource of interpolating functions added to the already enriched Lagrangian functions. Such differences can cause significant impact in results, reinforcing the need to compare results from both methods. Due to the modifications needed, new enriching functions can be defined as (Babuška and Banerjee, 2012):

${\tilde{h}}_{jk}=h_{jk}-I_{{\omega }_j}\left(h_{jk}\right);I_{{\omega }_j}={\displaystyle \sum _{j=1}^4{N}_j}{h}_{jk}$ (13)

in this case, N_j and h_jk are defined exactly as done for Eq. (9) and (10). ${\mathbf{I}}_{{\omega }_j}$ is the interpolating function defined by N_j and h_jk. Thus, a new enrichment function h̃_jk is defined and substituted in Eq. (9), creating a different shape function Θ̃^m_n (Babuška and Banerjee, 2012):

${\tilde{\Theta }}_n^m=\left\{{\left\{N_j\right\}}_{j=1}^n\cup {\left\{N_j{\tilde{h}}_{jk}\right\}}_{j=1}^n;j=\mathrm{1,}\cdots ,n;k=\mathrm{1,}\cdots ,I\left(j\right)\right\}$ (14)

Finally, the strategies presented along this Section, GFEM and SGFEM, allow the development of an analysis technique using coarser meshes compared to the same analysis using classical FEM.

4 NUMERICAL IMPLEMENTATIONS

The computational implementation of the TOM using GFEM formulation (TOM-GFEM) and SGFEM formulation (TOM-SGFEM) are built in an analogous manner, since both methods differ only in one stage. Initially, a mesh of the model is generated, followed by the definition of enrichment functions, global stiffness matrices and load vectors for a given domain. The computational code was implemented by using MATLAB^® language and any subroutine for filtering is not included. Nevertheless, the algorithm can be adapted to solve problems in three dimensions. To get a copy of the complete program, contact the authors via e-mail.

Since the enrichment process is quite similar in GFEM and SGFEM, both algorithms share most of the routines and subroutines (Figure 4). Initially the TOM-GFEM algorithm was written loosely based on an already existing MATLAB^® code for classical TOM (Sigmund, 2001). To adapt the SGFEM formulation into the TOM, changes made were in the main routine (TOP_GFEM.m that became TOP_SGFEM.m) and in Quad_NBmat.m subroutine, in which enrichment is applied to the PoUs. This subroutine defines the PoUs and also creates the B matrix (described in Eq. (8)), where enriched functions are stored for future calculations.

Figure 4
hierarchic flowchart of both MATLAB® codes.

Figure 5 shows a generic flowchart for both computational codes, where dotted lines represent the different paths the code may follow, depending on what code is running. As it can be checked, the first step consists in reading the input file information to allocate the necessary memory and global variables, such as the number of finite elements (mesh) and optimization parameters, including the penalty factor and desired volume, for example. In general, the code works with rectangular mesh, but it can be easily adapted to work with non-rectangular domains. The following steps consist in determining the enriched displacement field and update the design variables by checking the optimality criteria. The process is then repeated until a convergence criterion is satisfied. In this work, the optimization process is repeated if the change in objective function is greater than a certain value, to guarantee that every iteration has a significant amount of change the topology.

Figure 5
Generic code flowchart for both new algorithms.

5 NUMERICAL EXAMPLES

Two optimization problems are studied in this work. One of them is the Messerschmitt-Bölkow-Blohm beam (MBB Beam), and the other one is Fixed Ends Beam, design domains which are presented in Figure 6. To save computational time, only half of the structure domain model is employed by using symmetry conditions, as shown in Figure 6.

Figure 6
Benchmark models: MBB beam (a), and Fixed-Ends beam (b) (Arruda, 2015).

All the simulations are developed by assuming a rectangular domain discretized into 60x20 square finite elements. A linear material with elastic modulus equal 1.0 Pa, Poisson’s ratio equal 0.3 are adopted. Meanwhile, the volume fraction is adopted as f=0.5 (leading to topologies that present 50% of the volume of the design domain when considering all elements as solid phase) and the penalization factor (p) is changed from 1.0 to 3.0, to check any relation between the method of enrichment and the presence of the checkerboard pattern in the solution obtained.

The complete second order polynomial is adopted for nodal enrichment in TOM-GFEM – given by the following polynomial (x+y+x²+y²+xy), I(j)=5 - see Eq. (9); since it produces better results, as presented by Arruda (2015). Meanwhile, in TOM-SGFEM, due to the use of an interpolant the nodal enrichment adopted is x²+y², I(j)=2, as seen on Eq. (14).

5.1 MBB Beam

Figure 7 shows the finite element mesh adopted to simulate the MBB beam, as well as the applied force and boundary conditions (red arrow and blue triangles, respectively). A structured mesh is adopted for a better application/control of the enriched nodes. Additionally, a preliminary simulation is performed by using a classical TOM formulation that uses conventional finite element, in order to define the workbench output for the MBB Beam. It is important to notice that the present work does not aim to compare run time from the classical TOM computational code (Sigmund, 2001) and the proposed codes with GFEM and SGFEM formulations. Additionally, the result obtained preliminarily by the classical TOM helps to define the best layout for enriched nodes in the new proposed formulations.

Figure 7
Rectangular 60x20 mesh applied at the MBB Beam (Arruda, 2015).

It is shown, in Figure 8, that when sensitivity filters are removed, the checkerboard pattern can be noticed in varying degrees according to the penalization (p) defined for each simulation. Moreover, a direct relation between the amount of checkerboard pattern and each penalty value can be seen, since an increase in the value of p causes a decrease in the area covered by the checkerboard pattern. With this in mind, Figure 8 sets the benchmark topologies obtained by classical TOM code.

Figure 8
Preliminary topology results obtained by classical TOM code, without filtering technique, and used as benchmark.

Besides the topology results, other parameters of interest can be analyzed, such as number of computer iterations and the objective function value (see Table 1). These numbers vary greatly, since the classical TOM code was meant to be used with an active sensitivity filter (Sigmund, 2001). Also, time is not a variable of interest since, due to a larger data set, processing TOM-GFEM and TOM-SGFEM is expected to take considerably longer than the classical TOM code.

Table 1
Benchmark data for MBB Beam topology.

Penalization	Volume Fraction	Iteration	Objective Function
1.0	0.500	20	164.4165
2.0	0.500	97	185.4375
2.5	0.500	47	191.0886
3.0	0.500	65	203.2379

For the sake of simplicity, we chose a uniform distribution of enriched nodes, aiming to cover most of the area obtained in the mesh from Figure 7. This mesh, called Mesh 1 in this work (see Figure 9), has the objective to exemplify the importance of the placement of enriched nodes during simulation. Sparse placement of the modified functions in the mesh can lead to underperforming results in terms of topology, although reaching higher efficiency in terms of processing time and iterations. Moreover, uniform meshes are much simpler to write in a code, justifying the attempt to use them.

Figure 9
MBB Beam uniform nodal enrichment distribution (Mesh 1).

Analyzing the results from each code (TOM, TOM-GFEM, and TOM-SGFEM) and comparing them side by side, we could observe modest improvement in solving the checkerboard problem. For penalty factor p ≤ 2.0 no significant difference in topology is found. However, as the value of p grows, more significant changes can be seen with the best result being TOM-GFEM with p = 3.0, as seen on Figure 10. Nevertheless, in every simulation the checkerboard problem is still not corrected or avoided, indicating that a better distribution of enriched nodes could be provided to improve the results.

Figure 10
Comparison between MBB Beam results from all codes regarding Mesh 1.

Moreover, the value of the objective function also shows minor change when compared to the classical TOM code, especially when p = 3.0, corroborating with the small amount of change in topology. Even considering the difference in number of iterations, change was too little in order to analyze efficiency with a uniform distribution of enriched nodes. Table 2 shows data from all codes side by side for comparison.

Table 2
Data compilation for MBB beam - Mesh 1.

Penalization	TOM		TOM-GFEM		TOM-SGFEM
	Iteration	Obj. Function	Iteration	Obj. Function	Iteration	Obj. Function
2.0	97	185.4375	115	188.1589	78	188.0354
2.5	47	191.0886	55	195.4185	54	194.8524
3.0	65	203.2379	47	203.7820	35	203.7459

Carefully observing the results obtained with classical TOM, and also taking into account the results obtained with the enrichment pattern in Mesh 1, we could arbitrarily choose a non-uniform distribution of enriched nodes for the structural domain (also called Mesh 2), as seen in Figure 9. This layout was chosen based on how severely the checkerboard pattern affected the different regions of the domain during classical TOM simulations and also aiming to cover areas in which the checkerboard pattern remained after simulations made by using Mesh 1 (higher density of enriched nodes in center-right patch of the domain, as seen in Figures 8 and 10).

Here, it is important to notice that the simulations performed are restricted to a maximum number of enriched nodes, since highly enriched meshes lead to numerical stability problems, as demonstrated by Góis and Proença (2012).

With this in mind, we adopted a mesh of 280 enriched nodes (<0.25n_elements), shown in Figure 11, which is chosen specifically for the solution of this problem. This number was chosen based on several tests involving different meshes, and varied regions of nodal enrichment, as well as the number of enriched nodes applied to the domain (Arruda, 2015). The layout of the enriched nodes was chosen in a way to cover most of the areas affected by the checkerboard problem in TOM simulation without sensitivity filters, creating areas with different densities of modified nodes.

Figure 11
MBB Beam non uniform nodal enrichment distribution (Mesh 2).

After running the simulation for both GFEM and SGFEM, it became possible to see considerable difference in the severity of the remaining checkerboard pattern problem. In both cases the amount of instability also decreases as penalization grows. Regarding TOM-GFEM, the results show a similar layout to classical TOM results with considerably less checkerboard pattern in every penalty value assessed, as seen in Figure 12. The topology results for TOM-SGFEM, shown in Figure 13, although present comparable results overall, show that the code was able to remove more of the checkerboard pattern compared to TOM-GFEM.

Figure 12
MBB Beam topology obtained by using TOM-GFEM and Mesh 2.

Figure 13
MBB Beam topology obtained by using TOM-SGFEM and Mesh 2.

For better comparison of the results, Figure 14 shows topologies obtained in each of the three different codes used in this work, organized by penalization. Since results for penalization equal to 1.0 present little visual difference, only topology results obtained from penalizations 2.0, 2.5 and 3.0 are shown. However, it is visible in Table 2 that there are differences in number of iterations and objective function value, which are directly caused by the application of the GFEM and SGFEM formulation, since it leads to an increase of the stiffness matrix and the displacement vector of the problem.

Figure 14
Comparison between MBB Beam results from all codes regarding Mesh 2.

To complement the analysis of the images in Figure 14, the data presented in Table 3 show the values of interest already mentioned. Besides the topology results, it is worth noting that in most cases both TOM-GFEM and TOM-SGFEM have a smaller number of iterations with few exceptions. Objective function values were always greatest for SGFEM; however, followed closely by the GFEM code. Both proposed codes found better results regarding iterations and objective functions.

Table 3
Data compilation for MBB beam - Mesh 2.

Penalization	TOM		TOM-GFEM		TOM-SGFEM
	Iteration	Obj. Function	Iteration	Obj. Function	Iteration	Obj. Function
1.0	20	164.4165	19	164.4956	19	164.4968
2.0	97	185.4375	114	190.5152	87	190.9066
2.5	47	191.0886	46	195.7316	51	196.1413
3.0	65	203.2379	35	206.0362	48	206.3735

In addition to numerical data, a close comparison between the three different results, obtained by each code makes, it evident the impact of the non-conventional FEM formulations in TO problem. Looking specifically to the last row in Figure 14, it is possible to see significant difference in the checkerboard pattern distribution between TOM and TOM-GFEM and also there is virtually no checkerboard pattern in TOM-SGFEM.

5.2 Fixed-Ends Beam

The second problem studied in this work presents the simulation and obtained results for the Fixed Ends Beam. In this model (Figure 15), the blue squares and triangles represent the boundary conditions (fixed nodes) of the problem, while the red arrow indicates the applied force.

Figure 15
Fixed Ends Beam finite element model.

The benchmark results for the Fixed Ends Beam obtained by the classical TOM c code are shown in Figure 16. As previously performed, those results are presented for the absence of the sensitivity filter. Analyzing the obtained topology results, it is possible to observe that a small value for the penalization factor leads to a uniform distribution of the checkerboard pattern, which is observed in the left region of the structure (p=2.0), while increasing the penalization value, more uniform topology and less regions containing the checkerboard pattern are obtained. Thus, much like the TOM results presented for the MBB beam, a direct relation between penalization and the amount of checkerboard pattern problem can be seen in Figure 16.

Figure 16
Preliminary topology results obtained by classical TOM code, without filtering technique, and used as benchmark.

The benchmark data for this problem can be found in Table 4. Much like the results from MBB beam, the greatest number of iterations comes with p = 2.0, and the higher value for the objective function is a direct consequence of higher penalization values, due to the decrease in severity on the checkerboard pattern as penalization grows.

Table 4
Benchmark data for Fixed-Ends Beam topology.

Penalization	Volume Fraction	Iteration	Objective Function
1.0	0.500	13	58.9521
2.0	0.500	101	67.8013
2.5	0.500	65	71.7995
3.0	0.500	41	74.5923

In order to perform the first simulation with nodal enrichment for a Fixed-Ends Beam, a similar mesh to that used in the first simulation for the MBB beam was chosen, with the only difference being the fixed nodes at extreme right, as shown in Figure 17. Besides that, the density and position of enriched nodes is exactly the same as Mesh 1, justifying carrying the name over to this problem.

Figure 17
Fixed-Ends Beam uniform nodal enrichment distribution - Mesh 1.

In this application of Mesh 1, results show little to no improvement in resolving the checkerboard pattern (Figure 18). Similarly, to the MBB beam problem, there is small change in topology for lower values of penalization, including a worsening in the checkerboard pattern for TOM-GFEM when p = 2.0. Besides that, results for penalization values of 2.5 and 3.0, for both new formulations, resemble very closely in appearance and start to show some improvement on the proposed problem, although not close to eliminating it yet.

Figure 18
Comparison between MBB Beam results from all codes regarding Mesh 1.

Regarding the data of interest, there is no significant change in any of the values other than number of iterations for p = 2.5. All results can be seen in Table 5. Even though the value of the objective function showed small increase in both TOM-GFEM and TOM-SGFEM, these were not indication of complete checkerboard pattern removal, especially considering that TOM-GFEM had the worst topology results for penalization equal to 2.0. It is worth noting that the number of iterations shows a pattern of increase for this same penalization (p=2.0), creating longer simulations when compared with other penalization values.

Table 5
Data compilation for Fixed-Ends beam - Mesh 1.

Penalization	TOM		TOM-GFEM		TOM-SGFEM
	Iteration	Obj. Function	Iteration	Obj. Function	Iteration	Obj. Function
2.0	101	67.8013	116	69.8454	96	69.3969
2.5	65	71.7995	46	72.5065	38	72.3141
3.0	41	74.5923	35	75.7611	42	75.6333

Thus, using the same observation method proposed for the MBB beam, a new and improved mesh (called Mesh 3) can be designed specifically for this problem, with a higher density of modified nodes closer to the center of the domain and a sparser distribution of nodes closer to extremities, as shown in Figure 19. For a better evaluation, the number of enriched nodes is the same as Mesh 2 with the same vertical disposition, aiming to cover the center of the domain, where the checkerboard pattern is more severe.

Figure 19
Fixed-Ends Beam non uniform nodal enrichment distribution (Mesh 3).

Results from both TOM-GFEM and TOM-SGFEM show an improvement when compared to results obtained with Mesh 1, and even better when compared directly to the non-enriched mesh. TOM-GFEM had interesting results regarding the positioning of the remaining checkerboard pattern. Each penalty value changed the position of the problem to a different focus, as seen in Figure 20. When worked with the SGFEM, the code was able to improve the final result, especially with p = 3.0 (Figure 21).

Figure 20
Fixed-Ends Beam topology obtained with TOM-GFEM and Mesh 3.

Figure 21
Fixed-Ends Beam topology obtained with TOM-SGFEM and Mesh 3.

When comparing results side by side, the differences can be seen more clearly. Figure 22 shows a compilation of results, with each row representing a penalty factor value and each column showing results from a different code. Both TOM-GFEM and TOM-SGFEM had similar results compared to each other; however, they showed excellent improvement when compared to classical TOM. For penalization equal to 2.5, both new codes had considerably better results with the SGFEM formulation being able to avoid the checkerboard problem more efficiently in the lower right corner and far left end of the beam. Regarding penalization equal to 3.0, the SGFEM formulation has shown the best result, being able to correct more of the problem in the upper left corner of the domain.

Figure 22
Comparison between Fixed-Ends Beam results for all codes regarding Mesh 3.

The data regarding the number of iterations and objective function values can be seen in Table 6. For penalization values equal or greater than two (p ≥ 2.0) there is significant increase in the value of the objective function with both new proposed methods. Taking into account the number of iterations, GFEM has better performance for lower penalizations, compared to classical TOM, and SGFEM needs more iterations than both other codes when p = 2.0.

Table 6
Data compilation for Fixed-Ends beam - Mesh 3.

Penalization	TOM		TOM-GFEM		TOM-SGFEM
	Iteration	Obj. Function	Iteration	Obj. Function	Iteration	Obj. Function
1.0	13	58.9521	12	59.0276	12	59.0272
2.0	101	67.8013	87	71.6633	142	71.5407
2.5	65	71.7995	72	73.4259	53	73.3993
3.0	41	74.5923	42	76.2234	39	76.1382

Overall, both proposed formulations for the TOM could drastically reduce the checkerboard pattern problem, when combined to a more carefully designed enrichment mesh rather than a uniform distribution of enriched nodes, indicating that positioning of enriched nodes leads to better results.

6 CONCLUSIONS

This paper studies the possibility of removing the checkerboard pattern problem for topology optimization design, using two unconventional formulations of the Finite Elements Method, GFEM and SGFEM. By analyzing the results obtained from classical TOM code, we could carefully design a distribution of enriched nodes for each contemplated problem. In this case, finite element meshes have a higher density of enriched nodes in areas most affected by the checkerboard pattern problem and enrichment of some remaining nodes along the rest of the domain.

Initially, a uniform distribution of enriched nodes was evaluated, as a base for comparison. These results for both TOM-GFEM and TOM-SGFEM in both problems were inferior to the results obtained for the specifically designed meshes. Although the code with the uniform distribution of enriched nodes was able to remove some of the checkerboard pattern, it did not cause significant impact, with a considerable area still covered by the checkerboard pattern in both examples. Even considering numerical data such as number of iterations and value of the objective function, changes are minimal when compared to the more conventional formulation of TOM.

Regarding the non-uniform meshes, results are promising, with the topology obtained using TOM-SGFEM and penalization p=3.0 being the best in terms of solving the checkerboard problem in all proposed examples. In the MBB beam simulations, both new codes performed significantly better than classical TOM and the uniform meshes, more specifically when penalization numbers were higher. TOM-SGFEM showed virtually no checkerboard pattern. Considering the Fixed-Ends beam, results were similar to the first example, even though both new formulations showed a small amount of checkerboard pattern in the upper left corner of the final topology. It is worth noting that TOM-SGFEM provided less computational cost when compared to TOM-GFEM, since the former used I(j)=2, while the latter has I(j)=5.

These results indicate that both new methodologies are highly sensible to the density and localization of enriched nodes. With this in mind it is possible to assume that better results could be obtained with even more specialized meshes designed for each example. It is also worth noting that, as expected, the new algorithms took significantly longer to complete the optimization loops, due to increased number of variables and creation of new matrices for enriched functions and placement of the modified nodes. Even then, results for number of iterations and value of objective function showed improvement, most of the times, when TOM was combined with both new formulations. However, these codes had to be written according to visual parameters observed in topologies obtained from classical TOM in each example.

In summary, the new proposed algorithms (TOM-GFEM and TOM-SGFEM) have been able to improve very significantly the topology results in both examples when compared to results obtained by classical TOM without sensibility filters.

References

• Arruda, L. S. (2015). Application of the Generalized Finite Element Method (Mefg) in Topological Optimization. 158f. Dissertation (MSc in Mechanical Engineering) in Portuguese, Federal University of ABC, Santo André.
• Bendsøe, M. P. (1989). Optimal shape design as a material distribution problem. Structural optimization, 1(4), 193-202. https://doi.org/10.1007/BF01650949
  » https://doi.org/10.1007/BF01650949
• Babuška, I., Banerjee, U. (2012). Stable generalized finite element method (SGFEM). Computer methods in applied mechanics and engineering, 201, 91-111. https://doi.org/10.1016/j.cma.2011.09.012
  » https://doi.org/10.1016/j.cma.2011.09.012
• Bathe, K. J. (1996). Finite Element Procedures. 2^nd Ed. New Jersey: Prentice Hall.
• Bendsøe, M. P., & Kikuchi, N. (1988). Generating optimal topologies in structural design using a homogenization method. Computer methods in applied mechanics and engineering, 71(2), 197-224. https://doi.org/10.1016/0045-7825(88)90086-2
  » https://doi.org/10.1016/0045-7825(88)90086-2
• Bendsøe, M. P., & Sigmund, O. (2003). Topology optimization: theory, methods, and applications. Springer Science & Business Media.
• Christensen, P. W., & Klarbring, A. (2009). An introduction to structural optimization (Vol. 153). Springer Science & Business Media.
• Diaz, A., & Sigmund, O. (1995). Checkerboard patterns in layout optimization. Structural optimization, 10(1), 40-45. https://doi.org/10.1007/BF01743693
  » https://doi.org/10.1007/BF01743693
• 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. https://doi.org/10.1016/S0045-7949(99)00211-4
  » https://doi.org/10.1016/S0045-7949(99)00211-4
• Duarte, C. A., & Oden, J. T. (1995). A new meshless method to solve Boundary-Value Problems. In Proceedings of the XVI CILAMCE-Iberian Latin American Conference on Computational methods for engineering, Curitiba, Brazil (pp. 90-99).
• Góis, W., & Proença, S. P. B. (2012). Generalized Finite Element Method on Nonconventional Hybrid-Mixed Formulation. International Journal of Computational Methods, 9(03), 1250038. http://dx.doi.org/10.1142/S0219876212500387
  » http://dx.doi.org/10.1142/S0219876212500387
• Gupta, V., Duarte, C. A., Babuška, I., & Banerjee, U. (2013). A stable and optimally convergent generalized FEM (SGFEM) for linear elastic fracture mechanics. Computer methods in applied mechanics and engineering, 266, 23-39. https://doi.org/10.1016/j.cma.2013.07.010
  » https://doi.org/10.1016/j.cma.2013.07.010
• He, Z. C., Zhang, G. Y., Deng, L., Li, E., & Liu, G. R. (2015). Topology optimization using node-based smoothed finite element method. International Journal of Applied Mechanics, 7(06), 1550085.https://doi.org/10.1142/S1758825115500854
  » https://doi.org/10.1142/S1758825115500854
• Mlejnek, H. P. (1992). Some aspects of the genesis of structures. Structural optimization, 5(1), 64-69. https://doi.org/10.1007/BF01744697
  » https://doi.org/10.1007/BF01744697
• Navarro, L. A. P., Góis W. (2017). A MATLAB® implementation of Hybrid-Mixed Stress Formulation on Topology Optimization. In: XXXVIII IberianLatin American Congress on Computational Methods in Engineering 2017.
• Rozvany, G. I., Zhou, M., & Birker, T. (1992). Generalized shape optimization without homogenization. Structural optimization, 4(3-4), 250-252. https://doi.org/10.1007/BF01742754
  » https://doi.org/10.1007/BF01742754
• Shobeiri, V. (2016). Structural topology optimization based on the smoothed finite element method. Latin American Journal of Solids and Structures, 13, 378-390. https://doi.org/10.1590/1679-78252243
  » https://doi.org/10.1590/1679-78252243
• Sigmund, O. (2001). A 99 line topology optimization code written in Matlab. Structural and multidisciplinary optimization, 21(2), 120-127. https://doi.org/10.1007/s001580050176
  » https://doi.org/10.1007/s001580050176
• Sigmund, O., & Petersson, J. (1998). Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural optimization, 16(1), 68-75. https://doi.org/10.1007/BF01214002
  » https://doi.org/10.1007/BF01214002
• van den Boom, S. J., Zhang, J., van Keulen, F., & Aragón, A. M. (2021). An interface-enriched generalized finite element method for level set-based topology optimization. Structural and Multidisciplinary Optimization, 63(1), 1-20. https://doi.org/10.1007/s00158-020-02682-5
  » https://doi.org/10.1007/s00158-020-02682-5

Edited by

Editor: Pablo Andrés Muñoz Rojas

Publication Dates

• Publication in this collection
  23 May 2022
• Date of issue
  2022

History

• Received
  08 Nov 2021
• Reviewed
  15 Mar 2022
• Accepted
  03 May 2022