Multimaterial Topology Optimization of Variational Inequalities

. The paper is concerned with the analysis and the numerical solution of the multimaterial topology optimization problems for bodies in unilateral contact. The contact phenomenon with Tresca friction is governed by the elliptic variational inequality. The structural optimization problem consists in ﬁnding such topology of the domain occupied by the body that the normal contact stress along the boundary of the body is minimized. The original cost functional is regularized using the multi-phase volume constrained Ginzburg-Landau energy functional. The ﬁrst order necessary optimality condition is formulated. The optimal topology is obtained as the steady state of the phase transition governed by the generalized Allen-Cahn equation. The optimization problem is solved numerically using the operator splitting approach combined with the projection gradient method. Numerical examples are provided and discussed.


Introduction
Multimaterial topology optimization aims to find the optimal distribution of several elastic materials in a given design domain to minimize a criterion describing the mechanical or the thermal properties of the structure or its cost under constraints imposed on the volume or the mass of the structure [1].In recent years multiple phases topology optimization problems have become subject of the growing interest [1,5,15].The use of multiple number of phases during design of engineering structures opens a new opportunities in the design of smart and advanced structures in material science and/or industry.In contrast to single material design the use of multiple number of materials extends the design space and may lead to better design solutions.
Analytical and numerical aspects of the multimaterial structural optimization are subject of intensive research (see references in [1,15]).Many methods including the homogenization method [2], the Solid Isotropic Material Penalization (SIMP) method [3] or different methods based on the level set approach [8,11,12], successful in single material optimization, have been extended to deal with the multimaterial optimization.The extension of these methods faces several challenges.A crucial issue in the solution of the multimaterial optimization problems is the lack of physically based parametrization of the phases mixture [1,15].Although in the literature are proposed different material interpolation schemes, in general, they may influence the optimization path in terms of the computational efficiency and the final design.The level set methods can eliminate the need of the material interpolation schemes provided that interphase interfaces are actually tracked explicitly [1].Among others, in [1] a multimaterial topology optimization problem for the elliptic equation has been solved using the level set method.The elasticity tensor has been smeared out using the signed distance function.In [15] similar optimization problem has been solved numerically using a generalized Allen-Cahn equation.
The paper is concerned with the structural topology optimization of systems governed by the variational inequalities.The class of such systems includes among others unilateral contact phenomenon [9] between the surfaces of the elastic bodies.This optimization problem consists in finding such topology of the domain occupied by the body that the normal contact stress along the boundary of the body is minimized.In literature [11] this problem usually is considered as two-phase material optimization problem with voids treated as one of the materials.In the paper the domain occupied by the body is assumed to consist from several elastic materials rather than two materials.Material fraction function is a variable subject to optimization.The regularization of the objective functional by the multiphase volume constrained Ginzburg-Landau energy functional is used.The derivative formula of the cost functional with respect to the material fraction function is calculated and is employed to formulate a necessary optimality condition for the topology optimization problem.The cost functional derivative is also used to formulate a gradient flow equation for this functional in the form of the generalized Allen-Cahn equation governing the evolution of the material phases.Two step operator splitting approach [15] is used to solve this gradient flow equation.The optimal topology is obtained as a steady state solution to this equation.Finite difference and finite element methods are used as the approximation methods.Numerical examples are reported and discussed.

