# Multi-material topology optimization of structural load-bearing capacity using limit analysis

Leyla Mourad, Jeremy Bleyer, Romain Mesnil, Joanna Nseir, Karam Sab,

Wassim Raphael

► To cite this version:

Leyla Mourad, Jeremy Bleyer, Romain Mesnil, Joanna Nseir, Karam Sab, et al.. Multi-material topology optimization of structural load-bearing capacity using limit analysis. 2024. hal-03528107v3

HAL Id: hal-03528107

https://enpc.hal.science/hal-03528107v3

Preprint submitted on 5 Feb 2024

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
<figure><img src="image_1.png" /></figure>

HAL
open science

# Multi-material topology optimization of structural load-bearing capacity using limit analysis

Leyla Mourad, Jeremy Bleyer, Romain Mesnil, Joanna Nseir, Karam Sab, Wassim Raphael

► To cite this version:

Leyla Mourad, Jeremy Bleyer, Romain Mesnil, Joanna Nseir, Karam Sab, et al. Multi-material topology optimization of structural load-bearing capacity using limit analysis. Journal of Theoretical, Computational and Applied Mechanics, 2025, pp.1-24. 10.46298/jtcam.12427. hal-03528107v4

HAL Id: hal-03528107

https://enpc.hal.science/hal-03528107v4

Submitted on 12 Feb 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multi-material topology optimization of structural load-bearing capacity using limit analysis

Leyla Mourad<fnref n="1" />,<fnref n="2" />, Jeremy Bleyer<fnref n="1" />, Romain Mesnil<fnref n="1" />, Joanna Nseir<fnref n="2" />, Karam Sab<fnref n="1" />, and Wassim Raphael<fnref n="2" />

We extend the problem of finding an optimal structure with maximum load-bearing capacity to the case of multiple materials. We first consider a reinforcement optimization case where the structure consists of a fixed background matrix material with given strength properties and optimize the reinforcement topology within this material. We discuss the use of various isotropic and anisotropic strength criteria to model the reinforcing phase, including reinforcements with discrete orientations. In a second time, we investigate a bi-material formulation where we optimize the topology of two material phases simultaneously. Various choices for the material strength conditions are proposed and we apply this formulation to the optimization of pure tensile and compressive phases of a single material. In all cases, two optimization variants are proposed using concepts of convex optimization and limit analysis theory, namely maximizing the load-bearing capacity under a fixed volume constraint or minimizing the volume under a fixed loading. Both problems are convex and a penalization procedure is proposed. The underlying problems can be solved using conic programming solvers. Illustrative applications demonstrate the versatility of the proposed formulation, including the influence of the selected strength criteria, the possibility to obtain structures with members of fixed orientation or structures with different importance granted to tensile and compressive regions. Finally, we also draw a parallel with the generation of strut-and-tie models for the analysis of reinforced concrete structures.

Keywords topology optimization, multi-material optimization, limit analysis, homogenization, bearing capacity, conic programming, reinforced concrete

## 1 Introduction

Although topology optimization methods have become increasingly popular, many works concentrate on compliance optimization of elastic materials. Plastic design optimization has been much less investigated although fundamental theoretical works have dealt with both elastic and plastic design (Prager 1974; Strang et al. 1986; Rozvany et al. 1986). Despite nonlinear constitutive relations having been considered in the topology optimization process (Swan et al. 1997; Maute et al. 1998; Amir and Bogomolny 2012; Wallin, Jönsson, et al. 2016; Amir 2017), the numerical cost is very high and there are still many difficulties to be tackled when including local stress constraints (Pedersen 1998; Duysinx et al. 1998). In contrast, plastic design of truss structures (Dorn et al. 1964; Gilbert et al. 2003; He et al. 2015) shows less problems and is now able to tackle large scale problems with the development of efficient linear programming solvers. Interestingly, the extension of the concepts of plastic truss design to continuum topology optimization has been largely ignored until recently when strength-based topology optimization of von Mises plastic materials were proposed (Kammoun et al. 2014; Fin et al. 2018; Herfelt et al. 2018). In a more general manner using the concepts of limit analysis, we formulated in Mourad et al. (2021) the problem of maximizing a structure's load-bearing capacity subject to given material strength properties and a material volume constraint. This problem is also closely linked to the problem of minimizing the total volume under the constraint of carrying a fixed loading as well as to Michell's theory of optimal trusses (Michell 1904). In particular, this approach, which relies on convex optimization solvers, enables to treat the case of materials with different strength properties in tension and compression which are commonly encountered in civil engineering applications.

<footnote n="1">¹ Laboratoire Navier, Ecole des Ponts ParisTech, Univ Gustave Eiffel, CNRS, France</footnote>

<footnote n="2">² Université Saint Joseph, Faculté des sciences, Mar Roukos-Dekwaneh, Lebanon</footnote>
**Identifiers**

DOI 10.46298/jtcam.12427

HAL hal-03528107v4

**History**

Received Oct, 17, 2023

Accepted Apr, 24, 2024

Published Feb 07, 2025

**Associate Editor**

Alexander POPP

**Reviewers**

Fabian DUDDECK

Anonymous

**Open Review**

OAI hal-04693259

**Supplementary Material**

See article's addendum

**Licence**

CC BY 4.0

©The Authors

# Multi-material topology optimization of structural load-bearing capacity using limit analysis

Leyla MOURAD¹,², Jérémy BLEYER¹, Romain MESNIL¹, Joanna NSEIR², Karam SAB¹, and Wassim RAPHAEL²

¹ Laboratoire Navier, École des Ponts ParisTech, Université Gustave Eiffel, CNRS, France

² Université Saint Joseph, Faculté des sciences, Mar Roukos-Dekwaneh, Lebanon

We extend the problem of finding an optimal structure with maximum load-bearing capacity to the case of multiple materials. We first consider a reinforcement optimization case where the structure consists of a fixed background matrix material with given strength properties and optimize the reinforcement topology within this material. We discuss the use of various isotropic and anisotropic strength criteria to model the reinforcing phase, including reinforcements with discrete orientations. In a second time, we investigate a bi-material formulation where we optimize the topology of two material phases simultaneously. Various choices for the material strength conditions are proposed and we apply this formulation to the optimization of pure tensile and compressive phases of a single material. In all cases, two optimization variants are proposed using concepts of convex optimization and limit analysis theory, namely maximizing the load-bearing capacity under a fixed volume constraint or minimizing the volume under a fixed loading. Both problems are convex and a penalization procedure is proposed. The underlying problems can be solved using conic programming solvers. Illustrative applications demonstrate the versatility of the proposed formulation, including the influence of the selected strength criteria, the possibility to obtain structures with members of fixed orientation or structures with different importance granted to tensile and compressive regions. Finally, we also draw a parallel with the generation of strut-and-tie models for the analysis of reinforced concrete structures.

**Keywords:** topology optimization, multi-material optimization, limit analysis, homogenization, bearing capacity, conic programming, reinforced concrete

## 1 Introduction

Although topology optimization methods have become increasingly popular, many works concentrate on compliance optimization of elastic materials. Plastic design optimization has been much less investigated although fundamental theoretical works have dealt with both elastic and plastic design (Prager 1974; Strang and Kohn 1986; Rozvany and Ong 1986). Despite nonlinear constitutive relations having been considered in the topology optimization process (Swan and Kosaka 1997; Maute et al. 1998; Bogomolny and Amir 2012; Wallin et al. 2016; Amir 2017), the numerical cost is very high and there are still many difficulties to be tackled when including local stress constraints (Podersen 1998; Duysinx and Bendsøe 1998). In contrast, plastic design of truss structures (Dorn et al. 1964; Gilbert and Tyas 2003; He and Gilbert 2015) shows less problems and is now able to tackle large scale problems with the development of efficient linear programming solvers. Interestingly, the extension of the concepts of plastic truss design to continuum topology optimization has been largely ignored until recently when strength-based topology optimization of von Mises plastic materials were proposed (Kammoun and Smaoui 2014; Fin et al. 2018; Herfelt et al. 2018). In a more general manner using the concepts of limit analysis, we formulated in (Mourad et al. 2021) the problem of maximizing a structure's load-bearing capacity subject to given material strength properties and a material volume constraint. This problem is also closely linked to the problem of minimizing the total volume under the constraint of carrying a fixed loading as well as to Michell's theory of optimal trusses (Michell 1904). In
Apart from topology optimization, interesting mechanical performances can be achieved when mixing materials with different mechanical properties, motivating the extension of topology optimization approaches towards the multi-material case. For instance, the density-based SIMP methodology was extended to multiple material phases in various works (Sigmund et al. 1997; Gibiansky et al. 2000; Guest 2009; Zuo et al. 2017) whereas others used a phase-field approach based on Allen-Cahn or Cahn-Hillard models (Zhou et al. 2007; Tavakoli 2014; Wallin, Ivarsson, et al. 2015) or level-set methods (Allaire et al. 2014; Liu et al. 2018). Only few works considered multi-material topology optimization with material non-linearity such as (X. S. Zhang et al. 2018; Rostami et al. 2020; X. Zhang et al. 2020). The extension of continuous plastic-based (or limit-analysis based) topology optimization to multi-materials, which is the main purpose of the present work, has never been done.

As it will be illustrated later, our work also bears conceptual similarities with the so-called strut-and-tie method frequently used for the ultimate design of reinforced concrete (RC) structures using a truss analogy (Ritter 1899; Mörsch et al. 1909; Schlaich et al. 1987). Compressed concrete regions are idealized as compressive struts and steel rebars as tensile ties, both of them experiencing uniaxial stress states and connecting to each other at nodes which sustain multi-axial stress states. The Strut-and-Tie (ST) method is used in common engineering practice as a hand-based procedure to verify a RC structure bearing capacity once the ST model has been established. It is also used, at an earlier stage of the design, to find an efficient layout of steel rebars by considering different ST models. Optimizing for the steel rebars position is difficult to automate since it also strongly depends on the compressive struts' layout. Previous works proposed to identify a ST model based on the flow of elastic stresses of the continuous structure (Marti 1980; Collins et al. 1980) or by optimizing for the structure compliance using topology optimization concepts (Gaynor et al. 2013; Bruggi 2009), arguing that a stiff structure exhibits the best load-deformation behavior. Only few references such as Querin et al. (2010); Victoria Nicolás et al. (2011); Bruggi (2016); Smarslik et al. (2019) generated ST models using topology optimization by considering different tensile and compressive properties for steel and concrete phases. We can also mention references based on a ground-structure layout optimization approach (Chavez 2018; X. S. Zhang et al. 2018; Smarslik et al. 2019). Overall, the majority of such works are based on compliance minimization.

Our contribution therefore aims at bridging the gap between multi-material topology optimization, plastic or limit analysis-based topology optimization and the generation of strut-and-tie like models using materials with different tensile and compressive strengths. Doing so, we aim at formulating a generic framework for the optimization of the load-bearing capacity of reinforced structures, either by considering a fixed background matrix material and optimizing for the reinforcing phase only or by simultaneously optimizing for the topology of two different phases (in addition to a potential void phase). As it will be seen, this second possibility will offer a natural framework for the generation of strut-and-tie models by considering tensile-only and compression-only phases. Following the same spirit as our previous work (Mourad et al. 2021), our approach will be strongly based on the concepts of convex optimization, especially regarding the formulation of isotropic and anisotropic strength conditions promoting specific features such as uniaxial or pure tension/compression stress states. As a result, our numerical implementation uses convex optimization solvers and solves simultaneously for the mechanical state and the optimal density field(s), thereby avoiding the need to compute sensitivities for instance.

The present manuscript is organized as follows: Section 2 first recalls concepts and notations of the limit analysis-based topology optimization methodology proposed in Mourad et al. (2021), Section 3 then discusses its extension to reinforcement optimization and the choice of a strength condition for the reinforcement phase, Section 4 presents a bimaterial formulation optimizing two material phases in addition to a void phase and similarly discusses different possible choices of strength criteria for both phases and Section 5 finally discusses the numerical aspects and presents various illustrative applications of the proposed approach.
particular, this approach, which relies on convex optimization solvers, enables to treat the case of materials with different strength properties in tension and compression which are commonly encountered in civil engineering applications.

Apart from topology optimization, interesting mechanical performances can be achieved when mixing materials with different mechanical properties, motivating the extension of topology optimization approaches towards the multi-material case. For instance, the density-based SIMP methodology was extended to multiple material phases in various works (Sigmund and Torquato 1997; Gibiansky and Sigmund 2000; Guest 2009; Zuo and Saitou 2017) whereas others used a phase-field approach based on Allen-Cahn or Cahn-Hillard models (Zhou and Wang 2007; Tavakoli 2014; Wallin et al. 2015) or level-set methods (Allaire et al. 2014; Liu and Ma 2018). Only few works considered multi-material topology optimization with material non-linearity such as (Zhang et al. 2018; Rostami and Marzbanrad 2020; Zhang et al. 2020). The extension of continuous plastic-based (or limit-analysis based) topology optimization to multi-materials, which is the main purpose of the present work, has never been done.

As it will be illustrated later, our work also bears conceptual similarities with the so-called strut-and-tie method frequently used for the ultimate design of reinforced concrete (RC) structures using a truss analogy (Ritter 1899; Mörsch and Goodrich 1909; Schlaich et al. 1987). Compressed concrete regions are idealized as compressive struts and steel rebars as tensile ties, both of them experiencing uniaxial stress states and connecting to each other at nodes which sustain multi-axial stress states. The Strut-and-Tie (ST) method is used in common engineering practice as a hand-based procedure to verify a RC structure bearing capacity once the ST model has been established. It is also used, at an earlier stage of the design, to find an efficient layout of steel rebars by considering different ST models. Optimizing for the steel rebars position is difficult to automate since it also strongly depends on the compressive struts' layout. Previous works proposed to identify a ST model based on the flow of elastic stresses of the continuous structure (Marti 1980; Collins and Mitchell 1980) or by optimizing for the structure compliance using topology optimization concepts (Gaynor et al. 2013; Bruggi 2009), arguing that a stiff structure exhibits the best load-deformation behavior. Only few references such as Querin et al. (2010); Victoria et al. (2011); Bruggi (2016); Smarslik et al. (2019) generated ST models using topology optimization by considering different tensile and compressive properties for steel and concrete phases. We can also mention references based on a ground-structure layout optimization approach (Chavez 2018; Zhang et al. 2018; Smarslik et al. 2019). Overall, the majority of such works are based on compliance minimization.

Our contribution therefore aims at bridging the gap between multi-material topology optimization, plastic or limit analysis-based topology optimization and the generation of strut-and-tie like models using materials with different tensile and compressive strengths. Doing so, we aim at formulating a generic framework for the optimization of the load-bearing capacity of reinforced structures, either by considering a fixed background matrix material and optimizing for the reinforcing phase only or by simultaneously optimizing for the topology of two different phases (in addition to a potential void phase). As it will be seen, this second possibility will offer a natural framework for the generation of strut-and-tie models by considering tensile-only and compression-only phases. Following the same spirit as our previous work (Mourad et al. 2021), our approach will be strongly based on the concepts of convex optimization, especially regarding the formulation of isotropic and anisotropic strength conditions promoting specific features such as uniaxial or pure tension/compression stress states. As a result, our numerical implementation uses convex optimization solvers and solves simultaneously for the mechanical state and the optimal density field(s), thereby avoiding the need to compute sensitivities for instance.

The present manuscript is organized as follows: Section 2 first recalls concepts and notations of the limit analysis-based topology optimization methodology proposed in (Mourad et al. 2021), Section 3 then discusses its extension to reinforcement optimization and the choice of a strength condition for the reinforcement phase, Section 4 presents a bimaterial formulation optimizing two material phases in addition to a void phase and similarly discusses different possible choices of strength criteria for both phases and Section 5 finally discusses the numerical aspects and presents various illustrative applications of the proposed approach.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Two variants of multi-material optimization. Left: reinforcement optimization where the background "matrix" material is fixed and we optimize a single "reinforcement" phase, right: general bi-material optimization where we optimize over two materials and a void phase</figcaption></figure>

## 2 A review of limit analysis-based topology optimization

In this section, we recall the general concepts of limit analysis-based topology optimization. We follow the exposition introduced by Mourad et al. (2021) which presented two different formulations: a load maximization and a volume minimization formulation for generic strength criteria. We refer to (Kammoun et al. 2014), and later (Herfelt et al. 2018; Fin et al. 2018), for volume minimization formulations based on von Mises plasticity.

First, limit analysis aims at finding the maximum load amplification factor $\lambda$ for which there exists an internal stress field $\sigma$ which can balance the loading and still comply with the material strength condition $\sigma \in G$ at every point in the domain $\Omega$. Here, $G$ is the strength criterion, or the plastic yield criterion when refering to a perfectly plastic material. It is assumed to be a convex set containing 0. According to such a characterization, the limit load $\Lambda^+$ can be found as the solution to the following convex optimization problem:

$$ 
\begin{array}{ll}
\Lambda^+ = \max_{\lambda, \sigma} & \lambda \\
\text{s.t.} & \operatorname{div} \sigma = 0 \quad \text{in } \Omega \\
& \sigma \cdot n = \lambda T \quad \text{on } \partial\Omega_T \\
& \sigma \in G \quad \text{in } \Omega
\end{array}
\qquad (1) $$

in which we neglect body forces for simplicity and look for the maximal value of the reference surface tractions $T$ acting on some part $\partial\Omega_T$ of the boundary.

### 2.1 A load-maximization problem

Considering a computational domain $\mathcal{D}$, we proposed in Mourad et al. (2021) an extension of the concepts of limit analysis to the determination of an optimized structure $\Omega \subseteq \mathcal{D}$ which would have the maximum load-bearing capacity for a given material volume fraction constraint $|\Omega| \le \eta|\mathcal{D}|$. Formulating this problem as a non-convex binary optimization problem and considering its convexification through a continuous pseudo-density $\rho(x) \in [0; 1]$, we define the following load-maximization problem:

$$ 
\begin{array}{ll}
\lambda^+ = \max_{\lambda, \sigma, \rho} & \lambda \\
\text{s.t.} & \operatorname{div} \sigma = 0 \quad \text{in } \mathcal{D} \\
& \sigma \cdot n = \lambda T \quad \text{on } \partial\mathcal{D}_T \\
& \sigma \in \rho G \quad \text{in } \mathcal{D} \\
& \int_{\mathcal{D}} \rho \, dx \le \eta |\mathcal{D}| \\
& 0 \le \rho \le 1
\end{array}
\qquad (\text{LOAD-MAX}) $$

where the main difference with Equation (1) comes from the density-dependent strength criterion $\rho G$. Clearly, when $\rho = 0$, we have $\sigma = 0$ and when $\rho = 1$, $\sigma \in G$ which correspond to the initial material strength condition. Finally, one key property of the above problem is that it is convex and can be solved using dedicated conic programming solvers for many usual strength criteria.
## 2 A review of limit analysis-based topology optimization

In this section, we recall the general concepts of limit analysis-based topology optimization. We follow the exposition introduced by Mourad et al. (2021) which presented two different formulations: a load maximization and a volume minimization formulation for generic strength criteria. We refer to (Kammoun and Smaoui 2014), and later (Herfelt et al. 2018; Fin et al. 2018), for volume minimization formulations based on von Mises plasticity.

