Gradient Based Biobjective Shape Optimization to Improve Reliability and Cost of Ceramic Components

TL;DR

This paper presents a gradient-based biobjective shape optimization approach for ceramic components, balancing reliability and cost using PDE constraints and adjoint methods, demonstrated through 2D test cases.

Contribution

It introduces two gradient-based optimization methods for biobjective shape optimization with PDE constraints, utilizing efficient adjoint gradient computations and demonstrating trade-offs in 2D ceramic component design.

Findings

Pareto fronts effectively illustrate reliability-cost trade-offs

Innovative shapes emerge balancing conflicting objectives

Methods outperform traditional approaches in efficiency

Abstract

We consider the simultaneous optimization of the reliability and the cost of a ceramic component in a biobjective PDE constrained shape optimization problem. A probabilistic Weibull-type model is used to assess the probability of failure of the component under tensile load, while the cost is assumed to be proportional to the volume of the component. Two different gradient-based optimization methods are suggested and compared at 2D test cases. The numerical implementation is based on a first discretize then optimize strategy and benefits from efficient gradient computations using adjoint equations. The resulting approximations of the Pareto front nicely exhibit the trade-off between reliability and cost and give rise to innovative shapes that compromise between these conflicting objectives.

Equations

$$\partial \Omega =\text{cl}({\partial \Omega }_D)\cup \text{cl}({\partial \Omega }_{N_{\text{fixed}}})\cup \text{cl}({\partial \Omega }_{N_{\text{free}}}).$$

$$\forall x\in \partial \hat{\Omega }\exists {\zeta }_x\in {\mathbb{R}}^p,\Vert {\zeta }_x\Vert =1:y+C({\zeta }_x,\theta ,l)\subset \hat{\Omega }\forall y\in B(x,r)\cap \hat{\Omega },$$

$${\mathcal{O}}^{\text{ad}}:=\left\{\Omega \subseteq \hat{\Omega }:\partial {\Omega }_D\subseteq \partial \Omega ,\partial {\Omega }_{N_{\text{fixed}}}\subseteq \partial \Omega ,\hat{\Omega }\text{ and }\Omega \text{ satisfy the cone property}\right\}.$$

$$\begin{array}[]{rcll}-\text{div}(\sigma(u(x)))&=&\tilde{f}(x)&\text{for}\;x\in\Omega\\ u(x)&=&0&\text{for}\;x\in\partial\Omega_{D}\\ \sigma(u(x))\hat{n}(x)&=&\tilde{g}(x)&\text{for}\;x\in\partial\Omega_{N_{\text{fixed}}}\\ \sigma(u(x))\hat{n}(x)&=&0&\text{for}\;x\in\partial\Omega_{N_{\text{free}}}\end{array}$$

$$\displaystyle\sigma(x)=\sigma(r,\phi)=\frac{1}{\sqrt{2\pi r}}\Bigl{\{}K_{\text{I}}\tilde{\sigma}^{\text{I}}(\phi)+K_{\text{II}}\tilde{\sigma}^{\text{II}}(\phi)+K_{\text{III}}\tilde{\sigma}^{\text{III}}(\phi)\Bigr{\}}+R(r,\phi).$$

$$(x,a,n)\in \mathcal{C}:=\Omega \times (0,\infty )\times S^{p-1},$$

$$K_{\text{I}}=K_{\text{I}}(a,{\sigma }_n(Du(x)))=\frac{2}{\pi }{\sigma }_n(Du(x))\sqrt{\pi a},$$

$${\sigma }_n(Du(x)):=\max \{n^{\top }\sigma (Du(x))n,0\}.$$

$$a>a_c:=\frac{\pi }{4}{\left(\frac{K_{\text{I}c}}{{\sigma }_n(Du(x))}\right)}^2$$

$$A_c:=A_c(\Omega ,Du)=\{(x,a,n)\in \mathcal{C}:K_{\text{I}}(a,{\sigma }_n(Du(x)))>K_{\text{I}c}\}$$

$$p_s(\Omega |Du)=P(N(A_c(\Omega ,Du))=0)=\exp \{-\upsilon (A_c(\Omega ,Du))\}.$$

$$\upsilon (A_c(\Omega ,Du))=\frac{\Gamma (\frac{p}{2})}{2{\pi }^{\frac{p}{2}}}\underset{\Omega }{\int }\underset{S^{p-1}}{\int }\underset{a_c}{\overset{\infty }{\int }}\text{d}{\upsilon }_a(a)\text{d}n\text{d}x$$

$$\upsilon (A_c(\Omega ,Du))=\frac{\Gamma (\frac{p}{2})}{2{\pi }^{\frac{p}{2}}}\underset{\Omega }{\int }\underset{S^{p-1}}{\int }{\left(\frac{{\sigma }_n(Du(x))}{{\sigma }_0}\right)}^m\text{d}n\text{d}x,$$

$$f_1(\Omega ):=\upsilon (A_c(\Omega ,Du))$$

$$f_2(\Omega ):=\underset{\Omega }{\int }\text{d}x.$$

$$\begin{array}{rl}\underset{\Omega \in {\mathcal{O}}^{\text{ad}}}{\min }& \text{leavevmode}f(\Omega ):=(f_1(\Omega ),f_2(\Omega ))\\ \text{s.t. }& u\in H^1(\Omega ,{\mathbb{R}}^p)\text{ solves the state equation }(\text{ref}stateequation),\end{array}$$

$$\begin{array}[]{ccl}z^{1}\leqq z^{2}&\iff&z^{1}_{j}\leq z^{2}_{j},\;j=1,2\\[5.69046pt] z^{1}\leqslant z^{2}&\iff&z^{1}\leqq z^{2}\text{ and }z^{1}\neq z^{2}\\[5.69046pt] z^{1}<z^{2}&\iff&z^{1}_{j}<z^{2}_{j},\;j=1,2.\end{array}$$

$${\mathbb{R}}_{⩾}^2:=\{z\in {\mathbb{R}}^2:z⩾(0,0{)}^{\top }\}\text{and}\bar{z}+{\mathbb{R}}_{⩾}^2:=\{z\in {\mathbb{R}}^2:z⩾\bar{z}\}\text{forleavevmodenobreak }\bar{z}\in {\mathbb{R}}^2.$$

$$\begin{array}{rl}\underset{\Omega \in {\mathcal{O}}^{\text{ad}}}{\min }& f_{\omega }(\Omega ):=\omega f_1(\Omega )+(1-\omega )f_2(\Omega )\\ \text{s.t. }& u\in H^1(\Omega ,{\mathbb{R}}^p)\text{ solves the state equation }(\text{ref}stateequation).\end{array}$$

$$f_{\omega }(\Omega )$$

$$X_{i,j}:=\left(x_i,{\varrho }_i^{\text{ml}}+\frac{{\varrho }_i^{\text{th}}}{n_y-1}\left(j-\frac{n_y+1}{2}\right)\right)\in {\mathbb{R}}^2,i=1,\dots ,n_x,j=1,\dots ,n_y.$$

$${\hat{\varrho }}^{\text{ml}}(x):=\sum _{j=1}^{n_B}{\gamma }_j^{\text{ml}}{\vartheta }_j(x)\text{and}{\hat{\varrho }}^{\text{th}}(x):=\sum _{j=1}^{n_B}{\gamma }_j^{\text{th}}{\vartheta }_j(x),x\in \mathbb{R}.$$

$${\varrho }_i^{\text{ml}}:={\hat{\varrho }}^{\text{ml}}(x_i)\text{and}{\varrho }_i^{\text{th}}:={\hat{\varrho }}^{\text{th}}(x_i),i=1,\dots ,n_x.$$

$$\frac{\partial f_j}{\partial {\gamma }^{\text{ml}}}=\frac{\partial f_j}{\partial X}\frac{\partial X}{\partial {\varrho }^{\text{ml}}}\frac{\partial {\varrho }^{\text{ml}}}{\partial {\gamma }^{\text{ml}}}\text{and}\frac{\partial f_j}{\partial {\gamma }^{\text{th}}}=\frac{\partial f_j}{\partial X}\frac{\partial X}{\partial {\varrho }^{\text{th}}}\frac{\partial {\varrho }^{\text{th}}}{\partial {\gamma }^{\text{th}}},j=1,2.$$

