Sequential method of topological optimization in multi-component systems

Abstract Topology optimization research has focused on structures of a single domain or component. The single component configuration fails to capture the complexity of real multi-component structures. It is necessary to develop new methods and numerical strategies to solve multi-component systems. In this work, we propose a new approach considering a sequential method of topological optimization in multi-component systems. The proposed algorithm for topology optimization of multi-component systems is sequential. It optimizes the first component, the result found is included in the analysis of the optimization of the next component and thus it continues analyzing consecutively all the components until the last one. This methodology is the only one that performs sequential optimization using open-source tools. The implementation of sequential method is developed in four processes: development of the multi-component mesh, numerical structural analysis through the FEM, sensitivity analysis and a final optimization. A comparison is made with an optimization of multi-component systems in a normal way using three commonly seen examples. Graphical Abstract


INTRODUCTION
Many structures are composed of more than one structural component, known as multi-component.Examples include trusses, automotive systems, machinery, aircraft components, and many others as illustrated in Figure 1.Real structures are usually composed of several components that are interconnected with the aim of improving transportation or manufacturing, these are generally composed by the "assembly" of several individual components or multicomponent structures.However, optimizing the topology of multi-component elements can be a challenging, since the division of structures into their individual components may not be sufficient to consider the interaction between components.Thus, topology optimization in multi-component structures is still an extensive field to be explored.In general, the conventional process of analysis of multi-component systems involves the consideration of the components individually, in order to simplify the process of modeling and designing the project (LI, STEVEN, et al., 2001).Thus, considering this individual analysis of each component, the connection elements are treated as a load transfer system, and take into account a variety of simplifications of kinematic constraints that may not reflect the real conditions.This strategy can limit the accuracy of the analysis and may not capture all interconnected effects between components.The load transfer system that is imposed on predefined connection locations between components may require that each component be sized taking into account the constraints and constraints of the component set.However, evaluating loading and boundary conditions for all multi-component systems may not be clear or sufficient for the effective development of multi-component systems.To ensure the desirable performance of multi-component systems, it is necessary to incorporate the layout of the various connected elements and their components into topology optimization.This will allow the interactions between components to be more accurately considered and optimize the performance of the system as a whole.
A direct approach to solve this problem of optimizing the topology of multi-component elements is proposed by Chickermane and Gea (1997).In this approach the joints between the components are modeled as a grid of spring elements that transfer force potentially between the nodes of two parts.However, this straightforward approach may not capture all physical interactions between elements.
Latin American Journal of Solids and Structures, 2023, 20(6), e501 3/18 According to recent research carried out by Thomas, et al., (2020), several studies have addressed multi-component optimization on periodic elements, but using single component optimization with element decomposition.However, other works also used the direct decomposition of elements, as is the case of Menassa and Devries (1991), which inserted fixed positions between elements, and (CHIREHDAST, JIANG, 1996), which established fixed weld points.In Jiang and Chirehdast (1997), adhesive bonding elements were considered, which generated a lot of assumptions about the dynamics of movement between the components.Even in works such as Li, et al., (2001) which considered fixed elements with discrete joints in three dimensions, load transfers were still performed through theoretical assumptions between components.
There are optimization methods that combine topology and layout, using component location variables as part of the design variables to allow flexible placement of connecting elements.In the work of ZHU, et al., (2017) is incorporated this approach into a final topology optimization framework, and the method was refined by Ambrozkiewicz And Kriegesmann (2021), who included predefined fixed-form connector elements to model force transfer between components and the structure support.These element fasteners are connected to the structure using coupling equations to avoid the need for meshing when components move.More recently, Rakotondrainibe, et al., (2020) modeled rigid supports together using movable boundary condition zones, whose positions were part of the optimization in a defined level set topology optimization framework.
The sequential methodology of topological optimization in multi-component systems or sequential method, is a new approach that is based on the work carried out by Ferro and Pavanello (2022).In this work, a complete demonstration of structural topological optimization in multi-component elements is made.It shows the advantages of including the mesh insertion of a 3D environment, then the possibilities of reading the mesh by the environment to solve problems using FEM, how to use the methodology of the adjoint solution and finally the optimization.An advantage of this approach is, in fact, the ease of integration between the creation of the 3D mesh and the consideration of all components and their connections, including direct contact between their elements, with no simplification between domains.Also, another advantage is the fact that it is possible to choose which domain(s) will be considered as a volume restriction in the optimization and this falls directly into the creation of the new sequential method approach.And this even considering the entire numerical process in an open-source environment.
In the new approach of the sequential method, considering the advantages of the normal approach, it is possible to insert results obtained individually from each domain and later insert them to calculate the optimization of the next domain sequentially.And so, verifying the gains of the new approach, where the best results for Compliance are graphically demonstrated, with even smaller values obtained, the optimization process is more stable and the processing time is not necessarily much longer.Which can be very feasible for situations with many components being considered.And finally, in the sequential method it is also possible to consider the volume constraints individually in each component and thus showing the robustness and diversity of possibilities in the sequential method.