First, limit analysis aims at finding the maximum load amplification factor $\lambda$ for which there exists an internal stress field $\sigma$ which can balance the loading and still comply with the material strength condition $\sigma \in G$ at every point in the domain $\Omega$. Here, $G$ is the strength criterion, or the plastic yield criterion when referring to a perfectly plastic material. It is assumed to be a convex set containing 0. According to such a characterization, the limit load $\Lambda^+$ can be found as the solution to the convex optimization problem

$$ 
\begin{array}{lll}
\Lambda^+ = \max_{\lambda, \sigma} & \lambda & \\
\text{s.t.} & \operatorname{div} \sigma = 0 & \text{in } \Omega \\
& \sigma \cdot n = \lambda T & \text{on } \partial \Omega_T \\
& \sigma \in G & \text{in } \Omega
\end{array}
\qquad (1) $$

in which we neglect body forces for simplicity and look for the maximal value of the reference surface tractions $T$ acting on some part $\partial \Omega_T$ of the boundary.

### 2.1 A load-maximization problem

Considering a computational domain $\mathcal{D}$, we proposed in (Mourad et al. 2021) an extension of the concepts of limit analysis to the determination of an optimized structure $\Omega \subseteq \mathcal{D}$ which would have the maximum load-bearing capacity for a given material volume fraction constraint $|\Omega| \le \eta|\mathcal{D}|$. Formulating this problem as a non-convex binary optimization problem and considering its convexification through a continuous pseudo-density $\rho(x) \in [0; 1]$, we define the load-maximization problem

$$ 
\begin{array}{lll}
\lambda^+ = \max_{\lambda, \sigma, \rho} & \lambda & \\
\text{s.t.} & \operatorname{div} \sigma = 0 & \text{in } \mathcal{D} \\
& \sigma \cdot n = \lambda T & \text{on } \partial \mathcal{D}_T \\
& \sigma \in \rho G & \text{in } \mathcal{D} \\
& \int_{\mathcal{D}} \rho \, dx \le \eta |\mathcal{D}| & \\
& 0 \le \rho \le 1 &
\end{array}
\qquad (\text{LOAD-MAX}) $$

where the main difference with Equation (1) comes from the density-dependent strength criterion $\rho G$. Clearly, when $\rho = 0$, we have $\sigma = 0$ and when $\rho = 1$, $\sigma \in G$ which correspond to the initial material strength condition. Finally, one key property of the above problem is that it is convex and can be solved using dedicated conic programming solvers for many usual strength criteria.

### 2.2 A related volume-minimization problem

Finally, for a fixed loading level $\lambda$, one can also find the structure sustaining the loading with the minimum volume. A similar volume-minimization problem therefore reads as (Kammoun and Smaoui 2014; Herfelt et al. 2018)

$$ 
\begin{array}{lll}
\eta^- = \min_{\sigma, \rho} & \displaystyle\frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} \rho \, dx & \\
\text{s.t.} & \operatorname{div} \sigma = 0 & \text{in } \mathcal{D} \\
& \sigma \cdot n = \lambda T & \text{on } \partial \mathcal{D}_T \\
& \sigma \in \rho G & \text{in } \mathcal{D} \\
& 0 \le \rho \le 1 &
\end{array}
\qquad (\text{VOL-MIN}) $$
## 2.2 A related volume-minimization problem

Finally, for a fixed loading level $\lambda$, one can also find the structure sustaining the loading with the minimum volume. A similar volume-minimization problem therefore reads as (Kammoun et al. 2014; Herfelt et al. 2018):

$$ 
\begin{array}{ll}
\eta^{-} = \min_{\sigma, \rho} & \displaystyle \frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} \rho \, dx \\
\text{s.t.} & \displaystyle \text{div } \sigma = 0 \quad \text{in } \mathcal{D} \\
& \displaystyle \sigma \cdot n = \lambda T \quad \text{on } \partial \mathcal{D}_T \\
& \displaystyle \sigma \in \rho G \quad \text{in } \mathcal{D} \\
& 0 \le \rho \le 1
\end{array}
\qquad (\text{VOL-MIN})
$$

which is again a convex problem.

Note that, contrary to an elastic problem with stress constraints, both formulations Equation (LOAD-MAX) and Equation (VOL-MIN) do not presuppose any constitutive relation and are able to account for stress redistributions until plastic collapse, following the principles of limit analysis.

## 2.3 Choice of the material strength criterion

The choice of the material strength criterion $G$ will have an important impact on the optimal design. For instance, criteria with anisotropic strength properties will induce anisotropic features in the optimized design. Regarding isotropic materials, the optimized structure load-capacity depends essentially upon the material uniaxial tensile (resp. compressive) strength $f_t$ (resp. $f_c$). Indeed, optimized topologies usually have a truss-like layout with many members subjected to uniaxial stress states. In Mourad et al. (2021), the so-called $L_1$-Rankine criterion $G^{L_1}$-Rankine was introduced by the following strength condition:

$$ \sigma \in G^{L_1}$-Rankine \Leftrightarrow g(\sigma_I) + g(\sigma_{II}) \le 1 \qquad (2) $$

written here in 2D with $\sigma_I$, $\sigma_{II}$ being the principal stresses and where:

$$ g(\sigma) = \max \left\{ -\frac{\sigma}{f_c}; \frac{\sigma}{f_t} \right\} \qquad (3) $$

It was advocated in Mourad et al. (2021) that the $L_1$-Rankine criterion was a good choice for promoting uniaxial stress states in the optimal solutions. This argument will be revisited in Section 3.4. Moreover, with such a bounded criterion, we easily see that when $\rho \to 0$, the density dependent criterion $\rho G^{L_1}$-Rankine $\to \{0\}$ which correctly enforces a zero stress state in the limit of a void material. This is in contrast with unbounded criteria such as a 3D or plane-strain von Mises criterion for instance. In the limit $\rho \to 0$, the latter still sustain non-zero spherical stress states $\sigma_I$, with $I$ denoting the identity tensor. In this case, we would recommend to use a strength criterion with a large but finite hydrostatic strength to avoid any residual strength.

# 3 Reinforcement material optimization

In this section, we investigate the situation in which we optimize a single phase with a fixed background phase (see Figure 1-left). This situation essentially applied to the optimization of a reinforcing phase in a composite material, the background fixed phase then corresponding to the matrix material. We will respectively denote by $G^m$ and $G^r$ the strength criteria of the matrix and reinforcement bulk materials.

## 3.1 Reinforced material strength conditions

In general, reinforcements are present in a small volume fraction $\phi \ll 1$. In such a case, the strength criterion of the composite material $G^{\text{comp}}$ can be approximated by the following dilute
which is again a convex problem.

Note that, contrary to an elastic problem with stress constraints, both formulations Equation (LOAD-MAX) and Equation (VOL-MIN) do not presuppose any constitutive relation and are able to account for stress redistributions until plastic collapse, following the principles of limit analysis.

## 2.3 Choice of the material strength criterion

The choice of the material strength criterion G will have an important impact on the optimal design. For instance, criteria with anisotropic strength properties will induce anisotropic features in the optimized design. Regarding isotropic materials, the optimized structure load-capacity depends essentially upon the material uniaxial tensile (resp. compressive) strength $f_t$ (resp. $f_c$). Indeed, optimized topologies usually have a truss-like layout with many members subjected to uniaxial stress states. in (Mourad et al. 2021), the so-called L₁-Rankine criterion $G^{L_1-Rankine}$ was introduced by the strength condition

$$ \sigma \in G^{L_1-Rankine} \Leftrightarrow g(\sigma_I) + g(\sigma_{II}) \le 1 \quad (2) $$

written here in 2D with $\sigma_I$, $\sigma_{II}$ being the principal stresses and where $g(\sigma) = \max\{-\sigma/f_c; \sigma/f_t\}$. It was advocated in (Mourad et al. 2021) that the L₁-Rankine criterion was a good choice for promoting uniaxial stress states in the optimal solutions. This argument will be revisited in Section 3.4. Moreover, with such a bounded criterion, we easily see that when $\rho \to 0$, the density dependent criterion $\rho G^{L_1-Rankine} \to \{0\}$ which correctly enforces a zero stress state in the limit of a void material. This is in contrast with unbounded criteria such as a 3D or plane-strain von Mises criterion for instance. In the limit $\rho \to 0$, the latter still sustain non-zero spherical stress states $\sigma_I$, with I denoting the identity tensor. In this case, we would recommend to use a strength criterion with a large but finite hydrostatic strength to avoid any residual strength.

# 3 Reinforcement material optimization

In this section, we investigate the situation in which we optimize a single phase with a fixed background phase, see Figure 1-left. This situation essentially applied to the optimization of a reinforcing phase in a composite material, the background fixed phase then corresponding to the matrix material. We will respectively denote by $G^m$ and $G^r$ the strength criteria of the matrix and reinforcement bulk materials.

<figure><img src="image_3.png" /><figcaption>**Figure 1** Two variants of multi-material optimization. Left: reinforcement optimization where the background “matrix” material is fixed and we optimize a single “reinforcement” phase, right: general bi-material optimization where we optimize over two materials and a void phase.</figcaption></figure>

## 3.1 Reinforced material strength conditions

In general, reinforcements are present in a small volume fraction $\phi \ll 1$. In such a case, the strength criterion of the composite material $G^{comp}$ can be approximated by the dilute estimation (de Buhan et al. 2017)

$$ \sigma \in G^{comp} \Leftrightarrow \exists \sigma^m, \sigma^r \text{ s.t. } \begin{cases} \sigma = \sigma^m + \sigma^r \\ \sigma^m \in G^m \\ \sigma^r \in G^{r,eff} \end{cases} \quad (3) $$
estimation (de Buhan, Bleyer, et al. 2017):

$$ \sigma \in G^{\text{comp}} \Leftrightarrow \exists \sigma^m, \sigma^r \text{ s.t. } \begin{cases} \sigma = \sigma^m + \sigma^r \\ \sigma^m \in G^m \\ \sigma^r \in G^{r, \text{eff}} \end{cases} \qquad (4) $$

where $G^{r,eff} = \phi G^r$. This formulation states that the composite strength condition corresponds to a stress state which is the sum of a stress $\sigma^m$ satisfying the bulk matrix condition and a contribution of the effective strength condition of the reinforcement material, the strength property of which being that of the bulk reinforcement material scaled by the reinforcement volume fraction.

## 3.2 Convex problems formulation

We assume here that the matrix material is fixed and occupies the whole computational domain $\mathcal{D}$ and we aim at optimizing over the reinforcement phase only. The approach of Mourad et al. (2021) can be easily extended by considering the a density-dependent composite strength criterion $G^{\text{comp}}(\rho)$, defined by the following strength condition:

$$ \sigma \in G^{\text{comp}}(\rho) = G^m \oplus \rho G^{r, \text{eff}} \Leftrightarrow \exists \sigma^m, \sigma^r \text{ s.t. } \begin{cases} \sigma = \sigma^m + \sigma^r \\ \sigma^m \in G^m \\ \sigma^r \in \rho G^{r, \text{eff}} \end{cases} \qquad (5) $$

where $\oplus$ denotes the Minkowski sum between two sets<fnref n="1" />. Note that we still have $G^{\text{comp}}(\rho = 1) = G^{\text{comp}}$. The main difference with the single material optimization for which $G(0) = \{0\}$ was representing void is that we now have $G^{\text{comp}}(\rho = 0) = G^m$. As a result, the above strength criterion will interpolate for $\rho \in [0; 1]$ between the pure matrix strength condition $G^m$ and the reinforced composite criterion $G^{\text{comp}}$.

The two corresponding load-maximization Equation (LOAD-MAX) and volume minimization Equation (VOL-MIN) problems are therefore respectively given by:

$$ \begin{array}{l} \lambda^+ = \max_{\lambda, \sigma^m, \sigma^r, \rho} \lambda \\ \text{s.t.} \quad \begin{array}{ll} \text{div}(\sigma^m + \sigma^r) = 0 & \text{in } \mathcal{D} \\ (\sigma^m + \sigma^r) \cdot n = \lambda T & \text{on } \partial \mathcal{D}_T \\ \sigma^m \in G^m & \text{in } \mathcal{D} \\ \sigma^r \in \rho G^{r, \text{eff}} & \text{in } \mathcal{D} \\ \int_{\mathcal{D}} \rho \, dx \le \eta |\mathcal{D}| & \\ 0 \le \rho \le 1 & \end{array} \end{array} \qquad (\text{REINF-LOAD-MAX}) $$

and

$$ \begin{array}{l} \eta^- = \min_{\sigma^m, \sigma^r, \rho} \frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} \rho \, dx \\ \text{s.t.} \quad \begin{array}{ll} \text{div}(\sigma^m + \sigma^r) = 0 & \text{in } \mathcal{D} \\ (\sigma^m + \sigma^r) \cdot n = \lambda T & \text{on } \partial \mathcal{D}_T \\ \sigma^m \in G^m & \text{in } \mathcal{D} \\ \sigma^r \in \rho G^{r, \text{eff}} & \text{in } \mathcal{D} \\ 0 \le \rho \le 1 & \end{array} \end{array} \qquad (\text{REINF-VOL-MIN}) $$

where the total stress $\sigma$ has been replaced by the sum of the two partial stresses arising in the definition Equation (5). As a result, both problems are very similar to the single material counterpart, which can be obtained in the particular case $\sigma^m = 0$. The numerical implementation,

<footnote n="1">¹ $A \oplus B = \{a+b \text{ s.t. } a \in A, b \in B\}$</footnote>
where $G^{\text{r,eff}} = \phi G^r$. This formulation states that the composite strength condition corresponds to a stress state which is the sum of a stress $\sigma^m$ satisfying the bulk matrix condition and a contribution of the effective strength condition of the reinforcement material, the strength property of which being that of the bulk reinforcement material scaled by the reinforcement volume fraction.

## 3.2 Convex problems formulation

We assume here that the matrix material is fixed and occupies the whole computational domain $D$ and we aim at optimizing over the reinforcement phase only. The approach of Mourad et al. (2021) can be easily extended by considering a density-dependent composite strength criterion $G^{\text{comp}}(\rho)$, defined by the strength condition

$$ \sigma \in G^{\text{comp}}(\rho) = G^m \oplus \rho G^{\text{r,eff}} \Leftrightarrow \exists \sigma^m, \sigma^r \text{ s.t. } \begin{cases} \sigma = \sigma^m + \sigma^r \\ \sigma^m \in G^m \\ \sigma^r \in \rho G^{\text{r,eff}} \end{cases} \quad (4) $$

where $\oplus$ denotes the Minkowski sum between two sets<fnref n="1" />. Note that we still have $G^{\text{comp}}(\rho = 1) = G^{\text{comp}}$. The main difference with the single material optimization for which $G(0) = \{0\}$ was representing void is that we now have $G^{\text{comp}}(\rho = 0) = G^m$. As a result, the above strength criterion will interpolate for $\rho \in [0; 1]$ between the pure matrix strength condition $G^m$ and the reinforced composite criterion $G^{\text{comp}}$.

The two corresponding load-maximization (LOAD-MAX) and volume minimization (VOL-MIN) problems are therefore respectively given by

$$ \begin{array}{ll} \lambda^+ = \max_{\lambda, \sigma^m, \sigma^r, \rho} & \lambda \\ \text{s.t.} & \text{div}(\sigma^m + \sigma^r) = 0 \quad \text{in } D \\ & (\sigma^m + \sigma^r) \cdot n = \lambda T \quad \text{on } \partial D_T \\ & \sigma^m \in G^m \quad \text{in } D \\ & \sigma^r \in \rho G^{\text{r,eff}} \quad \text{in } D \\ & \int_D \rho \, dx \le \eta |D| \\ & 0 \le \rho \le 1 \end{array} \qquad (\text{REINF-LOAD-MAX}) $$

and

$$ \begin{array}{ll} \eta^- = \min_{\sigma^m, \sigma^r, \rho} & \frac{1}{|D|} \int_D \rho \, dx \\ \text{s.t.} & \text{div}(\sigma^m + \sigma^r) = 0 \quad \text{in } D \\ & (\sigma^m + \sigma^r) \cdot n = \lambda T \quad \text{on } \partial D_T \\ & \sigma^m \in G^m \quad \text{in } D \\ & \sigma^r \in \rho G^{\text{r,eff}} \quad \text{in } D \\ & 0 \le \rho \le 1 \end{array} \qquad (\text{REINF-VOL-MIN}) $$

where the total stress $\sigma$ has been replaced by the sum of the two partial stresses arising in the definition Equation (4). As a result, both problems are very similar to the single material counterpart, which can be obtained in the particular case $\sigma^m = 0$. The numerical implementation, including the penalty procedure, will therefore be a straightforward adaptation to that described in (Mourad et al. 2021). Note in particular that only a global equilibrium condition on the total stress $\sigma = \sigma^m + \sigma^r$ needs to be considered so that standard finite-element discretization can be used to enforce this equation weakly.

Before finishing this section, let us now discuss some particular choices for $G^{\text{r,eff}}$ in the case of uniaxial reinforcements.

<footnote n="1">¹ $A \oplus B = \{a+b \text{ s.t. } a \in A, b \in B\}$</footnote>
including the penalty procedure, will therefore be a straightforward adaptation to that described in Mourad et al. (2021). Note in particular that only a global equilibrium condition on the total stress $\sigma = \sigma^m + \sigma^r$ needs to be considered so that standard finite-element discretization can be used to enforce this equation weakly.

Before finishing this section, let us now discuss some particular choices for $G^{r,eff}$ in the case of uniaxial reinforcements.

### 3.3 The particular case of uniaxial reinforcements

When considering uniaxial reinforcements aligned along a direction $e_\alpha$ in small volume fraction, the composite homogenized strength criterion is in fact exactly given by the dilute estimation Equation (4) (de Buhan and Taliercio 1991; de Buhan, Bleyer, et al. 2017) which reads here:

$$ \sigma \in G^{\text{comp}} \Leftrightarrow \exists \sigma^m, \sigma^r \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma = \sigma^m + \sigma^r e_\alpha \otimes e_\alpha \\ \sigma^m \in G^m \\ -f_c^{r,\text{eff}} \le \sigma^r \le f_t^{r,\text{eff}} \end{array} \right. \qquad (6) $$

where $f_t^{r,eff}$ (resp. $f_c^{r,eff}$) denotes the reinforcement uniaxial effective tensile (resp. compressive) strength (per unit of transverse area). As a result, the case of uniaxial reinforcements is similar to Equation (4) except that the effective reinforcement strength criterion $G^{r,eff}$ is now given by the corresponding uniaxial strength condition:

$$ G^{r,eff} = \{ \sigma^r e_\alpha \otimes e_\alpha \text{ s.t. } -f_c^{r,eff} \le \sigma^r \le f_t^{r,eff} \} \qquad (7) $$

Finally, the homogenized strength criterion Equation (6) and Equation (7) easily generalize to a reinforcing material made of multiple reinforcement directions by summing the corresponding uniaxial stress contributions. For instance, an important practical case of interest is that of orthogonal reinforcements aligned with the global x, y directions (and possibly z in 3D). In this case, the effective strength criterion Equation (7) for two reinforcement directions $e_x$ and $e_y$ generalizes to:

