A multi-point constraints based layout design of multi-component systems

– The primary purpose of this paper was to present an efﬁcient method to implement the layout design of multi-component systems. In this study, two kinds of design variables: topology variables and geometry variables are simultaneously optimized to maximize the structure stiffness. The multi-point constraints (MPC) are used to simulate the connection behavior between these movable components and supporting structure. So during the optimization iterations, mesh regeneration caused by moving components is avoided compared to the linkage of nodal coincidence. At the same time, the precise geometry shapes and properties of components are maintained to reduce the model building error. At the last of the paper, the numerical example is presented to show the validity of this method.


Introduction
Topology optimization for continuum structures came out since the original idea of homogenization-based design method was proposed in the work of Bendsøe and Kikuchi [1]. More efforts including Solid Isotropic Material with Penalization (SIMP), Rational Approximation of Material Properties (RAMP) and level-set method were made by Bendsøe and Sigmund [2], Stolpe and Svanberg [3] and Michael et al. [4]. Up to now, topology optimization has received much success both in follow-up theoretical study and practical applications. Now topology optimization is still developing rapidly to solve more complicated problems with the help of computing technology. Among those problems, the layout design of multi-component systems is to find not only the configuration of supporting structure but also the proper mounting positions of several components with predefined geometrical and physical properties. The idea can be expressed in Figure 1.
Qian and Ananthasuresh [5], Xia et al. [6] and Zhang et al. [7] introduced different material interpolations including Heaviside function, level-set function and X-FEM to describe the components. Their methods hold the advantages of fixed mesh and analytical sensitivity analysis. However, the precise geometry properties cannot be ensured. Zhu et al. [8,9,10] developed the techniques of density point and embedded mesh to update the pseudo-density variables. At the same time, the finite circle method (FCM) is introduced to avoid the overlapping problem caused components and design domain boundaries. This method can keep the precise geometry properties of com-ponents but leads to element remeshing and difficulties in sensitivity analysis.
In our paper, we use the multi-point constraints to build the connection between components and structure, and its theory and advantages will be introduced below.
2 Formulation of layout design of multi-component systems

Formulation of MPC
Firstly, the design domain and components should be discretized into finite element, respectively, as shown in Figure 2. Assume P 1 is one of the connecting nodes of the components, and is projected inside the design domain element e 1 . Then we enable the following MPC equation.
where u P1 denote the displacement vectors of node P 1 . u e1 and N e1 are the displacement vectors of the element e 1 and the shape function coefficient matrix at the position where P 1 locates, respectively. The equation above indicates the displacements of node P 1 is consistent with the point where it locates in the structure element. When it involves several connecting nodes, the number of MPC equations will increase as a result. Actually, these MPC equations can be regarded as a linear combination of node displacements, so we can write it as the following form.
where H is the coefficient matrix which is determined by shape functions of structure elements and coordinates of connecting nodes. u denotes the displacement vector of the structure system with components.
Considering the above displacement constraints, we use the Lagrange multiplier method to revise the overall potential energy of the global system where K is the global stiffness matrix. F is the nodal loads. k is the Lagrange multiplier vector. Then we have the Euler conditions for a stationary point The displacement vector u and k can be finally obtained by solving the above equations.
These constraint equations show much convenience since remeshing of structure elements around components is not needed and global stiffness matrix re-assembling process can be avoided when components move to different positions during optimization iterations. That brings high computing efficiency during optimization.