OPEN-SOURCE SOFTWARES
With the advancement of technology and the spread of numerical methods, open-source software has become more improved and widely disseminated.One such software is FEniCS (ALNAES, BLECHTA, et al., 2015, LOGG, WELLS, et al., 2011), which offers excellent resources for partial differential equations using the finite element method.Dolfin-Adjoint (MITUSCH, FUNKE, et al., 2019), a complementary package to FEniCS, is used to calculate the model sensitivities, i.e., the derivatives of objective functions and volume constraints.To solve the constrained minimization problem, we can use the Ipopt (Interior Point OPTimizer) (WÄCHTER, BIEGLER, 2006).The optimization solver is a largescale nonlinear optimization package that finds local solutions.The Gmsh software (REMACLE, GEUZAINE, 2007) is used to generate finite element meshes in three dimensions, with an integrated CAD engine and a post-processor.
Thus, all work is done using the following sequence: Mesh generation with Gmsh, Finite Element Analysis with FEniCS, Sensitivity analysis with Dolfin Adjoint and Optimization with Ipopt.

Mesh generation: Gmsh software
Gmsh is a software that has a combination of three-dimensional mesh generator, CAD and post-processor (REMACLE, GEUZAINE, 2007).In this work , we will only use the three-dimensional mesh generation functionality, which will be later converted for reading in FEniCS.
To define a model in Gmsh, its boundary representation (BRep) is used, where a volume is bounded by a set of surfaces, a surface is bounded by a series of curves, and a curve is bounded by two endpoints.The model is composed of topological entities that deal only with the adjacencies in the model, and are implemented as a set of abstract Latin American Journal of Solids and Structures, 2023, 20(6), e501 4/18 topological classes.BRep can be extended by defining internal or embedded entities in the model: internal points, edges and surfaces can be embedded in volumes; and interior points and curves can be embedded in surfaces.
One of the great advantages of using Gmsh is that its geometry structure allows defining each volume of each domain, including the connection elements, and converting them for reading in FEniCS according to the mathematical process of the topological optimization method.And according to the types of meshes that can be read by FEniCS, the finite element models used in the work use tetrahedral linear finite elements.

Finite element method: FEniCS
The utilization of the FEniCS platform was crucial for the implementation of the finite element method, following a functional analysis-based approach (LANGTANGEN, MARDAL, 2019).This software is capable of solving Partial Differential Equations (PDEs) and simplifying the implementation process by utilizing the abstract formulation of the method, which is considered one of the two primary approaches of the finite element method, with the other being the engineering structural analysis formulation (LANGTANGEN, LOGG, 2016).FEniCS is one of the tools based on the first approach, as well as other works by authors such as Kirby ( 2007) and Larson and Bengzon (2013).
Using the GALERKIN-type weighted residuals method ("CG" or Continuous Galerkin), it is possible to write the formulation of the PDEs solution considering a nodal polynomial approximation by subdomains.The test functions used in the FEniCS programs are defined according to the spaces of functions that specify the properties of the approximations of the numerical values adopted (LANGTANGEN, LOGG, 2016).To reduce the number of degrees of freedom of the models, only Lagrange CG1 tetrahedral elements are used both for the analysis of the displacements and for the filtering system.

Sensitivity analysis: dolfin adjoint
The adjoint method is a mathematical technique used in optimization to efficiently calculate the gradient of an objective function with respect to a set of design variables.This method is used in dolfin-adjoint software, which is a tool for solving PDEs and optimizing their solutions.
The adjoint system is the same idea as the reverse algorithmic or automatic differencing mode.Where the term adjunct is also used in the literature as dual.
Thus, as seen in Farrell, et al., (2013) and Mitusch et al., (2019) one starts by fixing the choice of parameter , and then solves the Jacobian solution d d ⁄ associated with choosing .Then the inner product is made with a term   ⁄ and thus the gradient d ̂d ⁄ is calculated.According to Farrell, et al., (2013) and Mitusch, et al.,( 2019), the gradient calculation process in an optimization problem with PDE constraint starts with the fixed choice of the parameter , followed by the resolution of the Jacobian solution d d ⁄ corresponding to that choice.Then, the inner product is made with a term   ⁄ , in order to calculate the gradient d ̂d ⁄ .