$$\begin{array}{rl}\underset{\gamma \in \Gamma }{\min }& (f_1(\gamma ),f_2(\gamma ))\\ \text{s.t. }& u(X(\gamma ))\text{ solves the discretized state equation }(\text{ref}stateequation).\end{array}$$

$$\displaystyle\Bigl{\{}d\in\mathbb{R}^{2n_{B}}\;:\;\nabla f_{j}(\gamma)^{\top}\,d<0,\;j=1,2\Bigr{\}}=\varnothing,$$

$$\begin{array}{rl}\underset{\gamma \in \Gamma }{\min }& f_w(\gamma ):=\omega f_1(\gamma )+(1-\omega )f_2(\gamma )\\ \text{s.t. }& u(X(\gamma ))\text{ solves the discretized state equation }(\text{ref}stateequation),\end{array}$$

$$\displaystyle\max\Bigl{\{}t\!=\!\frac{1}{2^{\ell}}\;:\;\ell\in\mathbb{N}_{0},\,f_{\omega}(\gamma^{(k)}\!+\!td^{(k)})\leq f_{\omega}(\gamma^{(k)})\!+\!\beta\,t\nabla f_{\omega}(\gamma^{(k)})^{\top}d^{(k)}\Bigr{\}};$$

$$\underset{\rho \in \mathbb{R},d\in {\mathbb{R}}^{2n_B}}{\min }$$

$$\nabla f_j({\gamma }^{(k)}{)}^{\top }d\le \rho ,j=1,2.$$

$$f_j({\gamma }^{(k)}+t{d}^{(k)})\le f_j({\gamma }^{(k)})+\beta t\nabla f_j({\gamma }^{(k)}{)}^{\top }{d}^{(k)},j=1,2.$$

Full text

Key words: biobjective shape optimization, shape gradients, probability of failure, descent algorithms

**MSC (2010): ** 90B50, 49Q10, 65C50, 60G55

1 Introduction

The optimization of the design of mechanical structures is a central task in mechanical engineering. If the material for a component is chosen and the use cases are defined, implying in particular the mechanical loads, then the central task of engineering design is to define the shape of the component. Among all possible choices, those shapes are preferred that guarantee the desired functionality at minimal cost. The functionality, however, is only guaranteed if the mechanical integrity of the component is preserved. The fundamental design requirements of functional integrity and cost are almost always in conflict, which makes mechanical engineering an optimization problem with at least two objective functions to consider.

Mathematically, the task of choosing the shape of a structure is formulated by the theory of shape optimization, see e.g., [1, 9, 32, 49] for an introduction. We thus consider admissible shapes $\Omega\subseteq\mathbb{R}^{p}$, $p=2,3$, along with an objective function $f(\Omega)$ which returns lower values for better designs. The task then is to find an admissible shape $\Omega^{*}\in\arg\min f(\Omega)$. The existence of optimal shapes has been studied in [10, 26, 32] – for the specific objective function $f$ of compliance see [1]. On the algorithmic side, the adjoint approach to shape calculus has led to efficient strategies to calculate shape gradients, see e.g., [12, 16, 22, 23, 34, 46, 47, 49]. While theory and numerical algorithms of shape calculus are highly developed mathematically, most publications in the field neither deal with multiobjective optimization problems, nor directly consider mechanical integrity as one of the objective functions, see [32, 2, 20, 39] for some remarkable exceptions.

In mechanical engineering, mechanical integrity is one of the central objectives, see e.g., [3]. However, if objectives like the ultimate load that the structure can bear or the fatigue life of a component are formulated deterministically, then the objective function is in general non differentiable as it only depends on the point of maximal stress. In numerical shape optimization this would lead to shape gradients concentrated on a single node, resulting in highly instable optimization schemes. At the same time, as material properties are subject to considerable scatter, a deterministic approach is not realistic. To overcome these two shortcomings, an alternative probabilistic approach to mechanical integrity has been proposed by some of the authors and others [5, 6, 28, 29, 30, 43, 44, 45], which has a smoothing effect on the singularities that are typical for deterministic models. Note that the probabilistic description of the ultimate strength of ceramics has become a standard in material engineering since the ground breaking work of Weibull, see e.g. [3, 8, 38, 42, 51].

In practice, there usually is a trade-off between the mechanical integrity of a structure and its volume (cost), since an improved mechanical integrity usually comes at the cost of a larger volume. Instead of presetting a fixed bound on the allowable volume, the trade-off between these two conflicting goals can be analyzed in a biobjective model. Other relevant objective functions may be, for example, the minimal buckling load of a structure or its minimal natural frequency, see, for example, [32]. For a general introduction into the field of multiobjective optimization we refer to [21, 36]. In the context of shape optimization problems, two major solution approaches can be distinguished: Metaheuristic and, in particular, evolutionary algorithms are widely applicable solution paradigms that do not utilize the particular structure of a given problem [11, 14, 15, 52]. However, in combination with expensive numerical simulations such approaches tend to be inefficient. On the other hand, gradient-based algorithms [17, 24, 53] require efficient gradient computations and are often applied in the context of adjoint approaches and using weighted sum scalarizations of the objective functions. See [41] for a comparison.

In this paper, we suggest a biobjective PDE constrained shape optimization problem for the simultaneous optimization of the mechanical integrity and the cost of a ceramic component. Section 2 is devoted to a formal introduction of the problem, including a review of Weibull type models for the probability of failure and existence results for Pareto optimal shapes. The numerical implementation is based on a first discretize then optimize approach using Lagrangian finite elements, which is detailed in Section 3. Section 3 also contains an introduction to gradient-based optimization strategies for biobjective problems and some details on their efficient implementation. The approach is validated at 2D ceramic components in Section 4, and perspectives for future research are suggested in Section 5.

2 Biobjective Shape Optimization (of Ceramic Structures)

is not available in this case. This motivates the formulation of a biobjective shape optimization problem where mechanical integrity and volume (cost) are considered simultaneously as equitable objectives.

In this section, we first introduce a set of admissible (feasible) shapes and review the state equations that model the physical behavior of a shape under external forces according to the linear elasticity theory (Section 2.1). The considered objective functions, the intensity measure modelling the mechanical integrity, and the volume of the shape, are formally introduced in Sections 2.2 and 2.3, respectively. The overall problem is formulated as a biobjective optimization problem in Section 2.4, and the existence of Pareto optimal solutions is shown in Section 2.5.

2.1 Admissible Shapes and State Equation

We follow the description from [5, 6] and consider a compact body (also referred to as component or shape) $\Omega\subset\mathbb{R}^{p}$, $p=2,3$ with Lipschitz boundary that is filled with ceramic material. Furthermore, we assume that the boundary $\partial\Omega$ of $\Omega$ is subdivided into three parts with nonempty relative interior,

[TABLE]

$\partial\Omega_{D}$ describes the part of the boundary where the Dirichlet boundary condition holds, $\partial\Omega_{N_{\text{fixed}}}$ the part where surface forces may act on and $\partial\Omega_{N_{\text{free}}}$ the part of the boundary that can be modified in an optimization approach. It is assumed to be force free for technical reasons [6].

Since all feasible shapes have to coincide in $\Omega_{D}$ and in $\Omega_{N_{\text{fixed}}}$, it is natural to restrict the analysis to subsets of a sufficiently large bounded open set $\widehat{\Omega}\subset\mathbb{R}^{p}$ that satisfies $\partial\Omega_{D}\subseteq\partial\widehat{\Omega}$ and $\partial\Omega_{N_{\text{fixed}}}\subseteq\partial\widehat{\Omega}$ (see Figure 1). We additionally assume that $\widehat{\Omega}$ satisfies the cone property for a given angle $\theta\in(0,\pi/2)$ and radii $r,l>0$, $r\leq l/2$, i.e.,

[TABLE]

