# Multimaterial Topology Optimization of Variational Inequalities

Andrzej Myśliński

► To cite this version:

Andrzej Myśliński. Multimaterial Topology Optimization of Variational Inequalities. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.380-389, 10.1007/978-3-319-55795-3_36. hal-01626885

HAL Id: hal-01626885

https://inria.hal.science/hal-01626885

Submitted on 31 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multimaterial Topology Optimization of Variational Inequalities

Andrzej Myśliński

Systems Research Institute,
ul. Newelska 6, 01-447 Warsaw, Poland
{andrzej.myslinski}@ibspan.waw.pl

**Abstract.** 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 finding 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 first 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.

**Keywords:** topology optimization, unilateral contact problems, phase field regularization, operator splitting method

## 1 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 sev-
eral 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 elim-
inate the need of the material interpolation schemes provided that interphase
interfaces are actually tracked explicitly [1]. Among others, in [1] a multima-
terial 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 sys-
tems governed by the variational inequalities. The class of such systems includes
among others unilateral contact phenomenon [9] between the surfaces of the elas-
tic 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 func-
tion is a variable subject to optimization. The regularization of the objective
functional by the multiphase volume constrained Ginzburg-Landau energy func-
tional 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.

## 2 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₁(x), f₂(x)), x ∈ Ω. Moreover, the surface tractions p(x) =
(p₁(x), p₂(x)), x ∈ Γ, are applied to a portion Γ₁ of the boundary Γ. The body
is clamped along the portion Γ₀ of the boundary Γ and the contact conditions
are prescribed on the portion Γ₂. Parts Γ₀, Γ₁, Γ₂ of the boundary Γ satisfy:
Γᵢ ∩ Γⱼ = ∅, i ≠ j, i, j = 0, 1, 2, Γ = Γ̃₀ ∪ Γ̃₁ ∪ Γ̃₂.

The domain Ω is assumed to be occupied by s ≥ 2 distinct isotropic elas-
tic 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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Elastic body occupying domain $\Omega$ in unilateral contact with the foundation.</figcaption></figure>