Optimization: Ipopt
Ipopt is a free non-linear optimization software.It uses the interior point method to solve continuous optimization problems, subject to nonlinear and linear constraints (ALONSO, DE SÁ, et al., 2019, DE SOUZA, SILVA, 2020) .
The interior point method (ANDREI, 2022, MAAR, SCHULZ, 2000) is an optimization method that consists of iteratively approximating the solution of the problem in a series of steps.Each step consists of determining a search direction that improves the current solution.This direction is determined by solving an auxiliary problem, which is a quadratic approximation of the original problem.

STRUCTURAL TOPOLOGY OPTIMIZATION IN MULTI-COMPONENT ELEMENTS
To perform the structural topological optimization in multi-component elements, the problem of minimizing the structural compliance (mean sensitivity) is considered, through an objective function over the considered domains, according to Equation 1 (BENDSOE, SIGMUND, 2003, CHICKERMANE, GEA, 1997).
where  is the Cauchy stress tensor,   is the body forces,  represents the external normal direction in Γ  and ∇.  is the stress divergence.
And the following additional restrictions: Since the parameter () is the design variable (() = 1 means presence of material and () = 0 means absence of material),   , in Equation 4, is the desired final volume of each component.The strain displacement relations considering the small elastic deformations of a multi-component Ω  can be written as: (5) The constitutive equations can be written as: and the linear integral form of the equilibrium conditions can be written as: where u i is the displacement field of each domain,   () is the linear part of the Green strain stresses of each domain,  is the trace and  is the identity tensor.The Lamé constants are  and  according to: being in Equation 8,  the Young's modulus and  a Poisson coefficient.With this, the variational formulation is summarized on how to find , total displacements considering the influence of all domains, so that  ∈  where: The method presented in this study uses the SIMP technique -Isotropic Solid Material with Penalization (SIGMUND, 2001) to optimize the structural topology of multi-component elements simultaneously and individually.This method allows an ideal material distribution for each component, considering the loads imposed on the structure as a whole, without simplifications.In this way, it is possible to carry out a precise optimization of the distribution of materials in each component of the structure.
Thus, a continuous SIMP relaxation is performed, proposed by Bendsøe and Sigmund (2003) according to:

DENSITY AND PROJECTION FILTERS
The choice of the power law for the interpolation of materials, according to Equation 11, can lead to numerical problems of the "checkerboard" type with alternating solutions between () = 0 and () = 1 (DÍAZ, SIGMUND, 1995, JOG, HABER, 1996, SIGMUND, PETERSSON, 1998).To mitigate these problems, it is common to use a hearing density filter implicitly by solving a Helmholtz-type partial differential with Neumann boundary conditions (LAZAROV, SIGMUND, 2011) as described in Equation 12: The continuous representation of the unfiltered design variable is (), while the filtered design variable is  � (), and the parameter   has a role similar to   in classical approaches to filtering by the SIMP method.In the study by (Lazarov and Sigmund, 2011), an approximate relationship between the length scales for the classical filter and the Helmholtz approach is described as   =   2√3 ⁄ , where the parameter   is defined by the solution of a partial differential equation of the Helmholtz type with homogeneous Neumann boundary conditions, according to Equation 12.The solution of the equation can be expressed in the form of an integral convolution, which is equivalent to the classical filter.In terms of variational formulation, we have: being  � Hilbert subspace of  admissible functions.The density filtering method results in non-binary values between 0 and 1, and to obtain a smoother transition between these values, the Heaviside Projection is used.This technique consists of projecting the intermediate densities to their minimum and maximum limits.In (WANG, LAZAROV, et al., 2010), this projection is performed by calculating a new density field,  �(), as follows: The  parameter is used in the intermediate density projection and is responsible for weighting the intermediate values.
The  parameter is responsible for performing the projection.

NUMERICAL IMPLEMENTATION
The equilibrium equations described by Equation 9 are solved using the finite element method, through the FEniCS framework.This software converts models described by variational forms using high-level language and automatically generates the corresponding finite element code.The sensitivities are calculated using the adjoint method, using the Dolfin Adjoint package, which is capable of automatically deriving the adjoint and tangent discrete linear models from a direct model written in the Python interface for FEniCS.To solve the optimization problem presented in Equation 10, optimization routines from the Ipopt library are used, which is capable of dealing with large-scale nonlinear problems, through the implementation of an "interior point" and a "line search".The analysis and optimization methods employed are suitable for handling large problems, with up to millions of variables and constraints.Figure 2 presents a general scheme of the numerical implementation used for a normal optimization including all components or domains of a considered multi-component structure.In Figure 3, the new scheme of sequential optimization is presented, where the use of the previous result in the subsequent result of each component or domain is considered.The sequential optimization sequence follows in Figure 4.