where $C(\zeta_{x},\theta,l):=\{c\in\mathbb{R}^{p}\,:\,\|c\|<l,\,c^{\top}\zeta_{x}>\|c\|\cos(\theta)\}$ is a truncated circular cone oriented along $\zeta_{x}$ with height $l$ and opening angle $2\theta$, and $B(x,r)\subset\mathbb{R}^{p}$ is an open ball of radius $r$ centered at $x$. Now the set of admissible shapes $\mathcal{O}^{\text{ad}}\subset{\mathcal{P}}(\mathbb{R}^{p})$ can be defined as

[TABLE]

Ceramic components behave according to the linear elasticity theory [38]. The state equation can be described as an elliptic partial differential equation, see, e.g., [7]. More precisely, we get the state equation which describes the reaction of the ceramic component to external forces as a partial differential equation:

[TABLE]

Here, $\hat{n}(x)$ is the outward pointing normal at $x\in\partial\Omega$, which is defined almost everywhere on $\partial\Omega$ given that $\partial\Omega$ is piecewise differentiable. Furthermore, let $\tilde{f}\in L^{2}(\Omega,\mathbb{R}^{p})$ be the volume forces and $\tilde{g}\in L^{2}(\partial\Omega_{N_{\text{fixed}}},\mathbb{R}^{p})$ the forces acting on the surface $\partial\Omega_{N_{\text{fixed}}}$, e.g., the tensile load. The displacement caused by the acting forces is given by $u\in H^{1}(\Omega,\mathbb{R}^{p})$, where $H^{1}(\Omega,\mathbb{R}^{p})$ is the Sobolov space of $L^{2}(\Omega,\mathbb{R}^{p})$-functions with weak derivatives in $L^{2}(\Omega,\mathbb{R}^{p\times p})$. The linear strain tensor $\varepsilon\in L^{2}(\Omega,\mathbb{R}^{p\times p})$ is given by $\varepsilon(u(x)):=\frac{1}{2}(Du(x)+(Du(x))^{\top})$, where $Du$ is the Jacobi matrix of $u$. It follows for the stress tensor $\sigma\in L^{2}(\Omega,\mathbb{R}^{p\times p})$ that $\sigma(u(x))=\lambda\,\text{tr}(\varepsilon(u(x)))I+2\mu\varepsilon(u(x))$, where $\lambda,\mu>0$ are the Lamé constants derived from Young’s modulus $E$ and Poisson´s ratio $\nu$ as $\lambda=\frac{\nu E}{(1+\nu)(1-2\nu)}$ and $\mu=\frac{E}{2(1+\nu)}$.

From a numerical perspective, a variational formulation of the state equation (2.2) is usually preferred, see, e.g., [5, 6]. This still guarantees a unique weak solution $u$, see [19]. Thus, $u$ is uniquely defined by the shape $\Omega$ [19], and we will equivalently write $\sigma(Du(x)):=\sigma(u(x))$ for $x\in\Omega$ to highlight that $\sigma$ depends on the Jacobi matrix of $u$.

2.2 Probability of Failure

The primary objective function, the mechanical integrity of the ceramic component, is modelled based on the probability of failure analogous to [5, 6, 8, 51]. For the sake of completeness this is briefly summarized in the following.

We want to optimize the reliability of a ceramic body $\Omega$, i.e., its survival probability, by minimizing its probability of failure under tensile load. In that sense failure means that the ceramic body breaks under the tensile load due to cracks. Such cracks occur as a result of small faults in the material caused by the sintering process. To understand the mechanics of cracks, three types of crack opening are considered, see [31] and Figure 2a for an illustration. They are referred to as Modes I, II and III, respectively, and relate to different loads. Note that in the two-dimensional case, only Modes I and II can occur. We refer to [31] for a detailed introduction into this topic.

The stresses and strains close to a crack are represented by the crack-tip field which depends on the respective crack opening modes. It is described locally by a two-dimensional model, see Figure 2b for an illustration. With $K_{\text{I}},K_{\text{II}}$ and $K_{\text{III}}$ being the stress-intensity factors (also called $K$-factors) corresponding to Modes I, II, and III, respectively, one can describe the crack-tip field $\sigma$ locally according to linear fracture mechanics as

[TABLE]

Here, $r$ is the distance to the crack tip, and $\phi$ the angle w.r.t. the $x_{1}$-axis (aligned with the crack plane), see Figure 2b. The functions $\tilde{\sigma}^{\text{I,II,III}}(\phi)$ are known functions of the angle $\phi$, see again [31], and $R(r,\phi)$ is a regular function of the considered position in $x\in\Omega$ that is independent of the crack. Note that in the two-dimensional case, Mode III is omitted from (2.3) since it does not exist. Moreover, experimental evidence has shown that Mode I, which relates to tensile and compressive load, is the most relevant for the failure of ceramic structures [8], see [31] for approaches for multi-mode failure. We will thus focus on $K_{\text{I}}$ in the following as the driving parameter for crack development under tensile load.

In order to evaluate $K_{\text{I}}$ analogous to [6], we adopt the concept of equivalent circular discs to represent different crack shapes and crack sizes, and hence assume that the cracks are penny shaped. Then a particular crack can be identified by its configuration

[TABLE]

where $x\in{\Omega}$ is its location, $a\in(0,\infty)$ its radius, and $n\in S^{p-1}$ its orientation ($S^{p-1}$ denotes the unit sphere in $\mathbb{R}^{p}$). $\mathcal{C}$ is called the crack configuration space. Given a crack $(x,a,n)\in\mathcal{C}$, $K_{\text{I}}$ can be computed as a function of the radius $a$ and of the tensile load $\sigma_{n}(Du(x))$ in the normal direction $n$ of the stress plane at the crack location $x$ as

[TABLE]

see, e.g., Table 4.1 in [31]. Following [6] we set

[TABLE]

Note that negative values of $\sigma_{n}(x)$ correspond to compressive loads which can be ignored in the analysis of crack development, see Figure 2a above.

A crack $(x,a,n)\in\mathcal{C}$ becomes critical, i.e., a fracture occurs and the material fails, if $K_{\text{I}}$ exceeds a material-specific critical value $K_{\text{I}c}$ (the ultimate tensile strength of the material). Note that (2.4) implies that all cracks with radius

[TABLE]

are critical. We denote the set of critical configurations by

[TABLE]

and want to minimize the probability of finding a crack with configuration in $A_{c}$.

Following [5, 6], we assume that the parameters $(x,a,n)$ are random (i.e., they are not deterministically given by the sintering process), that the cracks are statistically homogeneously distributed in $\Omega$, and that their orientations are isotropic. Let $A\subseteq\mathcal{C}$ be a measurable subset of the configuration space. Then under quite general assumtions the random number $N(A)$ of cracks in $A$ is Poisson distributed (see [33, 50]), and hence $N(A)$ is a Poisson point process. It follows that $P(N(A)=k)=e^{-\upsilon(A)}\frac{\upsilon(A)^{k}}{k!}\sim Po(\upsilon(A))$, where $\upsilon$ is the (Radon) intensity measure of the process. Recall that a component fails if $N(A_{c})>0$. Given a displacement field $u\in H^{1}(\Omega,\mathbb{R}^{p})$, we can now write the survival probability of the component $\Omega$ as

[TABLE]

Hence, to maximize the survival probability of a component $\Omega$ we need to minimize the intensity measure $\upsilon$. Since only cracks $(x,a,n)$ with radius $a>a_{c}$ need to be considered (c.f. (2.5) above), [5, 6] determine the intensity measure as

[TABLE]

with $\text{d}x$ the Lebesgue measure on $\mathbb{R}^{p}$, $\text{d}n$ the surface measure on $S^{p-1}$, and $\text{d}\upsilon_{a}(a)=c\cdot a^{-\tilde{m}}\text{d}a$ being a positive Radon measure modelling the occurrence of cracks of radius $a$ in $\Omega$ ($c>0$ and $\tilde{m}\geq\frac{3}{2}$ are positive constants). Note that for $p=3$ the $\Gamma$-function takes the value $\Gamma(\frac{3}{2})=\frac{\sqrt{\pi}}{2}$ and for $p=2$ we obtain $\Gamma(1)=1$. With $m:=2(\tilde{m}-1)\geq 1$ and using again (2.5) the inner integral can be evaluated, yielding

[TABLE]