Optimization model
The goal of the optimization is to minimize the global strain energy of the structure system with a prescribed material volume fraction constraint. The design variables are the pseudodensity variables which control the material distribution and the geometry variables which describe the components' status (locations and orientations). Extra geometrical non-overlapping constraints have to be imposed to avoid the components' overlapping and keep all the components inside the design domain.
The optimization problem to minimize the strain energy of a multi-component system can be mathematically elaborated as: Find: g i ; i; ¼; 1; 2; . . . ; N d ; n jx ; n jy ; n jh À Á ; j; ¼; 1; 2; . . . ; N c ; where N d and N c are the numbers of the pseudo-density variables and the components, respectively. n jx , n jy and n jh are the geometry design variables i.e., the location and orientation related to the j th component. C is the strain energy of the structure system and K is the global stiffness matrix. V is the material volume fraction of the design domain with an upper limit of V U . X j , X j1 and X j are the areas of the jth, j1th and j2th components, respectively. X d denotes the area of the global design domain. Different from previous mentioned material interpolation scheme SIMP or RAMP, here we use a polynomial material interpolation in the work of Zhu et al. [9].
where g i denotes the ith pseudo-density. E(i) stands for the corresponding Young's modulus of one element, while E 0 (i) represents Young's modulus of a solid element. p is the penalty factor and a is polynomial coefficient factor. In this study, p and a are set to be 4 and 1/16. Component layout design without non-overlapping constraints may lead to an unreasonable solution; even fails in optimization process especially when the number of components is  very large. To deal with this problem, the finite circle method (FCM) is considered to be an effective mean. The main idea of FCM is illustrated in Figure 3: a series of circles are proposed to approximately describe the shape of two components and the rectangle design domain, so the overlap judgment of the two components or those of components and the structure design domain can be analytically evaluated by distance of the circles according to their center coordinates and radii.

Sensitivity analysis 3.1 Sensitivity analysis for pseudo density variables
Considering the static equation in equation (5), the differentiation with respect to the pseudo-density variable g i can be written as: Assume F are design independent loads, so oF/og i = 0, the equation below holds Then we have the derivative of the overall strain energy Substituting equations (8) and (9) into equation (10) and recalling oH T /og i = 0, we yield: As u T H T = 0, equation (11) can be simplified as: where the derivatives of the stiffness matrix K can be easily obtained according to the material interpolation model for the element stiffness.

Sensitivity analysis for geometry design variables
According to the static equation of equation (4) and noting oF/on j = 0, the differentiation with respect to one geometry variable n j can be written as: So the derivative of the strain energy with respect to the geometry design variable n j is similarly written as Suppose n j is a translational variable of the jth component, both the stiffness matrix of the supporting structures and the components will remain unchanged after a translational moving. So we yield If n j is a rotational variable, the stiffness matrix of the jth component after the rotation can be expressed as K Ã j is the initial stiffness matrix of the jth component and A j is the corresponding rotational transformation matrix. Suppose u j is the displacement vector of the jth component, then we can obtain The derivative of the rotational transformation matrix can be easily calculated according to the coordinate transformation matrix.
Consequently, equation (14) can be written as: Now the final geometry sensitivities can be concluded as: Note that H consists of the shape function coefficient matrix of the corresponding elements, which depend upon the connection positions and further the geometrical design variables describing the components' locations and orientations.

Numerical test
As displayed in Figure 4, a rectangle structure with a size of 1.8 m · 0.6 m · 0.02 m acts as the design domain and is divided into one layer of 90 · 30 shell elements. The loads and boundary conditions applied to the structure are also shown in Figure 4, and remain unchanged in the optimization iterations. Two same complex components are bonded onto the design domain elements by using MPC through the connection nodes locate on in the centers of the holes also as in Figure 4. The material properties of the components are defined as: E c = 2 · 10 11 MPa and m c = 0.3. While the material property set for the supporting structure is E s = 7 · 10 10 MPa and m s = 0.3.
With a volume fraction upper limit constraint of 0.5 and non-overlapping constraints, final design is achieved after 55 iterations using GCMMA (Globally converged method of moving asymptotes) method. As shown in Figure 5a-5f, the components work positively as an inseparable part of the  supporting structure to be subjected to loads. The iteration history of the objective function is shown in Figure 6. The stain energy value decreases from initial 0.00328J to final 0.000401J.

Conclusion
In this paper, a simultaneous design method for structure system with multi-components is proposed. The connections are automatically generated with the help of MPC between supporting structure and components. Based on that connection type, fixed mesh, analytical sensitivities and precise geometry of components can be achieved at the same time. The test under static loads is numerically implemented. Through the optimal results above, the method developed for optimization problems of multi-component systems is evidently proved valid and reliable. Some engineering applications of large scale can be expected in the near future using that method.