Biomimetic Approach to Compliance Optimization and Multiple Load Cases

The variational approach to shape optimization in linearized elasticity is used in order to improve convergence of a known heuristic algorithm. The speed method of shape optimization is applied to obtain necessary optimality conditions for representative test examples. The algorithm originates from the biomimetic approach to compliance optimization. The trabecular bone adapts its form to mechanical loads and is able to form structures that are lightweight and very stiff at the same time. In this sense, it is a problem pertaining to both the nature or living entities which is similar to structural optimization, especially topology optimization. The paper presents the biomimetic approach, based on the trabecular bone remodeling phenomenon, with the aim of minimizing the compliance in multiple load cases. The method employed aims at minimizing the energy and combines structural evolution inspired by trabecular bone remodeling and the shape gradient framework, with strict analysis based on functionals in the 3-dimensional elasticity model. The method is enhanced to handle the problem of structural optimization under multiple loads. The new biomimetic approach does not require volume constraints. Instead of imposing volume constraints, shapes are parameterized by the assumed strain energy density on the structural surface. The stiffest design is obtained by adding or removing material on the structural surface in virtual space. Structural evolution is based on shape gradient approximation by the speed method, and it is separated from the finite element method of the model solution. Numerical examples confirm that the heuristic algorithm for structural optimization is efficient.


Introduction
The variational approach to shape and topology optimization can be developed, e.g., within the speed and the topological derivative methods. In particular, the topological derivative method is used to this end [1][2][3]. Our aim in this paper is to improve the convergence of the biomimetic approach to compliance optimization and multiple load cases by an application of the speed method. This means that we introduce a variational problem for which the proposed optimality criterion of biomimetic approach becomes a necessary optimality condition. Such an approach is known in the theory of free boundary problems.
The field of structural optimization is still a relatively new area but with many applications in the automotive, aerospace, civil engineering, machine design and other engineering fields. The classical approach to structural optimization can be found, e.g., in the monographs [4,5]. Many researchers must face challenges related to structural optimization. One of the crucial ones is to develop structural optimization methods that are appropriate for multiple load issues.
From an engineer's point of view, single load case problems are rather rare. A more common problem, but also more challenging in terms of mechanical design, is the one of structural optimization under multiple loads. The biomimetic optimization method should be useful as, in real life, multiple independent loads are always present.
It is widely believed that the trabecular bone adapts its form to mechanical loads following Wolff's law [6,7]. A healthy tissue of the trabecular bone has a very sophisticated structure. The tissue forms a network of beams called trabeculae, and this structure is capable of handling a wide range of loads. The length of a trabecula can be 100 or 200 µm, whereas its diameter is approximately 50 µm. This structure is continually rebuilding itself so that the whole bone tissue is replaced within a few years. The process is called trabecular bone adaptation or remodeling.
Wolff's law is not a mechanical law as such but rather the hypothesis that the trabecular bone is a self-optimizing material. In addition, other researchers suggested [8,9] that the trabecular tissue can adapt its structural form in reaction to external loads. In this sense, it is a problem which resembles structural optimization, especially topology optimization.
The strain energy density on the structural surface also appears in the area of optimization research, which is distant from biomechanical studies [10,11]. The most important theorem concerning the energy density on the structural surface states that for the stiffest design, the energy density must be constant along the designed shape. The new formulation of biomimetic optimization based on the trabecular bone phenomenon used in [12] was established with the use of the shape derivative concept. The mathematical analysis of the so-called speed method in shape optimization is performed, e.g., in [13].

Structural Design with Multiple Load Cases
A multiple load problem arises when a real-life design has to be developed. In topology optimization, a common procedure uses minimization of a weighted sum of different loads. The compliance for each load case is computed, and then, the multi-load problem is replaced by minimization of the weighted sum of the compliances under each load, as first presented by Sigmund [14]. For optimization of the worst possible case, another approach is needed. One of the possible methods is the bound formulation approach [15][16][17]. The main feature of the approach is the introduction of an upper bound for other objectives and minimizing the upper bound. Another approach uses the aggregation technique [18,19]. The main advantage of the constraint aggregation method is that all load cases contribute throughout the optimization process, owing to which the process does not lead to local minima.
Relying on the approach presented in [12], using biomechanical models of the trabecular bone remodeling directly as a base for structural optimization, we will discuss a design with multiple load cases. In general, there is no solution to the problem of stiffest design (compliance minimization). If the volume of an object is increasing, the compliance is decreasing. Hence, in the standard approach to energy-based topology optimization, an additional constraint has to be added. Usually, the volume of material is limited. When using volume constraints, additional stress constraints should be introduced. As a result, topology optimization has to be separated from shape and size optimization within the numerical procedure [17]. In the case of the biomimetic approach, the role of an additional constraint is played by the strain energy density on the structural surface, and the volume or structural mass results from the optimization procedure. In the standard SIMP topology optimization method, the volume constraint is used. Also, in the stress-based topology optimization approach [20,21], the volume constraint is present. To start the optimization procedure, the volume constraint has to be imposed. In the biomimetic approach, instead of imposing a volume constraint, we parameterize shapes by the assumed energy density, which may be quite accurately predicted from yield criteria.