where $\sigma_{0}$ is an appropriately chosen positive constant. As highlighted in [5, 6], this is in accordance with the statistical model introduced by Weibull [51]. In this context, the parameter $m$ is referred to as Weibull module and typically assumes values between $5$ and $25$.

Summarizing the discussion above, we define our primary objective function $f_{1}:\mathcal{O}^{\text{ad}}\rightarrow\mathbb{R}$ as

[TABLE]

and refer to it as intensity measure, modelling the probability of failure (PoF) of the component $\Omega$. Recall that $u(\Omega)$ is uniquely defined by $\Omega$ and thus $f_{1}(\Omega)$ is completely defined by the shape $\Omega$ (given fixed boundary conditions $\tilde{f},\tilde{g}$).

2.3 Material Consumption

Improving the intensity measure $f_{1}$ of a ceramic component (and hence its PoF) usually comes at the price of an increased material consumption, which is directly correlated with the cost of the component. In order to avoid excessively expensive solutions, classical approaches thus set a predetermined bound on the allowable volume of the shape $\Omega$ (see, e.g., [5, 6]). We follow a more general approach in this manuscript and interpret the volume (and hence the cost) of the component as an equitable second objective function. This facilitates, in particular, the analysis of the trade-off between these two criteria and supports the engineer in finding a preferable design. We thus define $f_{2}:\mathcal{O}^{\text{ad}}\rightarrow\mathbb{R}$ as the volume of a shape $\Omega\in\mathcal{O}^{\text{ad}}$ given by

[TABLE]

2.4 Biobjective Optimization

When multiple conflicting goals are relevant in an optimization problem, a common approach is to use a weighted sum of the individual objectives as an overall objective function and then resort to classical optimization algorithms. The advantages and also the shortcomings of this so-called weighted sum scalarization are discussed in the following section, see also [21]. Particularly when choosing fixed weights, this method is of limited applicability. While fixed weights may represent the preferences of one decision maker, another decision maker may have other preferences, i.e., other weights. Moreover, the objective ranges and the scales of the objectives may be very different or even incomparable, which generally leads to numerical difficulties.

Another common approach to handle multiple conflicting goals is to select one “most important” objective function to minimize, e.g., the probability of failure, and set upper bounds on the acceptable objective function values of the other objective functions. In our case this would imply a constraint on the allowable material consumption, see, e.g., [5, 6]. This approach is referred to as $\varepsilon$-constraint scalarization, see again [21] for a general discussion of this topic. In addition to the numerical difficulties that may arise from adding potentially complicating constraints to the problem formulation, this approach has similar drawbacks as the weighted sum scalarization: The selection of meaningful upper bound values may be difficult, and trade-off information is ignored.

A more general approach is to formulate a multiobjective optimization problem, and hence to compute a set of relevant solution alternatives rather than one single “optimal” solution. By providing a set of solution alternatives the decision maker can not only choose a solution that aligns the most with his preferences, but he can also inspect the trade-off between alternative solutions and can adjust his preferences accordingly. A decision maker may, for example, prefer reliability over volume, but looking into the trade-off between solution alternatives there may be a solution that is some small percentage worse w.r.t. the reliability while it is a lot better regarding the volume. This may lead to a re-evaluation of the decision maker’s preferences.

With our two objective functions “intensity measure” ($f_{1}$, modeling the PoF) and “volume” ($f_{2}$), the following biobjective shape optimization problem arises:

[TABLE]

where $f_{1}$ and $f_{2}$ are defined according to Sections 2.2 and 2.3 above. Note that only $f_{1}$ depends on the displacement field $u(\Omega)$.

We call $f=(f_{1},f_{2}):{\mathcal{O}^{\text{ad}}}\longrightarrow\mathbb{R}^{2}$ the biobjective function vector and $\mathbb{R}^{2}$ the objective space. Let $Z:=f(\mathcal{O}^{\text{ad}})\subset\mathbb{R}^{2}$ denote the set of all feasible outcome vectors in the objective space, i.e., the set of all outcome vectors that are images of admissible shapes $\Omega\in{\mathcal{O}^{\text{ad}}}$. In contrast to single objective optimization we have to define optimality in the presence of two objectives, since there is no natural order on $\mathbb{R}^{2}$. For two shapes $\Omega_{1},\Omega_{2}\in{\mathcal{O}^{\text{ad}}}$, let $z^{1}=f(\Omega_{1})$ and $z^{2}=f(\Omega_{2})$ be the respective outcome vectors in $Z$. We write

[TABLE]

Note that $z^{1}\leqslant z^{2}$ implies that $z^{1}_{j}\leq z^{2}_{j}$ for $j=1,2$ with at least one strict inequality. We use the notation

[TABLE]

The notations $\mathbb{R}^{2}_{\geqq}$, $\mathbb{R}^{2}_{>}$, $\mathbb{R}^{2}_{\leqslant}$, $\mathbb{R}^{2}_{\leqq}$ and $\mathbb{R}^{2}_{<}$ are used accordingly.

We say that $z^{1}$ dominates $z^{2}$ if and only if $z^{1}\leqslant z^{2}$, i.e., if and only if $z^{1}\in z^{2}+\mathbb{R}^{2}_{\leqslant}$. An outcome vector $\bar{z}\in Z$ is called nondominated if there is no other outcome vector $z\in Z$ such that $z\leqslant\bar{z}$. Accordingly, an admissible shape ${\Omega}_{P}\in{\mathcal{O}^{\text{ad}}}$ is called Pareto optimal or efficient, if there is no other admissible shape $\Omega\in{\mathcal{O}^{\text{ad}}}$ such that $f(\Omega)\leqslant f({\Omega}_{P})$. We are mainly interested in Pareto optimal shapes since these are precisely those shapes that can not be improved in one objective without deterioration in the other objective. The set of all Pareto optimal shapes is called the Pareto front and denoted by $\mathcal{O}^{\text{ad}}_{P}$. Similarly, the set of all nondominated outcome vectors $Z_{N}:=f(\mathcal{O}^{\text{ad}}_{P})$ is referred to as the nondominated front in the objective space.

As in single-objective optimization, one often has to resort to local minima if the underlying optimization problem is nonconvex (and difficult). In the biobjective setting, an admissible shape ${\Omega}_{\ell P}\in{\mathcal{O}^{\text{ad}}}$ is called locally Pareto optimal or locally efficient, if there is a neighborhood $\mathcal{N}\subseteq\mathcal{O}^{\text{ad}}$ of $\Omega_{\ell P}$ such that there is no other admissible shape $\Omega\in\mathcal{N}$ with $f(\Omega)\leqslant f({\Omega}_{\ell P})$.

2.5 Existence of Pareto Optimal Shapes

In order to prove the existence of Pareto optimal shapes, we consider the weighted sum scalarization of problem (2.8) in which, given a weight $\omega\in(0,1)$, the two objective functions are combined into one single weighted sum objective:

[TABLE]

It is a well-known fact that every optimal solution of problem (2.9) is Pareto optimal for problem (2.8), see, e.g., [21].

Theorem 1**.**

If the crack size measure has the non decreasing stress hazard property (see [6] for a formal definition), then the set $\mathcal{O}^{\text{ad}}_{P}$ is non-empty.

Proof.

Suppose that $\omega\in(0,1)$ is chosen arbitrarily, but fixed. Then the weighted sum objective can be evaluated as

[TABLE]

Thus, the incorporation of $f_{2}$ into the scalarized objective function corresponds to the addition of a constant term in the shape integral of $f_{1}$. This does not affect the convergence analysis of [6], which is based on convexity of the integrand in $Du$, see [10, 26]. We can conclude that the weighted sum scalarization has an optimal solution for every $\omega\in(0,1)$. Since every such solution is Pareto optimal for (2.8), the result follows. ∎

3 Numerical Implementation

To actually compute locally Pareto optimal shapes, we adopt the finite element discretization implemented in [5] for two-dimensional instances (i.e., $p=2$). In this implementation, the shapes $\Omega\in\mathcal{O}^{\text{ad}}$, the state equation (2.2), the objective functions $f_{1}$ and $f_{2}$ and their gradients are discretized. Standard Lagrangian finite elements are used for the discretization of the state equation (2.2), and all integrals are calculated using numerical quadrature. The discretized shape gradients are obtained by an adjoint approach to reduce computational costs. We refer to [5] for a detailed description.