MODELS ANALYZED
To illustrate the application of the proposed approach, some common models in the literature will be generated and analyzed.One of the models is the Cantilever Beam with two domains, which is often used in structural topology optimization studies in multi-component elements, and presents the overlapping domains and their connection elements.Superimposed multi-component models, like the one shown in Figure 1 in 3D analysis, or models with even greater thickness, can lead to significant distortions in the final results due to the eccentricity of the load in relation to the components.However, this distortion is not as evident in 2D analyzes and is often overlooked in the cited studies, which generally focus on 2D analyses.
With the aim of making the models more realistic and circumventing the problems resulting from distortion in the final results caused by the eccentricity of the load in relation to the components, this work also proposes the inclusion of an additional connection element, such as a connection plate, to promote concentric loading between the multicomponents in the models.The numerical approach also allows choosing which domains can be optimized, only changing the need for volume constraint functions.To make the simulation even more realistic, the connection elements will be modeled as pins or hollow screws, reducing the rigidity of the connection and making it more flexible, improving the applicability in the analysis conditions in each individual element.The models of the individual components and their connections are considered with a material with modulus of elasticity  = 210. 10 9 and Poisson coefficient of  = 0.3.All individual components have the following dimensions: length: 3, height: 1 and width: 0.02 In all examples the differences between sequential optimization and normal optimization are shown.And due to the differentiated form of mesh creation presented here, with the use of external mesh generation and compatible with the solution to the problem in an FEM environment, there are few works for comparative purposes.In the work done by Ambrozkiewicz and Kriegesmann (2021), in its main example with two components and two connection elements, only the final result of Complince is shown, however, even so, it has a result similar to the one developed here.
In the sequential method it is also possible to consider the volume constraints individually in each component.As demonstrated by Ferro and Pavanello (2023) using a subclass of Dolfin-Adjoint, it is possible to consider the desired volume for each component individually according to Equation 4. And thus, a comparison considering variations between volumes is also shown at the end of each model.

Cantilever with two domains
The first example is the Cantilever with two domains, in which a multi-component with two overlapping beams and two fixed connecting elements was created -the left domain is fixed to the left face, while a load is applied to the middle part of the right face of the right domain.Volume reduction to 50%.The finite element model uses tetrahedral linear finite elements with 72.822 elements (FEniCS uses only tetrahedral mesh in 3D models).
In Figures 5 and 6 is the result for normal optimization, that is, being performed at the same time in all domains.And in Figures 7 and 8 the sequential optimization method is shown.An important analysis in structural optimization is the evaluation of the objective function or the Compliance, as a function of the iteration compared to the reduction of the volume of the structure.Thus, Figure 9 shows this result comparing the Compliance between a structure with normal optimization and a sequential optimization.After viewing the results obtained for both normal optimization and sequential optimization, it is observed that the final results are very close, presenting a variation of approximately 4% with better performance for sequential optimization.When comparing the time spent on numerical processing, the normal optimization is done in 151 seconds and the sequential one in 380 seconds in this example.Considering that for the sequential optimization there are ttwo sequential optimizations, this longer time is justified.And even with the additional time, even so, if a model with a larger mesh is considered, the final results will be even better for the sequential method, showing its ability to obtain better results.

Cantilever with two domains -different volume constraints
Now, the results will be shown considering that the volumes of each component can be individualized using the sequential method.
According to Figure 5 a) considering the Cantilever with two domains, the volumes of Components 1 and 2 will be changed.Figure 10 shows the result for the Volume of 30% for Component 1 and 50% for Component 2. And Figure 11 is showing the result for the Volume of 50% for Component 1 and 30% for Component 2.   After viewing the results obtained for both normal optimization and sequential optimization, in this model, it is also observed that the final results are very close, presenting a variation of approximately 4.4% with better performance for sequential optimization.When comparing the time spent on numerical processing, the normal optimization is done in Latin American Journal of Solids and Structures, 2023, 20(6), e501 13/18 131 seconds and the sequential one in 576 seconds in this example.In this case, sequential optimization is performed in three sequential optimizations, justifying the longer time.