value of Young modulus [2]. The materials distribution is described by a phase field vector $\rho = \{\rho_m\}_{m=1}^s$ where the local fraction field $\rho_m = \rho_m(x) : \Omega \to \mathbb{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 $\rho_m$

$$ \alpha_m \le \rho_m \le \beta_m, \text{ for } m = 1, ..., s, \text{ and } \sum_{m=1}^{s} \rho_m = 1, \quad (1) $$

where the constants $0 \le \alpha_m \le \beta_m \le 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

$$ \int_{\Omega} \rho_m(x) dx = w_m |\Omega|, \quad 0 \le w_m \le 1, \text{ for } m = 1, ..., s, \text{ and } \sum_{m=1}^{s} w_m = 1. \quad (2) $$

The parameters $w_m$ are user defined and $|\Omega|$ denotes the volume of the domain $\Omega$. From the equality (1) it results that $\rho_s = 1 - \sum_{m=1}^{s-1} \rho_m$ and the fraction $\rho_s$ may be removed from the set of the design functions. Therefore from now on the unknown phase field vector $\rho$ is redefined as $\rho = \{\rho_m\}_{m=1}^{s-1}$. Due to the simplicity and robustness the SIMP material interpolation model [3, 4] is used. Following this model the elastic tensor $A(\rho) = \{a_{ijkl}(\rho)\}_{i,j,k,l=1}^2$ of the material body is assumed to be a function depending on the fraction function $\rho$:

$$ A(\rho) = \sum_{m=1}^{s} g(\rho_m) A_m = \sum_{m=1}^{s-1} g(\rho_m) A_m + g(1 - \sum_{m=1}^{s-1} \rho_m) A_s, \quad (3) $$

with $g(\rho_m) = \rho_m^3$. The constant stiffness tensor $A_m = \{\tilde{a}_{ijkl}^m\}_{i,j,k,l=1}^2$ 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 $\tilde{a}_{ijkl}^m(x)$, $i, j, k, l = 1, 2, m = 1, ..., s$, of the elasticity tensors $\mathcal{A}$ and $\mathcal{A}_m$, respectively, satisfy [9] usual symmetry, boundedness and ellipticity conditions. Denote by $u = (u_1, u_2)$, $u = u(x)$, $x \in \Omega$, the displacement of the body and by $\sigma(x) = \{\sigma_{ij}(u(x))\}$, $i, j = 1, 2$, the stress field in the body. Consider elastic bodies obeying Hooke's law, i.e., for $x \in \Omega$ and $i, j, k, l = 1, 2$,

$$ \sigma_{ij}(u(x)) = a_{ijkl}(\rho)e_{kl}(u(x)), \quad e_{kl}(u(x)) \stackrel{\text{def}}{=} \frac{1}{2}(u_{k,l}(x) + u_{l,k}(x)), \quad (4) $$

where $u_{k,l}(x) = \frac{\partial u_k(x)}{\partial x_l}$. We use here and throughout the paper the summation convention over repeated indices [9]. The stress field $\sigma$ satisfies the system of equations in the domain $\Omega$ [9]

$$ -\sigma_{ij}(x),_j = f_i(x) \quad \sigma_{ij}(x),_j = \frac{\partial \sigma_{ij}(x)}{\partial x_j}, \quad x \in \Omega, i, j = 1, 2. \quad (5) $$

The following boundary conditions are imposed on the boundary $\partial\Omega$

$$ u_i(x) = 0 \quad \text{on } \Gamma_0, \quad \sigma_{ij}(x)n_j = p_i \quad \text{on } \Gamma_1, \quad i, j = 1, 2, \quad (6) $$

$$ (u_N + v) \le 0, \quad \sigma_N \le 0, \quad (u_N + v)\sigma_N = 0 \quad \text{on } \Gamma_2, \quad (7) $$

$$ |\sigma_T| \le 1, \quad u_T\sigma_T + |u_T| = 0 \quad \text{on } \Gamma_2, \quad (8) $$

where $n = (n_1, n_2)$ is the unit outward versor to the boundary $\Gamma$. Here $u_N = u_i n_i$ and $\sigma_N = \sigma_{ij} n_i n_j$, $i, j = 1, 2$, represent [9] the normal components of displacement $u$ and stress $\sigma$, respectively. The tangential components of displacement $u$ and stress $\sigma$ are given [9] by $(u_T)_i = u_i - u_N n_i$ and $(\sigma_T)_i = \sigma_{ij} n_j - \sigma_N n_i$, $i, j = 1, 2$, respectively. $|u_T|$ denotes the Euclidean norm in $\mathbb{R}^2$ of the tangent vector $u_T$. 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 $V_{sp} = \{z \in H^1(\Omega; \mathbb{R}^2) : z_i = 0 \text{ on } \Gamma_0, i = 1, 2\}$ and $K = \{z \in V_{sp} : z_N + v \le 0 \text{ on } \Gamma_2\}$ the space and the set of kinematically admissible displacements and by $\Lambda = \{\zeta \in L^2(\Gamma_2; \mathbb{R}^2) : |\zeta| \le 1\}$ the set of tangential tractions on $\Gamma_2$. Variational formulation of problem (5)-(8) has the form: for a given $(f, p, \rho) \in L^2(\Omega; \mathbb{R}^2) \times L^2(\Gamma_2; \mathbb{R}^2) \times L^\infty(\Omega; \mathbb{R}^{s-1}) \cap H^1(\Omega; \mathbb{R}^{s-1})$ find a pair $(u, \lambda) \in K \times \Lambda$ satisfying for $i, j, k, l = 1, 2$

$$ \int_{\Omega} a_{ijkl}(\rho) e_{ij}(u) e_{kl}(\varphi_i - u_i) dx - \int_{\Omega} f_i(\varphi_i - u_i) dx - \int_{\Gamma_1} p_i(\varphi_i - u_i) ds + \int_{\Gamma_2} \lambda_i(\varphi_{Ti} - u_{Ti}) ds \ge 0 \quad \forall \varphi_i \in K, \quad (9) $$

$$ \int_{\Gamma_2} (\zeta_i - \lambda_i) u_{Ti} ds \le 0 \quad \forall \zeta_i \in \Lambda. \quad (10) $$

The function $\lambda$ 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 $\Gamma_2$, i.e., $\lambda = -\sigma_{T|\Gamma_2}$ and belongs to space $H^{-1/2}(\Gamma_2; \mathbb{R}^2)$. Here following [9] function $\lambda$ is assumed to be more regular, i.e., $\lambda \in L^2(\Gamma_2; \mathbb{R}^2)$. Recall
from [9, 14] under assumptions imposed on the elasticity tensors by standard arguments for a given $(f, p, \rho) \in L^2(\Omega; R^2) \times L^2(\Gamma_2; R^2) \times L^\infty(\Omega; R^{s-1})$ there exists a unique solution $(u, \lambda) \in K \times \Lambda$ to the system (9)-(10).

### 3 Phase field based topology optimization problem

Before formulating a structural optimization problem for the system (9)-(10) let us introduce a set $U_{ad}^\rho$ of the admissible fraction functions:

$$ U_{ad}^{\rho} = \{ \rho \in L^{\infty}(\Omega; R^{s-1}) \cap H^1(\Omega; R^{s-1}) : 1 - \beta_s \le \sum_{m=1}^{s-1} \rho_m \le 1 - \alpha_s, \\ \alpha_m \le \rho_m \le \beta_m, \int_{\Omega} \rho_m dx = w_m |\Omega| \text{ for } m = 1, \dots, s-1 \}. \quad (11) $$

The set $U_{ad}^\rho$ 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_\eta : H^1(\Omega) \to R$ approximating the normal contact stress on the contact boundary $\Gamma_2$

$$ J_{\eta}(u(\rho)) = \int_{\Gamma_2} \sigma_N(u(\rho))\eta_N(x)ds, \quad (12) $$

depending on a given auxiliary bounded function $\eta(x) \in M^{st}$. The set $M^{st} = \{\eta = (\eta_1, \eta_2) \in H^1(\Omega; R^2) : \eta_i \le 0 \text{ on } \Omega, i = 1, 2, \|\eta\|_{H^1(\Omega; R^2)} \le 1\}$. Functions $\sigma_N$ and $\eta_N$ are the normal components of the stress field $\sigma$ and the function $\eta$, respectively. The optimization problem consisting in finding such $\rho \in U_{ad}^\rho$ to minimize the functional $J_\eta(u(\rho))$ 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(\rho) : U_{ad}^\rho \to R$ rather than the standard perimeter term [2, 7, 14]

$$ J(\rho, u(\rho)) = J_{\eta}(u(\rho)) + E(\rho). \quad (13) $$

The Ginzburg-Landau free energy functional $E(\rho)$ is expressed as [7, 15]

$$ E(\rho) = \sum_{m=1}^{s-1} \int_{\Omega} \psi(\rho_m) d\Omega, \quad \psi(\rho_m) = \frac{\gamma\epsilon}{2} |\nabla \rho_m|^2 + \frac{\gamma}{\epsilon} \psi_B(\rho_m), \quad (14) $$

where $\epsilon > 0$ is a real constant governing the width of the intrefaces, $\gamma > 0$ is a real parameter related to the interfacial energy density. Moreover $\nabla \rho_m \cdot n = 0$ on $\Gamma$ for each $m$. The function $\psi_B(\rho_m) = \rho_m^2(1 - \rho_m)^2$ 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 $\rho^* \in U_{ad}^\rho$ such that

$$J(\rho^*, u^*) = \min_{\rho \in U_{ad}^\rho} J(\rho, u(\rho)), \qquad (15)$$

where $u^* = u(\rho^*)$ denotes a solution to the state system (9)-(10) depending on $\rho^* \in L^\infty(\Omega; R^{s-1}) \cap H^1(\Omega; R^{s-1})$ and the set $U_{ad}^\rho$ is given by (11). The existence of an optimal solution $\rho^* \in U_{ad}^\rho$ to the problem (15) follows by classical arguments (see [5, 6, 14]).

## 4 Necessary optimality condition

Let us apply the Lagrangian approach to compute the derivative of the cost functional (13) with respect to the function $\rho$. The Lagrangian $L(\rho) = L(\rho, u, \lambda, p^a, q^a) : L^\infty(\Omega; R^{s-1}) \cap H^1(\Omega; R^{s-1}) \cap U_{ad}^\rho \times H^1(\Omega; R^2) \times L^2(\Gamma_2; R^2) \times K_1 \times \Lambda_1$ associated to the problem (15) is expressed for $i, j, k, l = 1, 2$ as

$$L(\rho, u, \lambda, p^a, q^a) = J_{\eta}(u(\rho)) + E(\rho) + \int_{\Omega} a_{ijkl}(\rho)e_{ij}(u)e_{kl}(p^a)dx - \\ \int_{\Omega} f_i p_i^a dx - \int_{\Gamma_1} p_i p_i^a ds + \int_{\Gamma_2} \lambda_i p_{Ti}^a ds + \int_{\Gamma_2} q_i^a u_{Ti} ds. \quad (16)$$

Using the same approach as in proof of [14, theorem 4.35, p. 219] an adjoint state $(p^a, q^a) \in K_1 \times \Lambda_1$ for $i, j, k, l = 1, 2$ is defined as the solution to the system

$$\int_{\Omega} a_{ijkl}(\rho)e_{ij}(\eta + p^a)e_{kl}(\varphi)dx + \int_{\Gamma_2} q_i^a \varphi_{Ti} ds = 0 \quad \forall \varphi_i \in K_1, \quad (17)$$

$$\int_{\Gamma_2} \zeta_i(p_{Ti}^a + \eta_{Ti})ds = 0 \quad \forall \zeta_i \in \Lambda_1. \quad (18)$$

The sets $K_1$ and $\Lambda_1$ are given by $K_1 = \{\xi \in V_{sp} : \xi_N = 0 \text{ on } A^{st}\}$ as well as $\Lambda_1 = \{\zeta \in \Lambda : \zeta(x) = 0 \text{ on } B_1 \cup B_2 \cup B_1^+ \cup B_2^+\}$, while the coincidence set $A^{st} = \{x \in \Gamma_2 : u_N + v = 0\}$. Moreover the other sets are determined as $B_1 = \{x \in \Gamma_2 : \lambda(x) = -1\}$, $B_2 = \{x \in \Gamma_2 : \lambda(x) = +1\}$, $\tilde{B}_i = \{x \in B_i : u_N(x) + v = 0\}$, $i = 1, 2$, $B_i^+ = B_i \setminus \tilde{B}_i$, $i = 1, 2$. From [12], [14, theorems 4.16, 4.27] it follows the mapping $\rho \to u(\rho)$ is Gâteaux differentiable. Using (9)-(10) and (16)-(18) the derivative of the Lagrangian $L$ with respect to $\rho$ is determined for all $\zeta \in H^1(\Omega; R^{s-1})$ and $i, j, k, l = 1, 2$ as

$$\int_{\Omega} \frac{\partial J}{\partial \rho}(\rho, u)\zeta dx = \int_{\Omega} \frac{\partial L}{\partial \rho}(\rho, u, \lambda, p^a, q^a)\zeta dx = \sum_{m=1}^{s-1} \int_{\Omega} [\gamma \epsilon \nabla \rho_m \cdot \nabla \zeta_m + \\ \frac{\gamma}{\epsilon} \psi'_B(\rho_m)\zeta_m] dx + \int_{\Omega} [a'_{ijkl}(\rho_m)e_{ij}(u)e_{kl}(p^a + \eta) - f_i(p_i^a + \eta_i)]\zeta_m dx. \quad (19)$$

From (1) and (3) it results the derivatives of the function $\psi_B(\rho_m)$ and the tensor element $a_{ijkl}(\rho_m)$ with respect to $\rho_m$ are equal to $\psi'_B(\rho_m) = 4\rho_m^3 - 6\rho_m^2 + 2\rho_m$ and $a'_{ijkl}(\rho) = 3\rho_m^2 \tilde{a}_{ijkl}^m - 3\rho_s^2 \tilde{a}_{ijkl}^s$, respectively. Using (19) the necessary optimality condition to the optimization problem (15) takes the form [2], [13, Lemma 2.21]:
**Theorem 1.** Let $U_{ad}^\rho$ be a nonempty closed convex subset of $H^1(\Omega; R^{s-1})$ and $\rho^* \in U_{ad}^\rho$ be an optimal solution to the structural optimization problem (15). Then

$$ \int_{\Omega} \frac{\partial J}{\partial \rho}(\rho^*, u^*)(\rho - \rho^*)dx \geq 0 \quad \forall \rho \in U_{ad}^\rho. \qquad (20) $$

The functions $(u^*, \lambda^*) \in K \times \Lambda$ and $(p^{a*}, q^{a*}) \in K_1 \times \Lambda_1$ in the derivative formula (19) denote the solutions to the systems (9)-(10) and (17)-(18) for $\rho = \rho^*$.

Using the orthogonal projection operator $P_{U_{ad}^\rho}: L^2(\Omega; R^{s-1}) \to U_{ad}^\rho$ from $L^2(\Omega; R^{s-1})$ on the set $U_{ad}^\rho$ condition (20) can be written [15] in the form: if $\rho^* \in U_{ad}^\rho$ is an optimal solution to the structural optimization problem (15), then for $\mu \in R$ and $\mu > 0$

$$ P_{U_{ad}^\rho}[\rho^* - \mu \frac{\partial J(\rho^*, u^*)}{\partial \rho}] - \rho^* = 0. \qquad (21) $$

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 $\rho$ depends not only on $x \in \Omega$ but also on time variable $t \in [0, T]$, $T > 0$ is a given constant, i.e., $\rho = \rho(x, t) = \{\rho_m(x, t)\}_{m=1}^{s-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 $\rho \in U_{ad}^\rho$ satisfying the initial boundary value problem:

$$ \frac{\partial \rho}{\partial t} = -P_{U_{ad}^\rho}[\rho - \mu \frac{\partial J(\rho, u)}{\partial \rho}] + \rho \quad \text{in } \Omega, \forall t \in [0, T), \qquad (22) $$

$$ \nabla \rho \cdot n = 0 \qquad \text{on } \partial\Omega, \forall t \in [0, T), \qquad (23) $$

$$ \rho(0, x) = \rho_0(x) \qquad \text{in } \Omega, t = 0, \qquad (24) $$

with $\rho_0(x) = \{\rho_{0m}(x)\}_{m=1}^{s-1}$ denoting a given $H^1(\Omega; R^{s-1})$ regular function. For $\rho_0 \in H^1(\Omega; R^{s-1})$ the system (22)-(24) possesses a solution $\rho \in L^\infty(0, T; H^1(\Omega)) \cap H^1(0, T; L^2(\Omega; 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 $\frac{\partial \rho}{\partial t} = 0$ the right hand side of the equation (22) vanishes and $\rho(x, t) = \rho^*(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(\rho, u) = J_1(\rho, u) + J_2(\rho)$ given by $J_1(\rho, u) = J_\eta(u(\rho)) + \sum_{m=1}^{s-1} \int_{\Omega} \frac{\gamma}{\epsilon} \psi_B(\rho_m) d\Omega$ and $J_2(\rho) = \sum_{m=1}^{s-1} \int_{\Omega} \frac{\gamma}{2} |\nabla \rho_m|^2 d\Omega$. The derivatives of these functionals result from formula (19). Assume the time interval $[0, T]$ is divided into N subintervals with stepsize $\Delta t = t_{k+1} - t_k$, $k = 1, ..., N$ and $\rho_k = \rho(t_k)$ is known. The design variable $\rho_{k+1}$ at the next time step $t_{k+1}$ is calculated in two substeps. First the trial value $\tilde{\rho}$ is calculated from the gradient flow equation (22) for the functional
$J_1$ only. Next this solution is updated to ensure its $H^1(\Omega)$ regularity [5] by solving the gradient flow equation (22) for the functional $J_2$ only with the boundary condition (23), i.e.,

$$ \frac{\partial \tilde{\rho}}{dt} = -P_{U_{ad}}^{\rho} [\tilde{\rho} - \frac{\partial J_1(\tilde{\rho}, u)}{\partial \rho}] + \tilde{\rho}, \quad \tilde{\rho}(t_k) = \rho_k, \quad t_k < t \le t_{k+1}. \quad (25) $$

$$ \frac{\partial \rho}{dt} = -\frac{\partial J_2(\rho)}{\partial \rho}, \quad \rho(t_k) = \tilde{\rho}_{k+1}, \quad t_k < t \le t_{k+1}. \quad (26) $$

## 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 $\rho_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

$$ \Omega = \{(x_1, x_2) \in \mathbb{R}^2 : 0 \le x_1 \le 8 \land 0 < v(x_1) \le x_2 \le 4\}, \quad (27) $$

is considered. The boundary $\Gamma$ of the domain $\Omega$ is divided into three disjoint pieces $\Gamma_0 = \{(x_1, x_2) \in \mathbb{R}^2 : x_1 = 0, 8 \land 0 < v(x_1) \le x_2 \le 4\}$, $\Gamma_1 = \{(x_1, x_2) \in \mathbb{R}^2 : 0 \le x_1 \le 8 \land x_2 = 4\}$, $\Gamma_2 = \{(x_1, x_2) \in \mathbb{R}^2 : 0 \le x_1 \le 8 \land v(x_1) = x_2\}$. The domain $\Omega$ and the boundary $\Gamma_2$ depend on the function $v(x_1) = 0.125 \cdot (x_1 - 4)^2$. Domain $\Omega$ is filled with $s = 3$ elastic materials. The Poisson's ratio of each material is $\nu = 0.3$. The Young's moduli of materials are: $E_1 = 6 \cdot E_0$, $E_2 = 3 \cdot E_0$ and $E_3 = E_0$, $E_0 = 2.1 \cdot 10^{11}$ Pa. The parameters $w_1, w_2, w_3$ are equal to .25, .5 and .25 respectively. As an initial design $\rho_0$ a feasible design with the uniform material distribution has been taken. The body is loaded by the boundary traction $p_1 = 0$, $p_2 = -5.6 \cdot 10^6$ N along the boundary $\Gamma_1$, the body forces $f_i = 0, i = 1, 2$. The auxiliary function $\eta$ is selected as a piecewise linear on $\Omega$ and is approximated by a piecewise linear function. The domain $\Omega$ is divided into 80 x 40 grid. The parameters $\varepsilon$ and $\gamma$ 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
<figure><img src="image_3.png" /><figcaption>Fig. 2. Optimal material distribution in domain $\Omega^*$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. Convergence of the cost functional value.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Initial and optimal normal contact stress.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Convergence of the cost functional gradient.</figcaption></figure>

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).

# 6 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.
References

1. Allaire, G., Dapogny, C., Delgado, G., Michailidis, G.: Multi-phase structural optimization via a level set method. ESAIM - Control Optimisation and Calculus of Variations 20, 576-611 (2014)

2. Allaire, G.: Shape optimization by the homogenization method. Springer, New York, (2001)

3. Bendsoe, M.P., Sigmund, O.: Topology Optimization: Theory, Methods, and Applications, Springer, Berlin (2004)

4. Blank, L., Butz, M., Garcke, H., Sarbu, L., Styles, V.: Allen-Cahn and Cahn-Hilliard variational inequalities solved with optimization techniques. In: G. Leugering, S. Engell, A. Griewank, M. Hinze, R. Rannacher, V. Schulz, M. Ulbrich, S. Ulbrich (eds.), Constrained Optimization and Optimal Control for Partial Differential Equations, International Series of Numerical Mathematics, vol. 160, pp. 21-35, Birkhäuser, Basel (2012)

5. Blank, L., Farshbaf-Shaker, M.H., M., Garcke, H., Rupprecht, C., Styles, V.: Multi-material Phase Field Approach to Structural Topology Optimization. In: Leugering, G., Benner, P., Engell, S., Griewank, A., Harbrecht, H., Hinze, M., Rannacher, R., Ulbrich, S.: Trends in PDE Constrained Optimization. International Series of Numerical Mathematics 165, pp. 231-246, Birkhäuser, Basel (2014)

6. Burger, M., Stainko, R.: Phase-field relaxation of topology optimization with local stress constraints. SIAM J. Control. Optim. 45, 1447-1466 (2006)

7. Dede, L., Boroden, M.J., Hughes, T.J.R.: Isogeometric analysis for topology optimization with a phase field model. Archives of Computational Methods in Engineering 19(3), 427-465 (2012)

8. van Dijk, N.P., Maute, K., Langlaar, M., van Keulen, F.: Level-set methods for structural topology optimization: a review. Structural and Multidisciplinary Optimization 48, 437-472 (2013)

9. Haslinger, J., Mäkinen, R.: Introduction to Shape Optimization. Theory, Approximation, and Computation. SIAM Publications, Philadelphia (2003)

10. Karali, G., Ricciardi, T.: On the convergence of a fourth order evolution equation to the Allen Cahn equation. Nonlinear Analysis 72, 4271-4281 (2010)

11. Myśliński, A.: Piecewise Constant Level Set Method for Topology Optimization of Unilateral Contact Problems. Advances in Engineering Software 80, 25-32 (2015)

12. Myśliński, A.: Level Set Method for Optimization of Contact Problems. Engineering Analysis with Boundary Elements 32, 986-994 (2008)

13. Tröltzsch, F.: Optimal control of partial differential equations: Theory, methods and applications, American Mathematical Society, Graduate Studies in Mathematics 112, Providence, Rhode Island (2010)

14. Sokołowski, J., Zolesio, J.P.: Introduction to Shape Optimization. Shape Sensitivity Analysis. Springer, Berlin (1992)

15. Tavakoli, R.: Multimaterial Topology Optimization by Volume Constrained Allen-Cahn System and Regularized Projected Steepest Descent Method. Comput. Meth. Appl. Mech. Eng. 276, 534-565 (2014)