3.1 Geometry Definition and Finite Element Mesh

The two-dimensional shapes $\Omega\in\mathcal{O}^{\text{ad}}\subset{\mathcal{P}}(\mathbb{R}^{2})$ are discretized by an $n_{x}\times n_{y}$ mesh $X:=X^{\Omega}=(X^{\Omega}_{ij})_{n_{x}\times n_{y}}$ (we write $X_{ij}:=X_{ij}^{\Omega}\in\mathbb{R}^{2}$ for short) using tetrahedrons, with $n_{x},n_{y}\in\mathbb{N}$ being the number of grid points in $x$ and $y$ direction, respectively. Given a shape $\Omega\in\mathcal{O}^{\text{ad}}$ and its discretization $X$, the objective function values $f_{1}(X)$ and $f_{2}(X)$ as well as their gradients $\nabla f_{1}(X)$ and $\nabla f_{2}(X)$ are computed using the implementation of [5].

For the optimization process, we fix the $x$-component of all grid points to equidistant values $x_{1},\dots,x_{n_{x}}$, and we only consider the $y$-components of those grid points that define the boundary of the shape to avoid deformation of the inner mesh structure. Note that this reformulation reduces the number of optimization variables from $2n_{x}n_{y}$ to $2n_{x}$. As a consequence, feasible shapes can alternatively be represented by a shape parameter $\varrho$ containing, for every relevant $x$-coordinate, the $y$-coordinate of the meanline $\varrho^{\text{ml}}_{i}\in\mathbb{R}$ of the shape, and the thickness $\varrho^{\text{th}}_{i}\in\mathbb{R}_{>}$ of the shape, $i=1,\dots,n_{x}$. Given a feasible shape represented by $\varrho:=(\varrho^{\text{ml}},\varrho^{\text{th}})\in\mathbb{R}^{2n_{x}}$ with $\varrho^{\text{th}}\in\mathbb{R}_{>}^{n_{x}}$, an associated mesh representation $X$ can be obtained using

[TABLE]

To further reduce the computational burden and to obtain smoother shapes, the shape parameters $\varrho^{\text{ml}}\in\mathbb{R}^{n_{x}}$ and $\varrho^{\text{th}}\in\mathbb{R}^{n_{x}}_{>}$ are modelled using B-splines. Let $n_{B}\in\mathbb{N}$, with $n_{B}<n_{x}$, be the number of B-spline basis functions, and let $\{\vartheta_{j}:\mathbb{R}\rightarrow\mathbb{R}_{\geq},\;j=1,\dots,n_{B}\}$ be a B-spline basis (see, e.g., [40]). Feasible shapes are then represented by B-spline coefficients $\gamma:=(\gamma^{\text{ml}},\gamma^{\text{th}})\in\mathbb{R}^{2n_{B}}$. The corresponding meanline and thickness values can be computed using the auxiliary functions

[TABLE]

These auxiliary meanline and thickness functions are then evaluated at the fixed $x$-coordinates of the gridpoints which yields

[TABLE]

Using the B-spline coefficients $\gamma=(\gamma^{\text{ml}},\gamma^{\text{th}})\in\mathbb{R}^{2n_{B}}$ as optimization variables yields a further reduction of the number of variables to $2n_{B}$. Moreover, the B-spline representation leads to an implicit regularization and smoothing of the represented shapes. In the following, we denote the set of feasible shape parametrizations by $\Gamma\subseteq\{(\gamma^{\text{ml}},\gamma^{\text{th}})\in\mathbb{R}^{2n_{B}}\}$.

To evaluate the objective functions $f_{j}(\gamma)$ and their gradients $\nabla f_{j}(\gamma)={\partial f_{j}}/{\partial\gamma}$, $j=1,2$, w.r.t. the new parametrization of shapes based on B-spline parameters $\gamma$, while still using the implementation of [5], we compute an associated grid $X$ using first (3.2) and then (3.1). While the resulting objective function values can be used immediately in the optimization process, the gradients computed w.r.t. the grid $X$ need to be translated to the space of B-spline coefficients, i.e.,

[TABLE]

The numerical computation of gradients of $f_{j}$, $j=1,2$, w.r.t. a B-spline representation $\gamma$ of a feasible shape $\Omega$ is thus based on a two-step projection of $\gamma$ onto the original grid $X$. The thus computed gradients of $f_{1}$ (the intensity measure) were validated, using finite differences, at the sample shape shown in Figure 4a. The validation is based on a grid $(X_{ij})_{41\times 7}$, i.e., $n_{x}=41$ and $n_{y}=7$. Consequently, for the corresponding meanline and thickness representation we have $\varrho=(\varrho^{\text{ml}},\varrho^{\text{th}})\in\mathbb{R}^{82}$, where $\varrho^{\text{th}}\in\mathbb{R}^{41}_{>}$. Moreover, we used a B-spline basis with five basis functions, i.e., $n_{B}=5$ and $\gamma=(\gamma^{\text{ml}},\gamma^{\text{th}})\in\mathbb{R}^{10}$. We computed all ten partial derivatives w.r.t. $\gamma$ via the respective transformations to the grid representation and compared them with finite differences. The results of this comparison, i.e., the absolute values of the differences between computed derivatives and finite differences, are shown in Figure 3a and 3b for the meanline and thickness parameters, respectively. The figures indicate in all cases that, when the finite differences are evaluated for decreasing values of the increment $\varepsilon$, then they correspond well to the computed gradients.

3.2 Pareto Critical Solutions

Given the parametrization of admissible shapes described in Section 3.1, the biobjective optimization problem (2.8) can now be restated as

[TABLE]

Recall that only $f_{1}$ depends on the displacement field $u(X)$.

Since derivative information is available, necessary optimality conditions can be formulated that generalize the concept of critical points from single-objective optimization. Towards this end, we omit the constraints implied by the parametric representation of admissible shapes to keep the exposition simple. All constraints will be handled implicitly in the numerical tests described in Section 4 below. Assuming that both objective functions are continuously differentiable a necessary condition for a solution $\gamma\in\mathbb{R}^{2n_{B}}$ to be locally Pareto optimal is that

[TABLE]

i.e., there does not exist a direction $d\in\mathbb{R}^{2n_{B}}$ that is a descent direction for both objectives. If $\gamma^{*}\in\mathbb{R}^{2n_{B}}$ satisfies this condition we call it a Pareto critical shape.

In this work, we aim at the efficient computation of Pareto critical shapes that, ideally, approximate the Pareto front. Since derivative information can be obtained for both objective functions, we select solution methods that efficiently utilize this information and that can be adopted such that a meaningful representation of a Pareto critical front is obtained. As two fundamental approaches in this category, a parametrized weighted sum method and a biobjective descent algorithm are chosen and explained in Sections 3.3 and 3.4, respectively. Their performance in the context of 2D shape optimization problems is compared in Section 4.

3.3 Weighted Sum Method

Maybe the easiest way to compute a representation of the Pareto front is to iteratively solve weighted sum scalarizations (2.9) with varying weights. The weighted sum scalarization of problem (3.4) can be restated as

[TABLE]

where $\omega\in(0,1)$ is the weight specifying the relative importance of $f_{1}$ and $f_{2}$, respectively. Recall that every solution of the weighted sum scalarization (3.6) is Pareto optimal for (2.8) [21]. A disadvantage of the weighted sum method is, however, that only solutions that map to the convex hull $\operatorname{conv}(Z)$ of the image set $Z=f(\Gamma)$ in the objective space can be found, and thus relevant compromise solutions in nonconvex areas of the nondominated front may be missed. Moreover, [13] showed at simple biobjective test instances that evenly distributed weights do in general not lead to well distributed outcome vectors in the objective space. This is particularly problematic if the considered objective function values are of largely different magnitude, which is the case here. In order to obtain solutions that are consistent with the preferences expressed by $\omega$, we thus normalize the objective functions by using appropriate scaling factors $c_{1},c_{2}>0$, and replace $f_{1}$ and $f_{2}$ in (3.6) by $c_{1}f_{1}$ and $c_{2}f_{2}$, respectively.