$$ G^{r,eff} = \{ \sigma^{r,x} e_x \otimes e_x + \sigma^{r,y} e_y \otimes e_y \text{ s.t. } -f_c^{r,eff} \le \sigma^{r,x}, \sigma^{r,y} \le f_t^{r,eff} \} \qquad (8) $$

### 3.4 The L₁-Rankine criterion for isotropically distributed reinforcements

Similarly to the concepts of the homogenization method in topology optimization, our ultimate goal is to find an optimal microstructure for the reinforcing phase at each material point. To do so, instead of considering that the reinforcing material is made of a fixed distribution of predefined orientations, we can consider a reinforcing material consisting of uniaxial reinforcements but with locally unknown orientation a priori, with the goal that the topology optimization process would naturally select the locally optimal orientation. In Mourad et al. (2021), the use of a L₁-Rankine was proposed in order to promote uniaxial stress fields. Let us now revisit this argument and exhibit the link with a material made of isotropically distributed uniaxial reinforcements.

Let us indeed consider that the reinforcing material is made of a distribution of uniaxial reinforcements belonging to a certain family $\mathcal{A}$ of orientations $\alpha$ and of similar effective tensile/compressive strengths $f_t^{r,eff}$, $f_c^{r,eff}$. In order to enforce that only one orientation is active at a given material point, we can write the following strength condition:

$$ \exists \sigma^{r,\alpha}, \zeta_{\alpha} \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma^r = \sum_{\alpha \in \mathcal{A}} \zeta_{\alpha} \sigma^{r,\alpha} e_{\alpha} \otimes e_{\alpha} \\ -f_c^{r,eff} \le \sigma^{r,\alpha} \le f_t^{r,eff} \quad \forall \alpha \in \mathcal{A} \\ \zeta_{\alpha} \in \{0; 1\} \quad \forall \alpha \in \mathcal{A} \\ \sum_{\alpha \in \mathcal{A}} \zeta_{\alpha} = 1 \end{array} \right. \qquad (9) $$
### 3.3 The particular case of uniaxial reinforcements

When considering uniaxial reinforcements aligned along a direction $e_\alpha$ in small volume fraction, the composite homogenized strength criterion is in fact exactly given by the dilute estimation Equation (3) (de Buhan and Taliercio 1991; de Buhan et al. 2017) which reads here

$$ \sigma \in G^{\text{comp}} \Leftrightarrow \exists \sigma^m, \sigma^r \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma = \sigma^m + \sigma^r e_\alpha \otimes e_\alpha \\ \sigma^m \in G^m \\ -f_c^{r,\text{eff}} \le \sigma^r \le f_t^{r,\text{eff}} \end{array} \right. \qquad (5) $$

where $f_t^{r,\text{eff}}$ (resp. $f_c^{r,\text{eff}}$) denotes the reinforcement uniaxial effective tensile (resp. compressive) strength (per unit of transverse area). As a result, the case of uniaxial reinforcements is similar to Equation (3) except that the effective reinforcement strength criterion $G^{r,\text{eff}}$ is now given by the corresponding uniaxial strength condition

$$ G^{r,\text{eff}} = \{ \sigma^r e_\alpha \otimes e_\alpha \text{ s.t. } -f_c^{r,\text{eff}} \le \sigma^r \le f_t^{r,\text{eff}} \}. \qquad (6) $$

Finally, the homogenized strength criterion Equation (5) and Equation (6) easily generalize to a reinforcing material made of multiple reinforcement directions by summing the corresponding uniaxial stress contributions. For instance, an important practical case of interest is that of orthogonal reinforcements aligned with the global $x, y$ directions (and possibly $z$ in 3D). In this case, the effective strength criterion Equation (6) for two reinforcement directions $e_x$ and $e_y$ generalizes to

$$ G^{r,\text{eff}} = \{ \sigma^{r,x} e_x \otimes e_x + \sigma^{r,y} e_y \otimes e_y \text{ s.t. } -f_c^{r,\text{eff}} \le \sigma^{r,x}, \sigma^{r,y} \le f_t^{r,\text{eff}} \}. \qquad (7) $$

### 3.4 The $L_1$-Rankine criterion for isotropically distributed reinforcements

Similarly to the concepts of the homogenization method in topology optimization, our ultimate goal is to find an optimal microstructure for the reinforcing phase at each material point. To do so, instead of considering that the reinforcing material is made of a fixed distribution of predefined orientations, we can consider a reinforcing material consisting of uniaxial reinforcements but with locally unknown orientation a priori, with the goal that the topology optimization process would naturally select the locally optimal orientation. in (Mourad et al. 2021), the use of a $L_1$-Rankine was proposed in order to promote uniaxial stress fields. Let us now revisit this argument and exhibit the link with a material made of isotropically distributed uniaxial reinforcements.

Let us indeed consider that the reinforcing material is made of a distribution of uniaxial reinforcements belonging to a certain family $\mathcal{A}$ of orientations $\alpha$ and of similar effective tensile/compressive strengths $f_t^{r,\text{eff}}$, $f_c^{r,\text{eff}}$. In order to enforce that only one orientation is active at a given material point, we can write the strength condition

$$ \exists \sigma^{r,\alpha}, \zeta_{\alpha} \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma^r = \sum_{\alpha \in \mathcal{A}} \zeta_{\alpha} \sigma^{r,\alpha} e_{\alpha} \otimes e_{\alpha} \\ -f_c^{r,\text{eff}} \le \sigma^{r,\alpha} \le f_t^{r,\text{eff}} \quad \forall \alpha \in \mathcal{A} \\ \zeta_{\alpha} \in \{0; 1\} \quad \forall \alpha \in \mathcal{A} \\ \sum_{\alpha \in \mathcal{A}} \zeta_{\alpha} = 1 \end{array} \right. \qquad (8) $$

where we introduced the binary variables $\zeta_\alpha$ which describe the activation or not of a specific orientation, the last constraint enforcing that only one such orientation can be active.

Clearly, criterion Equation (8) is non-convex due to the binary constraint on the $\zeta_\alpha$ variable which will result in the corresponding topology optimization problem being extremely difficult to solve. To alleviate this issue and follow the same kind of convex optimization methodology of our approach, a natural idea is to convexify $G^{r,\text{eff}}$ by relaxing the binary constraint. We therefore
where we introduced the binary variables $\zeta_\alpha$ which describe the activation or not of a specific orientation, the last constraint enforcing that only one such orientation can be active.

Clearly, criterion Equation (9) is non-convex due to the binary constraint on the $\zeta_\alpha$ variable which will result in the corresponding topology optimization problem being extremely difficult to solve. To alleviate this issue and follow the same kind of convex optimization methodology of our approach, a natural idea is to convexify $G^{\text{r,eff}}$ by relaxing the binary constraint. We therefore consider the following convexified formulation:

$$ \sigma^r \in G^{\text{r,eff}} \Leftrightarrow \exists \sigma^{r,\alpha}, \zeta_\alpha \text{ s.t. } \left\{ \begin{array}{ll} \sigma^r = \sum_{\alpha \in \mathcal{A}} \zeta_\alpha \sigma^{r,\alpha} \mathbf{e}_\alpha \otimes \mathbf{e}_\alpha & \\ -f_c^{\text{r,eff}} \le \sigma^{r,\alpha} \le f_t^{\text{r,eff}} & \forall \alpha \in \mathcal{A} \\ 0 \le \zeta_\alpha \le 1 & \forall \alpha \in \mathcal{A} \\ \sum_{\alpha \in \mathcal{A}} \zeta_\alpha = 1 & \end{array} \right. \tag{10} $$

which we recognize as the definition of the convex hull of the individual uniaxial strength conditions $G^\alpha = \{\sigma^{r,\alpha} \mathbf{e}_\alpha \otimes \mathbf{e}_\alpha \text{ s.t. } -f_c^{\text{r,eff}} \le \sigma^{r,\alpha} \le f_t^{\text{r,eff}}\}$ i.e.

$$ G^{\text{r,eff}} = \underset{\alpha \in \mathcal{A}}{\text{conv}}\{G^\alpha\} \tag{11} $$

which is indeed the tightest convexification of the union of all the $G^\alpha$. Finally, in the case where $\mathcal{A}$ spans all the possible directions in space, we can easily show that $G^{\text{r,eff}}$ is in fact equal to the $L_1$-Rankine criterion with tensile (resp. compressive) strength $f_t^{\text{r,eff}}$ (resp. $f_c^{\text{r,eff}}$) introduced in Mourad et al. (2021), see Appendix A.1 for the proof. This result justifies that the $L_1$-Rankine is the tightest convex criterion promoting uniaxial stress states in an isotropic fashion.

In the case of orthogonal reinforcements where $\alpha$ can only be aligned with either the $x$ or $y$ directions, the above convexified formulation Equation (10) reads:

$$ \sigma^r \in G^{\text{r,eff}} \Leftrightarrow \exists \sigma^{r,x}, \sigma^{r,y}, \zeta_x, \zeta_y \text{ s.t. } \left\{ \begin{array}{ll} \sigma^r = \zeta_x \sigma^{r,x} \mathbf{e}_x \otimes \mathbf{e}_x + \zeta_y \sigma^{r,y} \mathbf{e}_y \otimes \mathbf{e}_y & \\ -f_c^{\text{r,eff}} \le \sigma^{r,x}, \sigma^{r,y} \le f_t^{\text{r,eff}} & \\ 0 \le \zeta_x, \zeta_y \le 1 & \\ \zeta_x + \zeta_y = 1 & \end{array} \right. \tag{12} $$

Interestingly, the above criterion corresponds exactly to the $L_1$-Rankine criterion intersected with the plane $\sigma_{xy}^r = 0$.

The main interest of the above construction in the case of a fixed family of discrete orientations $\mathcal{A} = \{\alpha_1, \dots, \alpha_N\}$ is that the resulting strength condition will be anisotropic with larger strength in the corresponding directions. We therefore expect the strength-based topology optimization procedure to naturally result in designs locally oriented in one of these directions.

By way of illustration, Figure 2 displays the uniaxial tensile strength in direction $\mathbf{e}_\theta = \cos \theta \mathbf{e}_x + \sin \theta \mathbf{e}_y$ for a material consisting of such a family of discrete orientations. As expected, when the reinforcement material is made of only two reinforcement directions, the material possesses no shear strength so that the uniaxial strength is always zero except if $\theta$ is perfectly aligned with one of the two directions. For more than two directions, the material possesses a shear strength and, therefore, a non-zero tensile strength for any $\theta$. Again, we observe that the uniaxial strength is equal to $f_t$ when the loading direction is aligned with one of the reinforcement direction and is less than $f_t$ in-between. The resulting material therefore possesses anisotropic strength properties. In the limit of an isotropic continuous distribution of reinforcement direction, the uniaxial strength becomes a constant equal to $f_t$ since the resulting material strength properties are equivalent to an isotropic $L_1$-Rankine strength criterion.
consider the convexified formulation

$$ \sigma^r \in G^{r,\text{eff}} \Leftrightarrow \exists \sigma^{r,\alpha}, \zeta_\alpha \text{ s.t. } \left\{ \begin{array}{ll} \sigma^r = \sum_{\alpha \in \mathcal{A}} \zeta_\alpha \sigma^{r,\alpha} \mathbf{e}_\alpha \otimes \mathbf{e}_\alpha & \\ -f_c^{r,\text{eff}} \le \sigma^{r,\alpha} \le f_t^{r,\text{eff}} & \forall \alpha \in \mathcal{A} \\ 0 \le \zeta_\alpha \le 1 & \forall \alpha \in \mathcal{A} \\ \sum_{\alpha \in \mathcal{A}} \zeta_\alpha = 1 & \end{array} \right. \tag{9} $$

which we recognize as the definition of the convex hull of the individual uniaxial strength conditions $G^\alpha = \{\sigma^{r,\alpha} \mathbf{e}_\alpha \otimes \mathbf{e}_\alpha \text{ s.t. } -f_c^{r,\text{eff}} \le \sigma^{r,\alpha} \le f_t^{r,\text{eff}}\}$ i.e.

$$ G^{r,\text{eff}} = \text{conv}_{\alpha \in \mathcal{A}}\{G^\alpha\} \tag{10} $$

which is indeed the tightest convexification of the union of all the $G^\alpha$. Finally, in the case where $\mathcal{A}$ spans all the possible directions in space, we can easily show that $G^{r,\text{eff}}$ is in fact equal to the $L_1$-Rankine criterion with tensile (resp. compressive) strength $f_t^{r,\text{eff}}$ (resp. $f_c^{r,\text{eff}}$) introduced in (Mourad et al. 2021), see Appendix A for the proof. This result justifies that the $L_1$-Rankine is the tightest convex criterion promoting uniaxial stress states in an isotropic fashion.

In the case of orthogonal reinforcements where $\alpha$ can only be aligned with either the $x$ or $y$ directions, Equation (9) reads

$$ \sigma^r \in G^{r,\text{eff}} \Leftrightarrow \exists \sigma^{r,x}, \sigma^{r,y}, \zeta_x, \zeta_y \text{ s.t. } \left\{ \begin{array}{ll} \sigma^r = \zeta_x \sigma^{r,x} \mathbf{e}_x \otimes \mathbf{e}_x + \zeta_y \sigma^{r,y} \mathbf{e}_y \otimes \mathbf{e}_y & \\ -f_c^{r,\text{eff}} \le \sigma^{r,x}, \sigma^{r,y} \le f_t^{r,\text{eff}} & \\ 0 \le \zeta_x, \zeta_y \le 1 & \\ \zeta_x + \zeta_y = 1. & \end{array} \right. \tag{11} $$

Interestingly, the above criterion corresponds exactly to the $L_1$-Rankine criterion intersected with the plane $\sigma_{xy}^r = 0$.

The main interest of the above construction in the case of a fixed family of discrete orientations $\mathcal{A} = \{\alpha_1, \dots, \alpha_N\}$ is that the resulting strength condition will be anisotropic with larger strength in the corresponding directions. We therefore expect the strength-based topology optimization procedure to naturally result in designs locally oriented in one of these directions.

By way of illustration, Figure 2 displays the uniaxial tensile strength in direction $\mathbf{e}_\theta = \cos \theta \mathbf{e}_x + \sin \theta \mathbf{e}_y$ for a material consisting of such a family of discrete orientations. As expected, when the reinforcement material is made of only two reinforcement directions, the material possesses no shear strength so that the uniaxial strength is always zero except if $\theta$ is perfectly aligned with one of the two directions. For more than two directions, the material possesses a shear strength and, therefore, a non-zero tensile strength for any $\theta$. Again, we observe that the uniaxial strength is equal to $f_t$ when the loading direction is aligned with one of the reinforcement direction and is less than $f_t$ in-between. The resulting material therefore possesses anisotropic strength properties. In the limit of an isotropic continuous distribution of reinforcement direction, the uniaxial strength becomes a constant equal to $f_t$ since the resulting material strength properties are equivalent to an isotropic $L_1$-Rankine strength criterion.

# 4 Bi-material optimization

In this section, we investigate the concurrent optimization of multiple materials in addition to a void phase (see Figure 1-right). For simplicity, we restrict here to the case of two materials only, although the proposed procedure can be easily generalized to $n$ different materials. We denote the two materials by their phase index $i = 1, 2$ whereas void is associated with $i = 0$. We aim at enforcing that a given point $\mathbf{x}$ belongs to either phase 1, phase 2 or to the void. In terms of strength conditions, we would therefore have $\sigma \in G^1$, $\sigma \in G^2$ or $\sigma = 0$ which can be written as
<figure><img src="image_3.png" /><figcaption>Figure 2: Anisotropic uniaxial tensile strength for a material reinforced by a family of reinforcements of discrete orientations</figcaption></figure>

## 4 Bi-material optimization

In this section, we investigate the concurrent optimization of multiple materials in addition to a void phase (see Figure 1-right). For simplicity, we restrict here to the case of two materials only, although the proposed procedure can be easily generalized to n different materials. We denote the two materials by their phase index i = 1, 2 whereas void is associated with i = 0. We aim at enforcing that a given point x belongs to either phase 1, phase 2 or to the void. In terms of strength conditions, we would therefore have $\sigma \in G^1$, $\sigma \in G^2$ or $\sigma = 0$ which can be written as follows:

$$ \exists \tilde{\sigma}^0, \tilde{\sigma}^1, \tilde{\sigma}^2, \zeta_0, \zeta_1, \zeta_2 \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma = \zeta_0 \tilde{\sigma}^0 + \zeta_1 \tilde{\sigma}^1 + \zeta_2 \tilde{\sigma}^2 \\ \tilde{\sigma}^0 = 0 \\ \tilde{\sigma}^1 \in G^1 \\ \tilde{\sigma}^2 \in G^2 \\ \zeta_0 + \zeta_1 + \zeta_2 = 1 \\ \zeta_0, \zeta_1, \zeta_2 \in \{0, 1\} \end{array} \right. \qquad (13) $$

where the binary variables $\zeta_i$ indicate the membership to the corresponding phase, the constraint $\sum \zeta_i = 1$ enforcing that one and only one of the $\zeta_i = 1$ while the others are zero.

Obviously, due to the binary constraints, the above strength condition is not convex which will result in the corresponding topology optimization problem being extremely difficult to solve. To alleviate this issue and following the same ideas as in Section 3.4, we convexify the above condition with its tightest convex relaxation. To do so, we allow each $\zeta_i$ to take continuous values inside [0; 1] instead of being binary. To make a clear distinction, we will replace each $\zeta_i$ with $\rho_i$, interpreting these variables as the pseudo-density fields of topology optimization. Hence,
<figure><img src="image_4.png" /><figcaption>**Figure 2** Anisotropic uniaxial tensile strength for a material reinforced by a family of reinforcements of discrete orientations.</figcaption></figure>

follows:

$$ \exists \tilde{\sigma}^0, \tilde{\sigma}^1, \tilde{\sigma}^2, \zeta_0, \zeta_1, \zeta_2 \text{ s.t.} \left\{ \begin{array}{l} \sigma = \zeta_0 \tilde{\sigma}^0 + \zeta_1 \tilde{\sigma}^1 + \zeta_2 \tilde{\sigma}^2 \\ \tilde{\sigma}^0 = 0 \\ \tilde{\sigma}^1 \in G^1 \\ \tilde{\sigma}^2 \in G^2 \\ \zeta_0 + \zeta_1 + \zeta_2 = 1 \\ \zeta_0, \zeta_1, \zeta_2 \in \{0, 1\} \end{array} \right. \qquad (12) $$

where the binary variables $\zeta_i$ indicate the membership to the corresponding phase, the constraint $\sum \zeta_i = 1$ enforcing that one and only one of the $\zeta_i = 1$ while the others are zero.

Obviously, due to the binary constraints, the above strength condition is not convex which will result in the corresponding topology optimization problem being extremely difficult to solve. To alleviate this issue and following the same ideas as in Section 3.4, we convexify the above condition with its tightest convex relaxation. To do so, we allow each $\zeta_i$ to take continuous values inside $[0; 1]$ instead of being binary. To make a clear distinction, we will replace each $\zeta_i$ with $\rho_i$, interpreting these variables as the pseudo-density fields of topology optimization. Hence, we consider

$$ \exists \tilde{\sigma}^1, \tilde{\sigma}^2, \rho_0, \rho_1, \rho_2 \text{ s.t.} \left\{ \begin{array}{l} \sigma = \rho_1 \tilde{\sigma}^1 + \rho_2 \tilde{\sigma}^2 \\ \tilde{\sigma}^1 \in G^1 \\ \tilde{\sigma}^2 \in G^2 \\ \rho_0 + \rho_1 + \rho_2 = 1 \\ \rho_0, \rho_1, \rho_2 \in [0; 1] \end{array} \right. \qquad (13) $$

$$ \iff \exists \tilde{\sigma}^1, \tilde{\sigma}^2, \rho_1, \rho_2 \text{ s.t.} \left\{ \begin{array}{l} \sigma = \rho_1 \tilde{\sigma}^1 + \rho_2 \tilde{\sigma}^2 \\ \tilde{\sigma}^1 \in G^1 \\ \tilde{\sigma}^2 \in G^2 \\ \rho_1 + \rho_2 \le 1 \\ \rho_1, \rho_2 \in [0; 1] \end{array} \right. \qquad (14) $$

where we removed the void density $\rho_0$.

This motivates the introduction of a density-dependent strength condition $G(\rho_1, \rho_2)$, defined
we consider:

$$ \exists \tilde{\sigma}^1, \tilde{\sigma}^2, \rho_0, \rho_1, \rho_2 \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma = \rho_1 \tilde{\sigma}^1 + \rho_2 \tilde{\sigma}^2 \\ \tilde{\sigma}^1 \in G^1 \\ \tilde{\sigma}^2 \in G^2 \\ \rho_0 + \rho_1 + \rho_2 = 1 \\ \rho_0, \rho_1, \rho_2 \in [0; 1] \end{array} \right. \qquad (14) $$

$$ \iff \exists \tilde{\sigma}^1, \tilde{\sigma}^2, \rho_1, \rho_2 \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma = \rho_1 \tilde{\sigma}^1 + \rho_2 \tilde{\sigma}^2 \\ \tilde{\sigma}^1 \in G^1 \\ \tilde{\sigma}^2 \in G^2 \\ \rho_1 + \rho_2 \le 1 \\ \rho_1, \rho_2 \in [0; 1] \end{array} \right. \qquad (15) $$

where we removed the void density $\rho_0$.

This motivates the introduction of a density-dependent strength condition $G(\rho_1, \rho_2)$, defined by the following strength condition:

$$ \sigma \in G(\rho_1, \rho_2) \iff \exists \sigma^1, \sigma^2 \text{ s.t.} \quad \left\{ \begin{array}{l} \sigma = \sigma^1 + \sigma^2 \\ \sigma^1 \in \rho_1 G^1 \\ \sigma^2 \in \rho_2 G^2 \end{array} \right. \qquad (16) $$

in which we made the change of variable $\sigma^i = \rho_i \tilde{\sigma}^i$.

In particular, if both $\rho_1(x) = \rho_2(x) = 0$ at a given point $x$, we have $\sigma(x) = 0$ i.e. $x$ is in a void phase. If $\rho_1(x) = 1$, then $\rho_2(x) = 0$ and $\sigma(x) \in G^1$ i.e. $x$ belongs to material 1 and vice versa. Note that it is possible to find states where $\rho_1(x) \neq 0$ and $\rho_2(x) \neq 0$ which results in $x$ belonging to a fictitious material averaging the strength properties of both phases.

Finally, in the case where $G^1 = G^2 = G$, we have $G(\rho_1, \rho_2) = (\rho_1 + \rho_2)G = \rho G$ and we recover the single material formulation of Mourad et al. (2021). Moreover, let us consider phase 1 to be the matrix phase. If we fix $\rho_1 = 1$ everywhere and if we ignore the constraint $\rho_1 + \rho_2 \le 1$, we recover the reinforcement formulation of Section 3 where phase 2 corresponds to the reinforcement phase. In this context and contrary to the bi-material formulation, we see that we do not optimize the matrix phase density and that we do not consider any void phase, we only optimize for the reinforcement density.

## 4.1 Convex problem formulation

The two corresponding load-maximization Equation (LOAD-MAX) and volume minimization Equation (VOL-MIN) problems are therefore respectively given by:

$$ \begin{array}{lll} \lambda^+ = \max_{\lambda, \sigma^1, \sigma^2, \rho_1, \rho_2} & \lambda & \\ \text{s.t.} & \text{div}(\sigma^1 + \sigma^2) = 0 & \text{in } \mathcal{D} \\ & (\sigma^1 + \sigma^2) \cdot n = \lambda T & \text{on } \partial \mathcal{D}_T \\ & \sigma^1 \in \rho_1 G^1 & \text{in } \mathcal{D} \\ & \sigma^2 \in \rho_2 G^2 & \text{in } \mathcal{D} \\ & \int_{\mathcal{D}} c_{\omega}(\rho_1, \rho_2) \, dx \le \eta |\mathcal{D}| & \text{(BIMAT-LOAD-MAX)} \\ & 0 \le \rho_1 \le 1 \\ & 0 \le \rho_2 \le 1 \\ & \rho_1 + \rho_2 \le 1 \end{array} $$
by the following strength condition:

$$ \sigma \in G(\rho_1, \rho_2) \Leftrightarrow \exists \sigma^1, \sigma^2 \text{ s.t. } \begin{cases} \sigma = \sigma^1 + \sigma^2 \\ \sigma^1 \in \rho_1 G^1 \\ \sigma^2 \in \rho_2 G^2 \end{cases} \qquad (15) $$

in which we made the change of variable $\sigma^i = \rho_i \tilde{\sigma}^i$.

In particular, if both $\rho_1(\mathbf{x}) = \rho_2(\mathbf{x}) = 0$ at a given point $\mathbf{x}$, we have $\sigma(\mathbf{x}) = 0$ i.e. $\mathbf{x}$ is in a void phase. If $\rho_1(\mathbf{x}) = 1$, then $\rho_2(\mathbf{x}) = 0$ and $\sigma(\mathbf{x}) \in G^1$ i.e. $\mathbf{x}$ belongs to material 1 and vice versa. Note that it is possible to find states where $\rho_1(\mathbf{x}) \neq 0$ and $\rho_2(\mathbf{x}) \neq 0$ which results in $\mathbf{x}$ belonging to a fictitious material averaging the strength properties of both phases.

Finally, in the case where $G^1 = G^2 = G$, we have $G(\rho_1, \rho_2) = (\rho_1 + \rho_2)G = \rho G$ and we recover the single material formulation of Mourad et al. (2021). Moreover, let us consider phase 1 to be the matrix phase. If we fix $\rho_1 = 1$ everywhere and if we ignore the constraint $\rho_1 + \rho_2 \le 1$, we recover the reinforcement formulation of Section 3 where phase 2 corresponds to the reinforcement phase. In this context and contrary to the bi-material formulation, we see that we do not optimize the matrix phase density and that we do not consider any void phase, we only optimize for the reinforcement density.

## 4.1 Convex problem formulation

The two corresponding load-maximization (LOAD-MAX) and volume minimization (VOL-MIN) problems are therefore respectively given by

$$ \begin{array}{ll} \lambda^+ = \max_{\lambda, \sigma^1, \sigma^2, \rho_1, \rho_2} & \lambda \\ \text{s.t.} & \text{div}(\sigma^1 + \sigma^2) = 0 \quad \text{in } \mathcal{D} \\ & (\sigma^1 + \sigma^2) \cdot \mathbf{n} = \lambda T \quad \text{on } \partial \mathcal{D}_T \\ & \sigma^1 \in \rho_1 G^1 \quad \text{in } \mathcal{D} \\ & \sigma^2 \in \rho_2 G^2 \quad \text{in } \mathcal{D} \qquad (\text{BIMAT-LOAD-MAX}) \\ & \int_{\mathcal{D}} c_{\omega}(\rho_1, \rho_2) \, dx \le \eta |\mathcal{D}| \\ & 0 \le \rho_1 \le 1 \\ & 0 \le \rho_2 \le 1 \\ & \rho_1 + \rho_2 \le 1 \end{array} $$

and

$$ \begin{array}{ll} \eta^- = \min_{\sigma^1, \sigma^2, \rho_1, \rho_2} & \frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} c_{\omega}(\rho_1, \rho_2) \, dx \\ \text{s.t.} & \text{div}(\sigma^1 + \sigma^2) = 0 \quad \text{in } \mathcal{D} \\ & (\sigma^1 + \sigma^2) \cdot \mathbf{n} = \lambda T \quad \text{on } \partial \mathcal{D}_T \\ & \sigma^1 \in \rho_1 G^1 \quad \text{in } \mathcal{D} \qquad (\text{BIMAT-VOL-MIN}) \\ & \sigma^2 \in \rho_2 G^2 \quad \text{in } \mathcal{D} \\ & 0 \le \rho_1 \le 1 \\ & 0 \le \rho_2 \le 1 \\ & \rho_1 + \rho_2 \le 1 \end{array} $$

where we introduced $c_{\omega}(\rho_1, \rho_2) = 2\omega\rho_1 + 2(1 - \omega)\rho_2$ which is a weighted-average cost function measuring the amount of both materials. As it will be seen later, the introduction of the weighting factor $\omega \in [0; 1]$ gives us increased flexibility in obtaining various optimal design depending on the cost associated with the presence of material 1 over material 2. Note that the above choice gives $c_{1/2}(\rho_1, \rho_2) = \rho_1 + \rho_2 = \rho$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Splitting of a nominal strength criterion G (in black) into a purely compressive part G⁻ (in blue) and a purely tensile part G⁺ (in red) and the corresponding convex hull conv{G⁺, G⁻} (in green). Left: a Rankine strength criterion, right: a L₁-Rankine strength criterion in the plane of principal stresses</figcaption></figure>

and

$$
\begin{align*}
\eta^{-} &= \min_{\sigma^{1}, \sigma^{2}, \rho_{1}, \rho_{2}} \frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} c_{\omega}(\rho_{1}, \rho_{2}) \, dx \\
\text{s.t.} \quad & \operatorname{div}(\sigma^{1} + \sigma^{2}) = 0 \quad \text{in } \mathcal{D} \\
& (\sigma^{1} + \sigma^{2}) \cdot \mathbf{n} = \lambda T \quad \text{on } \partial\mathcal{D}_{T} \\
& \sigma^{1} \in \rho_{1}G^{1} \quad \text{in } \mathcal{D} \qquad (\text{BIMAT-VOL-MIN}) \\
& \sigma^{2} \in \rho_{2}G^{2} \quad \text{in } \mathcal{D} \\
& 0 \le \rho_{1} \le 1 \\
& 0 \le \rho_{2} \le 1 \\
& \rho_{1} + \rho_{2} \le 1
\end{align*}
$$

where we introduced $c_{\omega}(\rho_{1}, \rho_{2}) = 2\omega\rho_{1} + 2(1 - \omega)\rho_{2}$ which is a weighted-average cost function measuring the amount of both materials. As it will be seen later, the introduction of the weighting factor $\omega \in [0; 1]$ gives us increased flexibility in obtaining various optimal design depending on the cost associated with the presence of material 1 over material 2. Note that the above choice gives $c_{1/2}(\rho_{1}, \rho_{2}) = \rho_{1} + \rho_{2} = \rho$.

## 4.2 No-tension and no-compression materials

An important case of application of the previous bi-material formulation is concerned with the optimization of a no-tension and a no-compression phase. Practically, this could correspond to two different materials respectively possessing negligible tensile strength (e.g. concrete, rocks, masonry, etc.) and negligible compressive strength (e.g. thin membrane which would buckle under compression). Another possibility is to consider a single material for which we would like to distinguish members in tension from members in compression in the optimization process, for example in order to assign a different cost between the tensile and compressive "phase".

As regards this last point of view for a single material of nominal strength properties G, one could define the no-tension strength criterion $G^{1} = G^{-} = G \cap S^{-}$ and the no-compression strength criterion $G^{2} = G^{+} = G \cap S^{+}$ where:

$$S^{\pm} = \{\sigma \text{ s.t. } \pm \sigma \geq 0\} \qquad (17)$$

represents the cone of symmetric positive/negative stress tensors. In this case, since $G^{\pm} \subset G$ and G is convex, we have that conv{G⁺, G⁻} ⊆ G. Again, this formulation will tend to promote stress states either in pure tension or in pure compression. Figure 3 illustrates this construction in the case of a Rankine and L₁-Rankine criterion. Note that we have that conv{G⁺, G⁻} = G<sup>L₁-Rankine</sup> in this latter case.

Finally, as already discussed, the use of a L₁-Rankine strength criterion will even further promote uniaxial stress states. If the original material is isotropic and possesses a characteristic
<figure><img src="image_5.png" /><figcaption>**Figure 3** Splitting of a nominal strength criterion G (in black) into a purely compressive part G⁻ (in blue) and a purely tensile part G⁺ (in red) and the corresponding convex hull conv{G⁺, G⁻} (in green). Left: a Rankine strength criterion, right: a L₁-Rankine strength criterion in the plane of principal stresses.</figcaption></figure>

## 4.2 No-tension and no-compression materials

An important case of application of the previous bi-material formulation is concerned with the optimization of a no-tension and a no-compression phase. Practically, this could correspond to two different materials respectively possessing negligible tensile strength (e.g. concrete, rocks, masonry, etc.) and negligible compressive strength (e.g. thin membrane which would buckle under compression). Another possibility is to consider a single material for which we would like to distinguish members in tension from members in compression in the optimization process, for example in order to assign a different cost between the tensile and compressive "phase".

As regards this last point of view for a single material of nominal strength properties G, one could define the no-tension strength criterion $G^1 = G^- = G \cap S^-$ and the no-compression strength criterion $G^2 = G^+ = G \cap S^+$ where $S^\pm = \{\sigma \text{ s.t. } \pm \sigma \ge 0\}$ represents the cone of symmetric positive/negative stress tensors. In this case, since $G^\pm \subset G$ and G is convex, we have that conv{G⁺, G⁻} ⊆ G. Again, this formulation will tend to promote stress states either in pure tension or in pure compression. Figure 3 illustrates this construction in the case of a Rankine and L₁-Rankine criterion. Note that we have that conv{G⁺, G⁻} = $G^{L_1-Rankine}$ in this latter case.

Finally, as already discussed, the use of a L₁-Rankine strength criterion will even further promote uniaxial stress states. If the original material is isotropic and possesses a characteristic tensile strength $f_t$ and compressive strength $f_c$, a natural modeling strategy for obtaining truss-like designs when distinguishing the optimization of tensile and compressive members is therefore to consider

$$G^1 = G^{L_1-Rankine}(f_c, 0) \qquad (16)$$

$$G^2 = G^{L_1-Rankine}(0, f_t) \qquad (17)$$

where $G^{L_1-Rankine}(f_c, f_t)$ denotes the isotropic L₁-Rankine strength criterion of compressive (resp. tensile) strength $f_c$ (resp. $f_t$).

Finally, the above choice can also be adapted to the situations for which members are constrained to have fixed potential orientations by using criteria of the form Equation (9) instead of the isotropic L₁-Rankine criterion.

# 5 Illustrative applications

## 5.1 Numerical implementation and penalization procedure

Similarly to (Mourad et al. 2021), the corresponding discrete optimization problems are formulated using the fenics_optim package (Bleyer 2020a; Bleyer 2020b) which enables to couple the FEniCS finite-element software package (Alnæs et al. 2015; Logg et al. 2012) with the Mosek conic optimization solver (MOSEK ApS 2018). The equilibrium constraint is enforced weakly through the virtual work principle using a continuous P²-Lagrange (quadratic polynomial) interpolation for virtual displacement fields. The pseudo density fields are discretized using a continuous
tensile strength $f_t$ and compressive strength $f_c$, a natural modeling strategy for obtaining truss-like designs when distinguishing the optimization of tensile and compressive members is therefore to consider:

$$G^1 = G^{L_1\text{-Rankine}}(f_c, 0) \qquad (18)$$

$$G^2 = G^{L_1\text{-Rankine}}(0, f_t) \qquad (19)$$

where $G^{L_1\text{-Rankine}}(f_c, f_t)$ denotes the isotropic $L_1$-Rankine strength criterion of compressive (resp. tensile) strength $f_c$ (resp. $f_t$).

Finally, the above choice can also be adapted to the situations for which members are constrained to have fixed potential orientations by using criteria of the form Equation (10) instead of the isotropic $L_1$-Rankine criterion.

# 5 Illustrative applications

## 5.1 Numerical implementation and penalization procedure

Similarly to Mourad et al. (2021), the corresponding discrete optimization problems are formulated using the fenics_optim package (Bleyer 2020a; Bleyer 2020b) which enables to couple the FEniCS finite-element software package (Alnæs et al. 2015; Logg et al. 2012) with the Mosek conic optimization solver (MOSEK ApS 2018). The equilibrium constraint is enforced weakly through the virtual work principle using a continuous $\mathbb{P}^2$-Lagrange (quadratic polynomial) interpolation for virtual displacement fields. The pseudo density fields are discretized using a continuous $\mathbb{P}^1$-Lagrange (linear polynomial) interpolation, see Mourad et al. (2021) for more details. For all examples, we use unstructured meshes consisting of affine triangles. Mesh dependency issues are removed by adding a slope-control constraint (Petersson et al. 1998) for $\|\nabla\rho_i\|_2 \le 1/\ell$ for all density fields. This slope control avoids localization of thin members on bands of a single element width by controlling the gradient of the density field. Doing so, member width must be at least of thickness $2\ell$ which is chosen to be larger a few times than the mesh size. With this procedure, we also avoid checkerboard issues which can occur in topology optimization without spatial regularization.

Finally, we also extend the continuation procedure penalizing intermediate densities proposed in Mourad et al. (2021) to the present case. More precisely, each convex strength constraint of the form $\sigma^i = \rho_i G^i$ is replaced by a penalized power-law (non-convex) constraint $\sigma^i = (\rho_i)^p G^i$ following ideas of the SIMP method (Bendsøe et al. 2004). At each iteration of the penalization procedure, the power-law is linearized around the current density estimate $\rho_{i,n}$ and the exponent is progressively increased from 1 to a maximum value $p_{\max} > 1$, see again Mourad et al. (2021) for more details.

## 5.2 Numerical examples objectives

In the following, we will investigate three different examples which will have the common goal of assessing the versatility of the proposed methodology. Each of them will analyze some specific features which can be considered at the modeling stage, namely:

* MBB beam example:

- Analyze the "reinforcement optimization" formulation and assess the influence of the reinforcement strength criterion choice (in terms of overall shape or anisotropy) on the resulting design;

- Analyze the "bi-material optimization" formulation and assess the influence of the strength criterion of one of the phases on the resulting design;

- Assess the corresponding formulation in which one phase only sustains tension and the other phase only compression;

- Assess the corresponding formulation where the tensile phase is subjected to a prescribed anisotropy in terms of reinforcement directions and compare its efficiency against an isotropic distribution of reinforcements.

* Bridge example:

- Analyze the influence of asymmetric tensile and compressive strengths on the resulting design;
P¹-Lagrange (linear polynomial) interpolation, see (Mourad et al. 2021) for more details. For all examples, we use unstructured meshes consisting of affine triangles. Mesh dependency issues are removed by adding a slope-control constraint (Petersson and Sigmund 1998) for $||∇ρ_i||_2 \le 1/ℓ$ for all density fields. This slope control avoids localization of thin members on bands of a single element width by controlling the gradient of the density field. Doing so, member width must be at least of thickness $2ℓ$ which is chosen to be larger a few times than the mesh size. With this procedure, we also avoid checkerboard issues which can occur in topology optimization without spatial regularization.

Finally, we also extend the continuation procedure penalizing intermediate densities proposed in (Mourad et al. 2021) to the present case. More precisely, each convex strength constraint of the form $\boldsymbol{\sigma}^i = \rho_i G^i$ is replaced by a penalized power-law (non-convex) constraint $\boldsymbol{\sigma}^i = (\rho_i)^p G^i$ following ideas of the SIMP method (Bendsøe and Sigmund 2004). At each iteration of the penalization procedure, the power-law is linearized around the current density estimate $\rho_{i,n}$ and the exponent is progressively increased from 1 to a maximum value $p_{max} > 1$, see again (Mourad et al. 2021) for more details.

## 5.2 Numerical examples objectives

In the following, we will investigate three different examples which will have the common goal of assessing the versatility of the proposed methodology. Each of them will analyze some specific features which can be considered at the modeling stage, namely:

* MBB beam example:

- Analyze the “reinforcement optimization” formulation and assess the influence of the reinforcement strength criterion choice (in terms of overall shape or anisotropy) on the resulting design;

- Analyze the “bi-material optimization” formulation and assess the influence of the strength criterion of one of the phases on the resulting design;

- Assess the corresponding formulation in which one phase only sustains tension and the other phase only compression;

- Assess the corresponding formulation where the tensile phase is subjected to a prescribed anisotropy in terms of reinforcement directions and compare its efficiency against an isotropic distribution of reinforcements.

* Bridge example:

- Analyze the influence of asymmetric tensile and compressive strengths on the resulting design;

- Compare the obtained design with a topology optimization involving a single material (with asymmetric strength properties);

- Assess the influence of the cost function factor $\omega$ used to weight the cost of one phase with respect to another.

* Deep beam example:

- Consider the design of widely studied reinforced concrete structure with an opening and compare the result against a strut-and-tie model used in practice;

- Consider a variation in which the reinforcement phase can only be located in orthogonal directions;

- Compare the designs obtained with the reinforcement and bi-material formulations.

## 5.3 MBB beam

We first consider a MBB beam example (Figure 4) of length $l = 36$ cm and height $h = 6$ cm with simple supports on the left and roller supports on the right, a vertical force of reference intensity $P = 1$ is applied at the top. In the following, only one half of the model will be represented, taking symmetry into account. Both supports and force are distributed over a small distance $s = 0.5$ cm to mitigate stress concentrations. The mesh consists of approximately 40,000 elements.
<figure><img src="image_5.png" /><figcaption>Figure 4: A MBB beam example</figcaption></figure>

- Compare the obtained design with a topology optimization involving a single material (with asymmetric strength properties);

- Assess the influence of the cost function factor $\omega$ used to weight the cost of one phase with respect to another.

* Deep beam example:

- Consider the design of widely studied reinforced concrete structure with an opening and compare the result against a strut-and-tie model used in practice;

- Consider a variation in which the reinforcement phase can only be located in orthogonal directions;

- Compare the designs obtained with the reinforcement and bi-material formulations.

## 5.3 MBB beam

We first consider a MBB beam example (Figure 4) of length $l = 36$ cm and height $h = 6$ cm with simple supports on the left and roller supports on the right, a vertical force of reference intensity $P = 1$ is applied at the top. In the following, only one half of the model will be represented, taking symmetry into account. Both supports and force are distributed over a small distance $s = 0.5$ cm to mitigate stress concentrations. The mesh consists of approximately 40,000 elements.

### 5.3.1 Reinforcement optimization

We first consider the case of a fixed matrix material of strength condition $G^m$ being given by a Rankine criterion of compressive strength $f_c^m = 1$ and tensile strength $f_t^m = 0.05$. The reinforced material effective strength criterion $G^{r,eff}$ will either be:

* a plane stress von Mises criterion with uniaxial tension/compression effective strengths $f_c^{r,eff} = f_t^{r,eff} = 1$

* a $L_1$-Rankine criterion with the same uniaxial strengths

* a "no-compression" $L_1$-Rankine criterion with $f_c^{r,eff} = 0, f_t^{r,eff} = 1$

* an orthotropic $L_1$-Rankine criterion with $f_{tx}^{r,eff} = f_{ty}^{r,eff} = 1$ and $f_{cx}^{r,eff} = f_{cy}^{r,eff} = 0$

We consider the reinforcement load-maximization problem Equation (REINF-LOAD-MAX) with an imposed volume fraction $\eta = 0.2$.

Figure 5 displays the obtained optimized reinforcement density in black along with the principal compressive stress field in the matrix phase in blue. Owing to the fact that the matrix phase has a low tensile strength, reinforcements are primarily located on the bottom tensile face of the beam. The amount of reinforcement increases from the support to the beam mid-span where the bending moment is maximum. The choice of the reinforcement phase strength criterion mainly influences small details in the layout such as the precise location of the reinforcements. The use of a von Mises criterion (Figure 5(a)) for the reinforcement phase tends to favor biaxial stress states in the reinforcement region whereas the use of a $L_1$-Rankine criterion (Figure 5(b)) effectively promotes uniaxial stress states. Inclined reinforcements are also obtained to provide anchoring between the reinforcement tensile stress and the diffuse matrix compressive stresses. Moreover, when the reinforcement material possesses a non-zero compressive strength, we obtain a small amount of reinforcement on the top face at mid-span, thereby reinforcing the beam bending capacity. This region disappears in the case of a zero compressive strength (Figure 5(c)). Finally, when considering an orthotropic $L_1$-Rankine in Figure 5(d), we obtain a 90°-bend as
<figure><img src="image_6.png" /><figcaption>**Figure 4** A MBB beam example.</figcaption></figure>

### 5.3.1 Reinforcement optimization

We first consider the case of a fixed matrix material of strength condition $G^m$ being given by a Rankine criterion of compressive strength $f_c^m = 1$ and tensile strength $f_t^m = 0.05$. The reinforced material effective strength criterion $G^{r,eff}$ will either be:

* a plane stress von Mises criterion with uniaxial tension/compression effective strengths $f_c^{r,eff} = f_t^{r,eff} = 1$

* a $L_1$-Rankine criterion with the same uniaxial strengths

* a "no-compression" $L_1$-Rankine criterion with $f_c^{r,eff} = 0, f_t^{r,eff} = 1$

* an orthotropic $L_1$-Rankine criterion with $f_{tx}^{r,eff} = f_{ty}^{r,eff} = 1$ and $f_{cx}^{r,eff} = f_{cy}^{r,eff} = 0$

We consider the reinforcement load-maximization problem (RE INF-LOAD-MAX) with an imposed volume fraction $\eta = 0.2$.

Figure 5 displays the obtained optimized reinforcement density in black along with the principal compressive stress field in the matrix phase in blue. Owing to the fact that the matrix

<figure><img src="image_7.png" /><figcaption>**Figure 5** Reinforcement optimization of the MBB example for various reinforcement strength criteria: in black, reinforcement optimized density; in blue, principal compressive stress in the matrix phase.</figcaption></figure>

phase has a low tensile strength, reinforcements are primarily located on the bottom tensile face of the beam. The amount of reinforcement increases from the support to the beam mid-span where the bending moment is maximum. The choice of the reinforcement phase strength criterion mainly influences small details in the layout such as the precise location of the reinforcements. The use of a von Mises criterion (Figure 5(a)) for the reinforcement phase tends to favor biaxial stress states in the reinforcement region whereas the use of a $L_1$-Rankine criterion (Figure 5(b)) effectively promotes uniaxial stress states. Inclined reinforcements are also obtained to provide anchoring between the reinforcement tensile stress and the diffuse matrix compressive stresses. Moreover, when the reinforcement material possesses a non-zero compressive strength, we obtain a small amount of reinforcement on the top face at mid-span, thereby reinforcing the beam bending capacity. This region disappears in the case of a zero compressive strength (Figure 5(c)). Finally, when considering an orthotropic $L_1$-Rankine in Figure 5(d), we obtain a 90°-bend as commonly encountered in reinforced-concrete structures.
<figure><img src="image_6.png" /><figcaption>**Figure 5:** Reinforcement optimization of the MBB example for various reinforcement strength criteria: in black, reinforcement optimized density; in blue, principal compressive stress in the matrix phase.</figcaption></figure>

commonly encountered in reinforced-concrete structures.

### 5.3.2 Bi-material optimization

We now investigate the bimaterial load-maximization problem Equation (BIMAT-LOAD-MAX) with again $\eta = 0.2$ and $\omega = 1/2$.

As in Section 5.3.1, phase 1 corresponds to a Rankine criterion with $f_c = 1$, $f_t = 0.05$ and phase 2 with either a plane stress von Mises or a $L_1$-Rankine criterion of strengths $f_c = f_t = 1$. Note that we voluntarily used similar compressive strengths for both phases. Figure 6 represents the optimized two phase densities in both cases. We can see that phase 1, with low tensile strength, is mostly used where it is the most efficient i.e. in the top part of the structure subjected to compression. The second phase is used mainly in tensile regions, except for the loading region in Figure 6(a). This can be explained by the fact that this region is essentially under a biaxial compressive state, a stress state for which the von Mises criterion is larger than the Rankine criterion of phase 1. Let us also remark that the use of a $L_1$-Rankine criterion seems to yield a simpler design than that obtained with the von Mises material which can probably be attributed to the former promoting more efficiently uniaxial states.

Moreover, we can also observe in Figure 6 a transition between phase 1 and phase 2 materials along the inclined struts in compression on the right of the structure. This can be explained by the fact that we used the same uniaxial compressive strength $f_c$ for both phases. Both of them are therefore equally optimal for uniaxial compression. The precise location of the transition between both phases is probably dictated by the initial distribution of the phase density in the initial stages of the penalization procedure. In such stages, phase 1 material is essentially located on the top part of the beam and phase 2 on the bottom part. The penalization process being a continuation procedure, it will naturally converge to a uniaxial stress state corresponding to the phase which was active initially. Note that the use of a phase 1 material with a lower compressive strength, e.g. $f_c = 0.9$ would have resulted in a design involving only phase 2 material, since phase 1 would have been strictly weaker and, therefore, less efficient than phase 2. In this case, the formulation becomes equivalent to a single-material formulation. To conclude, this examples illustrates the efficiency of the convex relaxation procedure with respect to binary constraints. Indeed, since binary constraints are relaxed to continuous constraints and that both materials exhibit exactly the same strength properties in compression, we could have expected that both phases will be continuously mixed throughout the domain. However, this is not what we observe. Instead, the continuation penalization procedure succeeds in driving the design towards an almost binary distribution, non-binary values for $\rho_1$ and $\rho_2$ being only observed in a small transition region of order $\ell$.
### 5.3.2 Bi-material optimization

We now investigate the bimaterial load-maximization problem (BIMAT-LOAD-MAX) with again $η = 0.2$ and $ω = 1/2$.

As in Section 5.3.1, phase 1 corresponds to a Rankine criterion with $f_c = 1$, $f_t = 0.05$ and phase 2 with either a plane stress von Mises or a $L_1$-Rankine criterion of strengths $f_c = f_t = 1$. Note that we voluntarily used similar compressive strengths for both phases. Figure 6 represents the optimized two phase densities in both cases. We can see that phase 1, with low tensile strength, is mostly used where it is the most efficient i.e. in the top part of the structure subjected to compression. The second phase is used mainly in tensile regions, except for the loading region in Figure 6(a). This can be explained by the fact that this region is essentially under a biaxial compressive state, a stress state for which the von Mises criterion is larger than the Rankine criterion of phase 1. Let us also remark that the use of a $L_1$-Rankine criterion seems to yield a simpler design than that obtained with the von Mises material which can probably be attributed to the former promoting more efficiently uniaxial states.

Moreover, we can also observe in Figure 6 a transition between phase 1 and phase 2 materials along the inclined struts in compression on the right of the structure. This can be explained by the fact that we used the same uniaxial compressive strength $f_c$ for both phases. Both of them are therefore equally optimal for uniaxial compression. The precise location of the transition between both phases is probably dictated by the initial distribution of the phase density in the initial stages of the penalization procedure. In such stages, phase 1 material is essentially located on the top part of the beam and phase 2 on the bottom part. The penalization process being a continuation procedure, it will naturally converge to a uniaxial stress state corresponding to the phase which was active initially. Note that the use of a phase 1 material with a lower compressive strength, e.g. $f_c = 0.9$ would have resulted in a design involving only phase 2 material, since phase 1 would have been strictly weaker and, therefore, less efficient than phase 2. In this case, the formulation becomes equivalent to a single-material formulation. To conclude, this example illustrates the efficiency of the convex relaxation procedure with respect to binary constraints. Indeed, since binary constraints are relaxed to continuous constraints and that both materials exhibit exactly the same strength properties in compression, we could have expected that both phases will be continuously mixed throughout the domain. However, this is not what we observe. Instead, the continuation penalization procedure succeeds in driving the design towards an almost binary distribution, non-binary values for $ρ_1$ and $ρ_2$ being only observed in a small transition region of order $ℓ$.

<figure><img src="image_8.png" /><figcaption>**Figure 6** Bi-material load-maximization of the MBB example: in blue, phase 1 density; in red, phase 2 density.</figcaption></figure>

### 5.3.3 Splitting between tension and compression

We now investigate the formulation discussed in Section 4.2 where phase 1 (resp. phase 2) corresponds to a pure compression (resp. pure tension) phase of strength $f_c$ (resp. $f_t$). In practice, we use a $L_1$-Rankine strength criterion for both phases with a small residual tensile (resp. compression) strength in phase 1 (resp. phase 2) to avoid numerical instabilities. Figure 7 represents the corresponding optimized densities for both phases and various imposed volume fractions $η$. Comparing for instance Figure 7(c) with Figure 6(b), we can see that both final
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Bi-material load-maximization of the MBB example: in blue, phase 1 density; in red, phase 2 density</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Bi-material load-maximization of the MBB example with tension/compression splitting ($f_c = f_t = 1$) and various maximum volume fraction $\eta$</figcaption></figure>

### 5.3.3 Splitting between tension and compression

We now investigate the formulation discussed in Section 4.2 where phase 1 (resp. phase 2) corresponds to a pure compression (resp. pure tension) phase of strength $f_c$ (resp. $f_t$). In practice, we use a $L_1$-Rankine strength criterion for both phases with a small residual tensile (resp. compression) strength in phase 1 (resp. phase 2) to avoid numerical instabilities. Figure 7 represents the corresponding optimized densities for both phases and various imposed volume fractions $\eta$. Comparing for instance Figure 7(c) with Figure 6(b), we can see that both final designs are very similar but with a different repartition between both phases. In the present tension/compression splitting formulation Equation (19), each truss member belongs to a single phase, depending on its state of tension and compression whereas, in the previous formulation, some compression members could involve the phase 2 material or even both materials as already discussed.

### 5.3.4 Anisotropic strength properties for the tensile phase

We finish this example by considering again a tension/compression splitting formulation, except that the tensile phase (phase 2) now enjoys anisotropic strength properties of the form Equation (10) with here $f_c^{\text{r, eff}} = 0$ and $f_t^{\text{r, eff}} = 1$. We recall that if the set of allowed orientations $\mathcal{A}$ spans all directions, then the corresponding strength criterion is equivalent to a $L_1$-Rankine strength criterion so that we recover the results of Figure 7. In Figure 8, we report the results obtained when considering a tensile phase with allowed orientations of $\alpha = 0^\circ, \alpha \in \{0^\circ; \pm 30^\circ\}, \alpha \in \{0^\circ; \pm 45^\circ\}$ or $\alpha \in \{0^\circ; 90^\circ\}$. As expected, if tensile members can be aligned horizontally only (Figure 8(a)), the most efficient design is obtained with tensile members located at the bottom of the beam and with inclined compressive struts transmitting the load to the supports.
designs are very similar but with a different repartition between both phases. In the present tension/compression splitting formulation Equation (17), each truss member belongs to a single phase, depending on its state of tension and compression whereas, in the previous formulation, some compression members could involve the phase 2 material or even both materials as already discussed.

<figure><img src="image_9.png" /><figcaption>**Figure 7** Bi-material load-maximization of the MBB example with tension/compression splitting ($f_c = f_t = 1$) and various maximum volume fraction $\eta$.</figcaption></figure>

### 5.3.4 Anisotropic strength properties for the tensile phase

We finish this example by considering again a tension/compression splitting formulation, except that the tensile phase (phase 2) now enjoys anisotropic strength properties of the form Equation (9) with here $f_c^{\text{r,eff}} = 0$ and $f_t^{\text{r,eff}} = 1$. We recall that if the set of allowed orientations $\mathcal{A}$ spans all directions, then the corresponding strength criterion is equivalent to a $L_1$-Rankine strength criterion so that we recover the results of Figure 7. In Figure 8, we report the results obtained when considering a tensile phase with allowed orientations of $\alpha = 0^\circ, \alpha \in \{0^\circ; \pm 30^\circ\}, \alpha \in \{0^\circ; \pm 45^\circ\}$ or $\alpha \in \{0^\circ; 90^\circ\}$. As expected, if tensile members can be aligned horizontally only (Figure 8(a)), the most efficient design is obtained with tensile members located at the bottom of the beam and with inclined compressive struts transmitting the load to the supports. Interestingly, we obtain two individual tensile members in this case. When we further allow for inclined directions along $\pm 30^\circ$ or $\pm 45^\circ$, we obtain, in addition to a horizontal tensile member, secondary inclined members to which additional compressive struts can be connected. Note that the case $\alpha \in \{0^\circ; \pm 45^\circ\}$ (Figure 8(c)) is quite close to the isotropic case obtained with the $L_1$-Rankine criterion in Figure 7(c). Finally, the case with $\alpha \in \{0^\circ; 90^\circ\}$ (Figure 8(d)) indeed produces tensile members aligned either horizontally or vertically. Interestingly, despite the fact that tensile orientations are constrained, one still has freedom in the length and location of those members which enables to reach a design where compressive struts can more or less follow the same paths as in the isotropic case. Also note that the convexified formulation in Equation (9) authorizes in theory a superposition of the different orientations at a given point. Clearly, this is not observed since each tensile member is in a pure uniaxial stress state corresponding to a single well-defined orientation. It is only at points corresponding to junctions between tensile and/or compressive members that different orientations may coexist.

Figure 9 compares the corresponding load-bearing capacity of the previous discrete orientation designs with that of the isotropic case. As expected, the isotropic case is the most efficient in terms of load-bearing capacity, both at the beginning and at the end of the penalization procedure. The case with only horizontal tensile members is the less efficient whereas the case with $45^\circ$ shows only slightly lower bearing capacity than the isotropic case which can be explained by the fact that the corresponding design was already quite close to that of the isotropic case. Regarding the computational cost, each outer iteration of the continuation procedure requires solving a convex conic-constrained optimization problem using the Mosek solver. Resolution of each
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Bi-material load-maximization ($\eta = 0.2$) of the MBB example with discrete orientations for the tensile phase</figcaption></figure>

Interestingly, we obtain two individual tensile members in this case. When we further allow for inclined directions along ±30° or ±45°, we obtain, in addition to a horizontal tensile member, secondary inclined members to which additional compressive struts can be connected. Note that the case $\alpha \in \{0°; \pm 45°\}$ (Figure 8(c)) is quite close to the isotropic case obtained with the $L_1$-Rankine criterion in Figure 7(c). Finally, the case with $\alpha \in \{0°; 90°\}$ (Figure 8(d)) indeed produces tensile members aligned either horizontally or vertically. Interestingly, despite the fact that tensile orientations are constrained, one still has freedom in the length and location of those members which enables to reach a design where compressive struts can more or less follow the same paths as in the isotropic case. Also note that the convexified formulation Equation (10) authorizes in theory a superposition of the different orientations at a given point. Clearly, this is not observed since each tensile member is in a pure uniaxial stress state corresponding to a single well-defined orientation. It is only at points corresponding to junctions between tensile and/or compressive members that different orientations may coexist.

Figure 9 compares the corresponding load-bearing capacity of the previous discrete orientation designs with that of the isotropic case. As expected, the isotropic case is the most efficient in terms of load-bearing capacity, both at the beginning and at the end of the penalization procedure. The case with only horizontal tensile members is the less efficient whereas the case with 45° shows only slightly lower bearing capacity than the isotropic case which can be explained by the fact that the corresponding design was already quite close to that of the isotropic case. Regarding the computational cost, each iteration of the continuation procedure requires solving a convex conic-constrained optimization problem using the Mosek solver. Resolution of each inner convex problem is particularly robust throughout the overall procedure and does not pose any convergence issues in general. Depending on the selected tolerance criteria, each resolution requires roughly 20-30 iterations of the primal-dual interior point solver used in Mosek.

## 5.4 Bridge example

We now turn to a bridge-like problem as described in Figure 10. On this example, we investigate the influence of asymmetric tensile/compressive strengths. We again use load-maximization formulations with $\eta = 0.2$ for the imposed volume fraction and $\omega = 1/2$. In particular, Figure 11, Figure 12 and Figure 13 compare the obtained optimized design when using either a bi-material or a single material formulation in the cases where $f_c = f_t = 1$, $f_c = 1$, $f_t = 10$ and $f_c = 10$, $f_t = 1$ respectively. As expected, we observe relatively similar designs from the two formulations in terms of general layout of members in tension and compression. Nevertheless, small differences can be observed, especially on Figure 12. Such differences can be attributed to the fact that the bi-material formulation treats differently the behavior of connections between tension and compression members compared to the single material formulation. One can in
<figure><img src="image_10.png" /><figcaption>**Figure 8** Bi-material load-maximization ($\eta = 0.2$) of the MBB example with discrete orientations for the tensile phase.</figcaption></figure>

convex problem instance is particularly robust throughout the overall procedure and does not pose any convergence issues in general. Depending on the selected tolerance criteria, each outer resolution requires roughly 20-30 inner iterations of the primal-dual interior point solver used in Mosek.

<figure><img src="image_11.png" /><figcaption>**Figure 9** Evolution of the load-bearing capacity during the outer iterations of the penalization procedure for the MBB example with discrete tensile orientations.</figcaption></figure>

## 5.4 Bridge example

We now turn to a bridge-like problem as described in Figure 10. On this example, we investigate the influence of asymmetric tensile/compressive strengths. We again use load-maximization

<figure><img src="image_12.png" /><figcaption>**Figure 10** A bridge structure with a central uniformly distributed loading $T = -e_y$. Fixed supports are distributed over regions of length 0.1 at both extremities.</figcaption></figure>

formulations with $\eta = 0.2$ for the imposed volume fraction and $\omega = 1/2$. In particular, Figures 11 to 13 compare the obtained optimized design when using either a bi-material or a single material formulation in the cases where $f_c = f_t = 1$, $f_c = 1$, $f_t = 10$ and $f_c = 10$, $f_t = 1$ respectively. As expected, we observe relatively similar designs from the two formulations in terms of general layout of members in tension and compression. Nevertheless, small differences can be observed,
<figure><img src="image_10.png" /><figcaption>**Figure 9:** Evolution of the load-bearing capacity during the penalization procedure for the MBB example with discrete tensile orientations</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10:** A bridge structure with a central uniformly distributed loading $T = -e_y$. Fixed supports are distributed over regions of length 0.1 at both extremities.</figcaption></figure>

particular observe that such connections have a larger spatial extension in the bi-material case. Indeed, in locations where both tension and compression phases coexist, the bi-material formulation allows for reduced strength properties compared to the single material phase. For instance, if $ρ_1 = ρ_2$, then at most $ρ_1 = ρ_2 = 0.5$ such that the effective tension and compression strengths are here $f_c/2$ and $f_t/2$, that is half of that of the single material formulation which can reach full capacity. This modeling aspect reflects the fact that anchoring two different materials usually has a detrimental effect on the local strength compared to a single material possessing both tensile and compressive strengths.

Although, the single material formulation of Mourad et al. (2021) enables to obtain similar design as the bimaterial formulation, the latter has the advantage of offering additional modeling choices by using a different factor in the combined volume measure of both phases through the factor $ω$. For instance, Figure 14 represents the evolution of the obtained design when varying this weighting coefficient. We can see that when $ω < 0.5$, tensile (phase 2) material costs more than compression (phase 1) material. The number of tensile members therefore tends to decrease with decreasing $ω$. It must be noted that, although compressive members seem slightly thicker, the number of such members do not necessarily increase with decreasing tensile material since they still need to be connected with tensile members to be supported. Conversely, when $ω$ increases

<figure><img src="image_12.png" /><figcaption>**Figure 11:** Symmetric strengths $f_c = f_t = 1$</figcaption></figure>
especially in Figure 12. Such differences can be attributed to the fact that the bi-material formulation treats differently the behavior of connections between tension and compression members compared to the single material formulation. One can in particular observe that such connections have a larger spatial extension in the bi-material case. Indeed, in locations where both tension and compression phases coexist, the bi-material formulation allows for reduced strength properties compared to the single material phase. For instance, if $ρ_1 = ρ_2$, then at most $ρ_1 = ρ_2 = 0.5$ such that the effective tension and compression strengths are here $f_c/2$ and $f_t/2$, that is half of that of the single material formulation which can reach full capacity. This modeling aspect reflects the fact that anchoring two different materials usually has a detrimental effect on the local strength compared to a single material possessing both tensile and compressive strengths.

<figure><img src="image_13.png" /><figcaption>**Figure 11** Symmetric strengths $f_c = f_t = 1$.</figcaption></figure>

Although, the single material formulation of Mourad et al. (2021) enables to obtain similar design as the bimaterial formulation, the latter has the advantage of offering additional modeling choices by using a different factor in the combined volume measure of both phases through the factor $ω$. For instance, Figure 14 represents the evolution of the obtained design when varying this weighting coefficient. We can see that when $ω < 0.5$, tensile (phase 2) material costs more than compression (phase 1) material. The number of tensile members therefore tends to decrease with decreasing $ω$. It must be noted that, although compressive members seem slightly thicker, the number of such members do not necessarily increase with decreasing tensile material since they still need to be connected with tensile members to be supported. Conversely, when $ω$ increases above the value 0.5, tensile members cost less than compression members and can therefore be utilized more extensively. Figure 15 represents the evolution of the structure load-bearing capacity as a function of the penalization procedure iterations for the considered values of $ω$. Interestingly, the initial value corresponding to the convex problem (BIMAT-LOAD-MAX) does not strongly depend on the density cost parameter $ω$. However, during the penalization procedure, the choice of $ω$ leads to different optimized designs, some of which performing better than other in terms of bearing capacity. In particular, the case $ω = 0.9$ yields a more efficient structure which may be explained by the fact that more tensile material can be considered since it costs less in terms of weighted volume than compressive material.

## 5-5 Deep beam example

We now finish with a deep beam example classically considered when dealing with the design of massive reinforced concrete structures, see Figure 16. Such structures offer great challenges to engineers when aiming at proposing an efficient reinforcement steel layout. Many solutions are possible in practice, depending on the priority given to the simplicity of the design (number of members, orientations, etc.) or to its optimality in terms of steel consumption for instance. In this section, we show how our methodology can easily address both concerns.

With the proposed methodology we consider a total design load $Q = 2$ MN and formulate a bimaterial volume minimization problem. The first phase, representing concrete is modeled with

<figure><img src="image_14.png" /><figcaption>**Figure 12** Asymmetric strengths $f_c = 1, f_t = 10$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12:** Asymmetric strengths $f_c = 1, f_t = 10$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13:** Asymmetric strengths $f_c = 10, f_t = 1$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14:** Bridge example with varying density-cost parameter $\omega$ in the case $f_c = f_t = 1$</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13** Asymmetric strengths $f_c = 10$, $f_t = 1$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14** Bridge example with varying density-cost parameter $\omega$ in the case $f_c = f_t = 1$.</figcaption></figure>

a $L_1$-Rankine criterion with $f_c = 40$ MPa and $f_t = 0.1$ MPa, the second phase, representing steel rebar reinforcements, is modeled with a $L_1$-Rankine criterion with (almost) no compressive strength and an effective tensile strength $f_t = \chi f_y = 40$ MPa where $f_y = 400$ MPa is the steel tensile strength and $\chi = 0.1$ a strength reduction factor accounting for the fact that, in practice, steel rebars only occupy a small fraction (10% at maximum in the present case) of the total beam thickness. Both phases have identical volume cost ($\omega = 0.5$) and we use $l = 0.15$ m.

Figure 17(a) represents the obtained design for both phases. We can clearly see multiple compressive struts aiming at the beam supports, equilibrated by a bottom horizontal tensile reinforcement and a more complex curved reinforcement above the beam opening. Clearly the obtained design is quite similar to strut-and-tie models proposed in (Muttoni et al. 2015) based on elastic stress fields, see Figure 17(b). The global structural behavior in terms of strut and tie location is quite similar between both models. There are however some noticeable differences: since the ST model of Figure 17(b) is based on the interpretation of elastic stress fields, the two main compressive struts exhibit a characteristic “bottle” shape corresponding to the diffusion of elastic stress fields between singular points corresponding to point loads and supports. Such a diffusion induces the presence of small secondary transverse reinforcements which help preventing cracking of concrete in tension. Such features are typical of elastic computations and

<figure><img src="image_17.png" /><figcaption>**Figure 15** Evolution of the load-factor during the outer iterations of the penalization procedure for various cost parameters.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15:** Evolution of the load-factor during the penalization procedure for various cost parameters</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16:** Deep beam problem (25 mm thick) investigated in (Muttoni et al. 2015)</figcaption></figure>

above the value 0.5, tensile members cost less than compression members and can therefore be utilized more extensively. Figure 15 represents the evolution of the structure load-bearing capacity as a function of the penalization procedure iterations for the considered values of $\omega$. Interestingly, the initial value corresponding to the convex problem Equation (BIMAT-LOAD-MAX) does not strongly depend on the density cost parameter $\omega$. However, during the penalization procedure, the choice of $\omega$ leads to different optimized designs, some of which performing better than other in terms of bearing capacity. In particular, the case $\omega = 0.9$ yields a more efficient structure which may be explained by the fact that more tensile material can be considered since it costs less in terms of weighted volume than compressive material.

## 5.5 Deep beam example

We now finish with a deep beam example classically considered when dealing with the design of massive reinforced concrete structures, see Figure 16. Such structures offer great challenges to engineers when aiming at proposing an efficient reinforcement steel layout. Many solutions are possible in practice, depending on the priority given to the simplicity of the design (number of members, orientations, etc.) or to its optimality in terms of steel consumption for instance. In this section, we show how our methodology can easily address both concerns.

With the proposed methodology we consider a total design load $Q = 2$ MN and formulate a bimaterial volume minimization problem. The first phase, representing concrete is modeled with a $L_1$-Rankine criterion with $f_c = 40$ MPa and $f_t = 0.1$ MPa, the second phase, representing steel rebar reinforcements, is modeled with a $L_1$-Rankine criterion with (almost) no compressive strength and an effective tensile strength $f_t = \chi f_y = 40$ MPa where $f_y = 400$ MPa is the steel tensile strength and $\chi = 0.1$ a strength reduction factor accounting for the fact that, in practice,
<figure><img src="image_18.png" /><figcaption>**Figure 16** Deep beam problem (25 mm thick) investigated in (Muttoni et al. 2015).</figcaption></figure>

are therefore absent from limit-analysis based computation which yield very straight compressive struts with constant cross-section. Similarly, we can note that the “array” of struts emerging from the point load is much wider in the ST model than what is obtained with our procedure, for similar reasons.

<figure><img src="image_19.png" /><figcaption>(a) Optimized design from the bimaterial volume minimization.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Strut-and-tie model with inclined reinforcements based on elastic stress fields, taken from (Muttoni et al. 2015). Tensile ties are in solid lines and compressive struts in dotted lines.</figcaption></figure>

<caption>**Figure 17** Comparison of the proposed methodology with strut-and-tie models with inclined reinforcements.</caption>

As pointed out above, such a complex rebar layout might not be practical for on-site placing for instance. Engineers therefore aim at obtaining simpler solutions, although being less optimal, which might be easier to implement in practice. A typical simplifying choice in the design of reinforced concrete structures is to resort to orthogonal reinforcement layout. To achieve this, we simply modify the isotropic strength criterion of the tensile phase with an orthotropic criterion with admissible reinforcement orientations of 0° or 90°. The corresponding result is represented in Figure 18(a) which indeed results in steel rebars being either placed horizontally or vertically. Again, the obtained design is compared with a ST model including only horizontal or vertical reinforcements in Figure 18(b). We can again notice the striking similarity between both models with some differences in terms of steel rebar lengths for instance.

Similarly, the reinforcement optimization formulation where the concrete is not optimized has also been considered on the same example. Results are reported on Figure 19 for both the L₁-Rankine and orthotropic L₁-Rankine criterion for the reinforcement phase. Clearly, this formulation provides reinforcement layouts qualitatively similar to the bi-material case. We can
<figure><img src="image_18.png" /><figcaption>**Figure 17:** Comparison of the proposed methodology with strut-and-tie models with inclined reinforcements</figcaption></figure>

steel rebars only occupy a small fraction (10% at maximum in the present case) of the total beam thickness. Both phases have identical volume cost ($\omega = 0.5$) and we use $l = 0.15$ m.

Figure 17(a) represents the obtained design for both phases. We can clearly see multiple compressive struts aiming at the beam supports, equilibrated by a bottom horizontal tensile reinforcement and a more complex curved reinforcement above the beam opening. Clearly the obtained design is quite similar to strut-and-tie models proposed in (Muttoni et al. 2015) based on elastic stress fields, see Figure 17(b). The global structural behavior in terms of strut and tie location is quite similar between both models. There are however some noticeable differences: since the ST model of Figure 17(b) is based on the interpretation of elastic stress fields, the two main compressive struts exhibit a characteristic "bottle" shape corresponding to the diffusion of elastic stress fields between singular points corresponding to point loads and supports. Such a diffusion induces the presence of small secondary transverse reinforcements which help preventing cracking of concrete in tension. Such features are typical of elastic computations and are therefore absent from limit-analysis based computation which yield very straight compressive struts with constant cross-section. Similarly, we can note that the "array" of struts emerging from the point load is much wider in the ST model than what is obtained with our procedure, for similar reasons.

As pointed out above, such a complex rebar layout might not be practical for on-site placing for instance. Engineers therefore aim at obtaining simpler solutions, although being less optimal, which might be easier to implement in practice. A typical simplifying choice in the design of reinforced concrete structures is to resort to orthogonal reinforcement layout. To achieve this, we simply modify the isotropic strength criterion of the tensile phase with an orthotropic criterion with admissible reinforcement orientations of 0° or 90°. The corresponding result is represented in Figure 18(a) which indeed results in steel rebars being either placed horizontally or vertically. Again, the obtained design is compared with a ST model including only horizontal or vertical reinforcements in Figure 18(b). We can again notice the striking similarity between both models with some differences in terms of steel rebar lengths for instance.

Similarly, the reinforcement optimization formulation where the concrete is not optimized has also been considered on the same example. Results are reported on Figure 19 for both the $L_1$-Rankine and orthotropic $L_1$-Rankine criterion for the reinforcement phase. Clearly, this formulation provides reinforcement layouts qualitatively similar to the bi-material case. We can notice that the amount of reinforcement is however much less than in Figure 17(a) and Figure 18(a) since, in the present case, concrete strength can be mobilized in the whole computational domain,
<figure><img src="image_21.png" /><figcaption>(a) Optimized design from the bimaterial volume minimization with orthotropic reinforcement criterion.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) Strut-and-tie model with orthogonal reinforcements, taken from (Muttoni et al. 2015).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(a) $L_1$-Rankine criterion.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) orthotropic $L_1$-Rankine criterion.</figcaption></figure>

<caption>**Figure 19** Reinforcement optimization of the deep beam example: in black, reinforcement optimized density; in blue, principal compressive stress in the matrix phase.</caption>

notice that the amount of reinforcement is however much less than in Figures 17(a) and 18(a) since, in the present case, concrete strength can be mobilized in the whole computational domain, thereby requiring less reinforcement to sustain the same load. However, the concrete stress field is still well localized and compressive concrete struts can be observed. This observation is interesting since the reinforcement optimization formulation involves only a single density variable and is therefore much less computationally intensive than the bi-material formulation.

# 6 Conclusions, discussion and future work

This work proposed an extension of limit analysis-based topology optimization problems to the case of multiple materials. In particular, the structural load-bearing capacity can be maximized under a material cost constraint. Alternatively, the total material cost function can be minimized for a given imposed loading. Following similar numerical techniques (penalization procedure, filtering) as a previous contribution dedicated to a single material (Mourad et al. 2021), the main contribution of the present work is related to the way composite materials are represented.

We first proposed a reinforcement formulation, in which we aim at optimizing only a reinforcement phase (e.g. steel rebars, fibers, etc.) embedded in a fixed background matrix material (e.g. concrete, soil, resin, etc.). The chosen strength criterion for the reinforcement phase attempts at promoting uniaxial stress states. Reinforcements can be either isotropically distributed (using the $L_1$-Rankine criterion) or with fixed preferential orientation directions (anisotropic criterion in Equation (9)).

Second, we proposed a formulation in which two distinct material phases are simultaneously optimized, each of them possessing its own distinct strength criterion. The latter can either be an isotropic one or an anisotropic one such as Equation (9). Each phase can also represent either two intrinsically distinct materials or represent instead a fictitious phase with a specific stress state,
<figure><img src="image_19.png" /><figcaption>**Figure 18:** Comparison of the proposed methodology with strut-and-tie models with orthogonal reinforcements</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19:** Reinforcement optimization of the deep beam example: in black, reinforcement optimized density; in blue, principal compressive stress in the matrix phase</figcaption></figure>
e.g. pure tension and pure compression as in Section 4.2. The introduction of a material cost function $c_{\omega}(\rho_1, \rho_2)$ also enables to consider a different cost for each phase through the weighting coefficient $\omega$. All these possibilities make the proposed formulation extremely versatile in terms of modeling capabilities.

The question remains which of the reinforcement or bimaterial formulation is more appropriate. Obviously, there is no definitive answer as it depends on the targeted application and the potential need for exploring different topologies. Let us recall that the reinforcement formulation can be seen as a particular case of the bimaterial formulation in which one of the phase pseudo-density is fixed to a value of 1 and where the constraint $\rho_1 + \rho_2 \le 1$ is ignored. The bimaterial formulation is, in this respect, more general but is also much more computationally demanding as two phases must be optimized instead of one. Second, considering a fixed background phase occupying the whole domain produces much more diffuse matrix stress fields in some instances (see Figure 5). For some applications such as the strut-and-tie method for reinforced concrete structures, it might be more interesting to have a clear and simplified stress field in the concrete phase, despite the fact that concrete still occupies the whole structure. In such cases, the bi-material formulation seems particularly attractive. Finally, the latter approach also offers a richer way of exploring various topologies through the choice of the phase criteria or the weighting cost factor. Finally, one can also wonder about the optimality of the solutions found using the proposed penalization procedure. First, let us recall that we are solving a sequence of convex optimization problems so that global optimality of each solution is ensured, especially for the solution of the convexified problems at the first iteration. The question remains open regarding optimality with respect to the initial non-convex problems. We should point out that our penalization procedure follows similar ideas to those of the SIMP method used in elastic compliance minimization problems. In this context, it is known that SIMP solutions are sub-optimal since optimal solutions are given by multi-rank laminates. However, in practice they perform very well. Although we are not aware of any similar result in the limit analysis setting, it might be possible that the same holds in this setting as well.

As regards future works, the possibility of generating strut-and-tie models with the proposed formulation seems very promising. It remains to be checked that the obtained models qualitatively agree with engineering practice on a wider range of examples. One should also check that the corresponding internal forces also quantitatively agree with calculations based on Eurocode design norms. Finally, as the proposed bi-material formulation is computationally demanding, its extension to 3D would require developing dedicated numerical strategies in order to reduce its computational cost.

## A Proof that the L₁-Rankine criterion is the convex hull of isotropically distributed uniaxial strengths

Let us consider an infinite distribution of uniaxial reinforcements indexed by $\alpha$ with tensile and compressive strengths $f_t$ and $f_c$ (independent of $\alpha$). Let us denote by $G^{\alpha} = \{\sigma_{\alpha}e_{\alpha} \otimes e_{\alpha} \text{ s.t. } -f_c \le \sigma_{\alpha} \le f_t\}$ the uniaxial strength criterion of a given family of reinforcement. We show that:

$$G^{L_1\text{-Rankine}} = \text{conv}_{\alpha}\{G^{\alpha}\} \qquad (A.1)$$

where $G^{L_1\text{-Rankine}}$ is the L₁-Rankine criterion defined in Equation (2) which can also be equivalently written as:

$$\sigma \in G^{L_1\text{-Rankine}} \Leftrightarrow \left\{ \begin{array}{l} \sigma_I + \sigma_{II} \le f_t \\ -\sigma_I - \sigma_{II} \le f_c \\ \sigma_I/f_t - \sigma_{II}/f_c \le 1 \\ -\sigma_I/f_c + \sigma_{II}/f_t \le 1 \end{array} \right. \qquad (A.2)$$

Note that although everything is written in a 2D setting, the following proof also holds in 3D.

*Proof.* Let us first show that $G^{L_1\text{-Rankine}} \subseteq \text{conv}_{\alpha}\{G^{\alpha}\}$. For any given stress tensor, we have $\sigma = \sigma_Ie_I \otimes e_I + \sigma_{II}e_{II} \otimes e_{II}$. Then there exist two families of reinforcement $\alpha_I$ and $\alpha_{II}$ of respective
thereby requiring less reinforcement to sustain the same load. However, the concrete stress field is still well localized and compressive concrete struts can be observed. This observation is interesting since the reinforcement optimization formulation involves only a single density variable and is therefore much less computationally intensive than the bi-material formulation.

# 6 Conclusions, discussion and future work

This work proposed an extension of limit analysis-based topology optimization problems to the case of multiple materials. In particular, the structural load-bearing capacity can be maximized under a material cost constraint. Alternatively, the total material cost function can be minimized for a given imposed loading. Following similar numerical techniques (penalization procedure, filtering) as a previous contribution dedicated to a single material (Mourad et al. 2021), the main contribution of the present work is related to the way composite materials are represented.

We first proposed a reinforcement formulation, in which we aim at optimizing only a reinforcement phase (e.g. steel rebars, fibers, etc.) embedded in a fixed background matrix material (e.g. concrete, soil, resin, etc.). The chosen strength criterion for the reinforcement phase attempts at promoting uniaxial stress states. Reinforcements can be either isotropically distributed (using the $L_1$-Rankine criterion) or with fixed preferential orientation directions (anisotropic criterion Equation (10)).

Second, we proposed a formulation in which two distinct material phases are simultaneously optimized, each of them possessing its own distinct strength criterion. The latter can either be an isotropic one or an anisotropic one such as Equation (10). Each phase can also represent either two intrinsically distinct materials or represent instead a fictitious phase with a specific stress state, e.g. pure tension and pure compression as in Section 4.2. The introduction of a material cost function $c_{\omega}(\rho_1, \rho_2)$ also enables to consider a different cost for each phase through the weighting coefficient $\omega$. All these possibilities make the proposed formulation extremely versatile in terms of modeling capabilities.

The question remains which of the reinforcement or bimaterial formulation is more appropriate. Obviously, there is no definitive answer as it depends on the targeted application and the potential need for exploring different topologies. Let us recall that the reinforcement formulation can be seen as a particular case of the bimaterial formulation in which one of the phase pseudo-density is fixed to a value of 1 and where the constraint $\rho_1 + \rho_2 \le 1$ is ignored. The bimaterial formulation is, in this respect, more general but is also much more computationally demanding as two phases must be optimized instead of one. Second, considering a fixed background phase occupying the whole domain produces much more diffuse matrix stress fields in some instances (see Figure 5). For some applications such as the strut-and-tie method for reinforced concrete structures, it might be more interesting to have a clear and simplified stress field in the concrete phase, despite the fact that concrete still occupies the whole structure. In such cases, the bi-material formulation seems particularly attractive. Finally, the latter approach also offers a richer way of exploring various topologies through the choice of the phase criteria or the weighting cost factor. Finally, one can also wonder about the optimality of the solutions found using the proposed penalization procedure. First, let us recall that we are solving a sequence of convex optimization problems so that global optimality of each solution is ensured, especially for the solution of the convexified problems at the first iteration. The question remains open regarding optimality with respect to the initial non-convex problems. We should point out that our peanlization procedure follows similar ideas to those of the SIMP method used in elastic compliance minimization problems. In this context, it is known that SIMP solutions are sub-optimal since optimal solutions are given by multi-rank laminates. However, in practice they perform very well. Although we are not aware of any similar result in the limit analysis setting, it might be possible that the same holds in this setting as well.

As regards future works, the possibility of generating strut-and-tie models with the proposed formulation seems very promising. It remains to be checked that the obtained models qualitatively agree with engineering practice on a wider range of examples. One should also check that the
orientation **e**<sub>I</sub> and **e**<sub>II</sub>. Moreover, assuming that σ ∈ G<sup>L<sub>1</sub>-Rankine</sup> and introducing ζ<sub>I</sub> = g(σ<sub>I</sub>) and ζ<sub>II</sub> = g(σ<sub>II</sub>), we have:

$$-f_c \le \frac{\sigma_I}{\zeta_I} \le f_t \qquad (A.3)$$

$$-f_c \le \frac{\sigma_{II}}{\zeta_{II}} \le f_t \qquad (A.4)$$

and ζ<sub>I</sub> + ζ<sub>II</sub> ≤ 1. As a result, σ can indeed be written as Equation (9) with σ<sup>r,αI</sup> = σ<sub>I</sub>/ζ<sub>I</sub> and σ<sup>r,αII</sup> = σ<sub>II</sub>/ζ<sub>II</sub> so that G<sup>L<sub>1</sub>-Rankine</sup> ⊆ conv<sub>α</sub>{G<sup>α</sup>}.

Now let us show that conv<sub>α</sub>{G<sup>α</sup>} ⊆ G<sup>L<sub>1</sub>-Rankine</sup>. Let σ given by Equation (9). Let **e**<sub>I</sub> and **e**<sub>II</sub> be the principal stress directions of this stress state and introduce θ<sub>α</sub> such that **e**<sub>α</sub> = cos θ<sub>α</sub>**e**<sub>I</sub> + sin θ<sub>α</sub>**e**<sub>II</sub> for any α. Then:

$$\sigma_I = \sum_{\alpha} \zeta_{\alpha} \cos^2 \theta_{\alpha} \sigma^{r,\alpha} \qquad (A.5)$$

$$\sigma_{II} = \sum_{\alpha} \zeta_{\alpha} \sin^2 \theta_{\alpha} \sigma^{r,\alpha} \qquad (A.6)$$

Since ζ<sub>α</sub> ≥ 0 and Σ<sub>α</sub> ζ<sub>α</sub> = 1, we have:

$$\sigma_I + \sigma_{II} = \sum_{\alpha} \zeta_{\alpha} \sigma^{r,\alpha} \le f_t \sum_{\alpha} \zeta_{\alpha} = f_t$$

$$-\sigma_I - \sigma_{II} = \sum_{\alpha} \zeta_{\alpha} (-\sigma^{r,\alpha}) \le f_c \sum_{\alpha} \zeta_{\alpha} = f_c$$

$$\sigma_I/f_t - \sigma_{II}/f_c = \sum_{\alpha} \zeta_{\alpha} (\cos^2 \theta_{\alpha} \sigma^{r,\alpha} / f_t - \sin^2 \theta_{\alpha} \sigma^{r,\alpha} / f_c) \le \sum_{\alpha} \zeta_{\alpha} = 1$$

$$\sigma_I/f_t - \sigma_{II}/f_c = \sum_{\alpha} \zeta_{\alpha} (-\cos^2 \theta_{\alpha} \sigma^{r,\alpha} / f_c + \sin^2 \theta_{\alpha} \sigma^{r,\alpha} / f_t) \le \sum_{\alpha} \zeta_{\alpha} = 1$$

So that σ ∈ G<sup>L<sub>1</sub>-Rankine</sup>.

□

Finally, let us note that in the case where the consider family of reinforcement is discrete α ∈ {α<sub>1</sub>, ..., α<sub>N</sub>}, we only have that conv<sub>α</sub>{G<sup>α</sup>} ⊆ G<sup>L<sub>1</sub>-Rankine</sup>.

## References

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

Alnæs, M. S., J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells (2015). The FEniCS Project Version 1.5. Archive of Numerical Software 3(100):9–23. [DOI].

Amir, O. (2017). Stress-constrained continuum topology optimization: a new approach based on elasto-plasticity. Structural and Multidisciplinary Optimization 55(5):1797–1818. [DOI], [ARXIV].

Bendsøe, M. P. and O. Sigmund (2004). Topology Optimization: Theory, Methods, and Applications. 2nd ed. Engineering online library. Berlin: Springer. 370 pp. [DOI].

[SW] Bleyer, J., fenics_optim – Convex optimization interface in FEniCS. Version 2.0 2020. SWHID: <swh:1:dir:c2c2b16f0797cfbdd51c154d54c957c34a1ee207;origin=https://gitlab.enpc.fr/navier-fenics/fenics-optim;visit=swh:1:snp;bdeeb512425b5937ad9b305da5d6a1844f5bafe2;an chor=swh:1:rev:e9d61cb9c55fbddd8ea07d9ee1e7f9da5703513e>.

Bleyer, J. (2020b). Automating the formulation and resolution of convex variational problems: Applications from image processing to computational mechanics. ACM Transactions on Mathematical Software 46(3):1–33. [DOI], [HAL].

Bogomolny, M. and O. Amir (2012). Conceptual design of reinforced concrete structures using topology optimization with elastoplastic material modeling. International Journal for Numerical Methods in Engineering 90(13):1578–1597. [DOI], [HAL].
corresponding internal forces also quantitatively agree with calculations based on Eurocode
design norms. Finally, as the proposed bi-material formulation is computationally demanding, its
extension to 3D would require developing dedicated numerical strategies in order to reduce its
computational cost.

**Conflict of interest**

The authors declare that they have no conflict of interest.

Replication of results

The Python code for implementing the topology optimization is part of the fenics_optim
Python package (v2.0) (Bleyer 2020a), itself relying on the FEniCS finite-element software library
https://fenicsproject.org/ and the Mosek conic optimization solver https://www.mosek.com/.

References

Allaire, G., C. Dapogny, G. Delgado, and G. Michailidis (2014). “Multi-phase structural optimization via a level set method”. *ESAIM: control, optimisation and calculus of variations* 20.2, pp. 576–611.
DOI: 10.1051/cocv/2013076

Alnæs, M. S., J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells (2015). “The FEniCS Project Version 1.5”. *Archive of Numerical Software* 3.100. DOI: 10.11588/ans.2015.100.20553

Amir, O. (2017). “Stress-constrained continuum topology optimization: a new approach based on elasto-plasticity”. *Structural and Multidisciplinary Optimization* 55.5, pp. 1797–1818. DOI: 10.1007/s00158-016-1618-8

Amir, O. and M. Bogomolny (June 2012). “Conceptual design of reinforced concrete structures using topology optimization with elastoplastic material modeling”. *International Journal for Numerical Methods in Engineering* 90, pp. 1578–1597. DOI: 10.1002/nme.4253

Bendsøe, M. P. and O. Sigmund (2004). *Topology optimization: theory, methods, and applications*. 2. ed., corr. printing. Engineering online library. OCLC: 249399186. Berlin: Springer. 370 pp.

Bleyer, J. (Jan. 2020a). *fenics_optim – Convex optimization interface in FEniCS*. Version 2.0. Zenodo, https://doi.org/10.5281/zenodo.3604085. DOI: 10.5281/zenodo.3604085

Bleyer, J. (2020b). “Automating the Formulation and Resolution of Convex Variational Problems: Applications from Image Processing to Computational Mechanics”. *ACM Transactions on Mathematical Software* 0.ja. DOI: 10.1145/3393881

Bruggi, M. (Dec. 2009). “Generating strut-and-tie patterns for reinforced concrete structures using topology optimization”. *Computers & Structures* 87, pp. 1483–1495. DOI: 10.1016/j.compstruc.2009.06.003

Bruggi, M. (2016). “Analysis and design of reinforced concrete structures as a topology optimization problem”. *VII European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2016*, pp. 1–9

Chavez, G. M. (2018). “Strut-and-Tie models for the design of non-flexural elements: computational aided approach”. PhD thesis. Université Paris-Est

Collins, M. P. and D. Mitchell (1980). “Shear and torsion design of prestressed and non-prestressed concrete beams”. *Journal of the Prestressed Concrete Institute*, 25(5), 32–100. DOI: 10.15554/pcij.09011980.32.100

de Buhan, P. and A. Taliercio (1991). “A homogenization approach to the yield strength of composite materials”. *European Journal of Mechanics A-solids* 10, pp. 129–154. DOI: 10.1016/S0022-5096(97)00002-1

de Buhan, P., J. Bleyer, and G. Hassen (2017). *Elastic, Plastic and Yield Design of Reinforced Structures*. Elsevier

Dorn, W., R. E. Gomory, and H. G. Greenberg (1964). “Automatic design of optimal structures”. *Journal de Mécanique* 3, pp. 25–52
Bruggi, M. (2009). Generating strut-and-tie patterns for reinforced concrete structures using topology optimization. *Computers & Structures* 87(23-24):1483-1495. [DOI].

Bruggi, M. (2016). Analysis and design of reinforced concrete structures as a topology optimization problem. *VII European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS)*, pp 1-9. [DOI], [OA].

Chavez, G. M. (2018). Strut-and-Tie models for the design of non-flexural elements: computational aided approach. PhD thesis. France: Université Paris-Est. [OA].

Collins, M. P. and D. Mitchell (1980). Shear and torsion design of prestressed and non-prestressed concrete beams. *PCI Journal* 25(5):32-100. [DOI].

de Buhan, P. and A. Taliercio (1991). A homogenization approach to the yield strength of composite materials. *European Journal of Mechanics A-solids* 10:129-154.

de Buhan, P., J. Bleyer, and G. Hassen (2017). *Elastic, Plastic and Yield Design of Reinforced Structures*. Elsevier. ISBN: 978-1-78548-205-2.

Dorn, W. C., R. E. Gomory, and H. G. Greenberg (1964). Automatic design of optimal structures. *Journal de Mécanique* 3:25-52.

Duysinx, P. and M. P. Bendsøe (1998). Topology optimization of continuum structures with local stress constraints. *International Journal for Numerical Methods in Engineering* 43(8):1453-1478. [DOI].

Fin, J., L. A. Borges, and E. A. Fancello (2018). Structural topology optimization under limit analysis. *Structural and Multidisciplinary Optimization* 59(4):1355-1370. [DOI].

Gaynor, A. T., J. K. Guest, and C. D. Moen (2013). Reinforced concrete force visualization and design using bilinear truss-continuum topology optimization. *Journal of Structural Engineering* 139(4):607-618. [DOI].

Gibiansky, L. V. and O. Sigmund (2000). Multiphase composites with extremal bulk modulus. *Journal of the Mechanics and Physics of Solids* 48(3):461-498. [DOI].

Gilbert, M. and A. Tyas (2003). Layout optimization of large-scale pin-jointed frames. *Engineering Computations* 20(8):1044-1064. [DOI].

Guest, J. K. (2009). Topology optimization with multiple phase projection. *Computer Methods in Applied Mechanics and Engineering* 199(1-4):123-135. [DOI].

He, L. and M. Gilbert (2015). Rationalization of trusses generated via layout optimization. *Structural and Multidisciplinary Optimization* 52(4):677-694. [DOI].

Herfelt, M. A., P. N. Poulsen, and L. C. Hoang (2018). Strength-based topology optimisation of plastic isotropic von Mises materials. *Structural and Multidisciplinary Optimization* 59(3):893-906. [DOI].

Kammoun, Z. and H. Smaoui (2014). A direct approach for continuous topology optimization subject to admissible loading. *Comptes Rendus Mécanique* 342(9):520-531. [DOI], [OA].

Liu, J. and Y. Ma (2018). A new multi-material level set topology optimization method with the length scale control capability. *Computer Methods in Applied Mechanics and Engineering* 329:444-463. [DOI].

Logg, A., K.-A. Mardal, and G. Wells (2012). *Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book*. Vol. 84. Berlin Heidelberg: Springer. [DOI].

Marti (1980). *Zur Plastischen Berechnung von Stahlbeton (On plastic analysis of reinforced concrete)*. Tech. rep. Institut für Baustatik und Konstruktion (Institute of Structural Engineers), ETH, Zürich, Switzerland, No. 104. [OA].

Maute, K., S. Schwarz, and E. Ramm (1998). Adaptive topology optimization of elastoplastic structures. *Structural Optimization* 15(2):81-91. [DOI].

Michell, A. G. M. (1904). LVIII. The limits of economy of material in frame-structures. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science* 8(47):589-597.

Mörsch, E. and E. Goodrich (1909). *Concrete-steel Construction: (Der Eisenbetonbau)*. Engineering News Publishing Company. ISBN: 1245947451.

MOSEK ApS, (2018). *The MOSEK optimization API for Python 8.1.0*.

Mourad, L., J. Bleyer, R. Mesnil, J. Nseir, K. Sab, and W. Raphael (2021). Topology optimization of load-bearing capacity. *Structural and Multidisciplinary Optimization* 64 (3):1367-1383. [DOI] [HAL].

Muttoni, A., M. Fernández Ruiz, and F. Niketić (2015). Design versus Assessment of Concrete
Duysinx, P. and M. P. Bendsøe (1998). "Topology optimization of continuum structures with local stress constraints". *International journal for numerical methods in engineering* 43.8, pp. 1453-1478. doi: 10.1002/(SICI)1097-0207(19981230)43:8<1453::AID-NME480>3.0.CO;2-2

Fin, J., L. Borges, and E. Fancello (Nov. 2018). "Structural topology optimization under limit analysis". *Structural and Multidisciplinary Optimization* 59. doi: 10.1007/s00158-018-2132-y

Gaynor, A. T., J. K. Guest, and C. D. Moen (2013). "Reinforced Concrete Force Visualization and Design Using Bilinear Truss-Continuum Topology Optimization". *Journal of Structural Engineering* 139.4, pp. 607-618. doi: 10.1061/(ASCE)ST.1943-541X.0000692

Gibiansky, L. V. and O. Sigmund (2000). "Multiphase composites with extremal bulk modulus". *Journal of the Mechanics and Physics of Solids* 48.3, pp. 461-498. doi: 10.1016/S0022-5096(99)00043-5

Gilbert, M. and A. Tyas (2003). "Layout optimization of large-scale pin-jointed frames". *Engineering computations*. doi: 10.1108/02644400310503017

Guest, J. (Dec. 2009). "Topology optimization with multiple phase projection". *Computer Methods in Applied Mechanics and Engineering* 199, pp. 123-135. doi: 10.1016/j.cma.2009.09.023

He, L. and M. Gilbert (2015). "Rationalization of trusses generated via layout optimization". *Structural and Multidisciplinary Optimization* 52.4, pp. 677-694. doi: 10.1007/s00158-015-1260-x

Herfelt, M., P. Poulsen, and L. Hoang (Oct. 2018). "Strength-based topology optimisation of plastic isotropic von Mises materials". *Structural and Multidisciplinary Optimization*. doi: 10.1007/s00158-018-2108-y

Kammoun, Z. and H. Smaoui (2014). "A direct approach for continuous topology optimization subject to admissible loading". *Comptes Rendus Mécanique* 342.9, pp. 520-531. doi: 10.1016/j.crme.2014.06.003

Liu, J. and Y. Ma (2018). "A new multi-material level set topology optimization method with the length scale control capability". *Computer Methods in Applied Mechanics and Engineering* 329, pp. 444-463. doi: 10.1016/j.cma.2017.10.011

Logg, A., K.-A. Mardal, and G. Wells (2012). *Automated solution of differential equations by the finite element method: The FEniCS book*. Vol. 84. Springer Science & Business Media

Marti (1980). *On plastic analysis of reinforced concrete*. Institute of Structural Engineers, ETH, Zurich, Switzerland, Rep. No. 104

Maute, K., S. Schwarz, and E. Ramm (1998). "Adaptive topology optimization of elastoplastic structures". *Structural Optimization* 15.2, pp. 81-91. doi: 10.1007/BF01278493

Michell, A. G. M. (1904). "LVIII. The limits of economy of material in frame-structures". *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science* 8.47, pp. 589-597

Mörsch, E. and E. Goodrich (1909). *Concrete-steel Construction: (Der Eisenbetonbau)*. Engineering News Publishing Company

MOSEK ApS, (2018). *The MOSEK optimization API for Python* 8.1.0

Mourad, L., J. Bleyer, R. Mesnil, J. Nseir, K. Sab, and W. Raphael (2021). "Topology optimization of load-bearing capacity". *Structural and Multidisciplinary Optimization* 64 (3). doi: 10.1007/s00158-021-02923-1

Muttoni, A., M. Fernández Ruiz, and F. Niketić (2015). "Design versus Assessment of Concrete Structures Using Stress Fields and Strut-and-Tie Models". *ACI Structural Journal* 112.5. doi: 10.14359/51687710

Pedersen, P. (1998). "Some general optimal design results using anisotropic, power law nonlinear elasticity". *Structural Optimization* 15.2, pp. 73-80. doi: 10.1007/bf01278492

Petersson, J. and O. Sigmund (1998). "Slope constrained topology optimization". *International Journal for Numerical Methods in Engineering* 41.8, pp. 1417-1434. doi: 10.1002/(SICI)1097-0207(19980430)41:8<1417::AID-NME344>3.0.CO;2-N

Prager, W. (1974). *Introduction to structural optimization*. Springer

Querin, O. M., M. Victoria, and P. Martí (2010). "Topology optimization of truss-like continua with different material properties in tension and compression". *Structural and Multidisciplinary Optimization* 42.1, pp. 25-32. doi: 10.1007/s00158-009-0473-2

Ritter, W. (1899). "Die Bauweise Hennebique (the Hennebique system)". *Schweizerische Bauzeitung* 33.7, pp. 59-61
Structures Using Stress Fields and Strut-and-Tie Models. *ACI Structural Journal* 112(5). [DOI],
[OA].

Petersson, J. and O. Sigmund (1998). Slope constrained topology optimization. *International Journal for Numerical Methods in Engineering* 41(8):1417-1434. [DOI].

Podersen, P. (1998). Some general optimal design results using anisotropic, power law nonlinear elasticity. *Structural Optimization* 15(2):73-80. [DOI].

Prager, W. (1974). *Introduction to structural optimization*. Vienna: Springer. [DOI].

Querin, O. M., M. Victoria, and P. Martí (2010). Topology optimization of truss-like continua with different material properties in tension and compression. *Structural and Multidisciplinary Optimization* 42(1):25-32. [DOI].

Ritter, W. (1899). Die Bauweise Hennebique (the Hennebique system). *Schweizerische Bauzeitung* 33(7):59-61.

Rostami, P. and J. Marzbanrad (2020). Multi-material topology optimization of compliant mechanisms using regularized projected gradient approach. *Journal of the Brazilian Society of Mechanical Sciences and Engineering* 42(9). [DOI].

Rozvany, G. I. N. and T. G. Ong (1986). Optimal plastic design of plates, shells and shellgrids. *Inelastic Behaviour of Plates and Shells*. International Union of Theoretical and Applied Mechanics. Berlin, Heidelberg: Springer, pp 357-384. [DOI].

Schlaich, J., K. Schafer, and M. Jennewein (1987). Toward a consistent design of structural concrete. *PCI Journal* 32(3):74-150. [DOI].

Sigmund, O. and S. Torquato (1997). Design of materials with extreme thermal expansion using a three-phase topology optimization method. *Journal of the Mechanics and Physics of Solids* 45(6):1037-1067. [DOI].

Smarslik, M., M. A. Ahrens, and P. Mark (2019). Toward holistic tension-or compression-biased structural designs using topology optimization. *Engineering Structures* 199:109632. [DOI].

Strang, G. and R. V. Kohn (1986). Optimal design in elasticity and plasticity. *International Journal for Numerical Methods in Engineering* 22(1):183-188. [DOI].

Swan, C. C. and I. Kosaka (1997). Voigt-Reuss topology optimization for structures with nonlinear material behaviors. *International Journal for Numerical Methods in Engineering* 40(20):3785-3814. [DOI].

Tavakoli, R. (2014). Multimaterial topology optimization by volume constrained Allen-Cahn system and regularized projected steepest descent method. *Computer Methods in Applied Mechanics and Engineering* 276:534-565. [DOI].

Victoria, M., O. M. Querin, and P. Martí (2011). Generation of strut-and-tie models by topology design using different material properties in tension and compression. *Structural and Multidisciplinary Optimization* 44(2):247-258. [DOI].

Wallin, M., N. Ivarsson, and M. Ristinmaa (2015). Large strain phase-field-based multi-material topology optimization. *International Journal for Numerical Methods in Engineering* 104(9):887-904. [DOI].

Wallin, M., V. Jönsson, and E. Wingren (2016). Topology optimization based on finite strain plasticity. *Structural and Multidisciplinary Optimization* 54(4):783-793. [DOI].

Zhang, X. S., H. Chi, and G. H. Paulino (2020). Adaptive multi-material topology optimization with hyperelastic materials under large deformations: A virtual element approach. *Computer Methods in Applied Mechanics and Engineering* 370:112976. [DOI].

Zhang, X. S., G. H. Paulino, and A. S. Ramos (2018). Multi-material topology optimization with multiple volume constraints: a general approach applied to ground structures with material nonlinearity. *Structural and Multidisciplinary Optimization* 57(1):161-182. [DOI].

Zhou, S. and M. Y. Wang (2007). Multimaterial structural topology optimization with a generalized Cahn-Hilliard model of multiphase transition. *Structural and Multidisciplinary Optimization* 33(2):89-111. [DOI].

Zuo, W. and K. Saitou (2017). Multi-material topology optimization using ordered SIMP interpolation. *Structural and Multidisciplinary Optimization* 55(2):477-491. [DOI].
Rostami, P. and J. Marzbanrad (Sept. 2020). "Multi-material topology optimization of compliant mechanisms using regularized projected gradient approach". *Journal of the Brazilian Society of Mechanical Sciences and Engineering* 42. DOI: 10.1007/s40430-020-02549-2

Rozvany, G. and T. Ong (1986). "Optimal plastic design of plates, shells and shellgrids". *Inelastic behaviour of plates and shells*. Springer, pp. 357–384

Schlaich, J., K. Schäfer, and M. Jennewein (1987). "Toward a consistent design of structural concrete". *PCI journal* 32.3, pp. 74–150

Sigmund, O. and S. Torquato (1997). "Design of materials with extreme thermal expansion using a three-phase topology optimization method". *Journal of the Mechanics and Physics of Solids* 45.6, pp. 1037–1067. DOI: 10.1016/S0022-5096(96)00114-7

Smarslik, M., M. A. Ahrens, and P. Mark (2019). "Toward holistic tension-or compression-biased structural designs using topology optimization". *Engineering Structures* 199, p. 109632. DOI: 10.1016/j.engstruct.2019.109632

Strang, G. and R. V. Kohn (1986). "Optimal design in elasticity and plasticity". *International journal for numerical methods in engineering* 22.1, pp. 183–188. DOI: 10.1002/nme.1620220113

Swan, C. C. and I. Kosaka (1997). "Voigt-Reuss topology optimization for structures with nonlinear material behaviors". *International Journal for Numerical Methods in Engineering* 40.20, pp. 3785–3814. DOI: 10.1002/(SICI)1097-0207(19971030)40:20<3785::AID-NME240>3.0.CO;2-V

Tavakoli, R. (July 2014). "Multimaterial topology optimization by volume constrained Allen-Cahn system and regularized projected steepest descent method". *Computer Methods in Applied Mechanics and Engineering* 276. DOI: 10.1016/j.cma.2014.04.005

Victoria Nicolás, M., O. Querin, and P. Martí-Montrull (Aug. 2011). "Generation of strut-and-tie models by topology design using different material properties in tension and compression". *Structural and Multidisciplinary Optimization* 44, pp. 247–258. DOI: 10.1007/s00158-011-0633-z

Wallin, M., N. Ivarsson, and M. Ristinmaa (2015). "Large strain phase-field-based multi-material topology optimization". *International Journal for Numerical Methods in Engineering* 104.9, pp. 887–904. DOI: 10.1002/nme.4962

Wallin, M., V. Jönsson, and E. Wingren (2016). "Topology optimization based on finite strain plasticity". *Structural and multidisciplinary optimization* 54.4, pp. 783–793. DOI: 10.1007/s00158-016-1435-0

Zhang, X., H. Chi, and G. Paulino (Oct. 2020). "Adaptive multi-material topology optimization with hyperelastic materials under large deformations: A virtual element approach". *Computer Methods in Applied Mechanics and Engineering* 370, p. 112976. DOI: 10.1016/j.cma.2020.112976

Zhang, X. S., G. H. Paulino, and A. S. Ramos (2018). "Multi-material topology optimization with multiple volume constraints: a general approach applied to ground structures with material nonlinearity". *Structural and Multidisciplinary Optimization* 57.1, pp. 161–182. DOI: 10.1007/s00158-017-1768-3

Zhou, S. and M. Y. Wang (2007). "Multimaterial structural topology optimization with a generalized Cahn-Hilliard model of multiphase transition". *Structural and Multidisciplinary Optimization* 33.2, p. 89. DOI: 10.1007/s00158-006-0035-9

Zuo, W. and K. Saitou (Feb. 2017). "Multi-material topology optimization using ordered SIMP interpolation". *Structural and Multidisciplinary Optimization* 55. DOI: 10.1007/s00158-016-1513-3

## 7 Backmatter information

The following subsections will be added to the accepted version of the paper.

**Funding** This work is part of the PhD thesis of L. Mourad who is supported by Université Paris-Est and Université Saint-Joseph.

**Authors' contributions** Author 1 carried out most of the study, performed numerical simulations, and drafted the manuscript. Author 2 helped with implementation and numerical issues. All authors developed the methodology, conceived of the study, and participated in its design, coordination, and critical review of the manuscript. All authors read and approved the final manuscript.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** LM carried out most of the study, performed numerical simulations, and drafted the manuscript. JB helped with implementation and numerical issues. All authors developed the methodology, conceived the study, and participated in its design, coordination, and critical review of the manuscript. All authors read and approved the final manuscript.

**Supplementary Material** The Python code for implementing the topology optimization is part of the fenics_optim Python package (v2.0) (Bleyer 2020a), itself relying on the FEniCS finite-element software library https://fenicsproject.org/ and the Mosek conic optimization solver https://www.mosek.com/. Corresponding datasets are available at the permalink zenodo/14773973.

**Acknowledgements** The authors would like to thank Jean-Marc JAEGER, Michel BUÉ and Camille CAREME (Setec tpi) for fruitful discussions on this topic.

**Funding** This work is part of the PhD thesis of L. MOURAD who is supported by Université Paris-Est and Université Saint-Joseph.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.
**Ethics approval and consent to participate** If applicable.

**Consent for publication** If applicable.

**Competing interests** Competing interests are defined as financial and non-financial interests that could directly undermine, or be perceived to undermine the objectivity, integrity and value of a publication, through a potential influence on the judgments and actions of authors with regard to objective data presentation, analysis and interpretation. Examples:

* The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

* The authors declare that they have no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a full copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

# A Appendices

## A.1 Proof that the L₁-Rankine criterion is the convex hull of isotropically distributed uniaxial strengths

Let us consider an infinite distribution of uniaxial reinforcements indexed by $\alpha$ with tensile and compressive strengths $f_t$ and $f_c$ (independent of $\alpha$). Let us denote by $G^\alpha = \{\sigma_\alpha e_\alpha \otimes e_\alpha \text{ s.t. } -f_c \le \sigma_\alpha \le f_t\}$ the uniaxial strength criterion of a given family of reinforcement. We show that:

$$G^{L_1\text{-Rankine}} = \text{conv}_{\alpha}\{G^{\alpha}\} \qquad (20)$$

where $G^{L_1\text{-Rankine}}$ is the L₁-Rankine criterion defined in Equation (2) which can also be equivalently written as:

$$\sigma \in G^{L_1\text{-Rankine}} \Leftrightarrow \left\{ \begin{array}{l} \sigma_I + \sigma_{II} \le f_t \\ -\sigma_I - \sigma_{II} \le f_c \\ \sigma_I/f_t - \sigma_{II}/f_c \le 1 \\ -\sigma_I/f_c + \sigma_{II}/f_t \le 1 \end{array} \right. \qquad (21)$$

Note that although everything is written in a 2D setting, the following proof also holds in 3D.

*Proof.* Let us first show that $G^{L_1\text{-Rankine}} \subseteq \text{conv}_\alpha\{G^\alpha\}$. For any given stress tensor, we have $\sigma = \sigma_I e_I \otimes e_I + \sigma_{II} e_{II} \otimes e_{II}$. Then there exist two families of reinforcement $\alpha_I$ and $\alpha_{II}$ of respective orientation $e_I$ and $e_{II}$. Moreover, assuming that $\sigma \in G^{L_1\text{-Rankine}}$ and introducing $\zeta_I = g(\sigma_I)$ and $\zeta_{II} = g(\sigma_{II})$, we have:

$$-f_c \le \frac{\sigma_I}{\zeta_I} \le f_t \qquad (22)$$

$$-f_c \le \frac{\sigma_{II}}{\zeta_{II}} \le f_t \qquad (23)$$

and $\zeta_I + \zeta_{II} \le 1$. As a result, $\sigma$ can indeed be written as Equation (10) with $\sigma^{r,\alpha_I} = \sigma_I/\zeta_I$ and $\sigma^{r,\alpha_{II}} = \sigma_{II}/\zeta_{II}$ so that $G^{L_1\text{-Rankine}} \subseteq \text{conv}_\alpha\{G^\alpha\}$.

Now let us show that $\text{conv}_\alpha\{G^\alpha\} \subseteq G^{L_1\text{-Rankine}}$. Let $\sigma$ given by Equation (10). Let $e_I$ and $e_{II}$ be the principal stress directions of this stress state and introduce $\theta_\alpha$ such that $e_\alpha =$
$\cos \theta_{\alpha} e_{I} + \sin \theta_{\alpha} e_{II}$ for any $\alpha$. Then:

$$ \sigma_I = \sum_{\alpha} \zeta_{\alpha} \cos^2 \theta_{\alpha} \sigma^{r,\alpha} \qquad (24) $$

$$ \sigma_{II} = \sum_{\alpha} \zeta_{\alpha} \sin^2 \theta_{\alpha} \sigma^{r,\alpha} \qquad (25) $$

Since $\zeta_{\alpha} \ge 0$ and $\sum_{\alpha} \zeta_{\alpha} = 1$, we have:

$$ \sigma_I + \sigma_{II} = \sum_{\alpha} \zeta_{\alpha} \sigma^{r,\alpha} $$

$$ \le f_t \sum_{\alpha} \zeta_{\alpha} = f_t $$

$$ -\sigma_I - \sigma_{II} = \sum_{\alpha} \zeta_{\alpha} (-\sigma^{r,\alpha}) $$

$$ \le f_c \sum_{\alpha} \zeta_{\alpha} = f_c $$

$$ \sigma_I/f_t - \sigma_{II}/f_c = \sum_{\alpha} \zeta_{\alpha} (\cos^2 \theta_{\alpha} \sigma^{r,\alpha} / f_t - \sin^2 \theta_{\alpha} \sigma^{r,\alpha} / f_c) $$

$$ \le \sum_{\alpha} \zeta_{\alpha} = 1 $$

$$ \sigma_I/f_t - \sigma_{II}/f_c = \sum_{\alpha} \zeta_{\alpha} (-\cos^2 \theta_{\alpha} \sigma^{r,\alpha} / f_c + \sin^2 \theta_{\alpha} \sigma^{r,\alpha} / f_t) $$

$$ \le \sum_{\alpha} \zeta_{\alpha} = 1 $$

So that $\sigma \in G^{L_1-\text{Rankine}}$. $\square$

Finally, let us note that in the case where the consider family of reinforcement is discrete $\alpha \in \{\alpha_1, \dots, \alpha_N\}$, we only have that $\text{conv}_{\alpha}\{G^{\alpha}\} \subseteq G^{L_1-\text{Rankine}}$.