key: cord-0073935-z7tis882
authors: Han, Yongsheng; Xu, Bin; Duan, Zunyi; Huang, Xiaodong
title: Stress-based topology optimization of continuum structures for the elastic contact problems with friction
date: 2022-01-21
journal: Struct Multidiscipl Optim
DOI: 10.1007/s00158-022-03169-1
sha: 00e9118e60db16908ecc93537e16cd391510b24a
doc_id: 73935
cord_uid: z7tis882

Structural problems have various nonlinearities in the real world and these nonlinearities should be accommodated in structural topology optimization. This work proposes a topology optimization method for minimizing the maximum von Mises stress of elastic continuum structures with frictional contact under material usage constraint, using an extended Bi-directional Evolutionary Structural Optimization (BESO) method. Stresses are treated as global performance (objective) function, the global von Mises stress is measured by the p-norm stress aggregation approach, and the friction behavior is governed by the Coulomb friction law regularized in analogy with the perfect elasto-plastic theory. BESO method based on discrete variables which can avoid the well-known stress singularity and the numerical instability issue in frictional contact problems. The adjoint sensitivity analysis method is adopted to derive the sensitivity numbers. The effectiveness of the proposed method is validated through a series of comparison studies including elastic-rigid and elastic-elastic contact problems. The influence of varying friction coefficient on the optimized results and the stress distributions are investigated in comparison with the maximum stiffness design. The effect of different parameters including p-norm, volume fraction and mesh density on the optimized results are discussed. The optimized results, for elastic-rigid contact, indicate that the maximum stress can be reduced compared with elastic-elastic contact. The optimized stress decreases as the friction coefficient increases because the friction behavior resists the tangential deformation at the contact interface. The results also show that the proposed approach can achieve a reasonable design that effectively controls the stress level and reduces the stress concentration effect at the critical stress areas.

Topology optimization is a versatile and powerful tool to find the optimized material distribution in a specified design domain under certain constraints. From the seminal work by Bendsøe and Kikuchi (1988) , topology optimization has been used in structural conceptual design (Zhu et al. 2016) . The original contributions of topology optimization were focused on the stiffness maximization design under material usage constraint (Wei and Wang 2009; Li et al. 2017; Xia et al. 2018a, b, c) . Notably, the well-established contributions on maximizing the stiffness of the structures cannot ensure higher strength if the stress criterion is not taken into account (Duysinx and Bendsøe 1998) . However, considering the numerous underlying numerical issues, the stress-based topology optimization is a challenging problem.