Despite the difficulties mentioned above, the weighted sum method is usually well-suited to efficiently compute at least a rough approximation of the Pareto front. For this purpose, problem (3.6) is solved iteratively for varying weights (in our case, we choose $\omega\in\{0.2,0.25,0.3,\dots,0.9\}$ since numerical experiments showed that this yields meaningful trade-offs). Each single objective optimization problem (3.6) is then individually solved using a classical gradient descent algorithm with stepsizes determined according to the Armijo rule, see, for example, [4].

Under appropriate assumptions, the gradient descent algorithm in the inner loop of Algorithm 1 converges to a critical point of (3.6), see, e.g., [4]. In our implementation, the inner loop is also terminated when a prespecified maximum number of iterations is reached. However, in this case there is no guarantee that the final iterate is close to a Pareto critical solution.

Note that a critical point of the weighted sum scalarization (3.6) is necessarily Pareto critical for the biobjective shape optimization problem (3.4), while the converse is not true in general. This has some correspondence to the fact that global optimal solutions of a weighted sum scalarization (3.6) are always Pareto optimal, while nonconvex problems may have Pareto optimal solutions that are not optimal for any weighted sum scalarization (3.6), see, e.g., [21].

Note also that the search direction $d^{(k)}=-\nabla f_{\omega}(\gamma^{(k)})$ does not necessarily satisfy $\nabla f_{j}(\gamma^{(k)})^{\top}d^{(k)}<0$, $j=1,2$, in all iterations. In other words, one objective function may deteriorate during the optimization process if only the other objective function compensates for this.

3.4 Biobjective Descent Algorithm

Different from the weighted sum method described above, biobjective descent algorithms – as a natural generalization of single-objective gradient descent algorithms – are potentially capable of finding every Pareto optimal solution, if only the starting solution is chosen appropriately. While this is a rather theoretical advantage, biobjective descent algorithms are indeed highly efficient in finding (or approximating) one Pareto critical solution without the necessity to specify preferences. However, if a representation of the complete Pareto front is sought, they need to be combined with other search strategies.

We adopt the multiobjective descent algorithm proposed in [24] (see also [25]) for the biobjective optimization problem (2.8). Similar approaches have been suggested in [17, 18, 27].

Biobjective descent algorithms iteratively improve both objective functions simultaneously. This is based on the observation that, if a solution $\gamma\in\mathbb{R}^{2n_{B}}$ is not Pareto critical according to (3.5), then there exists a direction $d\in\mathbb{R}^{2n_{B}}$ which is a descent direction for both objectives. Thus, if in an iterative solution method the current iterate $\gamma^{(k)}\in\mathbb{R}^{2n_{B}}$ is not Pareto critical, a direction of steepest biobjective descent $d^{(k)}\in\mathbb{R}^{2n_{B}}$ can be defined according to [24] as a direction solving the auxiliary optimization problem

[TABLE]

Problem (3.7) is a convex quadratic optimization problem with linear inequality constraints. Note that the term $\frac{1}{2}\|d\|^{2}$ in the objective function ensures that the problem is bounded, and that the solution $\rho=0$, $d=0$ is always feasible. Note also that the optimal value $\rho^{*}$ is negative if and only if $d^{*}\neq 0$, i.e., if a direction of steepest biobjective descent exists.

When a direction of steepest biobjective descent $d^{(k)}\neq 0$ is found, then we move from $\gamma^{(k)}$ into the direction $d^{(k)}$ to a new point $\gamma^{(k+1)}:=\gamma^{(k)}+t_{k}d^{(k)}$. The step length $t_{k}>0$ is computed using an Armijo-like rule. Towards this end, let $\beta\in(0,1)$ be a prespecified constant. Then a step length $t$ is accepted if it guarantees a sufficient biobjective descent in the sense that

[TABLE]

In order to compute an acceptable step length $t$, we iteratively test the values $(\frac{1}{2})^{\ell}$, $\ell=0,1,2,\dots$ until condition (3.8) is satisfied. A proof for the finiteness of this procedure is given in [24]. The overall method is summarized in Algorithm 2.

If $f_{1}$ and $f_{2}$ are continuously differentiable and $\varepsilon=0$, then Algorithm 2 converges to a Pareto critical solution [24]. A natural stopping condition for practical implementations, motivated by (3.5), is that $\|t_{k}d^{(k)}\|\leq\varepsilon$, with $\varepsilon>0$ a prespecified small constant.

In practice, we also terminate the algorithm when a prespecified maximum number of iterations is reached. In this case, the final solution has to be used with caution since the optimization procedure has generally not yet converged.

The choice of the search direction using problem (3.7) together with condition (3.8) implies that the iterates of Algorithm 2 satisfy $f(\gamma^{(k+1)})<f(\gamma^{(k)})$ for all $k=1,2,\dots$. In other words, the objective vector $f(\gamma^{(k+1)})$ in iteration $k+1$ is bounded above by the objective vector $f(\gamma^{(k)})$ of the previous iteration $k$, i.e., $f(\gamma^{(k+1)})\in f(\gamma^{(k)})-\mathbb{R}^{2}_{>}$.

Several alternative Pareto critical solutions (and hence trade-off information between them) can be obtained, for example, by varying the starting solution. We follow a different approach in our implementation that is somewhat similar to the weighted sum method, and that is based on the observation that the optimal solution of problem (3.7) (i.e., the direction of steepest biobjective descent) depends on the scaling of the objective functions $f_{1}$ and $f_{2}$. Thus, Algorithm 2 is executed repeatedly, using different scalings of the objective functions. In our implementation, we use a scaling parameter $s:=\bar{\omega}r^{\max}>0$ and replace $f_{2}$ by $sf_{2}$ in the optimization process, where the parameter $r^{\max}>0$ is chosen as the largest ratio between partial derivatives of $f_{1}$ and $f_{2}$, evaluated at the starting solution $\gamma^{(1)}$. Note that the latter aims at the constraints in problem (3.7) in the sense that they should be comparable, i.e., both objective functions should equally contribute to active constraints and thus influence the choice of the search direction. By varying the parameter $\bar{\omega}\in\{0.5,0.6,\dots,2\}$, we can compute different solutions starting from the same initial shape. Note that the volume of the solutions can be expected to increase with larger values of $\bar{\omega}$.

Note also that the resulting parametric version of Algorithm 2 is fundamentally different from the weighted sum method in Algorithm 1 in the way the search directions are chosen and in the way the iterates converge to a Pareto critical solution.

3.5 Scalar Products and Gradients in Shape Optimization

The performance of Algorithms 1 and 2 depends largely on the choice of the search direction, which is computed based on the discretized gradients $\nabla f_{j}(\gamma)$, $j=1,2$. Michor and Mumford [35] showed that (continuous) shape gradients calculated with respect to the ordinary $L^{2}$-scalar product lead to an ill defined notion of the distance of two shapes, as the infimum over all deformation path lengths is zero. They suggest a modified scalar product given by

[TABLE]

and show that this indeed leads to a well defined Riemannian metric on the shape space. Here, $h,k$ are two vector fields in normal direction to the boundary of $\partial\Omega$, $\text{d}A$ is the induced surface measure, $\kappa$ is the scalar curvature of the surface, and $\xi>0$ is a regularization parameter. In practice, this corresponds to a transformation of function values on $\partial\Omega$ that, given some function $g:\partial\Omega\rightarrow\mathbb{R}^{2}$, can be described by $g_{\xi}(x)=\frac{g(x)}{1+\xi\kappa^{2}(x)}$ for $x\in\partial\Omega$.

We adopt a discretized version of this concept in the numerical implementation of shape gradients for both objectives $f_{j}$, $j=1,2$. More precisely, a discretized scalar curvature $\kappa$ is computed at grid points on the boundary $\partial\Omega$, which is represented by a polygonal approximation induced by the shape parameters $(\varrho^{\text{ml}},\varrho^{\text{th}})\in\mathbb{R}^{2n_{x}}$, $\varrho^{\text{th}}\in\mathbb{R}_{>}^{n_{x}}$. Since the upper and lower boundary of the shape $\Omega$ may have a different curvature at the same $x$-coordinate value $x_{i}$ ($i\in\{1,\dots,n_{x}\}$), we have to compute the curvature for upper and lower boundary points separately. For the upper boundary, this is realized by comparing the normals $n^{\text{u}}_{i}$ and $n^{\text{u}}_{i+1}$ on two consecutive facets of length $l^{\text{u}}_{i}$ and $l^{\text{u}}_{i+1}$, respectively. Similarly, for the lower boundary we use $n^{\text{l}}_{i}$, $n^{\text{l}}_{i+1}$ and $l^{\text{l}}_{i}$, $l^{\text{l}}_{i+1}$, and obtain

