# Binary Level Set Method for Topology Optimization of Variational Inequalities

Andrzej Myśliński

► To cite this version:

Andrzej Myśliński. Binary Level Set Method for Topology Optimization of Variational Inequalities. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.199-209, 10.1007/978-3-662-45504-3_19. hal-01286414

HAL Id: hal-01286414

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

Submitted on 10 Mar 2016

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
Binary level set method for topology
optimization of variational inequalities

Andrzej Myśliński

Systems Research Institute
ul. Newelska 6, 01-447 Warsaw, Poland
myslinsk@ibspan.waw.pl

**Abstract.** The paper is concerned with the topology optimization of the elliptic variational inequalities using the level set approach. The standard level set method is based on the description of the domain boundary as an isocountour of a scalar function of a higher dimensionality. The evolution of this boundary is governed by Hamilton-Jacobi equation. In the paper a binary level set method is used to represent sub-domains rather than the standard method. The binary level set function takes at convergence value 1 in each sub domain of a whole design domain and -1 outside this sub domain. The sub domains interfaces are represented by discontinuities of these functions. Using a two-phase approximation and a binary level set approach the original structural optimization problem is reformulated as an equivalent constrained optimization problem in terms of this level set function. Necessary optimality condition is formulated. Numerical examples are provided and discussed.

**Keywords:** topology optimization, unilateral contact problems, binary level set method, Uzawa method

# 1 Introduction

Topology optimization problem for an elliptic second order variational inequality is considered in the paper. This inequality governs unilateral contact between an elastic body and a rigid foundation. The results concerning the existence and the uniqueness of solutions to this inequality are provided in [9]. The topology optimization problem for the elastic body in unilateral contact consists in finding such material distribution within the domain occupied by the body in contact and/or the shape of its boundary that the normal contact stress along the boundary of the body is minimized. The volume of the body is bounded.

Topology optimization of continuum structures is widely investigated in literature [1, 5, 7, 11]. Among others the homogenization method and its simplified version solid isotropic material with penalization method (SIMP) as well as evolutionary structural optimization method have been proposed to solve these problems (see [5]). Recently, the level set approach [13] is employed in the numerical algorithms of structural optimization [7] for tracking the evolution of the domain boundary on a fixed mesh and finding an optimal solution to structural optimization problems. This approach, in classical form, is based on an
implicit representation of the boundaries of the optimized structure, i.e., the position of the boundary of the body is described as an isocountour of a scalar function of a higher dimensionality. The evolution of the domain boundary is governed by Hamilton-Jacobi equation. The solution of this equation requires reinitialization procedure to ensure that it is as close as possible to the signed distance function to the interface. Moreover this approach requires regularization of non-differentiable Heaviside and Dirac functions.

In order to avoid the drawbacks of the classical level set method an alternative piecewise constant level set method has been proposed, first in image processing area and next in structural optimization [14]. For a domain divided into $2^N$ sub-domains in classical level set approach is required $2^N$ level set functions to represent them. Piecewise constant level set method can identify an arbitrary number of sub-domains using only one discontinuous piecewise constant level set function. This function takes distinct constant values on each sub-domain. The interfaces between sub-domains are represented implicitly by the discontinuity of a set of characteristic functions of the sub-domains [14]. Comparing to the classical level set method, this method is free of the Hamilton-Jacobi equation and do not require the use of the signed distance function as the initial one. Binary level set method [4, 10, 15] is a special piecewise constant level set method where the function takes only two values either +1 or -1. Compared with general piecewise constant level set approach binary level set approach requires N level set functions to represent a structure of $2^N$ different material phases and is very close to the phase-field approach [12].

In the paper the original structural optimization problem is approximated by a two-phase material optimization problem. Using the binary level set method this approximated problem is reformulated as an equivalent constrained optimization problem in terms of the binary level set function only. Therefore neither shape nor topological complicated sensitivity analysis is required. During the evolution of the binary level set function small holes can be created without the use of the topological derivatives. Necessary optimality condition is formulated. This optimization problem is solved numerically using the augmented Lagrangian method. Numerical examples are provided and discussed.

## 2 Problem formulation