Problem Setting
Let us define the compliance functional for the standard elasticity system div σ (u) = 0 in Ω, Here, ∂Ω = Γ 0 ∪ Γ 1 ∪ Γ v and t is a constant vector of traction forces. This functional depends on the shape of Ω, and Γ v is a part of its boundary subject to modification. Now, suppose we want to take into account two possible loads (t 1 , t 2 ) acting on Γ 1 , the rest of the formulation remaining intact. This gives rise to two compliances and two solutions u 1 , u 2 . We wish to minimize somehow both compliances by modifying Γ v under the volume constraint We have described this setting in the classical way in order to start similarly to other considerations concerning multiple load problems. In the numerical procedure we propose, shape and topology optimization is performed without volume constraint and under the assumed energy density.
The assumption of constant strain energy density for the stiffest design [12] was justified in the same setting, but for a single load. The derivative of the Lagrange function gives (at the stationary point) the following condition for the stiffest design configuration: Similar results were presented previously by other researchers [11]. Our result is based on the shape derivative concept which does not need any additional assumptions. The value of λ is not known, but the assumed value of the strain energy density on the part of the boundary subject to modification could be related to material properties. Changing that value results in a change of the structural form-topology and volume. In this way, the final structural volume results from the optimization procedure. Instead of imposing a volume constraint, we parameterize shapes by energy density, which may be quite accurately predicted from yield criteria.

Insensitivity Zone Concept
The biomimetic approach based on trabecular bone remodeling allows comprising optimization of size, shape and topology in one numerical procedure, where the lazy zone concept is an important component of the algorithm. The lazy zone concept was originally proposed by Carter [22], as an enhancement of the mechanostat theory of Frost [23]. The concept was also exploited by other researchers [9,24]. The stiffest design is obtained by adding or removing material on the structural surface in virtual space, mimicking the natural process of trabecular bone remodeling [25]. Bone mass increases above a certain level of mechanical stimulation (measured by the strain energy density on the tissue surface) and decreases below a certain energy level. When the strain energy density on the structural surface is between these two levels, the bone mass is maintained, and this rate of bone mechanical stimulation is called adaptation or lazy zone because then the bone does not react to changes in the mechanical stimulation level. The processes of resorption and formation are coupled in basic multicellular units, regions smaller than the single trabecula where a sequence of successive bone tissue resorption and formation occurs, and this is a local change. However, bone formation depends on global mechanical stimulation of the whole bone [8,9].

Heuristic Algorithm
The numerical results of [26] show a similarity between trabecular bone remodeling and structural optimization. The results have no biological background, but have been rigorously proven mechanically and mathematically. Nonetheless, trabecular bone remodeling models could be used to realize the optimization procedure. Our heuristic algorithm is as follows: -it is assumed that the energy density σ (u) : ε(u) has a constant value λ on Γ v ; -if at a given point on Γ v this density is greater than λ + s, then the boundary is moved outside; -if at a given point on Γ v this density is smaller than λ − s, then the boundary is moved inside; -these steps are repeated until equilibrium is achieved; -the value of λ is modified if the final design is unsatisfactory.
Since apart from the boundary modification there is a need to take into account the increase or decrease in the surface area itself, an additional step is introduced to the algorithm. The κ curvature is incorporated into the procedure in the following way. The improvement of descent direction should result from the partially heuristic modification.
-If κ > 0 and F >F, then after a biomimetic modification the boundary is additionally moved inside Ω by 50% of the biomimetic step. -If κ < 0 and F >F, then after a biomimetic modification the boundary is additionally moved outside Ω by 50% of the biomimetic step.
Thus, the curvature term locally enhances or diminishes the purely biomimetic modification of the structure. Numerical experiments seem to confirm this claim, and the speed method furnishes a better descent direction compared to the already published results.
The problems will be transformed suitably and analyzed using the speed method.

