Multi-material structural topology optimization considering material interfacial stress constraints

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 363 (2020) 112887 www.elsevier.com/locate/cma

Multi-material structural topology optimization considering material interfacial stress constraints Pai Liua , Litao Shib , Zhan Kanga ,∗ a

State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China b Aerospace System Engineering Shanghai, Shanghai 201109, China Received 23 July 2019; received in revised form 22 January 2020; accepted 29 January 2020 Available online xxxx

Abstract For multi-material structures, ensuring material interface strength is particularly vital for their integrity and durability. In the present study, we incorporate material interfacial stress constraints into topology optimization of bi-material structures. A multi-material level set method, in conjunction with interface-conforming finite element meshes, is employed to describe the distribution of different material phases and to capture the evolution of the material interfaces. The use of interface-conforming meshes enables accurate analysis of both interfacial stresses and their design sensitivities. Noting that the interfacial strength failure is usually characterized by a tension/compression asymmetric mechanism, we adopt an equivalent interfacial stress (which was originally proposed for strength criterion concerning composite delamination) expressed by the interface tensile and tangential stresses in the considered strength criterion. To handle the local nature of such interfacial stress constraints, we propose a global stress measure, which is an approximation of the p-norm of the equivalent interfacial stress field. An adjoint sensitivity analysis scheme is derived by taking into account the interface transmission conditions. To treat multiple constraints (the volume constraints and the interfacial stress constraint) in level set-based topology optimization, the velocity field-level set method is employed. Numerical examples are presented to show effectiveness of the present method. It is also shown that the tension/compression asymmetric interfacial strength criteria may lead to asymmetric designs. c 2020 Elsevier B.V. All rights reserved. ⃝ Keywords: Topology optimization; Multi-material; Interfacial stress constraint; Interfacial transmission condition; Interface-conforming mesh; Level set

1. Introduction Since the seminal work of Bendsœand Kikuchi [1] on the homogenization method, structural topology optimization has experienced a rapid development during the last three decades. Currently, topology optimization has been applied in many fields to provide conceptual designs. Three major topology optimization frameworks have become popular, which are the solid isotropic material with penalization (SIMP) method [2,3], the level set method [4–6], and the Evolutionary Structural Optimization (ESO) method [7,8]. The level set method provides an implicit description of material boundaries and also enables easier extraction of their geometrical information ∗ Corresponding author.

E-mail address: [email protected] (Z. Kang). https://doi.org/10.1016/j.cma.2020.112887 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 1. Schematic illustration of interface failures in a multi-material structure.

(normal direction and curvature), and thus, it is particularly suitable for structural boundary/interface-related design problems [9–12]. Many variants of the level set method have been developed to fulfill different design requirements (see e.g. [13–18]). Here, we only mention two of them, whose main techniques are employed in this study. The multi-material level set method (MM-LS) [16] uses m level set functions to describe m +1 material phases, and thus avoids redundant combinations of the level set functions in the design space. The velocity field level set (VFLS) method proposed by Wang and Kang [17] constructs the normal evolution velocity field with the velocity design variables and uses it as input of the Hamilton–Jacobi equation to advance the level set function. This formulation facilitates the use of general mathematical programming methods as the optimizer, and thus more conveniently handles multiple constraints in the level set framework. Material interface behavior is an important issue concerning the structural performance and failure in topology optimization. This is because interface-related issues account for a majority of failure modes for some types of multi-material structures [19]. Liu et al. [9] and Liu and Kang [10] studied multi-material topology optimization and multi-component integrated topology optimization considering debonding of the cohesive material interfaces, respectively. Behrou et al. [11] considered the interface cohesion in the topology optimization of a specific material anchor structure. Da et al. [20] optimized the distribution of inclusions in composite materials to maximize the fracture resistance under interface damage. These studies employed the cohesive zone model as a constitutive model of the material interfaces, which can describe the initial elastic stage and the subsequent softening stage of the interfaces. The cohesive zone model describes the relationship between the surface traction and its separation during interface failure. This phenomenon is characterized by an increase of the surface traction (for a given imposed strain), followed then by a decrease with separation of the interface. There are several types of cohesive models (e.g. the exponential form, the bi-linear form), and they are widely used in fracture mechanics or in the description of bonded interfaces. For a comprehensive review of such models, the readers are referred to the review paper of Park and Paulino [21]. Using the cohesive model, nonlinear finite element analysis is usually required to predict the interfacial responses. A reasonable design requirement of multi-material structures is to ensure that the material interfaces do not enter the softening (damage) stage when the structure is in service (a schematic illustration of a multi-material structure undergoing interface failure is shown in Fig. 1, which endangers the multi-material structural stiffness and integrity). Moreover, as material interfaces usually have a short elastic range, they can be simply modeled as perfectly bonded ones before entering the softening stage. Thus, the damage-free requirement of interfaces can be achieved by explicitly restricting the stresses at the perfectly bonded material interfaces in the topology optimization. In this particular circumstance, the difficulty of resolving the interface softening stage or the strong discontinuity of the displacement field at the material interfaces can be circumvented. When considering perfectly bonded multi-material structures with discontinuous material properties across the interfaces, the interface transmission conditions become a necessary ingredient to derive the exact sensitivities, as discussed in Bernardi and Pironneau [22], and Hettlich and Rundell [23]. In [24], Pantz noted an issue in the article of Bernardi and Pironneau [22], and proposed a suitable treatment of the problem by means of the fast C´ea’s Lagrangian method. Allaire et al. first derived exact analytical sensitivities of multi-material structures for the damage and fracture propagation problem [25], and the compliance minimization problem [26]. Allaire et al. [27]

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

3

also proposed a mesh evolution method that can explicitly discretize the structure and interfaces, and accurately handle the interface transmission conditions in the sensitivity analysis. Stress-related topology optimization problems have been extensively studied [28–33]. One of the major issues that requires proper handling in stress-related topology optimization problems is the singularity of the stress constraints and the resulting degeneration of the feasible domain [34]. The singularity issue of stress constraints is associated with the finite stress values of bar members with infinitesimal cross-sectional areas or elements with zero densities, which may lead to degenerate feasible domains with zero measure [35]. This makes a gradient-based optimizer difficult to arrive at the optimum. Cheng and Guo [36], Bruggi [37], and Burger and Stainko [38] proposed various stress constraint relaxation schemes. Another issue is the large number of local stress constraints. Different global stress constraints on the basis of aggregation functions, such as the Kreisselmeier–Steinhauser (K–S) function [39] and the p-norm function, and grouping schemes to enhance the control capability of the global stress constraints have been proposed (e.g. [40–43]). Recently, compression/tension asymmetric failure criteria were also studied in structural topology optimization by Bruggi and Duysinx [44], and Luo and Kang [45]. Using the stress constraints, Luo et al. [46] proposed a topology optimization method to seek wrinkle-free designs of thin membrane structures under stretching. To the best of the authors’ knowledge, the stress constraints in the presence of multiple-phase material interfaces have been seldom studied in multi-material topology optimization. In the present study, we consider the stress constraints at the bonding interfaces between different materials, in order to avoid entering the interface softening stage and material debonding. It is noted that the mechanisms for material interfacial failures and those occurring within a material body are essentially different. For instance, compressive stresses in the normal direction of a material interface may not cause interface debonding failure. Therefore, it would be questionable to simply resort to usual material strength measures (e.g. the von Mises stress) to ensure integrity of multi-material interfaces. We thus employ an interfacial stress criterion [47], which was originally proposed for predicting the delamination in laminated composite structures. On the basis of the level set model for the material interfaces and boundaries, an interface-conforming (i.e. boundary and interface-fitted) mesh is generated to explicitly discretize the design domain, which facilitates an accurate evaluation of the interfacial stresses and their design sensitivities. As we use bodyfitted analysis meshes in the finite element discretization and thus avoid intermediate density elements, the stress constraint singularity issue is not encountered in the present study. A p-norm type global interfacial stress measure based on the considered interfacial stress criterion is first proposed and then approximated with a domain integral over a narrow band along the material interfaces. In the optimization model, the strain energy is to be minimized under a constraint on this global interface stress measure. We formulate the topology optimization problem using the velocity field level set method, which enables the multi-constraint optimization problem to be solved with a generic mathematical programming method. Analytical sensitivities are derived with the adjoint method, which take into account the interface transmission conditions. Three numerical examples are presented to demonstrate the capabilities of the proposed method in controlling the material interfacial stresses. The organization of the present article is as follows: Section 2 introduces the adopted level set model, the interface transmission conditions, and the proposed global interfacial stress measure. Section 3 discusses the velocity field level set method first, and then presents the formulation of the optimization problem. The sensitivity analysis is given in Section 4. In Section 5, several numerical examples are shown. Conclusions are drawn in the last section. 2. Interface modeling in bi-material topology optimization In this section, the multi-material level set model and generation procedures of body-fitted meshes in this topology description framework are first introduced. Then, the interface transmission conditions and the proposed stress constraint at the material interfaces are presented. 2.1. Level set modeling of material interfaces As we are mainly concerned with the stresses at material interfaces, a crisp geometrical description of the interfaces becomes necessary in this study. The level set method implicitly describes the structural boundary as the zero level set of a one dimension higher function, and thus can more naturally handle topological changes. Therefore, it can be suitably used here to define the material interface locations and describe the bi-material distribution. In