Consider deformations of an elastic body occupying two-dimensional domain $\Omega$ with the smooth boundary $\Gamma$ (see Fig.1). Assume $\Omega \subset D$ where D is a bounded smooth hold-all subset of $\mathbb{R}^2$. Let $E \subset \mathbb{R}^2$ and $D \subset \mathbb{R}^2$ denote given bounded domains. So-called hold-all domain D is assumed to possess a piecewise smooth boundary. Domain $\Omega$ is assumed to belong to the set $O_l$ defined as follows:

$$O_l = \{\Omega \subset \mathbb{R}^2 : \Omega \text{ is open, } E \subset \Omega \subset D, \#\Omega^c \le l\}, \quad (1)$$

where $\#\Omega^c$ denotes the number of connected components of the complement $\Omega^c$ of $\Omega$ with respect to D and $l \ge 1$ is a given integer. Moreover all perturbations $\delta\Omega$ of $\Omega$ are assumed to satisfy $\delta\Omega \in O_l$. The body is subject to body forces
$f(x) = (f_1(x), f_2(x))$, $x \in \Omega$. Moreover, surface tractions $p(x) = (p_1(x), p_2(x))$, $x \in \Gamma$, are applied to a portion $\Gamma_1$ of the boundary $\Gamma$. We assume, that the body is clamped along the portion $\Gamma_0$ of the boundary $\Gamma$, and that the contact conditions are prescribed on the portion $\Gamma_2$, where $\Gamma_i \cap \Gamma_j = \emptyset$, $i \neq j$, $i, j = 0, 1, 2$, $\Gamma = \bar{\Gamma}_0 \cup \bar{\Gamma}_1 \cup \bar{\Gamma}_2$. We denote by $u = (u_1, u_2)$, $u = u(x)$, $x \in \Omega$, the displacement

<figure><img src="image_2.png" /><figcaption>Fig. 1. Initial Domain $\Omega$.</figcaption></figure>

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}(x)e_{kl}(u(x)). \quad (2) $$

We use here and throughout the paper the summation convention over repeated indices [9]. The strain $e_{kl}(u(x))$, $k, l = 1, 2$, is defined by:

$$ e_{kl}(u(x)) = \frac{1}{2}(u_{k,l}(x) + u_{l,k}(x)), \quad (3) $$

where $u_{k,l}(x) = \frac{\partial u_k(x)}{\partial x_l}$. The stress field $\sigma$ satisfies the system of equations [9]

$$ -\sigma_{ij}(x),_j = f_i(x) \quad x \in \Omega, i, j = 1, 2, \quad (4) $$

where $\sigma_{ij}(x),_j = \frac{\partial \sigma_{ij}(x)}{\partial x_j}$, $i, j = 1, 2$. The following boundary conditions are imposed

$$ u_i(x) = 0 \quad \text{on } \Gamma_0, \quad i = 1, 2, \quad (5) $$

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

$$ u_N \le 0, \quad \sigma_N \le 0, \quad u_N \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 the normal components of the displacement $u$ and the stress $\sigma$, respectively. The tangential components of displacement $u$ and stress $\sigma$ are given 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 $R^2$ of the tangent vector $u_T$. The results concerning the existence and uniqueness of solutions to (2)-(8) can be found in [9].

## 2.1 Variational formulation of contact problem

Let us formulate contact problem (4)-(8) in variational form. Denote by $V_{sp}$ and $K$ the space and set of kinematically admissible displacements:

$$V_{sp} = \{z \in [H^1(\Omega)]^2 = H^1(\Omega) \times H^1(\Omega) : z_i = 0 \text{ on } \Gamma_0, i = 1, 2\}, \quad (9)$$

$$K = \{z \in V_{sp} : z_N \le 0 \text{ on } \Gamma_2\}. \quad (10)$$

Denote also by $\Lambda$ the set

$$\Lambda = \{\zeta \in L^2(\Gamma_2) : |\zeta| \le 1\}. \quad (11)$$

Variational formulation of problem (4)-(8) has the form: find a pair $(u, \lambda) \in K \times \Lambda$ satisfying

$$\int_{\Omega} a_{ijkl} e_{ij}(u) e_{kl}(\varphi - u) dx - \int_{\Omega} f_i(\varphi_i - u_i) dx - \int_{\Gamma_1} p_i(\varphi_i - u_i) ds + \int_{\Gamma_2} \lambda(\varphi_T - u_T) ds \ge 0 \quad \forall \varphi \in K, \quad (12)$$

