A duality principle and related computational method for a class of structural optimization problems in elasticity

TL;DR

This paper introduces a duality principle and a computational method for optimizing the topology of elastic structures to minimize internal energy, using convex analysis and duality theory, with numerical examples demonstrating its effectiveness.

Contribution

It develops a novel duality-based computational approach for structural topology optimization in elasticity, avoiding filtering steps to ensure critical point solutions.

Findings

The method successfully finds critical points for the optimization problem.

Numerical examples validate the theoretical duality approach.

The approach is applicable to various elastic structural problems.

Abstract

In this article we develop a duality principle and concerning computational method for a structural optimization problem in elasticity. We consider the problem of finding the optimal topology for an elastic solid which minimizes its structural inner energy resulting from the action of external loads to be specified. The main results are obtained through standard tools of convex analysis and duality theory. We emphasize our algorithm do not include a filter to process the results, so that the result obtained is indeed a critical point for the original optimization problem. Finally, we present some numerical examples concerning applications of the theoretical results established.

Equations

$$\hat{J}(u,t)=\frac{1}{2}⟨u_i,f_i{⟩}_{L^2(\Omega )}+\frac{1}{2}⟨u_i,{\hat{f}}_i{⟩}_{L^2({\Gamma }_t)},$$

$$\left\{\begin{array}[]{lll}(H_{ijkl}(t)e_{kl}(u))_{,j}+f_{i}=0\text{ in }\Omega,\\ \\ H_{ijkl}(t)e_{kl}(u)n_{j}-\hat{f}_{i}=0,\text{ on }\Gamma_{t},\;\forall i\in\{1,2,3\}.\end{array}\right.$$

$$U=\{u=(u_1,u_2,u_3)\in W^{1,2}(\Omega ;{\mathbb{R}}^3):u=(0,0,0)=0\text{ on }{\Gamma }_0\},$$

$$B=\left\{t:\Omega \to [0,1]\text{ measurable }:{\int }_{\Omega }t(x)dx=t_1|\Omega |\right\},$$

$$0<t_1<1$$

$${\sigma }_{ij}(u)=H_{ijkl}(t)e_{kl}(u),$$

$$e_{ij}(u)=\frac{1}{2}(u_{i,j}+u_{j,i}),\forall i,j\in \{1,2,3\}.$$

$$\{H_{ijkl}(t)\}=\{t{H}_{ijkl}^0+(1-t)H_{ijkl}^1\},$$

$$J_1(u,t)=-G(e(u),t)+⟨u_i,f_i{⟩}_{L^2(\Omega )}+⟨u_i,{\hat{f}}_i{⟩}_{L^2({\Gamma }_t)},$$

$$G(e(u),t)=\frac{1}{2}{\int }_{\Omega }{H}_{ijkl}(t)e_{ij}(u)e_{kl}(u)dx,$$

$$dx=d{x}_{1}d{x}_{2}d{x}_3.$$

$$J^{\ast }(u)$$

$$H_{ijkl}^0{z}_{ij}{z}_{kl}>c_0{z}_{ij}{z}_{ij}$$

$$H_{ijkl}^1{z}_{ij}{z}_{kl}>c_1{z}_{ij}{z}_{ij},\forall z=\{z_{ij}\}\in {\mathbb{R}}^{3\times 3},\text{ such that }z\ne 0.$$

$$J(u,t)=\hat{J}(u,t)+Ind(u,t),$$

$$Ind(u,t)=\left\{\begin{array}[]{ll}0,&\text{ if }(u,t)\in A^{*},\\ +\infty,&\text{ otherwise },\end{array}\right.$$

$$A_1=\{(u,t)\in U\times B:({\sigma }_{ij}(u){)}_{,j}+f_i=0,\text{ in }\Omega ,\forall i\in \{1,2,3\}\}$$

$$A_2=\{(u,t)\in U\times B:{\sigma }_{ij}(u)n_j-{\hat{f}}_i=0,\text{ on }{\Gamma }_t,\forall i\in \{1,2,3\}\}.$$

$$J(u_0,t_0)$$

$$G^{\ast }(\sigma ,t)$$

$$\{{\overline{H}}_{ijkl}(t)\}=\{H_{ijkl}(t){\}}^{-1}$$

$$C_1=\{\sigma \in Y^{\ast }:{\sigma }_{ij,j}+f_i=0,\text{ in }\Omega ,\forall i\in \{1,2,3\}\}$$

$$C_2=\{\sigma \in Y^{\ast }:{\sigma }_{ij}{n}_j-{\hat{f}}_i=0,\text{ on }{\Gamma }_t,\forall i\in \{1,2,3\}\}.$$

$$\underset{(u,t)\in U\times B}{\inf }J(u,t)$$

$$\frac{\partial J_1(u_0,t_0)}{\partial u}=0$$

$$(H_{ijkl}(t_0)e_{kl}(u_0){)}_{,j}+f_i=0\text{ in }\Omega ,$$

$$H_{ijkl}(t_0)e_{kl}(u_0)n_j-{\hat{f}}_i=0\text{ on }{\Gamma }_t,\forall i\in \{1,2,3\},$$

$$G(e(u_0))=\frac{1}{2}⟨(u_0{)}_i,f_i{⟩}_{L^2(\Omega )}+\frac{1}{2}⟨(u_0{)}_i,{\hat{f}}_i{⟩}_{L^2({\Gamma }_t)}.$$

$$J^{\ast }(u_0)$$

Full text

$\underline{\hskip 449.55362pt}$

1 Introduction

Consider an elastic solid which the volume corresponds to an open, bounded, connected set, denoted by $\Omega\subset\mathbb{R}^{3}$ with a regular (Lipschitzian) boundary denoted by $\partial\Omega=\Gamma_{0}\cup\Gamma_{t}$ where $\Gamma_{0}\cap\Gamma_{t}=\emptyset.$ Consider also the problem of minimizing the functional $\hat{J}:U\times B\rightarrow\mathbb{R}$ where

[TABLE]

subject to

[TABLE]

Here $\mathbf{n}=(n_{1},n_{2},n_{3})$ denotes the outward normal to $\partial\Omega$ and

[TABLE]

[TABLE]

where

[TABLE]

and $|\Omega|$ denotes the Lebesgue measure of $\Omega.$

Moreover $u=(u_{1},u_{2},u_{3})\in W^{1,2}(\Omega;\mathbb{R}^{3})$ is the field of displacements relating the cartesian system $(0,x_{1},x_{2},x_{3})$, resulting from the action of the external loads $f\in L^{2}(\Omega;\mathbb{R}^{3})$ and $\hat{f}\in L^{2}(\Gamma_{t};\mathbb{R}^{3}).$

We also define the stress tensor $\{\sigma_{ij}\}\in Y^{*}=Y=L^{2}(\Omega;\mathbb{R}^{3\times 3}),$ by

[TABLE]

and the strain tensor $e:U\rightarrow L^{2}(\Omega;\mathbb{R}^{3\times 3})$ by

[TABLE]

Finally,

[TABLE]

where $H^{0}$ corresponds to a strong material and $H^{1}$ to a very soft material, intending to simulate voids along the solid structure.

The variable $t$ is the design one, which the optimal distribution values along the structure are intended to minimize its inner work with a volume restriction indicated through the set $B$.

The duality principle obtained is developed inspired by the works in [10, 11]. Similar theoretical results have been developed in [3], however we believe the proof here presented, which is based on the min-max theorem is easier to follow (indeed we thank an anonymous referee for his suggestion about applying the min-max theorem to complete the proof). A theory for a topology optimization problem in elasticity is presented in [4], even though in our book [3] of 2014, we have developed a more general result with a proof based on the inverse function theorem. Also, dual methods for discrete structural optimization problems were used in [5] .

We highlight throughout this text we have used the standard Einstein sum convention of repeated indices. Related models, among others, are addressed in [9].

A Matlab code using a filter for the numerical computation of similar problems is presented in [7]. We emphasize in our algorithm we have not used a filter. The majority of topology optimization works use filtering to avoid the check-board effect [1, 12]. One of the proposals of this work is to increase discretization in the direction of the loads to avoid this problem.

Moreover, details on the Sobolev spaces addressed may be found in [8]. In addition, the primal variational development of the topology optimization problem has been described in [2].

One of the main contributions of this work is to present detailed theoretical developments for such a class of structural optimization problems, through duality theory and an application of the min-max theorem. We have also discovered that without the use of any filters, to avoid the up-surging of the check-board problem in some parts of the optimal structure, it is necessary to discretize more in the load direction, in which the displacements are much larger.

Finally, it is worth mentioning the numerical examples presented have been developed in a Finite Element (FE) context, based on the work of [7].

2 Mathematical formulation of the topology optimization problem

Our mathematical topology optimization problem is summarized by the following theorem.

Theorem 2.1**.**

Consider the statements and assumptions indicated in the last section, in particular those refereing to $\Omega$ and the functional $\hat{J}:U\times B\rightarrow\mathbb{R}.$

Define $J_{1}:U\times B\rightarrow\mathbb{R}$ by

[TABLE]

where

[TABLE]

and where

[TABLE]

Define also $J^{*}:U\rightarrow\mathbb{R}$ by

[TABLE]

Assume there exists $c_{0},c_{1}>0$ such that

[TABLE]

and

[TABLE]

Finally, define $J:U\times B\rightarrow\mathbb{R}\cup\{+\infty\}$ by

[TABLE]

where

[TABLE]

where $A^{*}=A_{1}\cap A_{2},$

[TABLE]

and

[TABLE]

Under such hypotheses, there exists $(u_{0},t_{0})\in U\times B$ such that

[TABLE]

where

[TABLE]

[TABLE]

and $C^{*}=C_{1}\cap C_{2},$ where

[TABLE]

and

[TABLE]

Proof.

Observe that

[TABLE]

Also, from this and the min-max theorem, there exist $(u_{0},t_{0})\in U\times B$ such that

[TABLE]

Finally, from the extremal necessary condition

[TABLE]

we obtain

[TABLE]

and

[TABLE]

so that

[TABLE]

Hence $(u_{0},t_{0})\in A^{*}$ so that $Ind(u_{0},t_{0})=0$ and $\sigma(u_{0})\in C^{*}.$

Moreover

[TABLE]

This completes the proof. ∎

3 About the computational method

The continuous topology optimization problem described in the previous section is discretized using the FE method, considering in plane deformations. The FE discretization is performed taking into account the bilinear isoparametric element as a master one, in similar way as in [6, 7].

To obtain computational results, we have defined the following algorithm.

Set $n=1$. 2. 2.

Set $t_{1}(x)=t_{1},\text{ in }\Omega.$ 3. 3.

Calculate $u_{n}\in U$ as the solution of equation

[TABLE]

that is

[TABLE] 4. 4.

Obtain $t_{n+1}$ by

[TABLE] 5. 5.

Set $n:=n+1$ and go to step 3 up to the satisfaction of an appropriate convergence criterion.

In the FE formulation, equations indicated in 9 stands for

[TABLE]

where $\mathbf{H}(t)$ is the global stiffness matrix, $\mathbf{U}$ is the global displacements vector and $\mathbf{f}$ is the global forces one.

Thus, for such a FE models ($N$ elements where $e\in\{1,...,N\}$), the primal optimization problem can be written in a matrix form as

[TABLE]

On the other hand, the dual problem may be expressed by

[TABLE]

and where $V_{e}$ is the area of element $e$.

Finally, the last minimization indicated corresponds to item 4 in the concerning algorithm. Indeed, such a procedure refers to minimize at each sub-iteration, through the Matlab Linprog routine (that is, in a sequentially linearized context), the function

[TABLE]

subject to $t\in B,$ where $p$ is a penalization parameter (typically, $p=3$).

4 Computational simulations and results

We present numerical results in an analogous two-dimensional context, more specifically for two-dimensional beams of dimensions $1\times l$ (units refer to the international system) represented by $\Omega=[0,1]\times[0,l]$, with $l=0.5$, $F=-10^{6}$ for the first case, $l=0.5$ and $F=-10^{7}$ for the second one, $l=0.6$ and $F=-10^{6}$ for the third case and, $l=1$ and $F=-10^{8}$ for the fourth one. $F$ is in the $y$-direction and corresponds to $f$ of the theoretical formulation presented above.

We consider the strain tensor as

[TABLE]

where $\mathbf{u}=(u,v)\in W^{1,2}(\Omega;\mathbb{R}^{2}),$ $e_{x}(\mathbf{u})=u_{x}$, $e_{y}(\mathbf{u})=v_{y}$ and $e_{xy}(\mathbf{u})=\frac{1}{2}(u_{y}+v_{x}).$

Moreover the stress tensor $\sigma(e(\mathbf{u}))$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

where $E_{0}=210*10^{9}$ (the modulus of Young) and $E_{1}\ll E_{0}.$ Moreover $\nu=0.33$.

As previously mentioned, we present four numerical simulations.

Case 1. For the first case see figure 1, on the left, for the concerning case, figure 1, in the middle, for the optimal topology for this case with no filter, figure 1, on the right, for the optimal topology for this first case with filter. For the objective function as function of number of iterations also for such a case with no filter, see figure 2, on the left, and for the objective function as function of number of iterations also for this first case with filter, see figure 2, on the right.

Case 2. For the second case see figure 3, on the left, for the concerning case, figure 3, in the middle, for the optimal topology for this case with no filter, figure 3, on the right, for the optimal topology for this second case with filter. For the objective function as function of number of iterations also for such a case with no filter, see figure 4, on the left, and for the objective function as function of number of iterations also for this second case with filter, see figure 4, on the right.

Case 3. For the third case see figure 5, on the left, for the concerning case, figure 5, in the middle, for the optimal topology for this case with no filter, figure 5, on the right, for the optimal topology for this third case with filter. For the objective function as function of number of iterations also for such a case with no filter, see figure 6, on the left, and for the objective function as function of number of iterations also for this third case with filter, see figure 6, on the right.

Case 4. For the fourth case see figure 7, on the left, for the concerning case, figure 7, in the middle, for the optimal topology for this case with no filter, figure 7, on the right, for the optimal topology for this fourth case with filter. For the objective function as function of number of iterations also for such a case with no filter, see figure 8, on the left, and the objective function as function of number of iterations also for this fourth case with filter, see figure 8, on the right.

We emphasize to have obtained in both optimized structures, without filter and with filter, robust topology from a structural point of view. One can note also in the figures that in all cases the objective functions, without filter and with filter, have similar final value, which indicates that the results obtained are consistent.

5 Final remarks and conclusions

In this article we have developed a duality principle and relating computational method for a class of structural optimization problems in elasticity. It is worth mentioning we have not used a filter to post-process the results, having obtained a solution $t:\Omega\rightarrow\{0,1\}$ (that is, $t(x,y)=0$ or $t(x,y)=1$ in $\Omega$), by finding a critical point $(u_{0},t_{0})\in U\times B$ for the functional $J_{1}:U\times B\rightarrow\mathbb{R}.$ This corresponds, in some sense, to solving the dual problem.

We address some final remarks and conclusions on the results obtained.

• •

For all examples, in a first step, we have obtained numerical results through our algorithm with a software which uses the Matlab-Linprog as optimizer at each iteration without any filter. In a second step, we obtain numerical results using the OC optimizer with filter, with a software developed based in the article [7] by Sigmund, 2001.

• •

We emphasize, to obtain good and consistent results, it is necessary to discretize more in the direction $y$, that is, the load direction, in which the displacements are much larger.

• •

If we do not discretize enough in the load direction, for the software with no filter, a check-board standard in the material distribution is obtained in some parts of the concerning struture.

• •

Summarizing, with no filter, the check-board problem is solved by increasing the discretization in the load direction.

• •

Moreover, with the OC optimizer with filter, the volume fraction of material is kept constant in 0.5 at each iteration during the optimization process, whereas for the case with no filter we start with a volume fraction of 0.95 which is gradually decreased to the value 0.5, using as the initial solution for a iteration with a specific volume fraction, the solution of the previous one.

• •

We also highlight the result obtained with no filter is indeed a critical point for the original optimization problem, whereas there is some heuristic in the procedure with filter.

• •

Once more we emphasize to have obtained more robust and consistent shapes by properly discretizing the approximate model in a FE context.

• •

Finally, it is also worth mentioning, we have obtained similar final objective function values without and with filter in all examples, even though without filter such values have been something smaller, as expected. The qualitative differences between the graphs without and with filter, for the objective function as function of the number of iterations, refer to the differences between the optimization processes, where in the case with filter the volume fraction is kept 0.5 and without filter it is gradually decreased from 0.95 to 0.5, as above described.

We highlight the results obtained may be applied to other problems, such other models of plates, shells and elasticity.