In this paper, a framework for topology optimization of assemblies considering unilateral contact- constraints is introduced. Here, an assembly consists of two or more components that are connected via joints, and the topology of all components of the assembly are optimized simultaneously. There- fore, linear-elastic material is considered and the augmented Lagrangian approach imposing Signorini’s contact condition is used. In doing so, the unilateral characteristic of contact is taken into account. Fur- thermore, the influence of different contact modeling techniques, i.e. augmented Lagrangian approach, a standard penalty approach and a unilateral penalty approach, on the optimized designs are dis- cussed. For topology optimization, the SIMP-approach is applied to all members of the assembly and sensitivity filtering is used. In consequence, the global design space of the optimization consists of the normalized densities of all members of the assembly. Using the proposed approach, one can optimize a topology for the whole assembly instead of optimizing all members of the assembly independently. This is of major importance, since changes in the topology of any member will lead to changes in the force transfer at the joints, which in consequence influences the optimization of all members of the assembly. Keywords: Topology Optimization, simultaneous optimization, Multicomponent optimization, contact 1 Introduction Topology optimization has been successfully applied to real world applications as well as many academic benchmarks. However, there are still limitations, because optimized designs are usually based on mechanical models and these models are limited by their assumptions. Especially in topol- ogy optimization, these assumptions should not significantly affect either the modeled mechanics or the optimization. This is a major challenge, as these models behave differently in optimizationcompared to mechanical analysis. Thus, assump- tions that are often made in mechanical models might lead to poor performance in optimization. One example of such an often neglected phe- nomenon is the unilateral nature of contact in joints. The necessity of modeling contact in Topology Optimization is demonstrated in several publica- tions, i.e. [10, 23, 2, 16, 11, 7, 12, 19], where different aspects of contact-constrained topology optimization are discussed. More precisely, the works of Fancello [10] and 1 Str¨ omberg and Klarbring [23] deal with linear- elastic materials and unilateral boundary con- straints. Here, a rigid contact obstacle is consid- ered, and the contact is modeled using unilat- eral contact-conditions such as unilateral penalty approaches, augmented Lagrangian approaches to satisfy Signorini’s Contact Conditions [22], as well as smooth approximations of Signorini’s Contact Conditions. Luo et al. [16] introduces hyperelastic materials in contact-constrained Topology Opti- mization. Furthermore, Andrade-Campos et al. [2] uses a node-on-segment approach in topology optimization, to model the contact between two discretized domains. Fernandez et al. [11] proposes the mortar segment- to-segment approach [17] for topology optimiza- tion, resulting in smooth contact forces and well- behaved convergence of the solution. In doing so, computing contact between multiple deformable bodies is possible, while the patch test is passed. The latter is not passed, if node-on-segment algo- rithms are used. Internal contact and/or self-contact are addressed in Bluhm et al. [7] and Frederiksen et al. [12]. Here, the third medium contact approach of Wriggers et al. [26] is used in Topology Optimization. Similar results can be achieved using the numerical stabi- lization method of Scherz et al. [19]. In doing so, void and low-density elements are stabilized for geometrical nonlinear topology optimization. Another challenging field is simultaneous mul- ticomponent topology optimization. In contrast to standard topology optimization, in which each component of an assembly is optimized inde- pendent of its surrounding, simultaneous topol- ogy optimization minimizes an objective that is defined for the whole assembly. Thus, if the design of any component of the assembly changes, the design of others might adapt, since its load- ing path is effected by the design modification. Therefore, it is not possible to decompose the mul- ticomponent optimization into several individual single component topology optimizations [8]. In the past, the following questions, among others, were examined in the context of assem- bly optimization. Where should the connections between the components be positioned? How many connection points are needed or how should the connection points be modeled? In this paper, the last questoin is addressed. The other questions are addressed for exampleby Chickermane and Gea [8], who proposes a straight forward approach to optimize the topol- ogy of the components as well as the layout of the fasteners. Here, a grid of spring elements is intro- duced to transfer the loads from one to the other component. In doing so, two design domains are introduced, one for the material distribution and one for the interconnections. The latter is similar to the RAMP approach [6]. Several modifications and extensions have been made in order to tackle challenges such as continuous positioning of joints [25, 1]. However, for the modeling of the connections between two components generally either linear spring-like connections (see Chickermane and Gea [8] or Ambrozkiewicz and Kriegesmann [1]) or even rigid constraints are imposed (see Rako- tondrainibe et al.[18]). Thus, the contact at the connecting joints is completely neglected in the assembly optimization process. This is particu- larly borderline, if bolt connections are considered. In consequence, the optimized designs are not optimized based on the actual loading paths at the connecting joints. This problem is either com- pletely neglected or workarounds such as artificial non-design spaces are introduced. This work extends the contact-constrained topology optimization framework by Str¨ omberg and Klarbring [23] to handle simultaneous opti- mization of assemblies. In doing so, the influence of contact-constraints using both spring-like con- nections and actual unilateral contact conditions are investigated. It is shown, that the unilateral characteristic must not be neglected in topology optimization and that the augmented Lagrangian appraoch is superior. The positions of the connecting joints are fixed and not part of the optimization, as the focus is on the unilateral contact behavior and its influence on the optimized assembly. Currently, only formula- tions based on linear spring models are developed to optimize the location of the joints simulta- neously. These are not compatible with unilat- eral contact formulations, see Ambrozkiewicz and Kriegesmann for more details [1]. Future work should revisit the optimization of joint locations in combination with unilateral contact-constraints. Next to multicomponent optimization, the pre- sented method is also an extension to the sin- gle component optimization. Using the proposed method, external loads can also be applied at the 2 continuous contact structure, which is not possible using the framework of Str¨ omberg and Klarbring [23] widening its applicability. The article is organized as follows. In section 2, the equllibrium equations for the assembly - connected via unilateral contact constraints - are developed. Next in the section 3, the compli- ance optimization problem is formulated, and the adjoint method is used to compute the deriva- tive of the objective efficiently. The optimization is performed to two different benchmark examples in section 4. Here also the influence of different contact laws in topology optimization is demon- strated. Last, conclusions are drawn in section 5. 2 Equilibrium equations of the contact problem As basis for the latter defined topology opti- mization problem, the equilibrium equations are introduced. In Sec. 2.1, the equilibrium for the individual linear-elastic members of the assembly are derived. Next, the connecting joint constraints are introduced in Sec. 2.2, whereas the contact constraints are derived in Sec. 2.3. In the last Section 2.4, the results of the previous sections are combined and the global equilibrium equations are obtained. In doing so, the unilateral contact- constraints will be introduced to implicitly couple adjacent members of the assembly. 2.1 Linear-elastic equilibrium The assembly consists of nComp components, and the equilibrium equations are based on the finite element analysis of each component j. Here for the topology optimization, a linear elastic constitutive law is considered and the modified SIMP approach [21] is applied. Hence, the Young’s modulus Ej iof thei-th finite element of the j-th component is a function of its normalized density φj i∈[0,1] Ej i=Ej min+ (φj i)p(Ej max−Ej min). (1) Thereby, φj iare the design variables in the topol- ogy optimization and φj i= 0 represents a void element, whereas φj i= 1 represents an element filled with material. For the j-th member of the assembly, the Young’s modulus of the void is Connecting jointBody I Body II qcjccj ccj ccjqIIqII dI,λIdII,λII Fig. 1 : Assembly connected via unilateral contact-constraints. defined as Ej minand the modulus of the material is Ej max. The penalization factor is set to p= 3. With the help of the physical Young’s modulus Ej i, one can assemble the stiffness matrix Kj=K(φj) of the j-th member of the assembly. Furthermore, the nodal displacement field dj, the nodal contact forces fj conand the external forces qjacting at thej-th member of the assembly are introduced. Thus, the equilibrium equations for each member of the assembly read K(φj)dj+fj con−qj=0. (2) Eq. (2) is amended by connecting joint constraints as well as contact conditions in order to compute the nodal contact forces fj con. 2.2 Connecting joint constraint The overall concept of connecting members of the assembly is visualized in Fig. 1, where an assem- bly of two components - namely body I and body II - is shown. Furthermore, an enlargement of the lower connecting joint is shown. The normal con- tact forces of body I are donated as λI, whereas the normal contact forces of body II are donated asλII. In order to transfer the occurring contact forces from body I to body II, a rigid continu- ous connecting joint is introduced. Note that the 3 following procedure is done for all ncjconnecting joints in the assembly. The overall idea is, that the displacement of the connecting joint ck cjis computed such that the sum of all contact forces λI/IIand possibly an additional external force acting on the joint qk cjis zero. Thus, the equilibrium equation for thek-th connecting joint is based on the con- tact forces λI/IIof adjacent components and its external forces qk cj. Thus, this formulation allows applying loads directly at the connecting joint and the equilibrium equation is Γk=nI conX AI=1 nI AIλI AI +nII conX AII=1 nII AIIλII AII +qk cj=0, (3) where the number of nodes at the contact surface of the body I and II are nI/II con. In this formula- tion, the number of nodes at the contact surface of both components can be different. Therefore, the indices AI/IIare introduced. The contact force of a specific contact node AIof body I is λI AI, whereas the AII-th contact force of body II is λII AII. Note that the following computation of the normal direction must be performed for all nodes at the connecting joint of all adjacent components. For readability, the body index is dropped and only A is used. The normal direction of A-th contact node nj A of the j-th component is shown in Fig. 2. The contact normal direction depends on the current nodal displacement dj Aof the contact nodes, as dj Anj Adj A,0 ck cj,0 ck cj Rk ηj A Fig. 2 : Contact kinematics of the A-th contact node at the k-th connecting joint.well as the joint’s displacement ck cj. Hence, the normal direction of the A-th contact node is nj A(dj A,ccj) =(ck cj,0+ck cj)−(dj A,0+dj A) ∥(ccj,0+ccj)−(dj A,0+dj A)∥,(4) where the undeformed positions of the contact nodedj A,0and the connection joint ccj,0are intro- duced. If only small deformations of the contact nodes and the center of the connecting joint are observed, it is sufficient to approximate the normal direction. In consequence, the normal direction does only depend on the positioning in the unde- formed configuration and is therefore constant. It reads nj A≈ck cj,0−dj A,0 ∥ck cj,0−dj A,0∥. (5) The normal directions are used to map the normal contact force vector λjinto the nodal domain. The nodal contact forces introduced in Eq. (2) are fj con= (Nj)Tλj, (6) where the mapping matrix Njholds the normal directions. Having the equilibrium for the mem- bers of the assembly introduced in Eq. (2) and the connecting joint constraints of Eq. (3), the unilateral contact law must be added in order to finally solve the non-linear and non-smooth contact problem. 2.3 Contact condition In order to compute the contact forces λjacting on each member of the assembly, it is necessary to distinguish between contact nodes, which are in contact and which are not in contact. Note that the contact conditions introduced in the following must be imposed for all adjacent members of the assembly. Looking at the example of Fig. 1, the contact conditions must be imposed connecting- joint wise for the contact nodes of body I and body II. This is done by introducing the contact gap ηA. Fig. 2 illustrates the deformed configuration of the A-th contact node at the connecting joint. 4 In doing so, the gap ηAcan be computed as ηj A= (nj A)T (ck cj,0+ck cj)−(dj A,0+dj A) −Rk, (7) where Rkis the radius of the connecting joint and nj Ais the normal direction of Eq. (4). Eq. (7) can be rewritten as ηj A= (nj A)T gj A,0+ck cj−dj A −Rk,(8) where the initial distance between the A-th con- tact node and the center of the k-th connecting joint is gj A,0=ck cj,0−dj A,0. (9) In order to link the contact gap ηj Awith the cor- responding contact force λj A, Signorini’s contact conditions [22] are introduced, which read ηj A≥0, λj A≥0, ηj Aλj A= 0. (10) The third condition of Eq. (10) imposes, that if there is a gap ηj A>0 between the A−thnode and the connecting joint, no contact force is trans- mitted. Thus, the contact force must be λj A= 0. Otherwise, the node is in contact with the con- necting joint, meaning that the gap is ηj A= 0 and therefore, a contact force λj A>0 is transmitted. For numerical implementations, Signorini’s con- tact conditions can be imposed by non-smooth functions such as Φj A=−λj A+ max(0 , λj A−ηj A) = 0 . (11) Next, Eq. (8) is inserted into Eq. (11) yielding to Φj A=−λj A+ max 0, λj A+rR− rh (nj A)T gj A,0+ck cj−dj Ai = 0,(12) where r >0 is a factor introduced to increase the numerical performance. If different contact laws, such as the smooth- ing approaches of Facchinei et al. [9] or penalty contact laws [27], are used, only the contact function given by Eq. (12) is replaced, wherebythe methodology remains the same. For contact- constrained topology optimization of a single com- ponent [23], it has been shown, that imposing Signorini’s contact conditions in the strong form of Eq. (12) performs slightly better than the smooth- ing approach of Facchinei et al. [9]. The latter is therefore not examined further in this work. The most commonly used contact models, are so-called penalty methods [27], where a penalty factor is introduced in order to penalize a devi- ation from the support. In almost all cases, this method neglects the one-sided characteristic of the contact. This both-sided contact laws have also been used in multicomponent topology opti- mization [8, 1]. In this case, the contact function is Φj A=−λj A−kηj A= 0, (13) where k > 0 is the penalty factor. Here, the both-sided penalty method is imposed using a Lagrange multiplier λj A, which is rather unusual. This formulation is equivalent to the standard for- mulation. It is chosen to keep the structure of the global equilibrium equations unchanged, even if different contact models are used. The global equi- libirum equations are introduced in section 2.4. If a unilateral characteristic of contact is taken into account, a unilateral penalty approach is used, its contact function reads Φj A=−λj A−min(0 , kηj A) = 0 . (14) Note that a negative gap ηj Ais equivalent to a pen- etration. The unilateral penalty contact law has also been successfully applied in stress-constrained topology optimization [13]. 2.4 Global equilibrium In the following, the derivation focuses only on Signorini’s contact conditions Eq. (12). If a penalty contact law is used, only the derivatives of the used contact law with respect to the state vector must be changed, which is straight forward for the penalty laws. More precisely, the Eqs. (22) to (24) must be adapted accordingly. Combining the contact function Eq. (12) with the equilibrium Eq. (2) and the connecting joint constraints of Eq. (3) allows computing all dis- placements dj, all contact forces λjand the 5 displacements of all kconnecting joints ck cj. There- fore, the overall set of equilibrium equations read hc(d,λ,ccj,φ) =   K(φj)dj+fcon(λj,c)−qj Φj(dj,λj,ccj) −Γk(d,λ,ck cj) =0,(15) where the curly brackets must be formulated for each member of the assembly, whereas the connecting joint constraints Γkmust be formu- lated for all connecting joints. In consequence, the unknown variables are the displacements dj and contact forces λjof each component, as well as the deformation of the center positions of all connecting joints ccj x= ( d λ)j ccj . (16) This set of nonlinear and non-smooth equations (15) can be efficiently solved using a line-search algorithm. In order to apply a line-search algo- rithm, its linearization is needed. The linearization is derived for the example of Fig. 1 to make it as legibly as possible. Thus, the number of components and number of connecting joints is nComp =ncj= 2. Furthermore, only small displacements are assumed and the normal direc- tions are assumed to be constants. The resulting linearization is ˆhc(x,φ) =Jcx−rc=   ∗ ∗ ∗ ∗I 0 0 ∂ΦI(x) ∂c1 cj  0 ∂ΦI(x) ∂c2 cj  0 ∗ ∗ ∗ ∗II 0 ∂ΦII(x) ∂c1 cj  0 ∂ΦII(x) ∂c2 cj   0−nI A 1 0−nII A 10 0  0−nI A 2 0−nII A 20 0  | {z } ˆJc ( d λ)I ( d λ)II c1 cj c2 cj  |{z} x−  q ∗I  q ∗II q1 cj q2 cj  |{z} rc=0. (17) Here, the curly brackets are  ∗ ∗ ∗ ∗j =  K(φj) (Nj)T ∂Φj(x) ∂dj∂Φj(x) ∂λj  j , (18) and  q ∗j =  qj 0|{z} no contact, r (nj A)Tgj A,0−Rj | {z } in contact  . (19) If unilateral contact laws are imposed, the current state of all contact node must be checked. Using Signorini’s contact condition, a node is in contact, if the argument of the max-function of the contact function Eq. (12) is active λj A−rηj A≥0. (20) If a unilateral penalty approach is used, the con- tact gap of Eq. (7) is used to determine the state of the contact node. Hence, a node is in contact, if the gap is ηj A≤0. (21) Last, the both-sided penalty contact law of Eq. (13) does not distinguish between nodes in contact and nodes that are not in contact. Next, the partial derivatives of the contact laws from Eq. (17) and Eq. (18) are introduced. Note that the contact law Φj Aof the A-th con- tact node and j-th component does only depend on its nodal displacement dj A, its contact force λj A, as well as the center of the considered connecting 6 jointck. The partial derivatives read ∂Φj A(dj A, λj A,ck) ∂dj A=( 0 , no contact r(nj A)T, in contact, (22) ∂Φj A(dj A, λj A,ck) ∂λj A=( −1 , no contact 0 , in contact,(23) ∂Φj A(dj A, λj A,ck) ∂ck=( 0 , no contact −r(nj A)T, in contact. (24) Finally, the non-linear equilibrium Eq. (15) and its linearization Eq. (17) are formulated. Given a density distribution ˜φof all members of the assembly and a guess for the state vector ˜x, a line-search algorithm is used to solve the non- linear equilibrium equations and therefore, com- pute the state vector x. In this work, the default line-search Newton method of the PETSc library is used [4, 3, 5]. This technique employs cubic backtracking [15]. Even though a discontinuous contact problem is considered, the line-search is applied successfully. A line-search based approach is also applied in Str¨ omberg and Klarbring [23]. The search direction of the line-search algorithm is s= (Jc(˜x,˜φ))−1h(˜x,˜φ). (25) If not very small deformations are observed, the normal direction of the A-th contact node Eq. (4) is a function of the nodal displacement dj A and the center of the considered connecting joint ck. In consequence, the derivatives with respect dj Aandckof the normal direction must be con- sidered in Eq. (17) as well. In 2D problems, the derivatives of the normal direction are symmetric 2x2 matrices ∂(nj A) ∂(dj A)= −(∆j A)2(∆j A)2 ∥∆j A∥3(∆j A)1(∆j A)2 ∥∆j A∥3 (∆j A)1(∆j A)2 ∥∆j A∥3−(∆j A)1(∆j A)1 ∥∆j A∥3  (26)and ∂(nj A) ∂(ck)=−∂(nj A) ∂(dj A). (27) The vector ∆j Ain Eq. (26) and Eq. (27) is the vec- tor between the deformed node Aand the center of the connecting joint ∆j A= (ccj,0+ccj)−(dj A,0+dj A), (28) and (∆j A)1/2describe the first or second entry of the vector. In 3D problems, the connecting joint has the form of a cylinder and its extrusion direction is nk cj. Therefore, the actual contact normal direction is obtained by projecting the cylinder surface onto a circle. In doing so, the contact nodes are free to deform tangential to the connecting joint. The resulting contact normal of the A-th contact node is  nj A(dj A,ccj) 3D=(∆j A)−(nk cj)T(∆j A)nk cj ∥(∆j A)−(nk cj)T(∆j A)nk cj∥, (29) and the actual gap is  ηj A 3D= nj A 3D (∗∗)−(nk cj)T(∗∗)nk cj −R. (30) The derivatives of the normal directions in 3D adapt accordingly. 3 Compliance optimization In this work, the compliance of the assembly is minimized. The compliance consists of the compli- ance of each linear-elastic member of the assembly as well as the compliance of each joint. The latter results from external forces at the joints qk cjand their center displacement. It reads c=nCompX j=1(qj)Tdj+ncjX k=1(qk cj)Tck cj =  q 0j qcj T  d λj ccj =rTx, (31) 7 where xis the solution of the contact problem of Eq. (15). The design variables of the optimization are the normalized densities φjof all elements of all members. They are summarized in the global design vector φ= [φ1,φ2, . . . ,φnComp]T. (32) Additionally, a volume constraint is imposed. One can either constrain the total volume of the assem- bly g(φ) =nCompX j=1 V(φj) −Vtotal, (33) or one constrains the volume of each member of the assembly individually gj(φ) =V(φj)−Vj 0, (34) where Vj 0is the maximum allowed volume for the j-th member of the assembly. The later results in nComp volume constraints. The resulting optimiza- tion problem reads minφrTx s.t.  hc(φ,x) =0, g(φ)≤0 org(φj)≤0, 0≤φj≤1.(35) In doing so, one imposes either the global vol- ume constraint of Eq. (33) or the local volume constraint of Eq. (34). Furthermore, the equilib- rium equations of the contact problem hcand the constraint on the design variables must be ful- filled. The optimization problem is solved using the method of moving asymptotes [24]. Therefore, the gradients with respect the design variables of the objective and the constraints are needed. 3.1 Gradient of the compliance The gradient of the compliance of the assembly can be derived similarly to the single component contact-constraint topology optimization frame- work [23]. The compliance from Eq. (31) consists of the constant external forces vector rand the state vector x. Thus, only the state vector xhold- ing the displacements dj, the contact forces λj and the displacements of the connecting jointsccjdepends implicitly on the design variables. In order to compute the derivative of the state vec- tor as efficient as possible, the adjoint vector ζis introduced. Then, the product of the adjoint vec- tor and the linearized equilibrium equation (17) is subtracted from the compliance. The linearized equilibrium equations are zero for the solution of the contact problem, so that the adjoint vector can be chosen freely, if the solution of the contact problem is considered. The extended compliance is ˜c(x) =rTx−ζThc(x,µz). (36) In consequence, the gradient of the compliance with respect to the e-th entry of the global density vector φis ∂˜c(x) ∂φe=∂r ∂φeT x+rT∂x ∂φe−ζT∂hc(x) ∂φe =rT∂x ∂φe−ζT∂Jc(x) ∂φex+Jc(x)∂x ∂φe −∂rc ∂φe =h rT−ζTJc(x)i∂x ∂φe −˜ζT∂Jc(x) ∂φex−∂rc ∂φe . (37) If the adjoint ζis chosen such that rT−ζTJc(x) =0, (38) the costly derivative∂x ∂φeis eliminated. Hence, the gradient simplifies to ∂c(x) ∂φe=∂˜c(x) ∂φe=−ζT∂Jc(x) ∂φex−∂rc ∂φe . (39) Furthermore, the vector rcis independent of φ, so that only the derivative of the Jacobian Jcwith respect to φmust be computed. Coming from the single component framework of Str¨ omberg and Klarbring [23], the derivative of the Jacobian is straight forward, since only the stiffness matrices Kj(φj) of each member of the assembly depend explicitly on the design variables. Note that each 8 entry of the design vector φdescribes the den- sity of one finite element of one particular member of the assembly. In consequence, a member wise approach is very efficient, since all entries of the Jacobian, that are not related to the considered member of the assembly must be zero. In doing so, the global design vector φis subdivided into the density vectors of each member φjand the gradient computation is done member by mem- ber. Following this approach, the gradient of the Jacobian is ∂Jc(x) ∂φj e=  ∂K(φ) ∂φe0 0 0  j . (40) Here the derivative of the stiffness matrix with respect to the e-th element is the derivative of the e-th element’s local stiffness matrix. In the deriva- tive given by Eq. (40), the SIMP approach Eq. (1) must be considered. In summary, the derivative of the Jacobian of the single component implemen- tation [23] must be done for each member of the assembly, in order to compute the gradient of the objective with respect to the design variables ∂c(x) ∂φe=∂c(x) ∂φj ˆe= (ζj)T  ∂K(φj) ∂φj ˆe0 0 0  j( d λ)j , where ζjis the part of the adjoint vector cor- responding to the nodal displacements and the contact forces of the j-th member. This member- wise concept is valid, since both the contact functions and the connecting joint constraints do not depend explicitly on the design variables. 3.2 Gradient of the volume constraints Next to the objective, the gradients of the volume constraints are needed. The constraints are intro- duced in Eq. (33) and Eq. (34). If the total volume of the assembly is constrained, the derivative with respect to the density of the e-th element reads ∂g(φ) ∂φe=Ve, (41) where the volume of the e-th undeformed element isVe. Otherwise, the volume of each member ofthe assembly is constrained. Thus, instead of only one volume constraint, nComp volume constraints are introduced. The gradient of each constraint is ∂gj(φ) ∂φe=( Vj e, φeis part of j-th component 0 ,else (42) For both volume constraints, the gradients with respect to the design variables are independent of the design variables and therefore constant. 4 Numerical Examples The developed simultaneous multicomponent optimization framework is applied to 2D and 3D examples. First, the single component bench- mark of Str¨ omberg and Klarbring [23] is modified, in order to demonstrate the differences between forces that are applied to the structure or to the contact support structure. In this context, the influence of penalty contact conditions and Sig- norini’s contact conditions is also investigated. Next, the L-beam example is analyzed for dif- ferent two component designs, where the L-beam is modeled as assembly. At last, the 3D L-beam is considered to proof the scalability of the pro- posed method. All shown numerical examples are computed on a standard workstation. 4.1 2D benchmark of Str¨ omberg and Klarbring The benchmark example of Str¨ omberg and Klar- bring [23] and the resulting optimized design is shown in Fig. 3. In this benchmark, the connecting joint is considered as rigid contact support struc- ture, which is fixed in space. The optimization is based on four load cases, where the applied load acts at the lower right corner of the design domain. The four load cases are q1= (1000 N ,0 N), q2= (−1000 N ,0 N), q3= (0 N ,1000 N) and q4= (0 N ,−1000 N). In the implementation [20], the design domain is discretized by 4444 finite elements, resulting in 50 contact nodes at the rigid contact support. The considered linear- elastic material is characterized by a Poisson’s ratio of 0 .3 and a Young’s modulus of 210 GPa. Further, the truncation radius of the used sensi- tivity filter [21] is 10 mm and the volume fraction 9 q x y (a) Design space. (b) Optimized design using Signorini’s contact condition [20]. Fig. 3 : Benchmark example of Str¨ omberg and Klarbring[23] for unilateral contact constrained topology optimization. is 0.5. A sensitivity filter is used, since the pub- lished results of Str¨ omberg and Klarbring [23] are also based on this filter. However, other filtering techniques can be applied without any problems. The original benchmark does not consider adapt- ing normal directions. Thus, the normal directions are considered to be constant and therefore only very small deformations are assumed. The nor- mal directions are computed based on Eq. (5). All presented optimizations are terminated after 100 iterations. The resulting compliance optimization prob- lem of the original benchmark reads minφn=4P i=11 4 qT idi s.t.  hc(φ,x)i=0, V(φ)−V0≤0, 0≤φ≤1.(43) 4.1.1 Modified benchmark using Signorini’s contact condition In contrast to the original benchmark of Str¨ omberg and Klarbring[23], the proposed assem- bly optimization approach also allows applying forces at the connecting joint (contact support). Thus, a force transfer at the connecting joint is modeled, which therefore widens the applicabil- ity of contact-constrained topology optimization also for the single component optimization. Notethat the very efficient framework for contact- constrained topology optimization of Str¨ omberg and Klarbring[23] can be seen as sub-problem of the proposed method, where the connecting joint cannot move. In doing so, the connecting joint constraint of Eq. (3) is replaced by Γk=ck cj=0. (44) The modified benchmark for the proposed method is shown in Fig. 4a. The same load cases are considered as in the original benchmark. How- ever, they are now applied at the connecting joint instead of the lower right corner of the design domain. This is possible, since the center of the connecting joint is not fixed in space anymore. As in the original benchmark, Signorini’s contact con- dition is used. Note that the modified as well as the initial Benchmark are not symmetric, so that non- symmetric optimized designs are expected and observed. The resulting optimization problem is now minφn=4P i=11 4(rixi) s.t.  hc(φ,x)i=0, V(φ)−V0≤0, 0≤φ≤1.(45) The obtained design using Signorini’s contact con- dition is shown in Fig. 4b. Next to this, the convergence plots of the original benchmark as qcjccj x y (a) Design space. (b) Optimized design using Signorini’s contact condition. Fig. 4 : Modified benchmark example using the multicomponent method. 10 0 20 40 60 80 10005001,000 Iterations [-]Objective [Nmm]OB Signorini MB Signorini MB both-sided penalty MB unilateral penalty Fig. 5 : Convergence of the performed optimiza- tions of Fig. 3, Fig. 4 and 6. (OB - Original Benchmark, MB - Modified Benchmark) well as the modified setup are shown in Fig. 5 and both methods converge as expected. Taking a closer look at Fig. 3b and Fig. 4b reveals, that the optimized design adapts to the modifications of the setup. Hence, it is crucial to differentiate between the existing contact- constrained topology optimization framework of Str¨ omberg and Klarbring[23] and the proposed method of this work also for the single component optimization. 4.1.2 Modified benchmark using penalty approaches In topology optimization as well as finite ele- ment analysis, contact is often modeled using penalty approaches. This is still nowadays the most common modeling approach, if contact is considered, and therefore, also successfully applied in the assembly optimization approaches of (see e.g. Chickermane and Gea[8], Thomas et al[25] or Ambrozkiewicz and Kriegesmann[1]). Especially in assembly optimization, it has the advantage that the locations of the connecting joints can be optimized easily, since the both-sided penalty contact law does not add any non-linearity to the problem and the unilateral characteristic of contact is neglected. However, the influence of penalty approaches on the final design is not often discussed. Therefore, the influence is inves- tigated using the modified benchmark of Fig. 4. The penalty factor is set to k= 2000 N /mm for the both-sided and unilateral contact law. (a) Assembly opt. using both-sided penalty con- tact law. (b) Assembly opt. using unilateral penalty contact law. Fig. 6 : Comparison of the optimized designs using penalty contact laws, where forces are applied at the connecting joints. The resulting optimized designs are visualized in Fig. 6 and the convergence history of both opti- mization is added to Fig. 5. The optimized design using the standard both-sided penalty approach is shown in Fig. 6a. If this design is compared to Fig. 4b (using Signorini), major differences are observed. In particular, significantly less material is provided in the vicinity of the contact struc- ture and some struts are split. Comparing the final design using the unilateral penalty approach (see Fig. 6b) to the both-sided penalty design (see Fig. 6a), significantly more material is distributed in the vicinity of the contact structure, while the struts are very similar. Nonetheless, the strut appearance is still very different in comparison to the design computed with Signorini’s contact conditions (see Fig. 4b). Furthermore, the compliance is reduced the most for the design using the Signorini’s contact condition. This is shown in Table 1, where the three final designs of the modified benchmark are evaluated twice. Once, it is evaluated using the Table 1 : Evaluation of the optimized designs from Fig. 4 and Fig. 6 based on the used contact law and Signorini’s contact conditions. final design Compl. in opt. Compl. using Signorini Signorini Fig. 4b 56.13 56.13 both-sided pen. Fig. 6a 73.79 58.22 unilateral pen. Fig. 6b 98.62 56.96 11 contact law used in optimization and once using Signorini’s contact condition for comparision. In comparison to the design of Fig. 4b, the perfor- mance of the design computed with the both-sided contact law is 3.7% worse. Thus, using both-sided penalty contact laws in topology optimizations results in different designs with noticeable losses in performance. Further, the performance of the design based on the unilateral contact law is a lot better. Nonetheless, the design is still 1.5% off compared to the design using Signorini’s con- tact condition. However, during optimization, the compliance is highest with the unilateral penalty approach. This is expected because, unlike the two-sided approach, only pressure regions trans- mit contact forces, leading to more penetration at the contact support. In consequence, a higher compliance is observed, if the penalty factor is the same for both penalty methods. These results imply that it is crucial to con- sider the unilateral characteristic of contact in topology optimization, in order to obtain good performing final designs. Note that one could also increase the penalty factor kin order to obtain better results for the penalty approaches, but the Jacobian (see Eq. (17)) will become more and more ill posted and the computation time will increase significantly. In this example, the whole optimization using the both-sided penalty method took 501 s, whereas the unilateral penalty opti- mization took 505 s. The optimization approach using Signorini’s contact condition takes approx- imately 4 times longer. Consequentially, the uni- lateral penalty approach comes at zero extra costs, in comparison to the both-sided penalty approach, whereas Signorini’s contact condition is the most expensive, while having the best performing design. 4.1.3 Adapting normal direction Until now, very small deformations were assumed, so that the normal direction of each contact node is not updated and therefore constant. Here, the modified benchmark is optimized using Signorini’s contact condition and changes of the contact nor- mal direction are considered. Thus, Eq. (4) is used to update the normal direction while solving the global equilibrium equations. The final design is shown in Fig. 7 and only minor differences to the Fig. 7 : Modified benchmark with adapting con- tact normal directions. design of Fig. 4b are observed. Thus, assuming small deformations is reasonable for this example. 4.2 2D L-beam example The considered L-beam example is visualized in Fig. 8. There, the design space as well as the single component optimized design without using con- tact are shown. The dimensions of the design space and boundary conditions are based on Holmberg et al. [14]. The material parameters are unchanged to Sec. 4.1, while the truncation radius of the sensitivity filter is set to 3.5. Dirichlet boundary conditions are imposed at the top and a vertical load of q= (0,−1500 N)Tis applied at the top of the tip. The volume fraction is set to 0.5. 4.2.1 Multicomponent optimization using Signorini’s contact condition In the following, the L-beam is optimized using the proposed multicomponent optimization approach. qBody I Body IIx 8080120120 y (a) Multicomponent design space. (b) Single component opt. design. Fig. 8 : Design space and single component refer- ence for the L-Beam. 12 (a) 2 connecting joints assembly. (b) 3 connecting joints assembly. (c) 4 connecting joints assembly. Fig. 9 : Optimized designs of the L-beam using Signorini’s contact condition. Therefore, the L-beam is decomposed into two beam-like components forming the assembly. In doing so, the red area of Fig. 8a is the over- lap of both components of the assembly. This is where connecting joints can be placed. The radius of each connecting joint is set to 5 mm. Here, three different setups with a different number of connecting joints are investigated and a global vol- ume constraint with a volume fraction of 0.5 is imposed. The obtained optimized designs are shown in Fig. 9 and the convergence plots of all opti- mized designs are shown in Fig. 10. It is worth mentioning, that the final compliance of the sin- gle component design (Fig. 8b) is the lowest. However, the design with 3 joints (Fig. 9b) and the 4 joint design (Fig. 9c) are close by. More precisely, the compliance of the final design of0 20 40 60 80 10000.511.5·104 Iterations [-]Objective [Nmm]SC Optimization MC 2 joints MC 3 joints MC 4 joints Fig. 10 : Convergence of the optimizations of the L-beam of Fig. 9. (SC - single component, MC - Multicomponent) the single component compliance optimization is 1935.0 Nmm, whereas the compliance of the 3 joint design is 2013.6 Nmm and the 4 joint design has a final compliance of 2040.3 Nmm. Only the design with 2 connecting joints performs significantly worse, with a final compliance of 2403.2 Nmm. This shows, that the number of joints as well as their positioning is challenging [1]. However, very good performance is observed, if the joints are placed inside struts of the standard single compo- nent compliance optimization results of Fig. 8b. In this example, this is done for the 3 joint design, which is the best performing assembly. Another noteworthy observation is the unused fourth contact support in Fig. 9c. All calculated contact forces at this support are about the same size as the relative tolerance of the convergence criteria used by the line-search algorithm solv- ing the contact problem. The relative tolerance is set to rtol= 10−8. Thus, no forces are transmit- ted at this joint. In these examples, the unilateral characteristic of contact is considered, so that the optimization algorithm has no interests in plac- ing material outside the pressure zones of the joints. This can be seen at any joint of Fig. 9a to Fig. 9c. As for the unused joint of Fig. 9c, the contact forces of contact nodes that belong to void elements are always in the order of magni- tude of the relative tolerance. Hence, forces are only transmitted at contact nodes of filled finite elements. 13 (a) 3 connecting joints assembly both-sided penalty. (b) 3 connecting joints assembly unilateral penalty. Fig. 11 : Optimized designs of the L-beam using penalty contact laws. 4.2.2 Multicomponent optimization using penalty approaches The 3 joint design is also optimized using the both- sided penalty and unilateral penalty approach. In these examples, the penalty factor is set again to k= 2000 N /mm. The obtained optimized designs are shown in Fig. 11. As in the former section 4.1, the both-sided penalty design results in a differ- ent design with a poorer performance, whereas considering the unilateral penalty law results in an overall similar design to Fig. 9b. Comparing the design using the unilateral penalty law (see Fig. 11b) and Signorini’s contact condition (see Fig. 9b), differences around the connecting joints are observed at the vertical component (body I), whereas the struts differ in size significantly at the horizontal member (body II) of the assembly. A closer look at the upper right joint of the hor- izontal component of Fig. 11a demonstrates the impact of neglecting the unilateral characteristic perfectly. In contrast to the unilateral designs, the diagonal struts do not go around this connecting joint, so that only tensile forces are transmitted. Thus, these struts are dysfunctional in practice. Furthermore, it is worth mentioning that once again the computation time of both penalty meth- ods is very similar, while using Signorini’s contactcondition is only 3 times slower. However, this changes drastically, if the four joint design is com- puted. If the both-sided penalty model is used, the line search algorithm suddenly needs a lot more iterations. In consequence, the numerical costs increase tremendously and exceed by magni- tudes the costs of the unilateral penalty law and Signorini’s contact condition. A closer look reveals that the progress of the line search algorithm shrinks drastically compared to the other two contact laws. In consequence, no optimization can be performed in reasonable time using the both-sided contact law under these conditions. (a) Signorini’s contact condition. (b) Both-sided penalty. (c) Unilateral penalty. Fig. 12 : L-beam with adapting contact normal directions. 14 4.2.3 Adapting normal direction The former L-beam optimizations are performed assuming small displacements, so that the con- tact normal directions can be approximated using Eq. (5). This assumption is removed and there- fore, Eq. (4) is used to compute the normal directions. Thus in this section, each contact nor- mal direction depends on the nodal displacement of the contact node as well as the displacement of the center of the connection joint. The opti- mized designs are shown in Fig. 12 and major differences are observed, in comparison to the solu- tion assuming small displacements of Fig. 9b and Fig. 11. Taking a closer look at the results imposing Signorini’s contact conditions, the vertical compo- nent (body I) of Fig. 9b and Fig. 12a are almost identical, while the horizontal component (body II) differs significantly inside the red box. The struts differ in size and more importantly, the strut at the lower connecting joint have a complete dif- ferent angle of attack. This happens, because the loading path inside the components changes, if the adapting normal directions are considered. In consequence, the optimized design adapts. One has to admit, that the design at the lower connecting joint of the horizontal com- ponent might fail in practice, since the struts might slip from the connecting joint. However, this can be easily solved by introducing a non- design space around the connecting joint, which is common practice in multicomponent optimization [1]. Using Signorini’s contact condition, the width of the non-design space is actually known, since the actual characteristic of contact mechanics are considered and therefore, the actual loading path inside the component is taken into account dur- ing optimization. In consequence, the width of the non-design space depends exclusively on the man- ufacturing constraints of the used fasteners. Using penalty contact laws, this assumptions hold only if extremely high penalty factors are used, since the absolute value of the penalty factor and the unilateral characteristic influence the loading path and therefore, the optimized design. The optimized design based on penalty contact laws ( k= 2000 N /mm) are shown in Fig. 12b and Fig. 12c. Also here, the optimized design consid- ering the adapting contact normal direction differ significantly from the previous designs assuming (a) Both-sided penalty. (b) Unilateral penalty. Fig. 13 : L-beam with adapting contact nor- mal directions and increased penalty factor k= 10000 N /mm. very small displacements of Fig. 11. If the both- sided penalty approach is used, the final design once again includes dysfunctional struts at the upper two connecting joints (see Fig. 12b). Here, the struts of the horizontal component transmit tension forces, which is nonphysical in the context of contact. Once again, one can introduce non- design spaces to overcome this problem. However using both-sided penalty approach, the width of the non-design spaces is not only constraint by the manufacturing constraints of the fasteners. The underlying model transmits tension forces at the joints. Consequentially, the underlying model will underestimate the needed material at the joints, which might lead to failure in practice. Thus, the width of the non-design space must be chosen very conservatively, which in turn has a signifi- cant impact on the performance of the optimized design. In Fig. 13, the penalty factor is increased to k= 10000 N /mm, in order to show the impact of the penalty factor of the both-sided and uni- lateral contact law. In doing so, the optimized design should become more and more similar to the design obtained by Signorini’s contact con- dition (see Fig. 12a). Especially, the unilateral 15 design transforms towards the design using Sig- norini’s contact condition. However, especially the appearance of the struts between the lower and left connection joint is still very different. At the same time, the computation times increase sig- nificantly for both penalty. Now, the both-sided penalty law is four times slower, than the opti- mizations using Signorini’s contact condition and the unilateral penalty law is two times slower, than using Signorini’s contact condition. Furthermore, its worth mentioning, that considering adapting contact normal directions reduces the computa- tion time by 20%, if Signorini’s contact conditions are imposed, while the computation times of both penalty methods increase in comparison to the constant contact normal directions. 4.3 3D L-beam example The presented simultaneous optimization approach is also applied to a 3D example of the L-beam. Therefore, the 3 joint design space is extruded into z-direction and the two components are placed next to eachother. The vertical beam is hold in place at it’s four upper corners and a load is distributed along the tip of the beam. Furthermore, symmetry is imposed so that the connecting joints cannot move along the z-axis. In doing so, a total of 242534 finite elements are considered, and each component has 1056 con- tact nodes. The results are shown in Fig. 14 and Fig. 15. The contact is modeled using Signorini’s contact conditions. Thus, no penalty approach is used here. The overall appearance of the vertical com- ponent of the 3D design is similar to the 2D design of Fig. 9b. However, there is also mate- rial distributed around the lower connecting joint and some differences are observed at the top. In contrast to the vertical component, the horizon- tal component differs significantly from the 2D design. Less struts are used and a ’I-profile’ like structure is found. Hereby, the width of the inner support differs over the height of the component. Interestingly, the appearance close to the connect- ing joint is similar to the 2D design. Due to the applied symmetry constraint, the actual appear- ance of the optimized assembly is shown in Fig. 15. (a) Frontside. (b) Backside. Fig. 14 : Optimized design of the L-Beam using Signorini’s contact condition. 5 Conclusion In this paper, the general framework for contact- constrained topology optimization of linear-elastic structures of Str¨ omberg and Klarbring [23] is extended to simultaneous optimization of assem- blies. Furthermore, external loads can be applied at the connecting joints extending the single com- ponent framework of Str¨ omberg and Klarbring [23] significantly. In doing so, different benchmark examples are shown and the necessity of mod- eling the unilateral characteristic of contact in optimization is demonstrated. 16 Fig. 15 : Final design including the symmetry. It is shown, that neglecting the unilateral characteristic in optimization, e.g. using the both- sided penalty approach, results in poorer perform- ing optimized design. One might argue, that one will always introduce a non-design space around the entire connecting joint. Nevertheless, unilat- eral contact laws must be considered, as otherwise both the compression and tension zones of the con- tact surface contribute to the stiffness of the joint. The optimizer will take advantage of this non- physical assumption, which in practice can lead to dysfunctional joints. When considering the unilateral characteristic, the actual contact mechanics are taken into account during the optimization. In consequence, topology optimization will place a sufficient amount of material in the pressure zones at the contact surface. Therefore, a non-design space around a connecting joint is only limited by man- ufacturing constraints, like the necessary support for fasteners (e.g., bolt connections). The size of these non-design spaces is known in advance. Furthermore, the augmented Lagrangian approach directly satisfying Signorini’s contact- conditions is compared to the unilateral penalty approach. Here, the penalty approach was suffi- cient, for the single component example, while in the multicomponent simulations, the computation time of the penalty approaches exceeded the com- putation time using the augmented Lagrangian approach significantly. Therefore, the augmented Lagrangian approach is recommended, since it is easier to handle in practice.As stated in the introduction, future work should also focus on optimizing the placement of the connecting joints. Here new coupling tech- niques must be developed, as the state of the art methods are not applicable using unilateral con- tact constraints. Furthermore, segment-to-segment contact formu- lations should be investigated, especially if stress- constrained topology optimization is of interest. In contrast to the node-on-segment approaches, the approach of Fernandez et. al. [11] results in a smoother contact pressure and passes the patch test. Declarations Funding. The current work is a part of the research training group “Simulation-Based Design Optimization of Dynamic Systems Under Uncer- tainties” (SENSUS) funded by the state of Ham- burg within the Landesforschungsf¨ orderung under project number LFF-GK11. Conflict of interest. On behalf of all authors, the corresponding author states that there is no conflict of interest. Replication of results. The used in-house Toolbox TOptiMuM cannot be published at this time. Nonetheless, all needed parameters to recompute the results are mentioned in the work.