Definition of the Speed Method
Let V(x) be a smooth vector field defined on IR 3 and T t : IR 3 → IR 3 the transformation defined by Let also Ω 0 ⊂ IR 3 be a certain domain and Ω t = T t (Ω 0 ). Note that T 0 (Ω 0 ) = Ω 0 . Next, we take a function u t (x) defined on Ω t and depending on it explicitly or, for example, as a solution of the elasticity boundary value problem posed in Ω t . Using this notation, we define in Ω 0 the shape derivative of u t with respect to the (fixed) vector field V(x) by Here, we need in fact an extension of u t to Ω 0 . The shape derivative is used to analyze the behavior of domain functionals depending on both Ω t and u t . Their general form is and, denoting Γ t = ∂Ω t , If we define the shape derivatives of these functionals in the direction of the field V(x) as then it can be proved that and Here, κ(x) denotes the average curvature of the surface Γ 0 at the point x, namely where R 1 (x) and R 2 (x) are the principal curvatures at this point.

Problem Transformation
Problem 1 may be transformed into an equivalent form:

Problem 1a
The corresponding Lagrange function has the form Problem 2 does not need transforming, and the corresponding Lagrange function has the form Next, one may formulate the Kuhn-Karush-Tucker conditions for both problems. To this end, we shall need shape derivatives of compliances.
Let us rewrite the state equation in the weak form Here, u, ϕ ∈ H 1 Γ 0 (Ω). Then, we take the shape derivative of the compliance using the speed method (denoted by (•) ): and the shape derivative of the weak formulation of the state equation After substituting ϕ := u in the weak state equation and ϕ := u in its shape derivative, we get The Kuhn-Karush-Tucker conditions for Problem 1a are Hence In the case of symmetry, when interchanging t 1 and t 2 does not change the problem, and we have λ 1 = λ 2 = 0.5, the optimality condition reads Obviously, we do not know the value of λ 0 . The same conditions for Problem 2 are obtained immediately:

Shape Modification Using Shape Derivative Without Volume Constraint
Let us define a function F(z) penalizing the deviation of z from 0 and taking into account the insensitivity zone: Using this function, we replace the heuristic strategy by minimization of the functional where for Problem 1a, while for Problem 2, For simplicity, we shall consider Problem 1a. After taking the shape derivative of the functional (observe that F = dF/dz), we have where Here, κ is the average curvature at the given point on Γ v and 1 , u 2 )). (28) In order to get rid of u 1 and u 2 , we introduce adjoint equations for p 1 and p 2 in the weak form where p 1 , p 2 , ϕ ∈ H 1 Γ 0 (Ω). Then, substituting again ϕ := p 1 or ϕ := p 2 in the shape derivative of the weak state equations and ϕ := u 1 or ϕ := u 2 in the adjoint equations, we obtain As a result, When the expression in square brackets is negative, we should add material (V·n > 0); otherwise, material has to be removed (V · n < 0). Nevertheless, in such a form the algorithm may have a tendency to liquidate the variable part of the boundary Γ v altogether. To control this unwanted behavior, we introduce the averaged versions of J λ and F: Observe that for |Γ v | = Γ v ds we have Hence, the shape derivative of the averaged functional may be expressed as As a result, the final formula has the form We assume that the heuristic algorithm for Problem 1a is realized as follows, taking into account the first three terms in (35): -it is assumed that the energy density σ (u) : ε(u) has a constant value λ 0 on Γ v ; -the energy densities σ (u 1 ) : ε(u 1 ) and σ (u 2 ) : ε(u 2 ) are computed for loads t 1 and t 2 acting on the same boundary Γ 1 ; -the average energy density is 1 2 for loads t 1 and t 2 ; -if at a given point on Γ v this density is greater than λ 0 + s, the boundary is moved outside; -if at a given point on Γ v this density is smaller than λ 0 − s, the boundary is moved inside; -these steps are repeated until equilibrium is achieved; -the value of λ 0 is modified if the final design is unsatisfactory.
The last term in (35), involving the curvature κ, is incorporated into the procedure in the following way: -If κ > 0 and F >F, then after a biomimetic modification the boundary is additionally moved inside Ω by 50% of the biomimetic step. -If κ < 0 and F >F, then after a biomimetic modification the boundary is additionally moved outside Ω by 50% of the biomimetic step.
The insensitivity range parametrized by s is needed to ensure existence of equilibrium. For Problem 2 the algorithm will be the same, with 1