[TABLE]

The upper and lower boundaries of the shape $\Omega$ are reconstructed from the meanline and thickness representation using the linear transformation $\varrho_{i}^{\text{u}}=\varrho_{i}^{\text{ml}}+\frac{1}{2}\varrho_{i}^{\text{th}}$ and $\varrho_{i}^{\text{l}}=\varrho_{i}^{\text{ml}}-\frac{1}{2}\varrho_{i}^{\text{th}}$, $i=1,\dots,n_{x}$. In other words, $(\varrho^{\text{u}},\varrho^{\text{l}})\in\mathbb{R}^{2n_{x}}$ is obtained from $(\varrho^{\text{ml}},\varrho^{\text{th}})\in\mathbb{R}^{2n_{x}}$, $\varrho^{\text{th}}\in\mathbb{R}_{>}^{n_{x}}$, as $(\varrho^{\text{u}},\varrho^{\text{l}})=M\,(\varrho^{\text{ml}},\varrho^{\text{th}})$, using an appropriate transformation matrix $M\in\mathbb{R}^{{2n_{x}}\times{2n_{x}}}$. This leads to a discretized representation of the respective boundaries by points $(x_{i},\varrho_{i}^{\text{u}})$ (upper boundary) and $(x_{i},\varrho_{i}^{\text{l}})$ (lower boundary), from which the $\kappa$ values can be computed according to (3.10).

Now (3.9) can be applied to the gradients of $f_{j}$ w.r.t. $(\varrho^{\text{u}},\varrho^{\text{l}})$, $j=1,2$, by multiplying the respective partial derivatives by

[TABLE]

Since we actually need the gradients of $f_{j}$ w.r.t. $\varrho=(\varrho^{\text{ml}},\varrho^{\text{th}})$, $j=1,2$, we additionally have to consider the linear tranformation $M$. Let $D_{\xi}=(d_{\xi,ij})_{2n_{x}\times 2n_{x}}\in\mathbb{R}^{{2n_{x}}\times{2n_{x}}}$ be a diagonal matrix with diagonal elements given by

[TABLE]

and set $\bar{D}_{\xi}:=M^{-1}\,D_{\xi}\,M$. Then we obtain the curvature adapted B-spline gradients as

[TABLE]

Note that for $\xi=0$ the matrix $\bar{D}_{0}$ is the identity matrix, and hence the $L^{2}$-gradient of $f_{j}$ w.r.t. $\gamma$, $j=1,2$, is recovered in this case, c.f. (3.3).

3.6 Control of Step Sizes

Large mesh deformations may cause numerical difficulties and thus have to be avoided. We thus limit the step size during the optimization procedure. Recall that the representation of feasible shapes, using meanline and thickness values $(\varrho^{\text{ml}}_{i},\varrho^{\text{th}}_{i})$ at fixed $x_{i}$ coordinates, $i=1,\dots,n_{x}$, implies that grid points can only move vertically. A natural choice for a maximum admissible step in one iteration of the optimization process is thus determined by the thickness of the shape, divided by the number $n_{y}$ of gridpoints in $y$-direction. Since in our case studies the shapes are fixed at the left boundary (i.e., at $x=x_{1}$) and hence their thickness is constant at $x_{1}$, we set

[TABLE]

i.e., to $80\%$ of the vertical distance between grid points on the left boundary of the initial shape. For a given search direction $d^{(k)}=(d^{\text{ml},(k)},d^{\text{th},(k)})\in\mathbb{R}^{2n_{B}}$ in iteration $k$ of the optimization algorithms, we check whether $\max_{i=1,\dots,2n_{B}}|d^{(k)}_{i}|\leq\delta^{\max}$. Otherwise, $d^{(k)}$ is scaled by a factor ${\delta^{\max}}/{\max_{i=1,\dots,2n_{B}}|d^{(k)}_{i}|}$. Then the step length $t\leq 1$ is computed according to the Armijo rule as indicated in Algorithms 1 and 2.

While $\delta^{\max}$ is derived from the mesh $X^{(1)}$, it still is a meaningful upper bound for a step $d^{(k)}$ in the B-spline representation. Indeed, if $\{\vartheta_{j},\,j=1,\dots,n_{B}\}$ is a B-spline basis and $\gamma^{(k)}=(\gamma^{\text{ml},(k)},\gamma^{\text{th},(k)})\in\Gamma$ is the current iterate, then the B-spline basis properties $\sum_{j=1}^{n_{B}}\vartheta_{j}(x)=1$ and $\vartheta_{j}(x)\geq 0$, $j=1,\dots,n_{B}$ (see, e.g., [40]) imply that, for all $i=1,\dots,n_{x}$,

[TABLE]

An analogous bound holds for the corresponding thickness parameters. Note that the above inequalities do in general not guarantee that all grid points of the corresponding mesh $X^{(k)}$ move by at most $80\%$, since this also depends on the current shape and the mutual movement of meanline and thickness values. In some situations it may thus be necessary to adapt this bound to a smaller value. However, this never occured in our numerical tests.

4 Case Studies

We consider 2D ceramic shapes made out of beryllium oxide (BeO) under tensile load. Therefore, we set Young’s modulus to $\texttt{E}=320\,\text{GPa}$ (see, e.g., [38]), Poisson’s ratio to $\nu=0.25$, and the ultimate tensile strength to $140\,\text{MPa}$, according to [48]. Weibull’s modulus is set to $m=5$, which is on the lower bound of industrial ceramics having $m$ between 5 and 30 as depending on the production process [37]. All considered shapes have a fixed length of $1.0\,\text{m}$ and a fixed height of $0.2\,\text{m}$ on the left and right boundaries. The shapes are fixed on the left boundary, where Dirichlet boundary conditions hold ($\partial\Omega_{D}$), and on the right boundary, where surface forces may act on and Neumann boundary conditions hold ($\partial\Omega_{N_{\text{fixed}}}$). The upper and lower boundaries are assumed to be force free ($\partial\Omega_{N_{\text{free}}}$). They can be modified within the optimization process. We set $\tilde{f}=0$ neglecting the gravity forces and $\tilde{g}=10^{7}\,\text{Pa}$, representing tensile load. Note that, in order to be consistent with 3D models, we define the force density w.r.t.  $\text{Pa}=\text{N}/\text{m}^{2}$ (and not w.r.t.  $\text{N}/\text{m}$). This is motivated by assuming a constant width of the 2D component of $1$ unit (i.e., $1\text{m}$). Then plane stresses and plane strains are obtained by neglecting Poisson effects in the third dimension.

The shapes are discretized by a $41\times 7$ grid (i.e., $n_{x}=41$ and $n_{y}=7$) using tetrahedrons as detailed in Section 3.1. The B-spline representation is based on $n_{B}=5$ basis functions. Moreover, the curvature regularization paramater is set to $\xi=10^{-4}$, see Section 3.5.

During the optimization process, we monitor the Euclidean norm of the update of the design variables in every iteration and stop when it is lower than $10^{-4}$. The implementation is realized in R version 3.4.2 and uses the adjoint finite element code of [5] as a subroutine.

4.1 A Straight Joint

In the first test case, a straight joint is sought that is fixed at the left side, while the tensile load acts on the right side. This is a particularly simple situation where the straight rod connecting from the left to the right can be expected to be optimal, with varying thickness depending on the trade-off between the intensity measure ($f_{1}$) and the volume ($f_{2}$). The optimization algorithms are challenged by providing a bended beam as a starting shape, which is clearly far from being optimal.