Problem Formulation
Consider deformations of an elastic body occupying two-dimensional bounded domain Ω with the smooth boundary Γ (see Fig. 1).The body is subject to body forces f (x) = (f 1 (x), f 2 (x)), x ∈ Ω.Moreover, the surface tractions p(x) = (p 1 (x), p 2 (x)), x ∈ Γ , are applied to a portion Γ 1 of the boundary Γ .The body is clamped along the portion Γ 0 of the boundary Γ and the contact conditions are prescribed on the portion Γ 2 .Parts Γ 0 , Γ 1 , Γ 2 of the boundary Γ satisfy: The domain Ω is assumed to be occupied by s ≥ 2 distinct isotropic elastic materials.Each material is characterized by Young modulus.The voids are considered as one of the phases, i.e., as a weak material characterized by low value of Young modulus [2].The materials distribution is described by a phase field vector ρ = {ρ m } s m=1 where the local fraction field ρ m = ρ m (x) : Ω → R, m = 1, ..., s, corresponds to the contributing phase.The phase field approach allows for a certain mixing between different materials.This mixing is restricted only to a small interfacial region.In order to ensure that the phase field vector describes the fractions the following pointwise bound constraints called in material science the Gibbs simplex [5,15] are imposed on every ρ m α m ≤ ρ m ≤ β m , for m = 1, ..., s, and where the constants 0 ≤ α m ≤ β m ≤ 1 are given and the summation operator is understood componentwise.The second condition in (1) ensures that no overlap and gap of fractions are allowed in the expected optimal domain.Moreover the total spatial amount of material fractions satisfies

, s, and
The parameters w m are user defined and | Ω | denotes the volume of the domain Ω.From the equality (1) it results that ρ s = 1 − s−1 m=1 ρ m and the fraction ρ s may be removed from the set of the design functions.Therefore from now on the unknown phase field vector ρ is redefined as ρ = {ρ m } s−1 m=1 .Due to the simplicity and robustness the SIMP material interpolation model [3,4] is used.Following this model the elastic tensor A(ρ) = {a ijkl (ρ)} 2 i,j,k,l=1 of the material body is assumed to be a function depending on the fraction function ρ: with g(ρ m ) = ρ 3 m .The constant stiffness tensor A m = {ã m ijkl } 2 i,j,k,l=1 characterizes the m-th elastic material of the body.For detailed discussion of the interpolation of the material elasticity tensor see [1,4,7].It is assumed, that elements a ijkl and ãm ijkl (x), i, j, k, l = 1, 2, m = 1, ..., s, of the elasticity tensors A and A m , respectively, satisfy [9] usual symmetry, boundedness and ellipticity conditions.Denote by u = (u 1 , u 2 ), u = u(x), x ∈ Ω, the displacement of the body and by σ(x) = {σ ij (u(x))}, i, j = 1, 2, the stress field in the body.Consider elastic bodies obeying Hooke's law, i.e., for x ∈ Ω and i, j, k, l = 1, 2, where u k,l (x) = ∂u k (x) ∂x l .We use here and throughout the paper the summation convention over repeated indices [9].The stress field σ satisfies the system of equations in the domain The following boundary conditions are imposed on the boundary ∂Ω where n = (n 1 , n 2 ) is the unit outward versor to the boundary Γ .Here [9] the normal components of displacement u and stress σ, respectively.The tangential components of displacement u and stress σ are given [9] by (u A gap between the bodies is described by a given function v. Let us formulate contact problem ( 5)-( 8) in the variational form.Denote by Γ2 The function λ is interpreted as a Lagrange multiplier corresponding to the term | u T | in equality constraint in (8) [9].This function is equal to tangent stress along the boundary Γ 2 , i.e., λ = −σ T |Γ 2 and belongs to space H −1/2 (Γ 2 ; R 2 ).Here following [9] function λ is assumed to be more regular, i.e., λ ∈ L 2 (Γ 2 ; R 2 ).Recall from [9,14] under assumptions imposed on the elasticity tensors by standard arguments for a given (f, p, ρ) ) there exists a unique solution (u, λ) ∈ K × Λ to the system ( 9)-( 10).