Cantilever with two domains and a connecting plate -different volume constraints
The results will be shown considering that the volumes of each component can be individualized using the sequential method.
According to Figure 12 a) considering the Cantilever with two domains and a connecting plate, the volumes of Components 1 and 2 and a connecting plate will be changed.Figure 17 shows the result for the Volume of 50% for Component 1 and 2 and 30% for a connecting plate.And Figure 18 shows the result for the Volume of 30% for Component 1 and 2 and 50% for a connecting plate.After viewing the results obtained for both normal optimization and sequential optimization, in this model, it is also observed that the final results are very close, presenting a variation of approximately 1% with better performance for sequential optimization.In this model, the sequential optimization is also performed in three sequential optimizations, justifying the longer numerical processing time, with the normal optimization performed in 126 seconds and the sequential optimization in 573 seconds.It is observed here that in the first iterations, there is a low variation between the normal and sequential optimization, but as the process continues, the improvement in the result in the sequential optimization is evidenced.

Cantilever with three domains -different volume constraints
The results will be shown considering that the volumes of each component can be individualized using the sequential method.
According to Figure 19 a) considering the Cantilever with three, the volumes of Components 1, 2 and 3 will be changed.Figure 24 shows the result for the Volume of 30% for Component 1 and 50% for a Components 2 and 3.And Figure 25 shows the result for the Volume of 30% for Components 1 and 2 and 50% for a Component 3. At first, the integration between the different platforms used may seem expensive, but the learning progress in this approach has a smooth curve over the works and models performed.Understanding the mesh generation process by Gmsh and how to make it available, within the multi-component elements approach, both for the finite element analysis performed by FEniCS and for the Dolfin-Adjoint sensitivity calculation, is a complex process and shown in this work integrated into all programs.
The great advantage of the numerical method developed in this work is the sequential use of results obtained for a component and inserted into the analysis to obtain the result of the next domain and so on until the last domain.
When comparing the results obtained here with well-known works such as those seen in Thomas, et al., (2020) that use commercial packages or even more recent works such as Ambrozkiewicz snd Kriegesmann (2021), we see that the results are very close, showing the viability of the codes generated here.However, this work is the first that makes a sequential analysis of its components individually.
Finally, the results obtained showed that the sequential method, even with a longer processing time, justified by the fact that there are a series of optimizations until reaching the final result, still proves to be advantageous in relation to the normal method, where in the final result, structures have the lowest compliance.Looking from the computational point of view, the increase in time is of the order of N components, therefore, for large amounts of components, or even models with larger meshes, the sequential method is even more advantageous if it is necessary to optimize components individually.

Figure 1
Figure 1 Multi-component structure 9) being  � the Hilbert subspace of  admissible functions.Isotropic elastic conditions are assumed in the examples of this work.The update material distribution () is done by updating the material stiffness through compliance minimization (Bendsøe and Sigmund, 2003) which can be written as follows: min.() s.a: (, ) = () ∀ ∈  � Latin American Journal of Solids and Structures, 2023, 20(6

Figure 2 Figure 3 Figure 4
Figure 2 Numerical implementation routine for Normal Optimization

Figure 9
Figure 9 Cantilever with two domains: Compliance/Volume Fraction curve versus Iteration

Figure 10 Figure 11
Figure 10 Cantilever with two domains (Sequential Optimization): Optimized Side Views -Volume of 30% for Component 1 and 50% for Component 2

Figure 12 18 Figure 14 Figure 15
Figure 12 Cantilever with two domains and a connecting plate (Normal Optimization): a) Base Geometry, b) Optimized Geometry

Figure 16
Figure 16 Cantilever with two domains and a connecting plate: Compliance/Volume Fraction curve versus Iteration

Figure 17 Figure 18
Figure 17 Cantilever with two domains and a connecting plate (Sequential Optimization): Optimized Side Views -Volume of 50% for Component 1 and 2 and 30% for a connecting plate

Figure 19 Figure 20 Figure 21 Figure 22
Figure 19 Cantilever with three domains (Normal Optimization): a) Base Geometry, b) Optimized Geometry

Figure 23
Figure 23 Cantilever with three domains: Compliance/Volume Fraction curve versus Iteration

Figure 24 18 Figure 25
Figure 24 Cantilever with three domains (Sequential Optimization): Optimized Side Views -Volume of 30% for Component 1 and 50% for a Components 2 and 3 (BENDSOE, SIGMUND, 2003)l forces of volume forces   and surface forces   applied to each domain,  represents the number of domains analyzed,  represents the subindex of each domain,  is the individual domains,  represents the test functions, and Γ   is the portion of the contour where the surface forces of each domain are applied.Equilibrium conditions(BENDSOE, SIGMUND, 2003)are considered as constraints of the problem and are given per: Latin American Journal ofSolids and Structures, 2023, 20(6), e5015/18where (