Numerical Examples
The heuristic biomimetic procedure is based directly on the biomechanical model of the trabecular bone remodeling phenomenon.
To compute the strain energy density on the structural surface, the finite element method is used. (The details of the numerical environment can be found in [12].) For The problem with square design domain with a central square rigid support; selected steps of the optimization procedure for the multiple load cases the multiple load case simulation, the problem with a square design domain with a central square rigid support was chosen. The reason was the existence of an exact analytical solution examined in many papers [27][28][29].
Two different load cases were examined, depicted in Fig. 1. The forces in each case are of the same magnitude. Material parameters were taken as follows: Young's modulus 2 × 10 11 Pa, Poisson's ratio 0.3 and the size of the insensitivity zone 80-180 MPa. Since the biomimetic method works in 3-D only, the design domain is in fact 3-dimensional, but due to very small thickness it can be treated as quasi 2-dimensional.
The results (selected steps) of the optimization procedure for the multiple load cases for force sets P and Q are presented in Fig. 2. The optimal configuration obtained is similar to the results presented previously by other researchers.
To compare our approach with other results derived by a more classical approach, a widely investigated T -bracket structure was chosen. The T -bracket problem model is presented in Fig. 3. Two load cases were analyzed. The forces for both load cases were of the same magnitude, and the material parameters were taken as in the previous example. For this example, 82 iterations were needed. The structural evolution was based on shape gradient approximation by the speed method, and it was separated from the finite element method computations. Each iteration included generation of a mesh (tetrahedral elements) and finite element calculations for each load case. The computational time for the initial configuration of 93,120 elements was 10 s (wall-clock time) for the mesh generation and 25 s for the finite element calculations for each load case. For the final configuration, there were 42,396 elements, 5 s for meshing and 20 s for structural calculation. The calculations were carried out with the use of our own mesh generation tool Cosmoprojector and the parallel finite element solver based on Message Passing Interface (MPI). The computations were run on 4 CPUs with the use of the cluster (16 GB RAM per node). The details of the mesh generation method can be found in [12]. Fig. 3 The T -bracket problem model. Left: Force located on the left part of the domain; right: force located on the right part of the domain. All loads are of the same magnitude Fig. 4 The T -bracket structure problem; selected steps of the optimization procedure for the multiple load cases  The results for the T -bracket structure problem analyzed in [30] with the level set method are very similar to the results presented in Fig. 4.
In order to test the procedure for 3-D problems, two different load cases were examined, and the design domains are presented in Fig. 5. It was assumed that all nodes of the back wall of the domain box could be fixed (clamped wall), and the bending force was applied to the middle of the opposite, front wall. The first load case was with a vertical bending force, and the second with the same definition of boundary conditions and a horizontal bending force. All loads were of the same magnitude. Material parameters and the size of the insensitivity zone were assumed to be as in the previous example.
Selected steps of the optimization procedure are presented in Fig. 6.

Conclusions
In this paper, a new formulation of biomimetic optimization based on trabecular bone remodeling applied to a multiple load problem is presented. The stiffest design is obtained by adding or removing material on the structural surface in virtual space. Three elements complete the formula driving the process of adding or removing material: non-local influence of boundary modification, spatial variability of the function penalizing deviation from the assumed strain energy density on the structural surface and measuring the increase or decrease in the area of the surface itself, by measuring the surface curvature. The structural evolution is based on shape gradient approximation by the speed method, and it is separated from finite element method computations. The finite element method is used to compute the distribution of the strain energy density only. An important advantage of our approach is that optimization results are in the form of a functional configuration. Accordingly, the assumed value of the strain energy density on the part of the boundary subject to modification can be related to material properties. Instead of imposing a volume constraint, we parameterize shapes by energy density, which may be quite accurately predicted from yield criteria. The functional configurations during the process of optimization allow including size, shape and topology optimization in one numerical procedure.
The 2-D result for the common benchmark is in line with previous results by other researchers [27][28][29]. The 3-D result shows the configuration for a real material. In view of the above, our 3-D results with two perpendicular bending forces as multiple loads could be a good benchmark for testing other approaches to multiple load optimization.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided 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.