Phase field based topology optimization problem
Before formulating a structural optimization problem for the system ( 9)-( 10) let us introduce a set U ρ ad of the admissible fraction functions: The set U ρ ad is assumed to be nonempty.The aim of the structural optimization of the bodies in contact is to reduce the normal contact stress responsible for wear, vibrations of the contacting surfaces or generated noise.The structural optimization problem with normal contact stress functional is difficult to analyze it and to solve it numerically.Therefore following [11] we shall use the cost functional J η : H 1 (Ω) → R approximating the normal contact stress on the contact boundary Γ 2 depending on a given auxiliary bounded function η(x) ∈ M st .The set Functions σ N and η N are the normal components of the stress field σ and the function η, respectively.The optimization problem consisting in finding such ρ ∈ U ρ ad to minimize the functional J η (u(ρ)) in general has no solutions [2,6,7,9,14,15].In order to ensure the existence of optimal solutions let us regularize the cost functional ( 12) by adding to it a regularizing term E(ρ) : U ρ ad → R rather than the standard perimeter term [2, 7, 14] The Ginzburg-Landau free energy functional E(ρ) is expressed as [7,15] where > 0 is a real constant governing the width of the intrefaces, γ > 0 is a real parameter related to the interfacial energy density.Moreover is a double-well potential [7] which characterizes the concentration of the material phases [15].The structural optimization problem for the system ( 9)-( 10) takes the form: find ρ ∈ U ρ ad such that J(ρ , u ) = min where u = u(ρ ) denotes a solution to the state system ( 9)-( 10) depending on ρ ∈ L ∞ (Ω; R s−1 ) ∩ H 1 (Ω; R s−1 ) and the set U ρ ad is given by (11).The existence of an optimal solution ρ ∈ U ρ ad to the problem (15) follows by classical arguments (see [5,6,14]).