Stress-based topology optimization has four prominent challenges (Le et al. 2010; Xia et al. 2018a) . First, the highly nonlinear stress behavior associated with local stress. Stress is extremely sensitive to the material redistribution, which results in the stress-based topology optimization is difficult to converge. The stress concentrations appear at the reentrant corners and zig-zag edges make convergence even worse. To stabilize the convergence, a global stress measure is adopted as a constraint function in the stress-based topology optimization formulation (Guo et al. 2011) . This is a possible approach followed by several works but local stress constraints are also used with success. Second, the well-known stress singularity which mainly arises in density-based methods with low density elements. The typical relaxation techniques for alleviating this difficulty include ε-relaxed formulation, qp approach (Bruggi and Duysinx 2012) . Stress singularity problem can be naturally avoided using discrete topology optimization methods, such as Evolutionary Structural Optimization (ESO) method (Xie and Steven 1993) and its improved Bi-directional ESO (BESO) method (Huang and Xie 2007) which have been widely adopted for both academic researches (Fritzen et al. 2016 ) and engineering applications (Xia et al. 2018b) . Third, stress is a local measure, which leads to enormous constraint functions and the optimization process is computationally intensive. To alleviate the dilemma, the active set strategy (Duysinx and Bendsøe 1998) has been proposed. The Augmented Lagrangian Method as another alternative to alleviate this issue. It was proposed by Pereira et al. (1995) and has been increasingly used with success since then (Silva et al. 2018 (Silva et al. , 2021 Emmendoerfer and Fancello 2014; Emmendoerfer et al. 2019) . Additionally, to reduce the computationally cost, an alternative way to handle the large number of local constraints is the adoption of an aggregation function, such as the p-norm and the Kreisselmeier-Steinhauser (K-S) aggregation functions (Duysinx and Bendsøe 1998; Xia et al. 2018a) . The fourth challenge is the accuracy of stress assessments. The inaccuracy in the evaluation of the stress field is related to the jagged nature of the optimized structure coming from the finite element mesh.

For decades, most of the contributions in the field of structural topology optimization are constructed based on a linear assumption. Mechanical contact exists extensively in structural designs concerning connection, fixation, load transmission, and so forth. The magnitude and distribution of contact stress have significant influences on structural performance, such as the wear process at the contact interface and the fatigue life of the involved mechanical structure. More seriously, crack initiation and subsequent premature fatigue failure in the presence of high contact stress concentration. Therefore, the structural optimization considering contact stress is of practical significance.

Design for contact problems was first considered by Benedict (1982) . Afterward, contact topology optimization was extended to continuum structures (Fancello 2006; Desmorat 2007; Luo et al. 2016 ). More recently, Ma et al. (2020) proposed an equivalent static displacements (ESD) method to solve the contact force optimization. Bluhm et al. (2021) proposed an extension of the third medium contact method for solving structural topology optimization problems that involve and exploit self-contact. A new regularization of the void region, which acts as the contact medium, makes the method suitable for cases with very large deformations. Lee and Park (2012) minimized the strain energy using the equivalent static loads (ESL) method for contact nonlinear static response structural optimization. Li et al. (2021) extended B-spline parameterization method to topology optimization of elastic contact problems for compliance minimization. Niu et al. (2020) presented a topology optimization method for the stiffness maximization design of elastic structures with frictional contact. Results shown that the optimized compliance generally decreases as the friction coefficient increases. A challenge in topology optimization with frictional contact is the numerical instability issue that is caused by the elements with weak materials along the contact boundary. It may result in nonconvergence of the frictional contact analysis procedure and even abnormal interruption of the optimization process. Niu et al. (2019) developed a topology optimization method for elastic continuum structures in frictionless contact to improve the uniformity of contact pressures including elastic-rigid and elastic-elastic frictionless contact problems.

However, existing works on contact topology optimization have been rather limited both in number and in breadth so far. Current literature on contact topology optimization of continua focus on stiffness maximization problems. As for the stress-based topology optimization of continuum structures considering contact nonlinearity, few work can be referenced. The work by Fancello (2006) is the first one in the literature that deals exactly with stress-constrained topology optimization and contact interfaces with successful results.

Novelly, the aim of this work is to accomplish stress minimization topology optimization for elastic continuum structures subjected to frictional contact and to study the influence of friction behavior on the stress-based topology optimization process and the optimized results. Using an extended BESO method, the well-known stress singularity and the numerical instability issue in frictional contact problems can be avoided. Stresses are treated as global performance (objective) function, the global von Mises stress is measured by the p-norm stress aggregation approach, and the friction behavior is governed by the Coulomb friction law. The sensitivity number is derived by adjoint method. Both sensitivity numbers and topological variables are filtered to overcome the highly nonlinear stress behavior and stabilize the optimization procedure. The effectiveness of the proposed method is validated through two comparison studies including elastic-rigid and elastic-elastic contact problems in comparison with the maximum stiffness design. The influence of varying friction coefficient on the optimized results and the stress distributions are investigated. The effect of different parameters including p-norm, volume fraction and mesh density on the optimized results are discussed. The methodology is generally applicable to any predetermined contact surface, as long as the two structures are connected at all nodes of the contact surface. In order to simplify, this paper only studies two numerical examples of horizontal contact.

The layout of this paper is organized as follows. In Sect. 2, the finite element formulation for node-to-node frictional contact analysis is described. Section 3 formulates the stress-based topology optimization model of continuum structures for elastic contact with friction. Section 4 derives the expression of sensitivity with respect to topological variables (elemental densities). The numerical techniques of the extended BESO are presented in Sect. 5. Section 6 validates the design method through a series of comparison design problems. Conclusions and future perspectives are drawn in Sect. 7.

The basic finite element formulation of node-to-node frictional contact analysis is briefly described in this section. The reader can refer to the monographs (Kikuchi and Oden 1988; Wriggers 2006; Yastrebov 2013) and articles (Niu et al. 2019 (Niu et al. , 2020 for detailed formulation. Without loss of generality, a contact problem between two planar elastic bodies of Ω I and Ω II is shown in Fig. 1 . The two elastic bodies after deformation are represented by Ω I and Ω II . f N denotes the contact and reaction forces on the contact surface. According to the theory of action and reaction force (Newton's third law), the two forces are equal in magnitude and opposite in direction. f T denotes the frictional force on the contact surface. Similarly, the two forces are equal in magnitude and opposite in direction. Γ I c and Γ II c denote the potential contact surfaces. The real contact surface Γ c , is usually unknown until the contact problem is solved. g N0 denotes the initial distance of the potential contact node pairs. ̃ I i and ̃ II i , respectively, denote the coordinate vectors of nodal points I i and II i in the deformed configuration. The node-to-node contact discretization along the potential contact surfaces, as shown in Fig. 2 . I i and II i denote the coordinate vectors representing a contact node pair indexed by i ( i = 1, 2, … , n p ) along the potential contact surfaces in the undeformed configuration. n p denotes the number of potential contact node pairs. i and i denote the normal and tangential unit vectors of contact node pair i , respectively. g N0i denotes the distance of the potential contact node pair i in undeformed configuration. Oxy represents the Cartesian coordinate system.

The normal contact gap g Ni of contact node pair i is calculated by projecting the relative displacement vector

where ̂ Ni is the normal kinematic transformation vector that contains the normal unit vectors related to contact node pair i (Zhang and Niu 2018) . = ( I ) T ( II ) T T and = ( I ) T ( II ) T T collect all the nodal coordinates and displacements of both discretized bodies I and II, respectively.

The tangential slip of contact node pair i can be obtained as Fig. 1 Contact problem of two planar elastic bodies (Niu et al. 2020 ). a Undeformed configuration. b Deformed configuration Fig. 2 Node-to-node contact discretization at undeformed configuration where ̂ Ti is the tangential kinematic transformation vector containing the tangential unit vectors of contact node pair i (Zhang and Niu 2018 ).

T h e n o r m a l c o n t a c t b e h av i o r o b e y s t h e Hertz-Signorini-Moreau conditions (Wriggers 2006) .

where f Ni is the normal contact force of contact node pair i.

The geometric nonpenetration condition g Ni ≥ 0 is regularized by allowing some small penetration (Kikuchi and Oden 1988; Wriggers 2006) . The normal contact force can be calculated as where N > 0 is the normal penalty parameter.

The friction force (Niu et al. 2020 ) of contact node pair i is computed as where T > 0 is the tangential penalty parameter and the constant parameter > 0 is the friction coefficient.

The finite element equation for the frictional contact problem shown in Fig. 1 can be obtained as where and denote the global elastic stiffness matrix and the displacement vector, respectively. , N and T are the global applied external load, the normal contact load and the friction load vectors, respectively. T can be further divided into the sticking part Tk and the sliding part Td . The finite element equation for the frictional contact nonlinear problem is generally solved using the incremental Newton-Raphson method (Wriggers 2006; Zhang and Niu 2018) . At each load increment n , the residual vectors can be obtained where n is the accumulated external load vector at the n th load increment.

Once the displacement vector n at load increment n obtained, N , Tk and Td can be expressed, respectively, as where N,n , Tk,n and Td,n represent the converged normal, sticking and sliding contact stiffness matrices detailed in (Zhang and Niu 2018) .

Thus, Eq. (7) can be rewritten as where is the converged tangential stiffness matrix of the n th load increment.

Particularly, for the frictionless contact problem, Eqs. (6-10) can be, respectively, simplified as

The design domain is discretized into N elements. The topological variables is denoted by = ( 1 , 2 , … , k , … , N ) T and assigned values of

where the values of one and min (= 10 -6 ) correspond to solid and void material, respectively.

The elemental constitutive matrix k of the k th element can be linked to the associated topological design variable k as

where 0 denote the constitutive matrix of solid material.

Following (Duysinx and Bendsøe 1998) , the overall stress criterion which aims to control the stress state at the microscopic level for rank-2 layered composites is considered. For a power-law type material interpolation model, according to (Xia et al. 2018a ), the consistent model of the effective stress vector k is given by where denotes the strain-displacement matrix. k and k are, respectively, the displacement and strain vectors of the k th element.

The p-norm global stress aggregation is adopted to approximate the maximum von Mises stress where vm,k denotes the von Mises stress of the k th element. p denotes the norm parameter.

The elemental von Mises stress is calculated as where and k denote the constant stress coefficient matrix and the stress vector of the k th element, respectively. The stress coefficient matrix can be written as

The von Mises stress can be obtained by substituting Eq. (21) into Eq. (20) as where xx,k , yy,k and xy,k denote the components of the stress vector k .

The topology optimization model for the minimization of the p-norm global stress under elastic contact with friction subject to volume constraint can be constructed as

where n and n denote the residual vectors. V req and V( ) are the required and total material volume, respectively. v k is the volume of the k th element.

The topology optimization problem as illustrated in Eqs. (23-25) is a stress minimization problem and not stressconstrained problem. While seeming similar, they are completely different in terms of complexity. In the stress minimization problem (present case), feasibility is easily achieved. In the stress-constrained problem, feasibility is not an easy task to be accomplished.

The sensitivities of the objective function with respect to the topological variables can be derived as 

Derivatives ( j ) of displacement vectors with respect to design variables j and derivation of adjoint vectors are shown in the Appendix. The sensitivity numbers can be further simplified as

Recalling the effective constitutive matrix Eq. (17), the derivatives of the stiffness matrices with respect to topological variables can be expressed as where N,k , Tk,k and Td,k are the normal, sticking and sliding contact stiffness matrices of the k th element, (0) N,j , (0)

Tk,j and (0) Td,j are the elemental normal, sticking and sliding contact stiffness corresponding to solid material. k is the stiffness matrix of the k th element and (0) j is the elemental stiffness corresponding to solid material, which can be, respectively, written as

where Ω j denotes the domain of the j th element.

Substituting Eqs. (32)-(35) into Eq. (31), the sensitivities can be eventually evaluated as

where j is the vector of the adjoint nodal values of the j th element. (0) tan is the elemental tangential stiffness corresponding to solid material, which can be written as Particularly, for the frictionless contact problem, the sensitivities Eq. (38) can be simplified as where j should be changed accordingly.

The sensitivities denote the relative ranking of elemental contribution of the design objective are used to guide material removal and addition. Based on Eq. (38), the elemental sensitivity number can be obtained as

The volume fraction V (l) at current optimization iteration ( lth) is obtained by where c er denotes the evolutionary ratio.

Sensitivity numbers are firstly smoothed by where w sen kj denotes weight factor where r sen and Δ(k, j) are the sensitivity filter radius and the element center-to-center distance, respectively. The sensitivities are averaged with their historical information

Two threshold parameters th del and th add obtained from a iterative algorithm (Huang and Xie 2007) are, respectively, used to update the topological design variables for material removal and addition It is known that due to the highly nonlinear stress behavior, it is more preferable to filter topological variables rather than sensitivities for stress-related topology optimization. However, it needs to be emphasized that in BESO method the sensitivity filtering scheme serves also for material recovery from void to solid by attributing filtered sensitivity numbers to design variables that are associated to voids. Therefore, an additional filtering on the topological variables using the same scheme as Eq. (43) is proposed where the linear weight factor w den kj is defined as with density filter radius r den . After the topological variables are filtered, the topological variables change from discrete ( min or 1) to continuous (from min to 1). For the BESO method, discrete variables ( min or 1) are adopted. Before moving on to the next topology optimization iteration, the resulting topology with intermediate densities is converted into the discrete design through an algorithm step (Xia et al. 2018a ) with a determined threshold value x th satisfying the current material volume usage constraint.

The optimization procedure is terminated if the convergence criterion is satisfied where cr and N itr denote the allowable convergence error and an integral number which are set to be 0.001 and 5 throughout this paper, respectively.

The topology optimization procedure is outlined in Fig. 3 , and it mainly consists of the following stages. 

The effectiveness of the method proposed in this paper is certified through two comparison design problems in this section, including elastic-rigid and elastic-elastic contact problems in comparison with the maximum stiffness design. The influence of varying friction coefficient on the optimized results and the stress distributions are investigated. The effect of different parameters including p-norm, volume fraction and mesh density on the optimized results are discussed. The methodology is generally applicable to any predetermined contact surface, as long as the two structures are connected at all nodes of the contact surface. In order to simplify, this paper only studies two numerical examples of horizontal contact.

The geometric configuration and boundary dimensions of two cantilever beams contact are illustrated in Fig. 4 . Two cantilever beams are fixed on the left and right sides, respectively. The two cantilever beams are of the same dimensions, 200 in length and 50 in width, respectively. The two cantilever beams are divided into 50 and 200 elements along the vertical and horizontal directions, respectively. The dimension of the initial potential contact region is 100 (half the length of the cantilever beam). In order to avoid stress concentration near the external load, the external load F = 500 acts on five nodes (red nodes) at the top side of the free end. The Young's modulus E and Poisson's ratio are set to be 3 × 10 6 and 0.3, respectively. The densities of void and solid elements are set as min = 10 −6 and 1, respectively. The evolutionary volume ratio c er and the maximum volume addition ratio AR max are prescribed to be c er = 0.01 and AR max = 0.005 , respectively. The volume fraction, filter radii for sensitivity numbers and topological variables are set to be 50%, r sen = 3 and r den = 4 , respectively. The p-norm parameter and the allowable convergence error are set to be p = 6 and cr = 0.001 . Friction coefficient = 0.3.

The evolution of the topological structures and von Mises stress distribution under frictional elastic-elastic contact for stress minimization design is shown in Fig. 5 . It can be seen from Fig. 5 that the global p-norm stress and the maximal 

Comparisons of the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization and stiffness maximization designs are shown in Fig. 6 . Compared to the stiffness maximization design, the maximal von Mises stress for stress minimization design has been reduced for 8.26% from 225.11 to 206.52. This shows that the proposed stress design method can well control the stress level of the continuum structure under elastic contact with friction. The final iteration step of the stress minimization design is 76 steps, and the speed of the convergence is better than the stiffness maximization design which the final iteration step is 116. Since the two continuum structures are discontinuous in the contact region, it can be seen from Fig. 6 that the stress distributions at the contact regions, for the both cantilever beams, are discontinuous which is different from the design for a holistic continuum structure with continuous stress distribution. The evolution histories of global p-norm stress, maximal von Mises stress and volume fraction of frictional elasticelastic contact for stress minimization and stiffness maximization designs are shown in Fig. 7a and b, respectively. The red, green, blue curves and coordinate axes, in Fig. 7 , are used to present volume fraction, the maximum von Mises stress and the global p-norm stress, respectively. The final iteration steps of the topological optimization designs are 76 Fig. 7 The evolution histories of global p-norm stress, maximal von Mises stress and volume fraction of frictional elastic-elastic contact: a Stress design. b Stiffness design and 116 steps for stress and stiffness designs, respectively. BESO method starts from the full design and gradually decreases the total volume of the structure, until the required material volume fraction is met. Due to the highly nonlinear behavior of stress and the action of contact nonlinearity, as can be seen from the curves shown in Fig. 7 , the topological optimization iteration processes present oscillation. However, the curves converge to constant values when the volume fraction remains constant at the final stage.

The optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact and frictional elastic-rigid contact for stress minimization design are shown in Fig. 8 . Compared to the stress minimization design of elastic-elastic contact, the global p-norm stress and the maximal von Mises stress of elastic-rigid contact have been effectively reduced for 13.56% from 392.04 to 338.89, and 21.74% from 206.52 to 161.62, respectively. Compared to elastic-elastic contact, the elastic-rigid contact can better resist the deformation of the active contact structure and improve the strength of the active contact structure.

To compare the influence of the friction coefficient on the results, the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization design are shown in Table 1 . It can be concluded from Table 1 that as the friction coefficient increases, the global p-norm stress and the maximal von Mises stress decrease, that is, the strength increases. That's because the friction behavior resists the tangential deformation at the contact interface. Table 1 Influence of the friction coefficient on the optimized topological structures Page 12 of 22

The influence of the p-norm parameter on the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization design are shown in Table 2 .

The evolution histories of global p-norm stress, maximal von Mises stress and volume fraction of frictional elasticelastic contact for stress minimization and stiffness maximization designs are shown in Fig. 7a and b, respectively. The red, green, blue curves and coordinate axes, in Fig. 7 , are used to present volume fraction, the maximum von Mises stress and the global p-norm stress, respectively. The final iteration steps of the topological optimization designs are 76 and 116 steps for stress and stiffness designs, respectively. BESO method starts from the full design and gradually decreases the total volume of the structure, until the required material volume fraction is met. Due to the highly nonlinear behavior of stress and the action of contact nonlinearity, as can be seen from the curves shown in Fig. 7 , the topological optimization iteration processes present oscillation. However, the curves converge to constant values when the volume fraction remains constant at the final stage. Table 2 Influence of the p-norm parameter on the optimized topological structures Table 3 Influence of the volume fraction on the optimized topological structures

To compare the influence of the volume fraction on the topological results, the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization design are shown in Table 3 .

It can be concluded from Table 3 that the maximum von Mises stress might not decrease monotonically as the volume increases. The stress distributions at the contact regions, for the both cantilever beams, are discontinuous.

The second numerical example considered here is the two L-shaped brackets contact as shown in Fig. 9 . The geometric configuration and boundary dimensions of the two L-shaped brackets contact are illustrated in Fig. 9 .

Two L-shaped brackets are fixed on the top and bottom sides, respectively. The L-shaped brackets are of the same dimensions, 100 in length and 100 in width, respectively. The L-shaped brackets are divided into 100 elements along the vertical and horizontal directions. The dimension of the initial potential contact region is 30. In order to avoid stress concentration near the external load, the external load F = 800 acts on eight nodes (red nodes) at the top side of the free end. The Young's modulus E and Poisson's ratio are set to be 30 × 10 6 and 0.3, respectively. The densities of void and solid elements are set as min = 10 −3 and 1, respectively. The evolutionary volume ratio c er and the maximum volume addition ratio AR max are prescribed to be c er = 0.01 and AR max = 0.005 , respectively. The volume fraction, filter radii for sensitivity numbers and topological variables are set to be 50%, r sen = 3 and r den = 3 , respectively. The p-norm parameter and the allowable convergence error are set to be p = 6 and cr = 0.001 . Friction coefficient = 0.3.

The evolution of the optimized topological structures together with the von Mises stress distribution under frictional elastic-elastic contact for stress minimization design are shown in Fig. 10 . It can be seen from Fig. 10 that the stress concentration phenomena have occurred at the orthogonal corners of the two L-shaped brackets during the initial iteration step. With the removal of material at the orthogonal corner during the second iterative step, and the maximal von Mises stress has been effectively reduced for 21.65% from 552.84 to 433.13. The global p-norm stress and the maximal von Mises stress have been effectively reduced for 15.00% from 700.05 to 595.04, and 27.55% from 522.84 to 378.78, respectively. The original right angles of the corbel structure are optimized into arcs. The proposed stress minimization design method can well control the stress level and solve the stress concentration phenomenon of the structure under elastic contact with friction. The final iteration step of the optimization is 69 steps, this shows that the proposed method has good convergence.

Comparisons of the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization and stiffness maximization designs are shown in Fig. 11 . Compared to the stiffness maximization design, the maximal von Mises stress for stress minimization design has been reduced for 16.51% from 453.71 to 378.78. The final iteration step of the stress minimization design is 69 steps, and similarly, the speed of the convergence is better than the stiffness maximization design which the final iteration step is 72. Since the two continuum structures are discontinuous in the contact region, it can be seen from Fig. 11 that the stress distributions at the contact regions, for the both L-shaped brackets, are discontinuous which is different from the design for a holistic continuum structure with continuous stress distribution.

Different from the stiffness maximization design, the optimized topological structure of stress minimization design is more robust. Stiffness design is to maximize the stiffness of the structure to resist the deformation of the optimized topological structure, while stress design is to minimize the maximum stress of the structure to improve the strength of the optimized topological structure, and then resist the failure (damage) of the structure.

The evolution histories of global p-norm stress, maximal von Mises stress and volume fraction of frictional elastic-elastic contact for stress minimization and stiffness Fig. 12 . The evolution histories of global p-norm stress, maximal von Mises stress and volume fraction of frictional elastic-elastic contact for stress minimization and stiffness maximization designs are shown in Fig. 12a and b, respectively. The final iteration steps of the topological optimization designs are 69 and 72 steps for stress minimization and stiffness maximization designs, respectively. Due to the highly nonlinear behavior of stress and the action of contact nonlinearity, as can be seen from the curves shown in Fig. 12 , the topological optimization iteration processes present oscillation. However, the curves converge to constant values when the volume fraction remains constant at the final stage.

The optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact and frictional elastic-rigid contact for stress minimization design are shown in Fig. 13 . Compared to the stress design of elastic-elastic contact, the global p-norm stress and the maximal von Mises stress of elastic-rigid contact have been effectively reduced for 40.17% from 595.04 to 356.01, and 57.83% from 378.78 to 240.00, respectively. Compared to elastic-elastic contact, the elastic-rigid contact can better resist the deformation of the active contact structure and improve the strength of the active contact structure.

The influence of the friction coefficient on the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization design are shown in Table 4 .

It can be concluded from Table 4 that as the friction coefficient increases, the global p-norm stress and the maximal von Mises stress decrease, that is, the strength increases. That's because the friction behavior resists the tangential deformation at the contact interface.

The influence of the p-norm parameter on the optimized topological structures and von Mises stress distributions Table 5 .

It can be concluded from Table 5 that as the p-norm parameter increases, the global p-norm stress and the maximal von Mises stress decrease, that is, the strength increases. Moreover, as the p-norm parameter increases, the difference between the maximum von Mises stress and the global p-norm stress becomes smaller. In other words, as the p-norm parameter increases, more precise approximation of the maximum stress can be obtained. However, the global p-norm stress becomes ill-conditioned when p increases.

The higher value p takes, the severe oscillation happens since only the peak stress is targeted by the algorithm.

The influence of the mesh density nelx × nely on the optimized topological structures and von Mises stress distributions under frictional elastic-elastic contact for stress minimization design are shown in Table 6 . nelx and nely represent the numbers of elements along the vertical and horizontal directions. The time consumption of different cases under different mesh densities are also investigated in this subsection. The computer processor is 8.00 GB RAM, Intel (R), Core (TM) i5-7400 CPU @ 3.00 GHz. The total time consumption are shown in the last column in Table 6 .

It can be concluded from Table 6 that as the mesh density increases, the global p-norm stress and the maximal von Mises stress increase. This is due to the increase of meshing density, the reduction of element size, and the element stress closer to the point stress. In other words, the denser the mesh, the smaller the element size and the greater the element stress. But as the mesh density increases, as shown in Table 6 , the computational time increases exponentially.

The main contribution of this work is to accomplish stress minimization topology optimization for elastic continuum structures subjected to frictional contact and to study the influence of friction behavior on the stress-based topology optimization process and the optimized results. The well-known stress singularity problem is avoided using an extended BESO method. Stresses are treated as global performance (objective) function, the maximal von Mises stress is approximately evaluated using a global stress measure based on the p-norm. The friction behavior is governed by the Coulomb friction law, and the sensitivities are derived by adjoint method. Both sensitivities and topological variables are filtered to overcome the highly nonlinear stress behavior and stabilize the optimization procedure. The effectiveness of the proposed method is validated through two comparison studies including elastic-rigid and elastic-elastic contact problems. The influence of varying friction coefficient on the optimized results and the stress distributions are highlighted in comparison with the maximum stiffness design. The effect of different parameters including p-norm, volume fraction and mesh density on the optimized results are discussed. The main conclusions of this research are as follows.

1. The proposed approach can achieve a reasonable contact design that effectively controls the stress level and reduces the stress concentration effect at the critical stress areas. 2. Compared to elastic-elastic contact, the elastic-rigid contact can better resist the deformation of the active contact structure and improve the strength of the active contact structure. 3. Due to the friction behavior resists the tangential deformation at the contact interface. The global p-norm stress and the maximal von Mises stress decrease (strength increases) with the friction coefficient increases. 4. The global p-norm stress and the maximal von Mises stress decrease (strength increases) with the p-norm parameter increases. However, the global p-norm stress becomes ill-conditioned when p-norm parameter increases, since only the peak stress is targeted by the algorithm. 5. Due to the increase of meshing density, the reduction of element size, and the element stress closer to the point stress. With the mesh density increases, the global p-norm stress and the maximal von Mises stress increase.

The work presented in this paper opens up a new avenue for integrating stress design (strength design) and the topology design of structures under elastic frictional contact nonlinearity. The methodology is generally applicable to any predetermined contact surface, as long as the two structures are connected at all nodes of the contact surface. In order to simplify, this paper only studies two numerical examples of horizontal contact. However, additional researches are needed to include other nonlinear problems, such as stress-constrained material or/and geometrically nonlinear structures design under elastic frictional contact nonlinearity. Moreover, the fracture-resistance problem under elastic frictional contact nonlinearity is also addressed urgently. In the future work, the authors consider extending the current framework to more complex cases, such as contact problems with considering the uncertainty of contact locations or the uncertainty of contact shapes. Finally, the proposed method is also expected to be extended and applied for the elastic frictional contact design of multi-material structures under cyclic loads to control the stress level. 

The final converged tangential stiffness matrix last tan and the global displacement vector last of the last load increment can be obtained by the modified Newton-Raphson iterative method. It is assumed that the prescribed global external load vector is independent with the topological structure. However, the normal contact load vector N and the friction load vector T are dependent on the topological change. Similarly, the final normal contact load vector last N and the friction load vector last T can also be obtained in the last Newton-Raphson iteration. Thus, by differentiating both sides of the equilibrium equation in the state of the final converged load increment One obtains 

Generating optimal topologies in structural design using a homogenization method

Maximum stiffness design for elastic bodies in contact

Internal contact modeling for finite strain topology optimization

Topology optimization for minimum weight with compliance and stress constraints

Topology optimization of continuum structures with stress constraints and uncertainties in loading

Local versus global stress constraint strategies in topology optimization: a comparative study

Structural rigidity optimization with frictionless unilateral contact

Topology optimization of continuum structures with local stress constraints

A level set approach for topology optimization with local stress constraints

Stress-constrained level set topology optimization for design-dependent pressure load problems

Topology optimization for minimum mass design considering local failure constraints and contact boundary conditions

Topology optimization of multiscale elastoviscoplastic structures

Stress-related topology optimization via level set approach

Convergent and mesh-independent solutions for bi-directional evolutionary structural optimization method

Contact problems in elasticity: A study of variational inequalities and finite element methods

Stress-based topology optimization for continua

Topology optimization for structures with nonlinear behavior using the equivalent static loads method

Eliminate localized eigenmodes in level set based topology optimization for the maximization of the first eigenfrequency of vibration

Topology optimization of elastic contact problems using B-spline parameterization

Topology optimization of hyperelastic structures with frictionless contact supports

Equivalent static displacements method for contact force optimization

Topology optimization of elastic contact problems with friction using efficient adjoint sensitivity analysis with load increment reduction

An algorithm for optimal control problems based on differential inclusions

Piecewise constant level set method for structural topology optimization

Computational contact mechanics

Topology optimization for maximizing the fracture resistance of quasi-brittle composites

Bi-directional evolutionary structural optimization on advanced structures and materials: a comprehensive review

Stress-based topology optimization using bi-directional evolutionary structural optimization method

A simple evolutionary procedure for structural optimization

Numerical methods in contact mechanics

A linear relaxation model for shape optimization of constrained contact force problem

Topology optimization in aircraft and aerospace structures design

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgements Yongsheng Han would like to thank Professor Bin Xu for his crucial guidance and critiques, as well as Qian Wang for valuable discussions and suggestions during the pursuit of this work. This paper was partly completed during the epidemic of Covid-19. The authors dedicate this paper to Chinese doctors and nurses due to their great efforts and even their lives against disease.Funding This work was financially supported by the National Natural Science Foundation of China (11872311), the Innovation Foundation for Doctor Dissertation of Northwestern Polytechnical University (CX2021014).

The authors declare that they have no conflicts of interest.

All the necessary data and information to reproduce the results reported in this work are provided in the manuscript. The interested reader may contact the corresponding author for further implementation details.