The starting shape is shown in Figure 4a, together with the $41\times 7$ tetrahedral discretization $X$. Its objective values are $f_{1}(X^{(1)})=0.769624$ (intensity measure) and $f_{2}(X^{(1)})=0.2$ (volume), respectively. The relatively high value for the intensity measure $f_{1}$ can be explained by the relatively high stresses that are illustrated in Figure 4b. Figure 4c shows that the B-spline representation based on only five basis functions leads to a rather inaccurate representation, particularly at the left and right boundary. This could be improved by fixing the slopes at the left and right boundary, however, at the price of a significantly reduced design space. Indeed, a majority of the Pareto critical shapes computed during our numerical tests do not have zero slopes at the left and right boundary, particularly in the case of the S-shaped joint considered in Section 4.2 below. Note that the smoothing induced by the B-spline representation in this case already leads to dominating objective values of $f_{1}(\gamma^{(1)})=0.453867$ and $f_{2}(\gamma^{(1)})=0.2$.

Results

Among all shapes with a fixed volume of $f_{2}(X)=0.2$, the straight rod shown in Figure 4d can be expected to have the mimimum possible intensity measure $f_{1}$. Indeed, the straight rod shown in Figure 4d achieves an objective value of $f_{1}(X)=0.00058$. Figures 4e and 4f show the results of the weighted sum method (Algorithm 1) with weight $\omega=0.8$ and of the biobjective descent algorithm (Algorithm 2) with scaling parameter $\bar{\omega}=1.8$. Both methods show a rather quick convergence (with the expected advantage for the biobjective descent algorithm) to solutions that are close to optimal. However, the solution of the biobjective descent algorithm seems to be a local solution with slightly higher stresses (and thus slightly higher objective value for $f_{1}$).

Figure 5 shows iteration histories of exemplary runs of the weighted sum method (Algorithm 1) and of the biobjective descent algorithm (Algorithm 2), respectively. It nicely illustrates that, in contrast to the biobjective descent algorithm, the weighted sum method permits iterations where one objective function deteriorates while the weighted sum objective is still decreasing. This may, in certain situations, help to overcome local Pareto critical solutions. On the other hand, the weighted sum method may get stuck in local minima as well. Indeed, independent of the chosen weight, the histories of the weighted sum method have a similar structure: First mainly the intensity measure (representing the PoF) is improved (since in early stages of the algorithm the gradient of $f_{1}$ is considerably larger than the gradient of $f_{2}$). Only at later stages of the algorithm, the volume is varied to a larger extent, depending on the given weight.

Note also that the final solution obtained with the biobjective descent algorithm largely depends on the starting solution, since the objective values can never deteriorate during the optimization process. Thus, when the starting solution has a volume of $f_{2}(X)=0.2$, then all Pareto critical shapes that can be computed with the biobjective descent algorithm have a volume of at most $0.2$, irrespective of the scaling.

Three shapes with progressively reduced volume (and hence lower cost) are shown in Figures 4g to 4i. As was to be expected, a lower cost comes at the price of a higher intensity measure (and hence higher PoF). A comparison between Figures 4h and 4g suggests that also for the low volume solutions, the weighted sum solutions slightly outperform the biobjective descent solutions.

Figure 6 summarizes the results of several optimization runs with varying weights (Algorithm 1) and varying scalings (Algorithm 2), respectively. The same starting solution was used in all cases, see Figure 4c. While the solution quality of the weighted sum method and of the biobjective descent algorithm is comparable, a clear advantage of the weighted sum method seams to be that it is not so much constrained by the (performance of the) starting solution. Indeed, the weighted sum solutions shown in Figure 6 span a large range of alternative objective values in the objective space and thus provide the decision maker with meaningful trade-off information and a variety of solution alternatives.

4.2 An S-Shaped Joint

A more complex situation is obtained when the left and right boundaries are not fixed at the same height, i.e., when an S-shaped joint is to be designed. In our tests, we fix the right boundary about $0.27\,\text{m}$ lower than the left boundary. The starting shape and its $41\times 7$ tetrahedral discretization $X$, that is used for all optimization runs, is shown in Figure 7a. Figure 7b highlights the stresses that are particularly strong towards the left boundary. The respective objective values are $f_{1}(X^{(1)})=1.520058$ (intensity measure) and $f_{2}(X^{(1)})=0.2$ (volume), respectively. As can be expected, the intensity measure (and hence also the PoF) is considerably higher than in the case of the straight joint discussed in Section 4.1. Despite the significant smoothing induced by the B-spline representation of the initial shape shown in Figure 7c, it has an even higher value of the intensity measure of $f_{1}(\gamma^{(1)})=1.910532$ (and hence a higher PoF value), while $f_{2}(\gamma^{(1)})=0.2$ remains constant.

Results

We observe that the resulting shapes resemble the profile of a whale. If we consider 1st principal stress of the stress tensor on the grid points of the initial shape resulting from tensile load, see Figure 7c, we observe an anti clockwise eddy in the left part of the joint. The hunch close to the left boundary of the optimized shapes gives room for the occurring stresses and therefore improves the intensity measure and, likewise, the PoF.

Note that, different from the case of the straight rod, we have no prior knowledge on the Pareto optimal shapes. For the solutions shown in Figures 7d and 7e, we can only guarantee that they are (approximately) Pareto critical, i.e., the respective optimization runs terminated due to the criticality test. Figure 7f shows a shape with a significantly higher volume of $f_{2}(X)=0.225906$, and with a largely improved intensity measure of $f_{1}(X)=0.196791$. This shape was obtained with the weighted sum method with weight $\omega=0.85$ after $150$ iterations. In this case, the algorithm terminated since it reached the maximum number of iterations and not due to convergence. We observed that all optimization runs of the weighted sum method with $\omega\geq 0.85$ were not converging in this setting. Thus in these cases it is not guaranteed, that the resulting solutions are Pareto critical. Note that, given a starting solution with a volume of $0.2$, this shape is not attainable with the biobjective descent algorithm.

However, there is no guarantee that the computed shapes are Pareto optimal. For example, the shape shown in Figure 7e obtained with the weighted sum method with weight $\omega=0.8$ achieves objective values of $f_{1}(X)=0.293853$ and $f_{2}(X)=0.188445$, and hence slightly dominates the shape shown in Figure 7d obtained with the biobjective descent algorithm with scaling parameter $\bar{\omega}=1.1$ that has objective values $f_{1}(X)=0.300996$ and $f_{2}(X)=0.188774$.

Figure 8 summarizes the results of several optimization runs of both Algorithms 1 and 2 in the objective space. Note that not all solutions of the weighted sum method lie on the convex hull of the computed points (and are thus not globally optimal for a weighted sum scalarization). In some cases, the biobjective descent algorithm also computes dominated points, while in other cases it found solutions that lie even below the convex hull of the weighted sum solutions (see, e.g., the result for $\bar{\omega}=0.5$ in Figure 8).

A larger range of alternative objective vectors is, as in the case of the straight rod, obtained with the weighted sum method. A cross-test between the two methods, where the final solution of Algorithm 1 was used as starting solution for Algorithm 2, confirms that local Pareto critical solutions were found for $\omega\leq 0.8$.

5 Conclusion and Outlook

We have developed a modelling and solution approach for biobjective PDE constrained shape optimization of ceramic components. The mechanical integrity of the component on one hand, and the cost of the component on the other hand, were considered as two pivotal optimization criteria. A probabilistic approach was used to assess the mechanical integrity (i.e., the reliability) of the component, which allows, in combination with a finite element discretization and an adjoint approach for gradient computations, the efficient calculation of derivative information. Approximations of the Pareto front were computed using two different approaches: (1) parametric weighted sum scalarizations in combination with a single objective gradient descent method, and (2) a biobjective descent algorithm with parametric scalings of the objective functions. Numerical results for 2D test cases visualize the trade-off between the reliability and the cost, and hence pave the way for an informed selection of a most preferred design. A generalization to 3D shapes seems possible and is the next natural step. Moreover, further optimization criteria like, for example, reliability w.r.t. other loading scenarios, minimal natural frequencies, and/or efficiency criteria, can be included into a general multiobjective shape optimization problem.

Acknowledgement

This work was supported by the federal ministry of research and education (BMBF, grant-no: 05M18PXA ) as a part of the GIVEN consortium.