Necessary optimality condition
Let us apply the Lagrangian approach to compute the derivative of the cost functional (13) with respect to the function ρ.The Lagrangian L(ρ) = L(ρ, u, λ, p a , q a ) : 15) is expressed for i, j, k, l = 1, 2 as Using the same approach as in proof of [14, theorem 4.35, p. 219] an adjoint state (p a , q a ) ∈ K 1 × Λ 1 for i, j, k, l = 1, 2 is defined as the solution to the system The sets K 1 and Λ 1 are given by K 1 = {ξ ∈ V sp : ξ N = 0 on A st } as well as Moreover the other sets are determined as From [12], [14, theorems 4.16, 4.27] it follows the mapping ρ → u(ρ) is Gâteaux differentiable.Using ( 9)-( 10) and ( 16)-( 18) the derivative of the Lagrangian L with respect to ρ is determined for all ζ ∈ H 1 (Ω; R s−1 ) and i, j, k, l = 1, 2 as From ( 1) and ( 3 Theorem 1.Let U ρ ad be a nonempty closed convex subset of H 1 (Ω; R s−1 ) and ρ ∈ U ρ ad be an optimal solution to the structural optimization problem (15).Then The functions (u , λ ) ∈ K ×Λ and (p a , q a ) ∈ K 1 ×Λ 1 in the derivative formula (19) denote the solutions to the systems ( 9)-( 10) and ( 17)-( 18) for ρ = ρ .
Using the orthogonal projection operator P U ρ ad : L 2 (Ω; R s−1 ) → U ρ ad from L 2 (Ω; R s−1 ) on the set U ρ ad condition (20) can be written [15] in the form: if ρ ∈ U ρ ad is an optimal solution to the structural optimization problem (15), then for µ ∈ R and µ > 0 Recall [7] the structural optimization problem (15) can be considered as a phase transition setting problem consisting in such evolution of the phases to minimize the cost functional ( 13) with respect to the initial configuration.In order to describe the evolution of phases in time let us assume that the phase field vector ρ depends not only on x ∈ Ω but also on time variable t ∈ [0, T ], T > 0 is a given constant, i.e., ρ = ρ(x, t) = {ρ m (x, t)} s−1 m=1 .The variable t may be interpreted as an artificial time or iteration number in the computational algorithm [7,15].Using the right hand side of (21) let us formulate the constrained gradient flow equation of Allen-Cahn type [5,7,15] for the cost functional ( 13): find function ρ ∈ U ρ ad satisfying the initial boundary value problem: with ρ 0 (x) = {ρ 0m (x)} s−1 m=1 denoting a given H 1 (Ω; R s−1 ) regular function.For ρ 0 ∈ H 1 (Ω; R s−1 ) the system (22)-(24) possesses a solution ρ ∈ L ∞ (0, T ; H 1 (Ω)) ∩H 1 (0, T ; L 2 (Ω; R s−1 )) (see [10]).The stationary solutions of ( 22)-(24) fulfill the first order necessary optimality conditions (20) or (21) for the problem (15) [5,7].For ∂ρ ∂t = 0 the right hand side of the equation ( 22) vanishes and ρ(x, t) = ρ (x, t) is an optimal solution to the problem (15).For the sake of numerical calculations we reformulate the initial boundary value problem ( 22)-(24) using the operator splitting approach [5,15].Remark, the cost functional (13) may be represented as a sum of two functionals, i.e., J(ρ, u) = J 1 (ρ, u) + J 2 (ρ) given by J The derivatives of these functionals result from formula (19).Assume the time interval [0, T ] is divided into N subintervals with stepsize ∆t = t k+1 − t k , k = 1, ..., N and ρ k = ρ(t k ) is known.The design variable ρ k+1 at the next time step t k+1 is calculated in two substeps.First the trial value ρ is calculated from the gradient flow equation ( 22) for the functional J 1 only.Next this solution is updated to ensure its H 1 (Ω) regularity [5] by solving the gradient flow equation ( 22) for the functional J 2 only with the boundary condition (23), i.e., 5 Numerical results The topology optimization problem (15) has been discretized and solved numerically.Time derivatives are approximated by the forward finite difference.Piecewise constant and piecewise linear finite element method is used as discretization method in space variables.The derivative of the double well potential is linearized with respect to ρ m .Primal-dual active set method has been used to solve the state and adjoint systems ( 5)-( 8) and ( 17)-( 18).The initial boundary value problem ( 22)-( 24) has been solved in two steps according to scheme (25)-( 26).The algorithms are programmed in Matlab environment.As an example a body occupying 2D domain is considered.The boundary Γ of the domain Ω is divided into three disjoint pieces The auxiliary function η is selected as a piecewise linear on Ω and is approximated by a piecewise linear function.The domain Ω is divided into 80 × 40 grid.The parameters ε and γ are equal to the mesh size and to 0.5, respectively.The total number of iterations k max in the optimization algorithm has been set to 90.It is approximately equivalent to final time T = 125 s.Fig. 2 presents the optimal topology domain obtained by solving structural optimization problem (15) using the necessary optimality condition ( 22)-( 24).The areas with the weak phases appear in the central part of the body and near the fixed edges.The areas with the strong phases appear close to the contact zone and along the edges.The rest of the domain is covered with the intermediate phase.The obtained normal contact stress for the optimal topology is almost constant along the contact boundary and has been significantly reduced comparing to the initial one (see Fig. 3).The convergence of the cost functional    value and its gradient as a function of the number of iterations are shown on Fig. 4 and 5, respectively.The cost functional value decreases almost monotonically when the number of iterations increases.At the beginning this decrease is significant and finally the cost functional value is almost steady.The regularized functional value at first is increasing and next rapidly decreasing (see Fig. 4).Similarly, after a few initial iterations the gradient of the cost functional also almost monotonically decreases to reach the steady state (Fig. 5).

Conclusions
The topology optimization problem for elastic contact problem with Tresca friction has been solved numerically in the paper.Obtained numerical results indicate that the optimal topologies are qualitatively comparable to the results reported in other phase-field topology optimization methods.Since the optimization problem is non-convex it has possibly many local solutions dependent on initial estimate.Gradient flow method employed in H 1 space is more regular and efficient than standard Allen-Cahn approach.

Fig. 1 .
Fig. 1.Elastic body occupying domain Ω in unilateral contact with the foundation.