$$\int_{\Gamma_2} (\zeta - \lambda) u_T ds \le 0 \quad \forall \zeta \in \Lambda, \quad (13)$$

$i, j, k, l = 1, 2$. Function $\lambda$ is interpreted as a Lagrange multiplier corresponding to term $|u_T|$ in equality constraint in (8) [9]. In general, function $\lambda$ belongs to the space $H^{-1/2}(\Gamma_2)$. Here following [9] function $\lambda$ is assumed to be more regular. The results concerning the existence and uniqueness of solutions to system (12)-(13) can be found, among others, in [9].

## 2.2 Structural optimization problem

Before formulating a structural optimization problem for the state system (12)-(13) let us introduce first the set $U_{ad}$ of admissible domains. Domain $\Omega$ is assumed to satisfy the volume constraint of the form

$$Vol(\Omega) - Vol^{giv} \le 0, \quad Vol(\Omega) \stackrel{def}{=} \int_{\Omega} dx, \quad (14)$$

where the constant $Vol^{giv} = const_0 > 0$ is given. Moreover this domain is assumed to satisfy the perimeter constraint [6]

$$Per(\Omega) \le const_1, \quad Per(\Omega) \stackrel{def}{=} \int_{\Gamma} dx. \quad (15)$$
The constant $const_1 > 0$ is given. The set $U_{ad}$ has the following form

$$ U_{ad} = \{ \Omega \in O_l : \begin{array}{l} \Omega \text{ is Lipschitz continuous,} \\ \Omega \text{ satisfies conditions (14) and (15)}. \end{array} \} \quad (16) $$

The set $U_{ad}$ is assumed to be nonempty. In order to define a cost functional we shall also need the following set $M^{st}$ of auxiliary functions

$$ M^{st} = \{ \eta = (\eta_1, \eta_2) \in [H^1(D)]^2 : \eta_i \le 0 \text{ on } D, i = 1, 2, \\ \qquad \qquad \qquad \qquad \qquad \qquad \| \eta \|_{[H^1(D)]^2} \le 1 \}, \quad (17) $$

where the norm $\| \eta \|_{[H^1(D)]^2} = (\sum_{i=1}^2 \| \eta_i \|_{H^1(D)}^2)^{1/2}$. Recall from [11] the cost functional approximating the normal contact stress on the contact boundary $\Gamma_2$

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

depending on the auxiliary given bounded function $\eta(x) \in M^{st}$. $\sigma_N$ and $\phi_N$ are the normal components of the stress field $\sigma$ corresponding to a solution $u$ satisfying system (12)-(13) and the function $\eta$, respectively.

Consider the following structural optimization problem: for a given function $\eta \in M^{st}$, find a domain $\Omega^* \in U_{ad}$ such that

$$ J_{\eta}(u(\Omega^*)) = \min_{\Omega \in U_{ad}} J_{\eta}(u(\Omega)) \quad (19) $$

**Lemma 1.** There exists an optimal domain $\Omega^* \in U_{ad}$ to the problem (19).

The proof follows from Šverák theorem and arguments provided in [3, Theorem 2]. Recall from [3] the class of domains $O_l$ determined by (1) is endowed with the complementary Hausdorff topology that guarantees the class itself to be compact. The admissibility condition $\#\Omega^c \le l$ is crucial to provide the necessary compactness property of $U_{ad}$ [3].

### 3 Level set approach

In [11] the standard level set method [13] is employed to solve numerically problem (19). Let $t > 0$ denote the time variable. Consider the evolution of a domain $\Omega$ under a velocity field $V = V(x, t)$. Under the mapping $T(t, V)$ we have

$$ \Omega_t = T(t, V)(\Omega) = (I + tV)(\Omega), \quad t > 0. $$

By $\Omega_t^-$ and $\Omega_t^+$ we denote the interior and the outside of the domain $\Omega_t$, respectively. This domain and its boundary $\partial\Omega_t$ are defined by a function $\phi = \phi(x, t) : \mathbb{R}^2 \times [0, t_0) \to \mathbb{R}$ satisfying the conditions:

$$ \begin{array}{l} \phi(x, t) = 0, \text{ if } x \in \partial\Omega_t, \quad \phi(x, t) < 0, \text{ if } x \in \Omega_t^-, \\ \phi(x, t) > 0, \text{ if } x \in \Omega_t^+. \end{array} \quad (20) $$
In the standard level set approach Heaviside function and Dirac function are used to transform integrals from domain $\Omega$ into domain $D$. Assume that velocity field $V$ is known for every point $x$ lying on the boundary $\partial\Omega_t$, i.e., such that $\phi(x, t) = 0$. Therefore the equation governing the evolution of the interface in $D \times [0, t_0]$, known as Hamilton-Jacobi equation, has the form [13]

$$ \frac{\partial\phi(x,t)}{\partial t} + V(x,t) \cdot \nabla_x\phi(x,t) = 0. \qquad (21) $$

Moreover $\phi(x, 0) = \phi_0$ where $\phi_0(x)$ is a given function close to the signed distance function [13].

## 3.1 Binary level set formulation

Recall from [4, 10] the notion of a binary level set function. Let $N$ be a given integer. Assume an open bounded domain $D$ in $R^2$ is partitioned into $2^N$ sub-domains $\{\Omega_j\}_{j=1}^{2^N}$ such that

$$ D = \bigcup_{j=1}^{2^N} (\Omega_j \cup \partial\Omega_j). \qquad (22) $$

$\partial\Omega_j$ denotes the boundary of the sub-domain $\Omega_j$. For $N = 2$ this function mapping $\phi: D \to R$, is defined as:

$$ \phi(x) = \begin{cases} +1, & \text{if } x \in \Omega_1, \\ -1, & \text{if } x \in \Omega_2 = D \setminus \Omega_1. \end{cases} \qquad (23) $$

The interface $\partial\Omega_1 = \partial\Omega_2$ is implicitly defined by the discontinuity of $\phi$, i.e., $\partial\Omega_1 = \{x \in D : \phi(x) = \kappa, \kappa \in (-1, 1)\}$. In order to ensure that for every $x \in D$ this function converges to values $+1$ and $-1$ it is supposed to satisfy:

$$ W(\phi) \stackrel{\text{def}}{=} (\phi^2 - 1) = 0. \qquad (24) $$

More generally, using $N$ binary level set functions $\phi_i$, $i = 1, 2, ..., N$, satisfying (24) we can represent $2^N$ sub-domains $\Omega_j$ of $D$. The characteristic functions $\chi_j$, $j = 1, 2, ..., 2^N$, of the sub-domains $\Omega_j$ in terms of binary level set functions $\phi_i$ are represented as

$$ \chi_j = \frac{(-1)^{s(j)}}{2^N} \prod_{i=1}^{N} (\phi_i + 1 - 2b_i^{j-1}) \quad \text{where} \quad s(j) = \sum_{i=1}^{N} b_i^{j-1}. \qquad (25) $$

For $j = 1, 2, ..., 2^N$ and $i = 1, ..., N$ numbers $b_i^{j-1} = 0 \lor 1$ denotes binary representation of $j-1$ th sub-domain. As long as each binary function satisfies (24) and $\chi_j(x)$ are defined by (25) then

$$ \text{supp}(\chi_j) = \Omega_j, \quad \chi_j = 1 \text{ in } \Omega_j \text{ and } \text{supp}(\chi_i) \cap \text{supp}(\chi_j) = \emptyset \text{ for } i \neq j, \qquad (26) $$

$$ \sum_j \text{supp}(\chi_j) = \Omega. \qquad (27) $$
Condition (26) ensures non overlapping of the phases while (27) prevents vacu-
ums. Basis functions $\chi_j$ are used to calculate the length of the boundary $\partial\Omega_j$ as
well as the area inside $\Omega_j$ using the integrals:

$$
| \partial \Omega_j | = \int_{\Omega_j} | \nabla \chi_j | dx \quad \text{and} \quad | \Omega_j | = \int_{\Omega_j} \chi_j dx. \qquad (28)
$$