4

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 2. A schematic illustration of bi-material distribution and the material interface.

the present study, the level set function values are associated with a fixed Cartesian grid and updated with the Hamilton–Jacobi equation. We employ a multi-material description model that is adapted from the MM-LS model [16]. Specifically, two level set functions are employed to describe the distribution of the two materials and void within the working domain. The characteristic functions are χ1 = H (ϕ1 ) (1 − H (ϕ2 )) , χ2 = H (ϕ2 ) , (1) χv = (1 − H (ϕ1 )) (1 − H (ϕ2 )) , in which the subscripts 1, 2, and v of χ respectively denote the two materials and the void phase; ϕ j ( j = 1, 2) are two level set functions, and H (·) is the Heaviside function: { ( ) 1 if ϕ j ≥ 0, H ϕj = (2) 0 if ϕ j < 0. It is noted that the employed MM-LS uses m LS level set functions to represent m LS + 1 material phase. As compared with the classical level set method for image segmentation [48] or the color level set method [15], which rely on four combinations of two level set functions for bi-material structures, the formulation of the characteristic functions in Eq. (1) effectively eliminates the need of artificial interpretation of the one redundant phase in the case of bi-material topology optimization [16]. A schematic illustration of this model is shown in Fig. 2. In the figure, Ω = Ω1 ∪ Ω2 represents the total material domain, D is the computational domain, and Σ is the interface between the two constituent materials. In particular, the material interfaces can be conveniently defined by ϕ1 > 0 and ϕ2 = 0, as shown by the red curve segment. The body-fitted mesh is constructed on the basis of the level set modeling. The procedures for generating such a mesh are as follows: we first bi-linearly seek the intersection points of the zero level sets ϕ j = 0 ( j = 1, 2) and the level set grid (as shown in Fig. 3(a)); then the interpolation curves for ϕ1 = 0 and ϕ2 = 0 are constructed respectively with these intersection points in the commercial software COMSOL; with the interpolation curves, we identify the portion of ϕ1 = 0 inside Ω2 (ϕ2 > 0) and delete them (see Fig. 3(b)); finally, we mesh the obtained geometry (excluding the void phase, where ϕ1 < 0 and ϕ2 < 0) with unstructured triangular elements (Fig. 3(c)). All these procedures are written in our Matlab code, which calls COMSOL during each iteration of optimization to generate the finite element mesh. Such a remeshing strategy is limited to 2D cases, as generating body-fitted meshes of complex 3D structures is still a challenging issue. As an alternative to the present meshing procedures, meshing algorithms that can generate body-fitted mesh on the basis of the level set functions (such as the method proposed by Allaire et al. [27]) can be adopted, especially for the 3D case. In the practical implementation, a generalized Hamilton–Jacobi equation with a diffusion term is employed to advance the level set functions. The diffusion term is applied as a regulation term, which penalizes the length of the structural boundaries represented by the level set functions ϕ1 and ϕ2 . The effect of this curvature term is similar to the perimeter constraint that controls the complexity of the optimized topology. The considered Hamilton–Jacobi

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

5

Fig. 3. A schematic illustration of the body-fitted meshing procedures: (a) bilinear search of the intersection points of the zero level sets and the level set grid, (b) creating the interpolation curves with the intersection points, (c) body-fitted meshing based on the created geometries.

Fig. 4. Illustration of level set function mapping between the finite element mesh and the level set grid.

equation is ⏐ ⏐ ⏐ ⏐ ∂ϕ j − V j ⏐∇ϕ j ⏐ − ccurv κ j ⏐∇ϕ j ⏐ = 0, (3) ∂t in which ccurv ⏐ is⏐)a constant for the curvature term [49], V j is the velocity for level set function ϕ j , and κ j = ( −∇ · ∇ϕ j / ⏐∇ϕ j ⏐ is the mean curvature of ϕ j . Re-initialization is implemented every five iteration steps to recover the signed distance property of the level set functions. As we have a body-fitted finite element mesh and a fixed Cartesian grid for the level set functions (see Fig. 4), we need to map the level set function values onto the finite element mesh in order to calculate the sensitivities (as

6

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

will be discussed in Section 4). For a point x in the body-fitted finite element mesh, its level set value is interpolated with the bi-linear shape functions. 2.2. Equilibrium equation and material interface transmission conditions In the present study, we consider linear elastic material behaviors under the small deformation assumption. The equilibrium equation reads ⎧ ( ) ⎨−div Aχ : ε (u) = f in Ω , u=0 on ΓD , (4) ⎩ Aχ : ε (u) · n = t on ΓN , ( ) in which ε (u) = ∇u + ∇uT /2 is the strain tensor, n is the unit outward normal vector to the structural boundary, f is the body force, t is the boundary traction, and Aχ is the elasticity tensor expressed by those of the candidate materials as Aχ = χ1 A1 + χ2 A2 ,

(5)

where the subscripts 1 and 2 denote material 1 and material 2, respectively. The transmission conditions of stresses must be satisfied at perfectly bonded material interfaces. These conditions require that the displacement u, the normal stresses σnn , σns and the tangent strain εss be continuous across the interface, while the normal strain εnn , εns and the tangent stresses σss be discontinuous. With the operator [·] = [·]1 − [·]2 expressing the jump of a quantity on the material interface (where the subscript denotes the side of the interface), we can express the interface transmission conditions for the stress field as [σnn ] = 0, [σns ] = 0,

(6)

[εss ] = 0. The quantities σnn , σns and εss on the material interface evaluated from either side can be computed with ( )⏐ σnn = nΣ · Aχ : ε (u) ⏐m · nΣ , m = 1 or 2, ( )⏐ σns = nΣ · Aχ : ε (u) ⏐m · sΣ , m = 1 or 2, εss = sΣ · ε (u)|m · sΣ ,

(7)

m = 1 or 2,

where m = 1, 2 denotes the material phases, nΣ (pointing from Ω2 to Ω1 ) and sΣ are the unit normal and tangential vectors of the material interface, which are expressed as −∇ϕ2 = {n Σ 1 n Σ 2 }T , (8) nΣ = |∇ϕ2 | sΣ = {n Σ 2 − n Σ 1 }T . (9) As stated in Allaire et al. [26], the interface transmission conditions are essential ingredients for a correct shapesensitivity analysis of the material interface. Therein, the authors derived the shape sensitivity in a sharp interface context, which includes the jump of the stress and strain on the material interface. These jumps thus need to be handled carefully, especially when using a fixed mesh. However, there is no trouble in computing these jumps in the present study, as we explicitly discretize the material interfaces and they can be easily evaluated from the elements on both sides of the interface. 2.3. Stress constraints at material interfaces A quadratic stress criterion was proposed by Brewer and Lagace [47] for the initiation of delamination in composites. We adopt this stress criterion for the considered bi-material interfaces. It is expressed by √ ( )2 1 t σint (x) = σnn + 2 σns 2 < cn , (10) β

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

7

t . Fig. 5. Illustration of the adopted interfacial stress criterion (a) the equivalent interfacial stress σint , (b) the expression of σnn

in which σint is the equivalent interfacial stress coupling the tensile stress and shear stress on the interface, β = cs /cn is the ratio between the maximum allowed interface shear strength cs and the interface normal tensile strength cn (in [47], cs and cn are referred to as “material system constants” and they can be obtained via material tests), and t σnn denotes the tensile stress at the interface, which is { σ if σnn ≥ 0, t σnn = nn (11) 0 if σnn < 0. The above interfacial stress criterion, which is tension/compression asymmetric, is schematically illustrated in Fig. 5. This stress criterion assumes that interface damage will occur when the equivalent interfacial stress σint exceeds cn . Eq. (10) is a point-wise constraint on the material interfaces. To handle the infinite number of such constraints, the p-norm or the K–S function can be used to aggregate the local constraints into a global one. Furthermore, t in Eq. (11), as σnn is not differentiable with respect to σnn at σnn = 0, we approximate its derivative at σnn = 0 with the average value of the derivatives on both sides when deriving the shape sensitivity. That is, its derivative is taken as 1/2 for σnn = 0. However, we point out here that the scenario of σnn = 0 (pure shear stress state at the interface) seldom occurs in practice. We consider a p-norm-type global stress measure. With a positive number P ≥ 1, we express the p-norm of the continuous equivalent interfacial stress field σint (x) on Σ as (∫ ) P1 ∥σint ∥ P ≡ |σint (x)| P dΓ . (12) Σ

To make the computation of sensitivity easier in the optimization process (as will be discussed in the sensitivity analysis section), considering the fact that the stress components used to compute the equivalent interfacial stress σint are continuous across the material interface, we define a global interfacial stress measure gσ in the form of a narrow band domain integral (see Fig. 6) to approximate ∥σint ∥ P (see Eq. (12)). The expression of gσ is obtained by first transforming Eq. (12) to a domain-integral form and then approximating the Dirac delta function with a smoothed function: (∫ ) P1 P gσ = H (ϕ1 ) δ (ϕ2 ) |∇ϕ2 | |σint | dΩ , (13) D

where δ (ϕ2 ) is a smoothed Dirac delta function (see Fig. 7): ( ⎧ 2) ⎨ 3 1 − ϕ2 , |ϕ | ≤ ∆, 2 δ (ϕ2 ) = 4∆ ∆2 ⎩ |ϕ2 | > ∆. 0,

(14)

where ∆ is half the width of the support domain of the smoothed Dirac delta function. The expression of this smoothed Dirac delta function is obtained as the derivative of a smooth polynomial approximation of the Heaviside function (see e.g. [5]).

8

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 6. Schematic illustration of the narrow band for stress constraint and the body-fitted mesh.

Fig. 7. The smoothed Dirac delta function.

Note that the narrow band-based expression Eq. (13) is only a formal approximation of the contour integralbased expression of the P-norm of the equivalent interfacial stress field. As gσ in Eq. (13) is computed in a narrow band of the material interface, the normal and tangential vectors (of non-zero level set contours) used to compute interfacial stress also need to be calculated in the neighborhood of the interface as n (x) = −∇ϕ2 (x) / |∇ϕ2 (x)| = [n 1 (x) n 2 (x)], s (x) = [n 2 (x) n 1 (x)]. ⏐ ⏐ Because the signed distance property (⏐∇ϕ j ⏐ = 1, j = 1, 2) of the level set functions is preserved with the re-initialization processes in our framework, the proposed global interfacial stress measure gσ in Eq. (13) can be further expressed as ) P1 (∫ . (15) gσ = H (ϕ1 ) δ (ϕ2 ) |σint | P dΩ D

3. Topology optimization considering interfacial stress constraints In the optimization model, we aim to minimize the strain of a bi-material structure. It can be expressed as ∫ 1 f = σ : εdΩ . (16) 2 Ω1 ∪Ω2

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Furthermore, we add the interface length to the objective function as a regularization term, which is ∫ lint = H (ϕ1 ) dΓ .

9

(17)

ϕ2 =0

Thus, the objective function is written as J =γ

f f

+ (1 − γ )

lint l int

,

(18)

in which γ is a weight coefficient, and f and l int are given reference values to scale the strain energy and the interface length. To handle multiple constraints in the level set framework, we use the VFLS method [17] as an alternative to the augmented Lagrange multiplier algorithm-based implementation of the level set method. Here we define the velocity design variables on the Cartesian level set grid, and bi-linearly interpolate the normal velocity fields Vn j of the zero level sets with the velocity design variables V ji defined on the level set grid points (i.e. the points of the Cartesian grid for the level set functions) as Vn j (x) =

N ∑

{ } Ni (x) V ji , x ∈ x|ϕ j (x) = 0 , j = 1, 2,

(19)

i=1

in which the subscript j denotes the corresponding level set function, N denotes the total number of the level set grid points, and Ni is the bi-linear shape function defined on the level set grid. The VFLS method seeks the optimal value of these velocity design variables with a general mathematical programming method (MMA here). In the level set framework, the topology optimization problem with multiple constraints can also be solved with other approaches. It is noted that Dunning and Kim [50] also proposed a method to handle multiple constraints in the level set framework, but using the sequential linear programming (SLP) concept and a different way to construct the boundary normal velocities. In [51,52], the level set method was implemented in conjunction with different customized optimizers. The considered optimization problem can thus be formulated as Find Vd = {V11 , . . . , V1N , V21 , . . . , V2N }T min

J

s.t.

( ( ) ) α u ϕ j , v = l (v) , ∀v ∈ Uadm , gσ ≤ cn ,

(20)

gV1 = f V1 − f V1 ≤ 0, gV2 = f V2 − f V2 ≤ 0, Vmin ≤ V ji ≤ Vmax , j = 1, 2; i = 1, . . . , N , (∫ ) (∫ ) in which Vd is the vector of design variables, f V1 = D χ1 dΩ /AD and f V2 = D χ2 dΩ /AD are the volume fractions of material 1 and 2, respectively (AD is the volume of the design domain), f V1 and f V2 are the imposed amount of each material, Vmin and Vmax are respectively the lower and upper limits of the design variables, which are determined by fulfilling the CFL condition of the Hamilton–Jacobi equation. In the variational form of the equilibrium state, Uadm is the space of kinematically admissible displacements. The energy bilinear form and the load linear form in the equilibrium equation are ∫ α (u, v) = σ (u) : ε (v) dΩ , Ω1 ∪Ω2

l (v) =

∫

∫ t · vdΓ + Γ1N ∪Γ2N

(21) f · vdΩ .

Ω1 ∪Ω2

In the optimization model, the bilinear shape functions defined on the level set grid are used to interpolate the boundary normal velocities with the design variables, as shown in Eq. (19).

10

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 8. Flowchart of the optimization procedures.

During the optimization process, the level set functions of the current design (the qth iteration) is approximated with the normal velocity Vd as ⏐ ⏐ q ϕ j ≈ q−1 ϕ j + Vn j ∆t ⏐∇ϕ j ⏐ , (22) Here, ∆t is a pseudo-time step. For the optimization formulation in Eq. (20), the design sensitivities of the objective function and interfacial stress constraint with respect to Vd can be derived with the adjoint method. With these sensitivities as inputs, we employ a general mathematical programming method (here, the MMA method [53]) to optimize the design variables Vd . Then, the normal evolution velocities on the structural boundaries and interfaces are constructed as in Eq. (19). As will be introduced in Sections 4 and 5, we first use the velocity design variables to interpolate the normal velocities at a narrow band around the structural boundary, and then compute the sensitivities with contour and domain integration. The adopted narrow band width for the proposed global interfacial stress constraint is relatively small as compared with the size of the finite elements. We need to extend the constructed boundary velocities to the entire design domain. Here, this velocity extrapolation process is implemented with the fast marching method [54]. A flowchart is given in Fig. 8 to show the optimization procedures. It is seen that the optimization procedure for the optimization formulation in Eq. (20) enables convenient handling of multiple constraints in the level set-based optimization framework.

11

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

4. Sensitivity analysis of the sharp-interface problem When deriving the sensitivities of a multi-phase material design under weak discontinuity of the displacement field (which means that the displacement itself is continuous but its gradient is discontinuous) across the interface, a vital issue needs to be carefully taken care of. That is, the state variable u, as a solution to the equilibrium equation in Eq. (4), is not shape differentiable with respect to the interface location (see [22,25]). In order to resolve this issue, we adopt the method suggested by Allaire et al. [25], who expressed the equilibrium state of a two-material solid as a two-phase system coupled by the transmission conditions. With similar notations as in Allaire et al. [25], we consider the restrictions u1 and u2 of u to Ω1 and Ω2 respectively, and define them as the displacement solutions to the following problem ⎧ −div (Am : ε (um )) = fm in Ωm ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ A : ε (um ) · nm = tm on ΓmN ⎪ ⎨ m um = 0 on ΓmD (23) (m = 1, 2) ⎪ ⎪ ⎪ ⎪u1 = u2 on Σ ⎪ ⎪ ⎪ ⎩ −A1 : ε (u1 ) · nΣ + A2 : ε (u2 ) · nΣ = 0 on Σ , where the subscript m indicates the corresponding quantities on Ωm , and ΓmN = ∂Ωm ∩ ΓN , ΓmD = ∂Ωm ∩ ΓD . From Eq. (23), for a general structural response J˜, one can construct the Lagrangian function as 2 ∫ ∑ q1m · (div (Am : ε (um )) + fm ) dΩ L = J˜ + m=1 Ωm

+

2 ∫ ∑

q2m · (Am : ε (um ) · nm − tm ) dΓ +

m=1 ΓmN

∫

q4 · (u1 − u2 ) dΓ +

+ Σ

2 ∫ ∑

q3m · um dΓ

(24)

m=1 ΓmD

∫

q5 · (−A1 : ε (u1 ) · nΣ + A2 : ε (u2 ) · nΣ ) dΓ ,

Σ

in which q1m , q2m , q3m (m = 1, 2), q4 , and q5 are the Lagrangian multipliers. They can be determined by ⟨∂ L/∂um , δum ⟩ = 0 (m = 1, 2). We now consider the p-norm type global interfacial stress measure ∥σint ∥ P given in Eq. (12), which is expressed with an integral on the material interfaces. The relationship ⟨∂ ∥σint ∥ P /∂um , δum ⟩ = 0 is used to derive the adjoint equations. The expression of ⟨∂ ∥σint ∥ P /∂um , δum ⟩ is ⟩ (∫ ) P1 −1 ∫ ( ) ⟨ 1 ∂ |σint | ∂ ∥σint ∥ P |σint | P dΓ , δum = P |σint | P−1 δ (∇um ) dΓ . (25) ∂um P ∂∇um Σ Σ The integral on the right-hand side of Eq. (25) is a boundary integral, and its integrand involves the gradient of um , which increases the difficulty to solve the resulting adjoint equation. To make the sensitivity computation more convenient, the proposed domain-integral form of global interfacial stress measure gσ (see Eq. (13)) is used. Although we do not use the contour-integral form interfacial stress measure (see Eq. (12)) in the present study, the sensitivity of this stress measure can be derived with the global interfacial stress measure, J˜ in Eq. (24) is replaced with (∑ the∫ method in [25]. With gσ )being 1/P 2 P , and the corresponding Lagrangian function is denoted by L σ . Since m=1 Ωm H (ϕ1 ) δ (ϕ2 ) |σint | m dΩ (∑ )1/P ∫ 2 P |σ | the sensitivities of H δ dΩ can be easily obtained with the chain rule using the (ϕ ) (ϕ ) 1 2 int m m=1 Ωm ∫ ∑ sensitivities of the term in the parenthesis, we only consider the case of J˜ = 2m=1 Ωm H (ϕ1 ) δ (ϕ2 ) |σint | P m dΩ in the following derivation. With ⟨∂ L σ /∂um , δum ⟩ = 0, we can derive the adjoint equation, resolve the relationships between the Lagrangian multipliers and express them with expressions of q1m . Thus the Lagrangian function L σ can be further expressed

12

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

as Lσ =

2 ∫ ∑

H (ϕ1 ) δ (ϕ2 ) |σint | P m dΩ

m=1 Ωm

−

2 ∫ ∑

ε (q1m ) : Am : ε (um ) dΩ +

m=1 Ωm

+

2 ∫ ∑

2 ∫ ∑

q1m · fm dΩ +

m=1 Ωm

(

2 ∫ ∑

q1m · tm dΓ

m=1 ΓmN

) q1m · Am : ε (um ) · nm + um · Am : ε (q1m ) · nm − H (ϕ1 ) δ (ϕ2 ) um · σˆ m · nm dΓ

m=1 ΓmD

(26)

∫ 1 (u1 − u2 ) · (A1 : ε (q11 ) + A2 : ε (q12 )) · nΣ dΓ 2 Σ ∫ 1 − (q11 − q12 ) · (A1 : ε (u1 ) + A2 : ε (u2 )) · nΣ dΓ 2 Σ ∫ ( ) 1 + (u1 − u2 ) · H (ϕ1 ) δ (ϕ2 ) σˆ 1 + H (ϕ1 ) δ (ϕ2 ) σˆ 2 · nΣ dΓ . 2 Σ ∫ ∑ The derivative of the term 2m=1 Ωm H (ϕ1 ) δ (ϕ2 ) |σint | P m dΩ in Eq. (26) with respect to um can be written as −

2 ∫ ∑

) ( ∂ H (ϕ1 ) δ (ϕ2 ) |σint | P m δum dΩ ∂um m=1 Ωm ( ) 2 ∫ ∑ ∂ |σint | P m = H (ϕ1 ) δ (ϕ2 ) : Am : δε (um ) dΩ ∂σm m=1 Ωm =

2 ∫ ∑ m=1 ∂Ωm

−

2 ∫ ∑

(27)

H (ϕ1 ) δ (ϕ2 ) σˆ m · nm · δum dΓ

( ) div H (ϕ1 ) δ (ϕ2 ) σˆ m · δum dΩ

m=1 Ωm

( ) Here, we define σˆ m ≡ Am : ∂ |σint | P m /∂σm . The adjoint equations obtained with ⟨∂ L σ /∂um , δum ⟩ = 0 (m = 1, 2) read ⎧ ) ( ⎪ in Ωm div (Am : ε (q1m )) − div H (ϕ1 ) δ (ϕ2 ) σˆ m = 0 ⎪ ⎪ ⎪ ⎪ ⎪ on ΓmD ⎪ ⎨q1m = 0 Am : ε (q1m ) · nm = H (ϕ1 ) δ (ϕ2 ) σˆ m · nm ⎪ ( ) ⎪ ⎪ ⎪ (A1 : ε (q11 ) − A2 : ε (q12 )) · nΣ = H (ϕ1 ) δ (ϕ2 ) σˆ 1 − σˆ 2 · nΣ ⎪ ⎪ ⎪ ⎩ q11 − q12 = 0

on ΓmN

(28)

on Σ on Σ

In the present study, we assume that the boundaries with prescribed displacements and the boundaries with traction forces are fixed in the design problem. Furthermore, in the sensitivity derivation, we adopt the approximation ⏐ ⏐ ψ j = Vn j ⏐∇ϕ j ⏐ in Eq. (22) (ψ j is the perturbation of ϕ j , and a unit pseudo-time step is used here). Thus the Fr´echet derivatives of L σ with respect to ϕ j ( j = 1, 2) have the following expressions ⟨ ⟩ ∫ ∂ Lσ , ψ1 = δ (ϕ2 ) |σint | P 1 Vn1 dΓ ∂ϕ1 ϕ1 =0 (29) ∫ + (−ε (q11 ) : A1 : ε (u1 )) (1 − H (ϕ2 )) Vn1 dΓ . ϕ1 =0

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

⟨

⟩ ∫ ) ( ′ ( ) ∂ Lσ , ψ2 = H (ϕ1 ) δ (ϕ2 ) |σint | P + δ (ϕ2 ) |σint | P ,n n2,ϕ2 Vn2 |∇ϕ2 | dΩ 2 ∂ϕ2 D ∫ + (ε (q11 ) : A1 : ε (u1 ) H (ϕ1 ) − ε (q12 ) : A2 : ε (u2 )) Vn2 dΓ ϕ2 =0 ∫ 1 ∂ (u1 − u2 ) − · (A1 : ε (q11 ) + A2 : ε (q12 )) · nΣ Vn2 dΓ 2 Σ ∂nΣ ∫ 1 ∂ (q1 − q2 ) − · (A1 : ε (u1 ) + A2 : ε (u2 )) · nΣ Vn2 dΓ 2 Σ ∂nΣ ∫ ) ∂ (u1 − u2 ) ( 1 + · H (ϕ1 ) δ (ϕ2 ) σˆ 1 + H (ϕ1 ) δ (ϕ2 ) σˆ 2 · nΣ Vn2 dΓ 2 Σ ∂nΣ

13

(30)

in which n2,ϕ2 is defined as the derivative of n2 with respect to ϕ2 , and the term n2,ϕ2 Vn2 |∇ϕ2 | can be expressed as ⟩ (( ) ) ⟨ ∂ Vn2 ∂ Vn2 ∂n2 , ψ2 = −∇s Vn2 = − e1 + e2 · s s. (31) ∂ϕ2 ∂x ∂y For the strain energy, the Lagrangian function becomes ∫ 2 ∑ 1 LC = ε (um ) : Am : ε (um ) dΩ 2 Ωm m=1 ) ∑ ∫ 2 ∫ 2 ( ∫ ∑ q1m · fm dΩ + − + (ε (q1m ) : Am : ε (um )) dΩ + +

m=1 2 ∫ ∑

Ωm

Ωm

q1m · tm dΓ

m=1 ΓmN

(q1m · Am : ε (um ) · nm + um · Am : ε (q1m ) · nm − um · Am : ε (um ) · nm ) dΓ

m=1 ΓmD

(32)

∫ 1 (u1 − u2 ) · (A1 : ε (q11 ) + A2 : ε (q12 )) · nΣ dΓ 2 Σ ∫ 1 − (q11 − q12 ) · (A1 : ε (u1 ) + A2 : ε (u2 )) · nΣ dΓ 2 Σ ∫ 1 + (u1 − u2 ) · (A1 : ε (u1 ) + A2 : ε (u2 )) · nΣ dΓ 2 Σ

−

From ⟨∂ L C /∂um , δum ⟩ = 0 (m = 1, 2), we obtain the following adjoint equations ⎧ ⎪ ⎪div (Am : ε (q1m )) − div (Am : ε (um )) = 0 in Ωm ⎪ ⎪ ⎪ ⎪ ⎪ on ΓmD ⎨q1m = 0 Am : ε (q1m ) · nm = Am : ε (um ) · nm on ΓmN ⎪ ⎪ ⎪ ⎪(A1 : ε (q11 ) − A2 : ε (q12 )) · nΣ = 0 on Σ ⎪ ⎪ ⎪ ⎩ q11 − q12 = 0 on Σ

(33)

It is seen that the problem is self-adjoint for the strain energy function, i.e. q1m = um . The Fr´echet derivative of L C with respect to ϕ j ( j = 1, 2) can be derived as ⟩ ⟨ ∫ 1 ∂ LC , ψ1 = − ε (u1 ) : A1 : ε (u1 ) (1 − H (ϕ2 )) Vn1 dΓ (34) ∂ϕ1 ϕ1 =0 2 ⟨ ⟩ ( ) ∫ ∂ LC 1 1 , ψ2 = ε (u1 ) : A1 : ε (u1 ) H (ϕ1 ) − ε (u2 ) : A2 : ε (u2 ) Vn2 dΓ ∂ϕ2 2 ϕ2 =0 2 (35) ∫ 1 ∂ (u1 − u2 ) − · (A1 : ε (u1 ) + A2 : ε (u2 )) · nΣ Vn2 dΓ 2 Σ ∂nΣ

14

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

For the material volume constraint and the interface length constraint, the sensitivities are directly derived as ⟨ ⟩ ∫ ∂ f V1 ϕ =0 (1 − H (ϕ2 )) Vn1 dΓ , ψ1 = 1 . (36) ∂ϕ1 AD ∫ ⟨ ⟩ − Σ Vn2 dΓ ∂ f V1 , ψ2 = . (37) ∂ϕ2 AD ⟨ ⟩ ∫ ∂ f V2 ϕ =0 Vn2 dΓ , ψ2 = 2 . (38) ∂ϕ2 AD ⟨ ⟩ ∫ ∂lint , ψ2 = κ2 Vn2 dΓ . (39) ∂ϕ2 Σ 5. Numerical examples This section presents numerical examples regarding a cantilever beam, an L-shaped beam, and a simply supported beam with consideration of interfacial stress constraints. Body-fitted meshes composed of triangular elements are employed here to discretize the structures. The maximum and minimum element edge lengths are 1 and 0.5, respectively. The quadratic elements are used for finite element analyses; while the fourth-order elements are used for solving the adjoint equation corresponding to the global interfacial stress measure, in order to resolve possibly existing large gradients of the adjoint variables in the narrow band around the material interface. The level set grid has unit distances between each pair of adjacent points, and is fixed during the design iterations. The CFL condition for solving the Hamilton–Jacobi equation is set such that the step length of level set function evolution is less than 0.5, and the coefficient of the curvature term in the Hamilton–Jacobi equation is 0.05. Unless otherwise stated, the narrow width for evaluating the global interfacial stress measure is 0.2. According to our numerical experience, a large value of P (e.g., P = 15) in the global interfacial stress measure usually leads to localized sensitivities, while a too small value of P (e.g. P ≤ 2) may lead to a weaker control of the peak stress value along the portion of interfaces that has relatively large stresses. We found P = 4 (which is also suggested by Duysinx and Sigmund [40]) to be a suitable choice for the considered problem. In the present study, the normal and shear stresses at the material interfaces are obtained by averaging the stresses computed at the edges of the adjacent elements using the gradient of the displacement field. 5.1. Design of a cantilever beam We first consider the design optimization of a cantilever beam, as shown in Fig. 9, under the interfacial stress constraint. The Young’s moduli of the materials are E 1 = 3 × 104 , E 2 = 6 × 104 , and the Poisson ratios are v1 = v2 = 0.3. The material volumes are restricted by f V1 = 0.25 and f V2 = 0.3. The weight coefficient in the objective function is γ = 0.8. The initial design is given in Fig. 10(a), in which the blue color represents the weaker material, and the red color represents the stronger material. The equivalent interfacial stress distribution, which is non-symmetric, for the initial design is shown in Fig. 10(b).

Fig. 9. The cantilever beam design problem.

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

15

Fig. 10. Initial design of the cantilever beam: (a) the material distribution, (b) the equivalent interfacial stress distribution. . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. The optimized designs of the cantilever beam under different constrained interfacial stress values: (a) cn = 25, (b) cn = 10, (c) cn = 7.5, (d) cn = 5, and (e) cn = 3.5.

Firstly, we set β = 2 in the interfacial stress constraint (which indicates that the interfacial shear strength is twice that of the normal tensile strength) and consider different threshold values cn = 25, 10, 7.5, 5 and 3.5 (for the definitions of β, cn and cs , see Eq. (10) and the explanations). For the cases of cn = 25, 10, 7.5 and 5, the optimized designs obtained after 600 iterations are shown in Fig. 11(a–d). For the case of cn = 3.5, as the constraint was more stringent, the optimization process was terminated after 800 iterations, yielding the optimized design shown in Fig. 11(e). In this case, the parameter P in the p-norm function was increased from 4 to 6 during the last 100 iterations to enhance better control of the localized peak interfacial stress. It is seen that the optimized

16

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887 Table 1 The strain energy, global interfacial stress measure, and maximum interfacial stress of the optimized cantilever beam designs for different constrained interfacial stress values. cn

f

gσ

Max. σint

25.0 10.0 7.5 5.0 3.5

5.79 5.76 6.40 5.76 5.76

25.10 10.47 7.78 5.08 3.54

17.43 9.42 5.06 3.14 2.18

Fig. 12. Equivalent interfacial stresses of the optimized cantilever beams under different constrained interfacial stress values: (a) cn = 25, (b) cn = 10, (c) cn = 7.5, (d) cn = 5, and (e) cn = 3.5.

design varies as the constraint value of interfacial stress decreases. In particular, when the constrained value is large enough, the optimized design exhibits a symmetric topology (see Fig. 11(a)), just as a conventional solution without interfacial strength consideration does. For a relatively stringent interfacial stress constraint, the optimized designs are asymmetric and the bi-material interfaces mainly locate in the compressive stress zones to avoid interfacial failure (see Fig. 11(c–e)). This is a major difference between multi-material topology optimization solutions with and without interface strength constraints. In this example, though the initial design and the loading condition are symmetric with respect to the center line, the distribution of the equivalent interfacial stresses defined by Eqs. (10) and (11) is non-symmetric. This explains why the optimized structures, except for the case of a very large allowable interfacial stress limit cn = 25, exhibit non-symmetric geometries.

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

17

Fig. 13. Optimized designs and material interfacial stresses of the cantilever beam obtained with different narrow band widths for evaluating the global interfacial stress: (a) 2∆ = 0.8; (b) 2∆ = 1.2; (c) 2∆ = 1.5.

The strain energy values, global interfacial stress measure values, and maximum interfacial stresses are listed in Table 1. As can be observed, for different constrained interfacial stress values cn , the strain energy has similar values except for the case of cn = 7.5, which may be caused by the occurrence of a local minimum and high nonlinearity of the global interfacial stress constraint. The global interfacial stress measure gσ is very close to its specified bound cn , and the maximum interfacial stress σint is well below cn , which shows effectiveness of the global interfacial stress constraint. Some slight violations of the global interfacial stress measure may be due to the relatively weak control capability of the p-norm type global interfacial stress measure on localized stress concentration along the interfaces and the high nonlinearity of this constraint. To verify the optimized designs, the equivalent interfacial stresses expressed in Eq. (10) are depicted in Fig. 12. It is seen that the proposed global interfacial stress constraint can effectively restrict the peak interfacial stress values for all cases of the constrained values. In order to explore the effects of the narrow band width, we reran the optimization process for the case cn = 10 using an increased narrow-band width 2∆ = 0.8, 1.2 and 1.5 (where ∆ is half the narrow band width). The optimized designs are depicted in Fig. 13. The strain energies for these three cases are very close to each other (5.76, 5.74 and 5.72), and the maximum equivalent interfacial stresses are respectively 9.45, 3.08 and 4.23. As seen in Fig. 13(a), the optimized design obtained with relatively smaller narrow-band width for evaluating the interfacial stress constraint is similar to the design shown in Fig. 11(b). When larger narrow-band width are used, the optimized designs only show some different local features, as shown in Fig. 13(b) and (c). For the three considered narrow band widths, it is observed that the optimized designs with larger narrow band widths have smaller maximum equivalent interfacial stress values than the threshold value. This may be caused by the occurrence of local minima due to the high nonlinearity of the proposed global interfacial stress constraint.

18

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 14. Iteration history of the cantilever beam design for the case of cn = 5: (a) strain energy and material volumes; (b) global interfacial stress measure.

The iteration histories for the case of constrained value cn = 5 is shown Fig. 14. As can be seen, the interfacial stress measure has some oscillations during the first three hundred iterations, and then achieves a steady convergence. The oscillations may be caused by the breakage of the bars and the formation of new material interfaces in the first stage of the optimization iterations, which may temporarily introduce large material interfacial stresses. As expected, the volume fraction constraints for all the optimized designs shown in Figs. 11 and 13 are active. Next, we change the value of β (the ratio between the shear and the normal interface strength) to examine its effect. In the first case, β = 10 is considered and the constrained interfacial stress is cn = 5. The optimized design and its equivalent interfacial stress are plotted in Fig. 15. It is observed that although the interfacial shear strength in this case is much larger (i.e. cs = β × cn = 50), the optimized design is still similar to the design obtained with β = 2 and cn = 5 (and thus cs = 10). This indicates that the optimized designs are insensitive to the interfacial shear strength when it is too large in comparison with the interfacial normal strength. In both cases of β = 2 and β = 10, the interfacial normal tensile behavior is the most critical one in the equivalent interfacial stress. Again, the volume fraction constraint is active in the optimized design. The strain energy of the optimized design in Fig. 15 is 5.73, and the global interfacial stress measure gσ is 3.38. The maximum interfacial stress is 1.92, which is again well restricted by the constraint.

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

19

Fig. 15. The optimized design of the cantilever beam obtained with β = 10 and cn = 5: (a) material distribution and (b) interfacial stress.

Fig. 16. The optimized design of the cantilever beam obtained with β = 0.5 and cn = 15: (a) material distribution and (b) interfacial stress.

Further, we study a case where β is smaller than 1 (i.e. the interface shear stress becomes more critical). Specifically, we set β = 0.5 and cn = 15 (i.e. cs = 7.5). The optimized design and its equivalent interfacial stresses are shown in Fig. 16. From this figure, it is obvious that the topology of the optimized design differs from those obtained when the interface normal stresses are more critical. The optimized design for this setting has a larger interface length, which is beneficial to decrease the interfacial shear stress. The strain energy of the optimized design is 5.99, and the global interfacial stress measure is 14.91. The maximum interfacial stress is 11.79, which is again well restricted. The volume fraction constraint is also active in the optimized design. 5.2. Design of an L-shaped beam In this example, we consider the design problem of an L-shaped beam as shown in Fig. 17. The material properties are E 1 = 3 × 104 , E 2 = 1.5 × 105 , and v1 = v2 = 0.3, in which the ratio of Young’s moduli between the

Fig. 17. The L-shaped beam design problem.

20

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 18. Initial design of the L-shaped beam.

Fig. 19. The optimized designs of the L-shaped beam (a) with interfacial stress constraint and (b) without interfacial stress constraint. . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 2 Comparison of the strain energy values, global interfacial stress measure and maximum interfacial stresses between the optimized L-shaped beam designs with and without the interfacial stress constraint.

f gσ Max. σint

With interfacial stress constraint

Without interfacial stress constraint

5.15 5.15 3.40

5.08 – 9.63

strong and the weak material is larger than that used in the previous example. The imposed amounts of material are f V1 = 0.14 and f V2 = 0.29. The weight coefficient in the objective function is γ = 0.98. The parameters of the global interfacial stress constraint are β = 2 and cn = 5. The initial design is shown in Fig. 18. The optimization process was terminated after 1000 iterations, yielding an optimized design as shown in Fig. 19(a). In this design problem, the parameter P in the p-norm is increased from 4 to 6 during the last 100 iterations to enhance the control over the peak interfacial stress. The iteration history is shown in Fig. 20. It is seen that the global interfacial stress measure has some oscillations during the optimization iterations, which again may be due to breakage of bar members, as well as formation and evolution of new material interfaces. For comparison, we also plot the design optimized solely for minimum strain energy in Fig. 19(b). The two designs have distinctly different topologies. In the former design, the weaker material (in blue color) tends to attach

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

21

Fig. 20. Iteration history for the L-shaped beam design with interfacial stress constraint: (a) strain energy and material volumes; (b) global interfacial stress measure.

Fig. 21. Equivalent interfacial stresses of the optimized L-shaped beams (a) with interfacial stress constraint and (b) without interfacial stress constraint.

22

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

to the stronger material members (in red color) as reinforcements, while in the latter design the weaker material forms independent load-bearing members. For minimum strain energy topology optimization of multi-material structures, the formation of independent structural members by the weaker material is also observed by Wang et al. [16], under the circumstance that the Young’s moduli of the stronger and weaker materials have comparable orders of magnitudes (e.g. steel and aluminum). The strain energy, the global interfacial stress measure and the maximum interfacial stress for the two designs are listed in Table 2. Compared with the minimum strain energy design without interfacial strength constraints, the optimized design obtained with the present method has a slightly lower stiffness, but a significantly smaller maximum interfacial stress (the equivalent interfacial stress of both designs are shown in Fig. 21). It is seen that the interfacial stress is effectively controlled in the design obtained with the present method (Fig. 21(a)), which highlights the significance of the interfacial stress constraint in this bi-material design problem. For the L-shaped beam obtained under the material interfacial stress constraint, the von Mises stress at the reentrant corner is predicted to be 210.40. Restricting the maximum stresses of the bulk material has been studied by many researchers, but is not concerned in this study. 5.3. Design of a simply supported beam The topological design of a simply supported beam shown in Fig. 22 is treated in this example. The material parameters are E 1 = 3 × 104 , E 2 = 1.5 × 105 , and v1 = v2 = 0.3. The given material volume ratios are f V1 = 0.05 and f V2 = 0.25. We set the weight coefficient to be γ = 0.5, and impose a global interfacial stress constraint with cn = 3.5 and β = 2. The initial design is depicted in Fig. 23. Here, we use a very thin layer of the weaker material (blue) attached to the stronger material (red) to approximately satisfy the relatively low volume fraction constraint of the weaker material (which is 0.05). The optimization process is terminated after 900 iterations. The iteration history is shown in Fig. 24. Again, in the optimized design, the volume fraction constraints are active.

Fig. 22. The simply supported beam design problem.

Fig. 23. Initial design of the simply supported beam and an enlargement of the material distribution. . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

23

Fig. 24. Iteration history for the simply supported beam design with interfacial stress constraint: (a) strain energy and material volumes; (b) global interfacial stress measure.

Table 3 Comparison of the strain energy values, global interfacial stress measure and maximum interfacial stresses between the optimized simply supported beam designs with and without the interfacial stress constraint.

f gσ Max. σint

With interfacial stress constraint

Without interfacial stress constraint

1.42 3.47 2.14

1.39 – 5.00

The optimized design, and the design optimized solely for minimum strain energy without the interfacial stress constraint, are both shown in Fig. 25. In the former, the weaker material acts as reinforcements to the main arc member, where the material interfaces mainly undergo shear stress. In the optimized design without the interfacial stress constraints, however, the weaker material forms two independent bars to directly transmit force from the loading point to the main arc, and the material interfaces obviously undergo tensile stress. The values of the strain energy, the global interfacial stress measure and the maximum interfacial stresses for both designs are given in Table 3. As expected, the design obtained with the proposed interfacial stress

24

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

Fig. 25. The optimized designs of the simply supported beam (a) with interfacial stress constraint and (b) without interfacial stress constraint.

Fig. 26. Equivalent interfacial stresses of the optimized simply supported beams (a) with interfacial stress constraint and (b) without interfacial stress constraint.

constraint yields a slightly larger strain energy (1.42 vs. 1.39), but a significantly smaller maximum interfacial stress (2.14 vs. 5.00) (the equivalent interfacial stress is depicted in Fig. 26), which again shows the effects of the proposed global interfacial stress constraint. 6. Conclusions In the present paper, we introduce a global material interfacial stress constraint into multi-material structural topology optimization for compliance minimization. A stress criterion taking into account both the shear and normal

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

25

tensile stresses on the material interface is adopted. To handle the local nature of the interfacial stress constraints, we aggregate the effective interfacial stresses with a p-norm type expression, which contains a material interface integral. Then we approximate this integral with a narrow-band domain integral around the material interface for ease of analytical derivation. The multi-material topology and material interface location are described with the level set method, and body-fitted meshes are employed in the interfacial stress field analysis. The velocity field level set method is adopted to tackle the multiple constraints. The sensitivities are derived with the adjoint method, and the material interface transmission conditions are observed to provide exact sensitivity. The numerical results show that the proposed method is able to control the effective material interfacial stress. Comparisons are made between the optimized designs and those obtained without consideration of interfacial stresses, showing the significance of including the interfacial strength constraints into multi-material topology optimization problems. It is also demonstrated that the velocity field level set method is effective in handling multiple constraints in an implicit level set framework. The presented numerical examples show that interfacial strength constrained multi-material topology optimization usually requires more iterations to achieve a converged design, which may significantly increase the computational cost. In theory, the proposed method can also be used for 3D problems and problems involving more than two materials. However, in such cases, the generation of body-fitted meshes for multi-material design is not an easy task. The difficulty of meshing can be circumvented by resorting to the extended finite element method, but the computation of stresses of the elements cut by multiple level set functions needs to be carefully handled. It is interesting to note that the weaker material tends to be attached to the main load-bearing components as a reinforcing phase in the optimized bi-material designs obtained in the considered numerical examples. This suggests that the weaker material should not be used as independent tensile force-bearing structural members to avoid interfacial strength failures under circumstances of relatively stringent interfacial strength constraints. For the same reason, the material interfaces should preferably be located at the compressive stresses zones. Also, the optimized solutions imply that forming reasonably larger interface lengths between the weaker and stronger materials is beneficial for improving the interfacial strength. These findings are in accordance with the engineering intuition, and may provide useful guidance to the design of bi-material composite structures in practice. Acknowledgments The authors acknowledge the support of the National Science Foundation of China (11902064, U1608256, 11872140) and the National Science Fund for Distinguished Young Scholars of China (11425207). The authors also thank the anonymous reviewers for their insightful comments and suggestions for improving the manuscript. References [1] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (1988) 197–224. [2] M.P. Bendsøe, Optimal shape design as a material distribution problem, Struct. Optim. 1 (1989) 193–202. [3] G.I. Rozvany, M. Zhou, T. Birker, Generalized shape optimization without homogenization, Struct. Optim. 4 (1992) 250–252. [4] J.A. Sethian, A. Wiegmann, Structural boundary design via level set and immersed interface methods, J. Comput. Phys. 163 (2000) 489–528. [5] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Engrg. 192 (2003) 227–246. [6] G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (2004) 363–393. [7] X. Huang, Y.-M. Xie, A further review of ESO type methods for topology optimization, Struct. Multidiscip. Optim. 41 (2010) 671–683. [8] K. Nabaki, J. Shen, X. Huang, Stress minimization of structures based on bidirectional evolutionary procedure, J. Struct. Eng. 145 (2018) 04018256. [9] P. Liu, Y. Luo, Z. Kang, Multi-material topology optimization considering interface behavior via XFEM and level set method, Comput. Methods Appl. Mech. Engrg. 308 (2016) 113–133. [10] P. Liu, Z. Kang, Integrated topology optimization of multi-component structures considering connecting interface behavior, Comput. Methods Appl. Mech. Engrg. 341 (2018) 851–887. [11] R. Behrou, M. Lawry, K. Maute, Level set topology optimization of structural problems with interface cohesion, Internat. J. Numer. Methods Engrg. 112 (2017) 990–1016. [12] Z. Kang, C. Wu, Y. Luo, M. Li, Robust topology optimization of multi-material structures considering uncertain graded interface, Compos. Struct. 208 (2019) 395–406.

26

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

[13] Z. Luo, M.Y. Wang, S. Wang, P. Wei, A level set-based parameterization method for structural shape and topology optimization, Internat. J. Numer. Methods Engrg. 76 (2008) 1–26. [14] G. Pingen, M. Waidmann, A. Evgrafov, K. Maute, A parametric level-set approach for topology optimization of flow domains, Struct. Multidiscip. Optim. 41 (2010) 117–131. [15] M.Y. Wang, X. Wang, Color level sets: a multi-phase method for structural topology optimization with multiple materials, Comput. Methods Appl. Mech. Engrg. 193 (2004) 469–496. [16] Y. Wang, Z. Luo, Z. Kang, N. Zhang, A multi-material level set-based topology and shape optimization method, Comput. Methods Appl. Mech. Engrg. 283 (2015) 1570–1586. [17] Y. Wang, Z. Kang, A velocity field level set method for shape and topology optimization, Internat. J. Numer. Methods Engrg. 115 (2018) 1315–1336. [18] J. Liu, Y. Ma, A new multi-material level set topology optimization method with the length scale control capability, Comput. Methods Appl. Mech. Engrg. 329 (2018) 444–463. [19] K.W. Schlichting, N. Padture, E. Jordan, M. Gell, Failure modes in plasma-sprayed thermal barrier coatings, Mater. Sci. Eng. A 342 (2003) 120–130. [20] D. Da, J. Yvonnet, L. Xia, G. Li, Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage, Internat. J. Numer. Methods Engrg. 115 (2018) 604–626. [21] K. Park, G.H. Paulino, Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces, Appl. Mech. Rev. 64 (2011) 060802. [22] C. Bernardi, O. Pironneau, Sensitivity of Darcy’s law to discontinuities, Chin. Ann. Math. 24 (2003) 205–214. [23] F. Hettlich, W. Rundell, The determination of a discontinuity in a conductivity from a single boundary measurement, Inverse Problems 14 (1998) 67. [24] O. Pantz, Sensitivity of the heat equation to jumps of conductivity, C. R. Math. 341 (2005) 333–337. [25] G. Allaire, F. Jouve, N. Van Goethem, Damage and fracture evolution in brittle materials by shape optimization methods, J. Comput. Phys. 230 (2011) 5010–5044. [26] G. Allaire, C. Dapogny, G. Delgado, G. Michailidis, Multi-phase structural optimization via a level set method, ESAIM Control Optim. Calc. Var. 20 (2014) 576–611. [27] G. Allaire, C. Dapogny, P. Frey, Shape optimization with a level set based mesh evolution method, Comput. Methods Appl. Mech. Engrg. 282 (2014) 22–53. [28] P. Duysinx, M.P. Bendsøe, Topology optimization of continuum structures with local stress constraints, Internat. J. Numer. Methods Engrg. 43 (1998) 1453–1478. [29] M. Bruggi, P. Venini, A mixed FEM approach to stress-constrained topology optimization, Internat. J. Numer. Methods Engrg. 73 (2008) 1693–1714. [30] G. Allaire, F. Jouve, Minimum stress optimal design with the level set method, Eng. Anal. Bound. Elem. 32 (2008) 909–918. [31] Q. Xia, T. Shi, S. Liu, M.Y. Wang, A level set solution to the stress-based structural shape and topology optimization, Comput. Struct. 90 (2012) 55–64. [32] Z. Fan, L. Xia, W. Lai, Q. Xia, T. Shi, Evolutionary topology optimization of continuum structures with stress constraints, Struct. Multidiscip. Optim. 59 (2019) 647–658. [33] A. Takezawa, G.H. Yoon, S.H. Jeong, M. Kobashi, M. Kitamura, Structural topology optimization with strength and heat conduction constraints, Comput. Methods Appl. Mech. Engrg. 276 (2014) 341–361. [34] G. Sved, Z. Ginos, Structural optimization under multiple loading, Int. J. Mech. Sci. 10 (1968) 803–805. [35] G. Cheng, Z. Jiang, Study on topology optimization with stress constraints, Eng. Optim. 20 (1992) 129–148. [36] G. Cheng, X. Guo, ε-Relaxed approach in structural topology optimization, Struct. Optim. 13 (1997) 258–266. [37] M. Bruggi, On an alternative approach to stress constraints relaxation in topology optimization, Struct. Multidiscip. Optim. 36 (2008) 125–141. [38] M. Burger, R. Stainko, Phase-field relaxation of topology optimization with local stress constraints, SIAM J. Control Optim. 45 (2006) 1447–1466. [39] G. Kreisselmeier, R. Steinhauser, Systematic control design by optimizing a vector performance index, in: Computer Aided Design of Control Systems, Elsevier, 1980, pp. 113–117. [40] P. Duysinx, O. Sigmund, New developments in handling stress constraints in optimal material distribution, in: 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, 1998, p. 4906. [41] C. Le, J. Norato, T. Bruns, C. Ha, D. Tortorelli, Stress-based topology optimization for continua, Struct. Multidiscip. Optim. 41 (2010) 605–620. [42] Y. Luo, M.Y. Wang, Z. Kang, An enhanced aggregation method for topology optimization with local stress constraints, Comput. Methods Appl. Mech. Engrg. 254 (2013) 31–41. [43] C. Kiyono, S. Vatanabe, E. Silva, J. Reddy, A new multi-p-norm formulation approach for stress-based topology optimization design, Compos. Struct. 156 (2016) 10–19. [44] M. Bruggi, P. Duysinx, Topology optimization for minimum weight with compliance and stress constraints, Struct. Multidiscip. Optim. 46 (2012) 369–384. [45] Y. Luo, Z. Kang, Topology optimization of continuum structures with Drucker–Prager yield stress constraints, Comput. Struct. 90 (2012) 65–75. [46] Y. Luo, J. Xing, Y. Niu, M. Li, Z. Kang, Wrinkle-free design of thin membrane structures using stress-based topology optimization, J. Mech. Phys. Solids 102 (2017) 277–293.

P. Liu, L. Shi and Z. Kang / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112887

27

[47] J.C. Brewer, P.A. Lagace, Quadratic stress criterion for initiation of delamination, J. Compos. Mater. 22 (1988) 1141–1155. [48] L.A. Vese, T.F. Chan, A multiphase level set framework for image segmentation using the mumford and shah model, Int. J. Comput. Vis. 50 (2002) 271–293. [49] N.P. van Dijk, K. Maute, M. Langelaar, F. van Keulen, Level-set methods for structural topology optimization: a review, Struct. Multidiscip. Optim. 48 (2013) 437–472. [50] P.D. Dunning, H.A. Kim, Robust topology optimization: minimization of expected and variance of compliance, AIAA J. 51 (2013) 2656–2664. [51] C. Barbarosie, S. Lopes, A.-M. Toader, An algorithm for constrained optimization with applications to the design of mechanical structures, in: International Conference on Engineering Optimization, Springer, 2018, pp. 272–284. [52] F. Feppon, G. Allaire, C. Dapogny, Null space gradient flows for constrained optimization with applications to shape optimization, 2019. [53] K. Svanberg, The method of moving asymptotes—a new method for structural optimization, Internat. J. Numer. Methods Engrg. 24 (1987) 359–373. [54] D. Adalsteinsson, J.A. Sethian, The fast construction of extension velocities in level set methods, J. Comput. Phys. 148 (1999) 2–22.