The length of the boundary $\partial\Omega_j$ of sub-domain $\Omega_j$ equals the total variation of
$\chi_j$ [2]. Consider piecewise constant density function $\rho = \rho(x) : D \to \mathbb{R}^2$ defined
as

$$
\rho(x) = \begin{cases} c_1 & \text{if } x \in D \setminus \bar{\Omega}, \\ c_2 & \text{if } x \in \Omega, \end{cases} \tag{29}
$$

where $0 < c_1 < c_2 < \infty$ denote two given material densities. This function can
be constructed as a weighted sum of the characteristic functions $\chi_j$. Denoting
by $\{c_j\}_{j=1}^{2N}$ a set of real scalars, we can represent a piecewise constant function
$\rho$ taking these $2N$ distinct constant values in sub-domains $\Omega_j$ by

$$
\rho(x) = \sum_{j=1}^{2^N} c_j \chi_j(x). \tag{30}
$$

We confine to consider a two-phase problem in the domain $D$, i.e., we set $N = 2$
and $c_1 = \epsilon$, $\epsilon > 0$ as well as $c_2 = 1$. Since $\Omega$ consists from two sub-domains
one binary level set function $\phi$ satisfying (24) will be used to describe these
sub-domains. Therefore

$$
\chi_1(x) = \phi(x) + 1 \quad \text{and} \quad \chi_2(x) = 1 - \phi(x), \tag{31}
$$

$$
\rho(x) = \sum_{j=1}^{2} c_j \rho_j = c_1 \chi_1(x) + c_2 \chi_2(x) = c_1(\phi(x) + 1) + c_2(1 - \phi(x)). \quad (32)
$$

Using (22) as well as (32) the structural optimization problem (19) can be trans-
formed into the following one: find $\phi \in U_{ad}^{\phi}$ such that

$$
\min_{\phi \in U_{ad}^{\phi}} J_{\eta}(\phi) = \int_{\Gamma_2} \rho(\phi) \sigma_N(u_{\epsilon}) \eta_N ds, \quad (33)
$$

where the set $U_{ad}^{\phi}$ of the admissible functions is given as

$$
U_{ad}^{\phi} = \{ \phi \in H^1(D) : \mathrm{Vol}(\phi) - \mathrm{Vol}^{giv} \le 0, W(\phi) = 0, \mathrm{Per}(\phi) \le \mathrm{const}_1 \}, (34)
$$

$$
\mathrm{Vol}(\phi) \stackrel{\mathrm{def}}{=} \int_{\Omega} \rho(\phi) dx, \quad \mathrm{Per}(\phi) \stackrel{\mathrm{def}}{=} \int_{\Omega} |\nabla \phi| dx. \qquad (35)
$$
The element $(u_\epsilon, \lambda_\epsilon) \in K \times \Lambda$ depending on $\epsilon$ satisfies the state system (12)-(13) in the domain $D$ rather than $\Omega$:

$$ \int_D \rho(\phi) a_{ijkl} e_{ij}(u_\epsilon) e_{kl}(\varphi - u_\epsilon) dx - \int_D \rho(\phi) f_i(\varphi_i - u_{\epsilon i}) dx - \\ \int_{\Gamma_1} p_i(\varphi_i - u_{\epsilon i}) ds + \int_{\Gamma_2} \lambda_\epsilon (\varphi_T - u_{\epsilon T}) ds \geq 0 \quad \forall \varphi \in K, \quad (36) $$

$$ \int_{\Gamma_2} (\zeta - \lambda_\epsilon) u_{\epsilon T} ds \leq 0 \quad \forall \zeta \in \Lambda. \quad (37) $$

**Lemma 2.** There exists an optimal solution $\phi \in H^1(D)$ to the optimization problem (33)-(37).

The proof follows from the lower semicontinuity in $L^1(D)$ of the regularization term in (34) see [2, Theorem 3.2.1, p. 75]).

## 3.2 Necessary optimality conditions

In order to formulate the necessary optimality condition for the optimization problem (33)-(37) we introduce the Lagrangian $L(\phi, \tilde{\lambda}) = L(\phi, u_\epsilon, \lambda_\epsilon, p^a, q^a, \tilde{\lambda})$:

$$ L(\phi, \tilde{\lambda}) = J_{\eta}(\phi) + \int_{D} \rho(\phi) a_{ijkl} e_{ij}(u_{\epsilon}) e_{kl}(p^{a}) dx - \int_{D} \rho(\phi) f_{i} p_{i}^{a} dx - \\ \int_{\Gamma_{1}} p_{i} p_{i}^{a} ds + \int_{\Gamma_{2}} \lambda_{\epsilon} p_{T}^{a} ds + \int_{\Gamma_{2}} q^{a} u_{\epsilon T} ds + \tilde{\lambda} d(\phi) + \sum_{i=1}^{3} \frac{1}{2\mu_{i}} d_{i}^{2}(\phi), \quad (38) $$

where $i, j, k, l = 1, 2$, $\tilde{\lambda} = \{\tilde{\lambda}_i\}_{i=1}^3$, $d(\phi) = \{d_i(\phi)\}_{i=1}^3 = [\mathrm{Vol}(\phi), \mathrm{W}(\phi), \mathrm{Per}(\phi)]^T$, $d^T(\phi)$ denotes a transpose of $d(\phi)$, $\mu_m > 0$, $m = 1, 2, 3$, is a given real. Element $(p^a, q^a) \in K_1 \times \Lambda_1$ denotes an adjoint state defined as follows:

$$ \int_D \rho(\phi) a_{ijkl} e_{ij}(\eta + p^a) e_{kl}(\varphi) dx + \int_{\Gamma_2} q^a \varphi_T ds = 0 \quad \forall \varphi \in K_1, \quad (39) $$

$$ \int_{\Gamma_2} \zeta(p_T^a + \eta_T) ds = 0 \quad \forall \zeta \in \Lambda_1. \quad (40) $$

The sets $K_1$ and $\Lambda_1$ are given by

$$ K_1 = \{ \xi \in V_{sp} : \xi_N = 0 \text{ on } A^{st} \}, \quad (41) $$

$$ \Lambda_1 = \{ \zeta \in \Lambda : \zeta(x) = 0 \text{ on } B_1 \cup B_2 \cup B_1^+ \cup B_2^+ \}, \quad (42) $$

while the coincidence set $A^{st} = \{x \in \Gamma_2 : u_N = 0\}$. Moreover $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) = 0\}$, $i = 1, 2$, $B_i^+ = B_i \setminus \tilde{B}_i$, $i = 1, 2$. The derivative of the Lagrangian $L$ with respect to $\phi$ has the form:

$$ \frac{\partial L}{\partial \phi}(\phi, \tilde{\lambda}) = \int_D \rho'(\phi) [a_{ijkl} e_{ij}(u_\epsilon) e_{kl}(p^a + \eta) - f(p^a + \eta)] dx + \\ \tilde{\lambda} d'(\phi) + \sum_{i=1}^{3} \frac{1}{\mu_i} d(\phi) d'(\phi), \quad (43) $$
where $\rho'(\phi) = c_1 - c_2 = 1 - \epsilon$, $d'(\phi) = [Vol'(\phi), W'(\phi), Per'(\phi)]$ and

$$Vol'(\phi) = 1, \quad W'(\phi) = 2\phi, \quad Per'(\phi) = \qquad (44)$$

$$\chi_{\partial\Omega=\text{const}_0} \max\{0, -\nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right)\} - \chi_{\partial\Omega>\text{const}_0} \nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right). \qquad (45)$$

Using (39)-(45) we can formulate the necessary optimality condition:

**Lemma 3.** If $\hat{\phi} \in U_{ad}^{\phi}$ is an optimal solution to the problem (33)-(37) than there exists Lagrange multiplier $\tilde{\lambda}^* = (\tilde{\lambda}_1^*, \tilde{\lambda}_2^*, \tilde{\lambda}_3^*) \in R^3$ such that $\tilde{\lambda}_1^*, \tilde{\lambda}_3^* \ge 0$ satisfying

$$L(\hat{\phi}, \tilde{\lambda}) \le L(\hat{\phi}, \tilde{\lambda}^*) \le L(\phi, \tilde{\lambda}^*) \quad \forall (\phi, \tilde{\lambda}) \in U_{ad}^{\phi} \times R^3. \qquad (46)$$

Proof follows from standard arguments [6, 8]. Recall [6, 9] condition (46) implies that for all $\phi \in U_{ad}^{\phi}$ and $\tilde{\lambda} \in R^3$

$$\frac{\partial L(\hat{\phi}, \tilde{\lambda})}{\partial \phi} \ge 0 \quad \text{and} \quad \frac{\partial L(\phi, \tilde{\lambda}^*)}{\partial \tilde{\lambda}} \le 0. \qquad (47)$$

# 4 Numerical experiments

The optimization problem (33)-(37) is discretized using the finite element method [8, 9]. The finite difference method is used to approximate interface evolution (gradient flow) equation [2]. The discretized structural optimization problem (33)-(37) is solved numerically. We employ Uzawa type algorithm to solve numerically optimization problem (33)-(37). The algorithm is 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\}, \qquad (48)$$

is considered. The boundary $\Gamma$ of the domain $\Omega$ is divided into three disjoint pieces

$$\begin{align}
\Gamma_0 &= \{(x_1, x_2) \in \mathbb{R}^2 : x_1 = 0, 8 \land 0 < v(x_1) \le x_2 \le 4\}, \notag \\
\Gamma_1 &= \{(x_1, x_2) \in \mathbb{R}^2 : 0 \le x_1 \le 8 \land x_2 = 4\}, \tag{49} \\
\Gamma_2 &= \{(x_1, x_2) \in \mathbb{R}^2 : 0 \le x_1 \le 8 \land v(x_1) = x_2\}. \notag
\end{align}$$

The domain $\Omega$ and the boundary $\Gamma_2$ depend on the function $v$ given as in Fig. 1. The obtained optimal domain is presented in Fig. 2. The areas with low values of density function appear in the central part of the body and near the fixed edges. The obtained normal contact stress is almost constant along the optimal shape boundary and has been significantly reduced comparing to the initial one.

# References

1. Allaire, G., Jouve, F., Toader, A.: Structural Optimization Using Sensitivity Analysis and a Level Set Method, Journal of Computational Physics, 194, 363–393 (2004)
<figure><img src="image_3.png" /><figcaption>Fig. 2. Optimal domain $\Omega^*$.</figcaption></figure>

2. Aubert, G., Kornprobst, P.: Mathematical Problems in Image Processing, Springer (2006)

3. Chambolle, A.: A density result in two-dimensional linearized elasticity and applications, Arch. Ration. Mech. Anal., 167, 211-233 (2003)

4. Dai, X., Tang, P., Cheng, X., Wu, M.: A Variational Binary Level Set Method for Structural Topology Optimization, Commun. Comput. Phys., 13(5), 1292-1308 (2013)

5. Deaton, J.D., Grandhi, R.V.: A survey of structural and multidisciplinary continuum topology optimization: post 2000, Struct. Multidisc. Optim., 49, 1-38, (2014)

6. Delfour, M.C., Zolesio, J.P.: Shapes and Geometries: Analysis, Differential Calculus and Optimization, SIAM Publications, Philadelphia (2001)

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

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

9. Hlaváček, I., Haslinger, J., Nečas, J., Lovíšek, J.: Solving of Variational Inequalities in Mechanics, Springer (1988)

10. Lie, J., Lysaker, M., Tai, X.C.: A Binary Level Set Model and Some Applications to Mumford Shah Image Segmentation, IEEE Transactions on Image Processing, 15(5), 1171-1181 (2006)

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

12. Yamada, T., Izui, K., Nishiwaki, S., Takezawa, A.: A topology optimization method based on the level set method incorporating a fictitious interface energy, Comput. Methods Appl. Mech. Engrg., 199(45-48), 2876-2891 (2010)

13. Osher, S., Fedkiw, R.: Level Set Methods and Dynamic Implicit Surfaces, Springer, New York, New York (2003)

14. Zhu, S., Wu, Q., Liu, C.: Shape and topology optimization for elliptic boundary value problems using a piecewise constant level set method, Applied Numerical Mathematics, 61, 752-767 (2011)

15. Zhu, S., Liu, C., Wu, Q.: Binary level set methods for topology and shape optimization of a two-density inhomogeneous drum, Computer Methods in Applied Mechanics and Engineering, 199(45-48), 2970-2986 (2010)