An evolutionary topology optimization approach with variationally controlled growth

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801 www.elsevier.com/locate/cma

An evolutionary topology optimization approach with variationally controlled growth Dustin R. Jantos a,∗ , Philipp Junker b , Klaus Hackl a a Institute of Mechanics of Materials, Ruhr-Universit¨at Bochum, Germany b Institute of Mechanics, Bergische Universit¨at Wuppertal, Germany

Received 16 March 2016; received in revised form 12 July 2016; accepted 18 July 2016 Available online 3 August 2016

Abstract Previous works of Junker and Hackl (2016) have presented a variational growth approach to topology optimization in which the problem of checkerboarding was suppressed by means of a discontinuous regularization scheme. This approach did not require additional filter techniques and also optimization algorithms were not needed any more. However, growth approaches to topology optimization demand some limitations in order to avoid a global and simultaneous generation of mass. The limitation has been achieved by a rather simple approach with restricted possibilities for controlling. In this contribution, we eliminate this drawback by introducing a Lagrange multiplier to control the total mass within the model space for each iteration step. This enables us to achieve directly controlled growth behavior and even find optimized structures for prescribed structure volumes. Furthermore, a modified growth approach, which we refer to as the Lagrange shift approach, results a numerically stable model that is easy to handle. After the derivation of the approach, we present numerical solutions for different boundary problems that demonstrate the potential of our model. c 2016 Elsevier B.V. All rights reserved. ⃝

Keywords: Topology optimization; Growth; Variational modeling; Constraint evolution; Discontinuous Galerkin approach

1. Introduction The objective of topology optimization is to find the topology of a mechanical structure that possess maximum stiffness for a given design space with the prescribed supports and loads. A popular method for topology optimization is the SIMP approach (Solid Isotropic Material with Penalization), in which the spatial distribution of structure mass is the design variable. This can be achieved e.g. by defining a density field within a model space. To prevent the algorithm from filling the whole design space with mass, which can be referred to as the trivial solution, constraints for the total mass or volume of the structure are usually introduced. These approaches use a finite element environment to determine the displacement field in combination with external optimization algorithms to find the spatial distribution of structure mass. Examples for SIMP approaches are given in [1–3]. A different method for topology optimization are growth approaches that find optimal solutions for all volume structures by increasing very same step by step. The ∗ Corresponding author.

E-mail address: [email protected] (D.R. Jantos). http://dx.doi.org/10.1016/j.cma.2016.07.022 c 2016 Elsevier B.V. All rights reserved. 0045-7825/⃝

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

781

simulation is stopped as soon as the prescribed volume or stiffness of the structure is exceeded [4–6]. Combinations of SIMP and growth approaches also exist [7–9]. Within the BESO method (Bidirectional Evolutionary Structural Optimization) [10,11] or methods based on genetic algorithms [12,13], structural parts are systematically added and removed within multiple iterations to find optimal structures. Optimization can be achieved by modeling structural fronts and boundaries within level-set methods [14–16] or phasefield methods [17,18]. Reviews and comparisons of different approaches are given in [19–21]. In our previous work, we applied variational principles known from material modeling to the problem of topology optimization. In material modeling, internal variables, denoted by υ, are introduced that describe the current microstructural state (e.g. plastic strains, volume fractions, damage parameters). The time-dependent behavior of the internal variables is modeling in terms of so-called evolution equations that have a structure given by υ˙ = f (ε, υ, . . .). There are various fundamental modeling schemes that enable a thermodynamically consistent formulation of the evolution equation, among which the Hamilton principle and the associated variational techniques are very well-known. For the problem of topology optimization, we also applied the Hamilton principle, which resulted in a variational growth approach. The density field was indirectly described by a compliance parameter, which served as internal variable. Similarly to material modeling, the time evolution of the compliance parameter was determined by evaluating the evolution equation within a solitary finite element environment, i.e. without any additional optimization algorithm. To prevent the trivial solution of a solid design space, mass generation was limited by defining an adaptive dissipation parameter that served directly as a threshold value for the driving forces, which determine locally the direction and rate of the unconstrained mass generation or removal. The final “decision” as to whether the compliance parameter is changing is given by a yield function that is similar e.g. to plasticity [22,23]. In this paper, we present a new approach by deriving an evolution equation in which the structure volume of the system over (pseudo-)time is directly controlled by a given growth function instead of utilizing a yield function. A Lagrange multiplier is introduced for this purpose. As we will show later, the Lagrange multiplier can be interpreted as a threshold value (i.e. yield limit) similar to our previous model. The new approach enables us to combine the advantages of a growth approach and approaches in which the volume structure is prescribed: solutions for different volume structures are found within one simulation and by keeping the volume structure constant for several simulation steps, the solution converges to a stiffer solution for these (prescribed) volume structures. By using a special growth function, which we refer to as the Lagrange shift approach, an easy handling of the global growth (velocity) is provided. The new approach also includes the discontinuous Galerkin approach for regularization from [23] to suppress intra-element checkerboarding. Discontinuous shape functions for the internal variable describing the density field enable evaluation of the evolution equations by consideration of its gradient locally within each finite element separately. Thus, the density field is not an additional field unknown within a Newton–Raphson iteration scheme (i.e. nodal unknown for the finite element treatment), and the numerical effort is thus greatly reduced. Starting with derivations of the variational equations for the new model in Section 2, the numerical treatment of very same is given in Section 3. In Section 3.1, we define a finite element approximation including discontinuous shape functions for the internal variable describing the density field. This makes it possible to find an analytical solution for the Lagrange multiplier in Section 3.2 and to derive the evolution equation for the compliance parameter. The numerical treatment closes with the presentation of adaptive numerical parameters in Section 3.3 and growth functions in Section 3.4. Numerical results are presented in Section 4, and a conclusion together with an outlook are given in Section 5. 2. Variational controlled growth model 2.1. Preliminaries The following relations in this section have already been presented in [22,23] and are repeated here for convenience. More details are given in the original publications. In contrast to SIMP, in which a relative mass density p is usually used as design variable, we define the material distribution in the present approach by an internal variable χ (x) ∈ [0, 1] that can be interpreted as a compliance parameter. A special interpolation for the internal variable will be introduced by ρ(χ ) = 1/ f (χ ), where ρ = ρ(x) ∈ [κ, 1] is the mass density given by p n in SIMP approaches with some numerical parameter n (see [24]). The governing equation for the key variable χ in our approach will be derived employing Hamilton’s principle. A direct formulation in terms of the density ρ, however, would yield a model describing damage processes (thus

782

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

reducing stiffness for maximum load). More details on this interesting aspect can be found in [22,23]. The quantity f (χ) is some mapping penalizing values of ρ ∈]κ, 1[ (see Eq. (6)). The mapping f (χ ) in our approach fulfills the same purpose as the power law p n in SIMP approaches. The parameter 0 < κ ≪ 1 is a small, but non-zero numerical factor to prevent singular tangent operator matrices given by vanishing elastic constants. According to Bendsøe [24], the problem of topology optimization for a linear-elastic material reads    1 σ : C(χ ) : σ dV min min E∈Ead σ ∈S 2 Ω with S = {σ | ∇ · σ + b = 0 in Ω , σ · n = t on ∂σ Ω }

(1)

where σ are the stresses, b are the body forces, and t are the given external tractions that act on the boundary ∂σ Ω . The design space is indicated by Ω and the material compliance is C. The set of admissible stiffness tensors Ead can be defined in various ways, see [24]. The stresses are defined as σ = E : ε.

(2)

The material quantities E and C = E−1 for the stiffness and compliance, respectively, depend on the spatial material distribution given by the internal variable χ (x), i.e. C(χ ) = f (χ )C0 ⇔ E(χ ) =

1 E0 = ρ(χ )E0 f (χ )

(3)

where C0 is the minimum material compliance and E0 is the maximum material stiffness. The stiffness E(χ ) is equal to the set of admissible stiffness tensors Ead . The structure volume V can be determined by evaluating   1 ρ(χ ) dV = V = dV. (4) f (χ ) Ω Ω To penalize ρ ̸∈ {κ, 1}, we choose a density interpolation function f (χ ) as a nonlinear function that penalizes values of χ ̸∈ {0, 1}. According to C(χ ) = χ 2 Cvoid + (1 − χ 2 )C0 C0 , the interpolation function f (χ ) follows as   1 − 1 χ2 f (χ ) = 1 + κ

with Cvoid =

(5)

1 κ

(6)

which ensures f (χ = 0) → full material (ρ = 1) → E(χ = 0) = E0 f (χ = 1) → “void” (ρ = κ) → E(χ = 1) = κ E0 .

(7)

Compared with power laws expressed as E = p n E0 where p ∈ {0, 1} and n > 1 in [1–3], the density change ′ = − ff 2(χ) , using the approach given in Eq. (6), is zero for both desired final configurations χ ∈ {0, 1} (see Figs. 1 (χ) and 2). A simple power law would cause a consecutive increase in structure volume for any full material configuration if no further restrictions are applied, which thus, increases the complexity for controlling the total mass in the system and granting a stable numerical convergence. Therefore, a simple power law is less convenient for a variational growth approach. By using a power law, the cardinality of the penalization of intermediate densities p ̸∈ {0, 1} can be controlled by the exponent n. In Eq. (6) this penalization can be increased by reducing the numerical parameter κ. The displacement field u and the spatial distribution of χ = χ (x) are the unknown quantities of the problem. The objective of the topology optimization is to determine these quantities so that the work of the external forces (often referenced as the compliance energy) will be minimized. Although there exists no physical time in topology optimization, the iteration steps can be interpreted as discretized (pseudo-)time steps so that topology optimization can also be performed using growth approaches. Evolution equations expressed as χ˙ = function(u, χ , x, . . .) are thus one possibility to close the system of equations for the displacement field and the spatial distribution of material. ∂ρ ∂χ

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

783

Fig. 1. Comparison of a simple power law (dashed) and the approach using Eq. (6) for density interpolation with different values for the exponent n and the numerical parameter κ.

Fig. 2. Slope of the density interpolation for the two different approaches.

2.2. Variational approach The Hamilton principle for dissipative continua reads  ∂D δG + δχ dV + cons = 0 Ω ∂ χ˙

(8)

where G is the Gibbs energy, D is the dissipation function, and cons are problem-specific constraints. For evaluation of Eq. (8), specifications are required for the Gibbs energy and the dissipation function. A deduced variational method for deriving material models is the so-called principle of the minimum of the dissipation potential [25]. Let us begin with the Gibbs energy, which we use in its standard form but expand it by an additional constraint    G= Ψ dV − b · u dV − t · u dA + λg(χ ). (9) Ω

Ω

∂Ω

The Lagrange multiplier λ ensures with  1 1 ! g(χ ) = dV − ϱ(t) = 0 VΩ Ω f (χ )

(10)

784

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

the adjustment of the desired structure volume V (t) = VΩ ϱ(t), where ϱ(t) ∈ [κ, 1] denotes the volume part of the evolving structure as a given growth function and VΩ is the model space volume. We define the Helmholtz free energy Ψ = Ψ (ε, χ ) as Ψ := Ψm + Ψr

(11)

where Ψm denotes the mechanical contribution of Ψ and is chosen similar to [22,23] as Ψm :=

1 1 1 ε: E0 : ε = σ : f (χ )C0 : σ 2 f (χ ) 2

(12)

and the energy for regularization similar to [23] as α |∇χ|2 (13) 2 where α > 0 is a numerical factor. The choice of a particular dissipation function D = D(χ˙ ) mainly influences the structure of the resulting evolution equation for the internal variable χ . We choose a rate-dependent approach analogous to a viscous material (in contrast to our previous model in [23] in which we chose an elasto-viscoplastic type). We define Ψr :=

1 r χ˙ 2 (14) 2 which results in evolution equations of the desired type (rate-dependent). The parameter r controls the evolution velocity and is thus referred to as viscosity of our approach. More details are given in [22]. As we see later, the viscosity r will be very important to ensure controlled growth of the structure volume in the system during numerical solution which also influences the numerical stability of the model. Since r is a numerical parameter in our approach and is not directly related to a “real” viscosity, we will present an adaptive approach for r as presented in Section 3.3. Finally, the constraints cons have to be specified, in this case as  γ δχ dV (15) cons := D=

Ω

where the Kuhn–Tucker parameter γ constrains the interval of the internal variable χ ∈ [0, 1] with  : χ˙ > 0 ∧ χ = 1  γ¯ γ = −γ¯ : χ˙ < 0 ∧ χ = 0  0 : else.

(16)

The constraint χ ∈ [0, 1], in turn, provided by Eq. (16) yields the admissible material stiffness E(x) ∈ [κE0 , E0 ] ∀ x ∈ Ω. The Gibbs energy depends on the displacement field u(x) and the spatial distribution of the internal variable χ (x). Thus, the variation has to be performed for both u and χ, which yields   ∂D ! δu G + δχ G + δχ dV + γ δχ dV = 0. (17) ∂ χ ˙ Ω Ω The variations δu and δχ can be evaluated independently which leads to   ∂D δu G = 0 ∧ δχ G + δχ dV + γ δχ dV = 0. Ω ∂ χ˙ Ω

(18)

The variation of G with respect to the displacement with g ̸= g(u) (for small deformations) yields the balance of linear momentum in its weak form    ∂Ψm δu G = : δε dV − b · δu dV − t · δu dA = 0 ∀ δu (19) Ω ∂ε Ω ∂Ω in which the energy-stress relation ∂Ψm /∂ε = σ = ε : E(χ ) can be inserted.

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

The variation of G with respect to the internal variable χ yields     ∂Ψr λ f ′ (χ ) ∂Ψm δχ dV + δχ G = − · ∇δχ dV = 0 2 ∂χ V ∂∇χ f (χ ) Ω Ω Ω

∀

Inserting this result into the Hamilton principle in Eq. (8) finally gives     ∂Ψr ∂Ψm λ f ′ (χ ) ∂D − + + γ δχ dV + · ∇δχ dV = 0 2 ∂χ VΩ f (χ ) ∂ χ˙ Ω Ω ∂∇χ

δχ .

∀

785

(20)

δχ .

(21)

Considering the presented specifications for the mechanical and regularization parts of Ψ , Eqs. (12) and (13), and the approach for D according to Eq. (14) allows us to derive the “balance equation” for χ as     λ f ′ (χ ) − pχ + − + r χ˙ + γ δχ dV + α ∇χ · ∇δχ dV = 0 ∀ δχ (22) VΩ f (χ )2 Ω Ω with the mechanical driving force, see (22, 23) for more details, pχ := −

∂Ψm 1 f ′ (χ ) ε : E0 : ε. =− ∂χ 2 f (χ )2

(23)

The driving forces accomplish a comparable task as the sensitivity analysis: both define the direction to a “more optimal” structure. In the presented approach, the Hamilton principle according to Eq. (8) can be interpreted as the stationarity condition for an appropriate functional. This functional, in turn, constitutes as the objective function in our approach. Consequently, the driving forces, which are necessary for the calculation of the stationarity condition in Eq. (8), fulfill the same purpose as variational (continuous) derivatives in “classical” approaches: both determine the direction of steepest decent [26,27]. The Lagrange multiplier λ in Eq. (22) is still unknown, and a direct formulation cannot be obtained due to the gradient of the internal variable χ. We will introduce some numerical formulations before we can find an analytical formulation for the Lagrange multiplier in Section 3.2. 3. Numerical treatment 3.1. Finite element approach For the numerical solution, we use different Galerkin approximations for the displacement field and the internal variable similar to [23]. For the displacement field, we use common trilinear shape functions that are collected in the vector Nu and enable the discretizations u(x) = Nu (x) · u, ˆ

δu(x) = Nu (x) · δ uˆ

(24)

where (ˆ·) denotes the nodal values of the corresponding property in vectorial form. For the internal variable, we use a discontinuous Galerkin approach χ (x) = Nχ (x) · χ˜ ,

χ˙ (x) = Nχ (x) · χ˙˜ ,

δχ (x) = Nχ (x) · δ χ˜

(25)

with  √     3 3 1 1 1 Nχ,i (ξ1 , ξ2 , ξ3 ) ∈ √ ± ξ1 √ ± ξ2 √ ± ξ3 8 3 3 3

(26)

where ξ1 , ξ2 , ξ3 define the natural coordinates of the isoparametric finite element. Consequently, the tilde (˜·) denotes the values at the Gauß points instead of the nodal quantities indicated by (ˆ·). For more details, we refer to [23]. The checkerboarding phenomenon renders the problem of topology optimization to be ill-posed when discretized with linear shape functions. In contrast to SIMP approaches, we do not assign one compliance parameter to an entire finite element. The compliance parameter may change from Gauß point to Gauß point. This “finer” discretization causes the checkerboarding phenomenon to become present at a lower level: in contrast to inter-element checkerboarding for SIMP approaches, in which entire elements oscillate between mass and no mass, checkerboarding takes

786

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

place in an intra-element way when discretized at the Gauß point level. Due to this intra-element checkerboarding, a discontinuous Galerkin approach is appropriate for a simplified numerical evaluation of the balance equation (22) (for more details see [23]). In accordance with ε = 21 (∇u + u∇), we define the B-operator matrix as εV = B · u, ˆ

δεV = B · δ uˆ

(27)

where the index V indicates properties in Voigt notation. Eq. (19) in its finite element approximation is then    1 ! R := BT · · E0,V · ε V dV − NuT · b dV − NuT · t dA = 0 f (χ ) Ω Ω ∂Ω

(28)

which can be solved by a Newton–Raphson scheme by evaluating ∂R ! · ∆uˆ = 0 (29) ∂ u˜ with (i + 1) and (i) as the later and former iteration steps, respectively. The tangent operator ∂R/∂ uˆ is given by    1 ∂R f ′ (χ ) ∂χ BT · = E0,V − 2 ⊗ E0,V · ε V · B dV. (30) ∂ u˜ f (χ ) f (χ ) ∂εV Ω R(i+1) = R(i) +

The missing derivative ∂χ /∂εV will be given in Section 3.2. It can be determined by evaluating the evolution equation (35). Inserting the discretizations according to Eq. (25) into δχ G in Eq. (21) yields     λ f ′ (χ ) ! ˙ − pχ − α(∇Nχ · χ˜ ) · (∇Nχ · δ χ˜ ) dV = 0 ∀ δ χ˜ (31) + r Nχ · χ˜ + γ Nχ · δ χ˜ dV + 2 VΩ f (χ ) Ω Ω which can be evaluated for arbitrary δ χ˜ as     λ f ′ (χ ) ! ˙ − pχ − α(∇Nχ )T · (∇Nχ · χ˜ ) dV = 0. + r N · χ ˜ + γ N dV + χ χ VΩ f (χ )2 Ω Ω By introducing the abbreviation    T   =: ∆χ GP α ∇Nχ,h · ∇Nχ,h · χ˜ dV Ωe

(32)

(33)

GP

we can reduce Eq. (32) to    λ f ′ (χ ) ∆χ GP  − pχ − + r χ ˙ + + γ  2  VΩ f (χ ) Ve

= 0.

(34)

x=xGP,e

The index h in Eq. (33) indicates that the gradient has to be evaluated only elementwise — there is no further information on the gradient of other elements needed: by aid of Eq. (34), the “balance” equation for the internal variable χ can now be evaluated locally at each Gauß point xGP,e of each element e (except for ∆χ GP and the element volume Ve , which needs to be evaluated beforehand by an integral over each finite element separately). Therefore, the later evolution equation can be also evaluated locally without solving the internal variable as an additional field of unknowns (analogously to the displacement field). The numerical effort for the solution can thus be greatly reduced. The iterative evaluation of the evolution equation (34) is similar to sequential mathematical programming solutions in “classical” topology optimization approaches[21,24,28]. For further information, please see [23]. 3.2. Determination of the Lagrange multiplier To determine a formulation for the Lagrange multiplier, we solve Eq. (34) for the rate of χ which is   ∆χ GP 1 λ f ′ (χ ) χ˙ = + pχ − −γ . r VΩ f (χ )2 Ve

(35)

787

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

Inserting Eq. (35) into the rate form of the condition for structure volume (see Eq. (10))  1 f ′ (χ ) ! χ˙ dV − ϱ(t) ˙ =0 g(χ ˙ )=− 2 V f (χ ) Ω Ω

(36)

leads to   λ = VΩ

Ω − pχ +

∆χ GP Ve

 ′ (χ) + γ ff (χ) ˙ 2 dV − r VΩ ϱ(t) .   2  f ′ (χ) dV 2 Ω f (χ)

(37)

∆χ

(χ) GP The Lagrange multiplier with VλΩ ff (χ) − γ , which allow 2 globally shifts the zero level of the driving forces pχ − Ve the sign of the rate χ˙ of the internal variable to change locally resulting in local mass generation or removal. The density will thus decrease in areas in which the driving forces are relatively low and will increase in areas with high driving forces. The local decrease and increase in density balance each other globally so that the rate of the relative volume of the total structure equals the given function ϱ(t). ˙ For all internal variables χ = 1 (“void”), the Lagrange multiplier can be interpreted as a threshold for growth: negative rates χ˙ of the internal variable cannot decrease the local density to negative values due to the Kuhn–Tucker condition, so that the density only increases in areas in which the driving forces exceed the zero level. ∂χ The missing derivative ∂ε in Eq. (30) can be calculated using Eqs. (35), (37), and (23): by applying an explicit V time integration, the evolution equation (35) can be discretized as    ∂ pχ ∂χ ∂χ (i+1) ∂ χ˙  ∆t 1 ∂λ f ′ (χ ) = = ∆t = + ∂εV ∂εV ∂εV χ =χ (i) r VΩ ∂εV f (χ )2 ∂εV χ =χ (i)   ′ (i) ∆t 1 ∂λ f (χ ) = − E0,V · εV r VΩ ∂εV f (χ (i) )2    ′ (i) 2 f (χ ) ′ (i) dV E · ε V 0,V ∆t  Ω f (χ (i) )2  f (χ ) = − E · ε   V 0,V   2  r f (χ (i) )2 f ′ (χ (i) ) dV Ω f (χ (i) )2     ′ (i) 2 f (χ ) ′ (i) E · B · u dV 0,V (i) 2 ∆t  Ω f (χ )  f (χ ) − E · B · u = . (38)   0,V   2  r f (χ (i) )2 f ′ (χ (i) ) dV (i) 2 Ω f (χ ) ′

Due to the explicit formulation in Eq. (38), ∆χ GP only depends on the known values of χ at the previous iteration step. Consequently, the only quantities that enter the derivative are ∂ pχ /∂εV and ∂λ/∂εV . 3.3. Adaptive viscosity and gradient penalization As discussed above, the dissipation parameter r fulfills the purpose of the viscosity in our rate-dependent approach and thus controls the velocity of the local growth. The maximal local growth rate in nature is a given value. However, we are not interested in modeling natural growth with accurate behavior in time, but try to find the optimal structure as fast as possible (meaning fast calculation time). Thus, we calculate the value for r so that the simulation time is reduced while simultaneously keeping the simulation numerically stable. Alternatively, we could have defined an adaptive time increment ∆t instead of an adaptive r . The Kuhn–Tucker condition, which provides the obedience of χ ∈ [0, 1], can be realized by a simple if-condition:  for χ (i+1) < 0 : χ (i+1) → 0 (39) (i+1) for χ >1: χ (i+1) → 1. Therefore, any interval transgression will be cut off, which leads to a deviation in the structure volume and can cause numerically unstable behavior. Since the evolution equation (35) contains r −1 as a prime factor, too small values

788

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

for the viscosity r increase the interval transgression and can cause numerically unstable behavior. Additionally the viscosity r appears in the Lagrange multiplier in Eq. (37) as a factor to the prescribed structure volume. However, by inserting (37) into (35), the viscosity r disappears as a factor for the prescribed structure volume. Thus, a given (global) growth rate ϱ(t) ˙ does not take into account the (local) growth velocity of the internal variable represented by the viscosity r . This leads to the following problem: if the values for the viscosity r and/or the growth rate ϱ(t) ˙ are relatively large compared to the driving forces, the Lagrange multiplier can become negative:    ∆χ GP f ′ (χ ) dV ⇒ λ < 0. (40) +γ r VΩ ϱ(t) ˙ > − pχ + Ve f (χ )2 Ω A negative Lagrange multiplier means a negative zero level for the driving forces. This causes growth in the whole model space and leads to the trivial solution since there is no interfering restriction for the driving forces. Even if the Lagrange multiplier becomes positive, due to the strict obedience of the global growth rate, high global growth rates ϱ(t) ˙ in combination with low local growth rates (large values for r ) can cause coarse structures and/or intermediate density configurations. The reason for this is that if each point in the model space has a low growth rate, large areas of the model space must be active to achieve the required global growth rate. Summing these facts up, small values for the viscosity r can cause unstable numerical behavior, whereas too large values cause intermediate configurations or can even lead to a trivial solution. Since the driving forces can change drastically within the simulation time and thus the optimal value for the viscosity, we introduce an adaptive viscosity r . (i) For each internal variable χ, a viscosity rχ will be calculated for which the specific χ would reach the interval border χ ∈ {0, 1} at the iteration step (i). Let us introduce some abbreviations (i−1)

∆χ GP F := − −γ Ve    ′ (i−1) (i−1) (i−1) f χ λn := −F 2 dV  Ω f χ (i−1)   ′  (i−1)  2 f χ (i−1) λd := dV  2 Ω f χ (i−1)   (i) λn − rχ VΩ ϱ˙ t (i) . λ = VΩ λd (i−1)

pχ(i−1)

Then, the evolution equation (35) can be written as    (i)   ′  (i−1)  2 (i−1) (i) f χ λ − r V ϱ ˙ t 1 χ n Ω χ˙ (i) = (i)  + F (i−1)  .  2 (i−1) rχ λd f χ (i−1)

(41) (42)

(43)

(44)

(45)

We define the individual rate of internal variable so that in each case the interval border for the internal variable is reached  1 − χ (i−1)   , χ˙ (i−1) > 0  (i) ∆t (i) ˙ χ˙ = χ˜ = (46) (i−1)    −χ , χ˙ (i−1) ≤ 0 ∆t where the tilde (˜·) was introduced for Gauß point values (see Section 3.1). This leads to rχ(i)

(i−1)

=

F (i−1) λd

(i) (i−1) χ˜˙ λd

 2 (i−1) ′  (i−1)  f χ (i−1) + λn f χ  2    . f χ (i−1) + VΩ ϱ˙ t (i) f ′ χ (i−1)

(47)

Eq. (47) needs to be evaluated for each (discretized) internal variable χ. If the global viscosity becomes the maximum (i) (i) of all rχ , the internal variable(s) corresponding to the maximal rχ will reach the interval border χ ∈ {0, 1}, whereas

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

789

all other internal variables are too “slow” to reach any interval border. Due to the explicit assumption in terms of the strains for the adaptive viscosity, interval transgression is still possible and numerical instabilities can still occur by (i) simply defining the maximal rχ as the global viscosity. Therefore, we introduce the safety factor rS > 1 that is a plain multiplier to the adaptive viscosity to slow down growth and prevent interval transgression. The global viscosity r (i) for the current iteration step is then   r (i) = rS max rχ(i) . (48) For the regularization, we choose the parameter α for the gradient penalization based on our previous work [23]:    Re · ue  ∀ Ω (49) α = Ve max e Ve  with the residuals Re , displacements ue , and volume Ve of the corresponding element containing Ωe . 3.4. Growth function Since we focus on a growth approach, we chose χ = 1 ∀ χ ∈ Ω (“void”) as the initial condition for the internal variable. By defining the rate of the growth function ϱ(t) ˙ in Eq. (37), it is possible to prescribe any arbitrary growth behavior for the structure with respect to time. Specifying the growth behavior ϱ(t) over time by a given function for the relative structure volume can easily be performed by defining   ϱ t (i+1)  V − V (i) Ω (i+1) ϱ˙ t = (50) ∆t  with the structure volume from the previous step V (i) = Ω f χ1(i) dV . Due to the independence of ϱ ̸= ϱ(χ ) in ( ) Eq. (10), the actual function ϱ(t) ˙ does not depend on the internal variable except for the known value in the previous step in Eq. (50)). With Eq. (50) it is also possible to define a constant target structure volume by defining  (as shown   ϱ t (i+1) as a constant value, more precisely than by defining ϱ˙ t (i+1) = 0. Eq. (50) can be used to compare the current and the desired (constant) structure volume in each iteration step. The rate ϱ(t ˙ (i+1) ) is then calculated so that the volume divergence due to rounding errors is reduced. It is possible to keep the structure volume quite accurately (∆ϱ < 0.1%), whereas the model converges to a “more optimal” solution with same mass. In this way, it is also possible to hold the structure volume constant during the whole simulation time by changing the initial condition to homogeneous distribution of the target structure mass within the model space. For the numerical evaluation of the model, we have considered two different growth functions. On the one hand, we introduce an exponential growth approach with an upper bound for the growth rate, and on the other hand, we define a growth function that depends on the non-restricted driving forces in each step. We refer to the later as the Lagrange shift approach. For the exponential growth we define   V (i) ϱ˙ t (i+1) := η VΩ

(51)

with η > 0 as a given parameter describing the relative structure volume increase  within  each  step.  In addition, we define an upper bound for the growth rate by prescribing a maximum increment ϱ t (i+1) − ϱ t (i) ≤ ∆ϱmax . This is recommended since a purely exponential growth would lead to a high growth rate in the later simulation steps, which would reduce the formation of fine structures and lead to an early trivial solution. There are two possible “extreme cases” for the growth that should be prevented: no growth (ϱ(t) ˙ = 0) and simultaneous growth in the entire model space leading to the trivial solution if the growth rate ϱ(t) ˙ is too high (see Section 3.3). Every growth rate ϱ(t) ˙ between these two cases leads to localized, i.e. non-simultaneous, growth granting non-trivial solutions. While the first case ϱ(t) ˙ = 0 can be easily prevented, the second case cannot be prevented directly, because the upper bound for the growth rate ϱ(t) ˙ = 0 is unknown a priori. However, to prevent the second case, we introduce the Lagrange shift approach, which enables us to “shift” continuously between these two cases by simply defining a single numerical parameter λS ∈ [0, 1] within the Lagrange multiplier instead of directly controlling

790

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

the growth rate ϱ(t) ˙ = 0. For the Lagrange shift approach, we define    1 1 − λS dV  ϱ(t) := VΩ Ω f (χ (t)) λ=0

(52)

where 0 < λS < 1 is a numerical parameter. The growth function in Eq. (52) states that the volume of the structure within each time step becomes a fraction (i.e. 1 − λS ) of the volume structure which would appear if there were no restrictions (i.e. λ = 0) to the growth and the mass would evolve freely according to the local driving force. Therefore, this particular growth function adapts to the growth defined by the actual driving forces, making this approach numerically quite stable since the driving forces are mainly influencing the structural evolution. Although the growth function in Eq. (52) depends on the internal variable, it can be easily shown that this particular function does not change δχ g of Eq. (10): ∂g(χ (λ)) ∂g ∂λ δχ = δχ ∂χ ∂λ ∂χ         ∂ 1 ∂ 1 − λS 1 1 ∂λ = dV δχ + dV  δχ ∂χ VΩ Ω f (χ) ∂λ VΩ Ω f (χ) λ=0 ∂χ   

δχ g =

=0

 1 f ′ (χ) = dV δχ (53) VΩ Ω f (χ)2 which is identical to δχ g in Eq. (20). Thus, all previous equations hold true. To find the Lagrange multiplier for the Lagrange shift approach, we first have to evaluate the time derivative of the growth function as    (1 − λS ) f ′ (χ ) ϱ(t) ˙ := . (54) − χ˙ dV  2 VΩ f (χ ) Ω λ=0 Inserting the evolution equation (37) yields       (1 − λS ) ∆χ GP f ′ (χ ) 1 λ f ′ (χ )  ϱ(t) ˙ := − + pχ − −γ dV  2 2  VΩ r V det J f (χ ) f (χ ) Ω Ω λ=0    (1 − λS ) f ′ (χ ) ∆χ GP − pχ + = − γ dV. 2 r VΩ det J Ω f (χ ) By inserting (55) into (37), we finally find  ′    ∆χ GP f (χ) S) + γ dV − r VΩ (1−λ − p + χ Ω Ω Ve r VΩ f (χ)2 λ = VΩ   f ′ (χ) 2 dV Ω f (χ)2  ′   ∆χ GP f (χ) Ω − pχ + Ve + γ f (χ)2 dV = λS VΩ   f ′ (χ) 2 dV Ω f (χ)2

f ′ (χ) f (χ)2

(55)

  ∆χ − pχ + VeGP − γ dV

(56)

as the Lagrange multiplier for the Lagrange shift approach which is equal to the solution for the Lagrange multiplier with ϱ˙ = 0 multiplied by the factor λS . If λS = 1, there is no growth and the model preserves the structure volume. Choosing λS = 0 would lead to a model without any restrictions for the structure growth which would obviously result in the trivial solution since the density would increase simultaneously in the entire model space. Therefore, the parameter λS shifts the zero level for the driving forces pχ from a state in which the structure volume does not change (ϱ(t) ˙ = 0, λS = 1) towards the solution in which no restrictions for the growth are given (λ = λS = 0). The benefit of the Lagrange shift approach is that only one bounded parameter λS ∈ ]0, 1[ must be chosen for the growth rate. In this case, any admissible 0 < λS < 1 will grant a positive Lagrange multiplier which solves the problem of the trivial solution stated above and in Section 3.3. The Lagrange shift approach is numerically more stable since the growth function now depends on the driving forces which can drastically change within the simulation

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

791

Fig. 3. Boundary conditions for the L-shaped cantilever.

progress causing unstable behavior. Furthermore, the Lagrange multiplier λ is independent of the viscosity r so that the choice of r does not influence the fineness of the structure. Thus, the viscosity r can be chosen so that the simulation is most stable without any disadvantages for the structural evolution. If the solution for a given structure volume is desired, the growth function can be switched according to Eq. (50) as soon as the structure volume exceeds the given target volume. 4. Numerical results For all simulations, we choose κ = 5.0 × 10−3 and the adaptive viscosity according to Eq. (47) with rS = 1.5 to ensure numerical stability. By using ∆t = 1.0, we unify the time step t with the number of iteration steps. As linear elastic material parameters, we choose E = 210 MPa for Young’s modulus and ν = 0.3 for the Poisson’s ratio. For the first and second example, we choose the Lagrange shift approach with λS = 0.1, whereas we will compare different values of λS and an exponential growth in our third example. We run the simulations on an Intel Core i7 processor with 2.8 GHz without parallelization (single core evaluation) on a 64-bit system. Let us introduce the quasi-2D1 problem of an L-shaped cantilever given by [24]. The boundary conditions are given in Fig. 3. In a first step, we simulate the free growth of the structure given by the Lagrange shift approach. Fig. 4 shows the evolution of the structure volume over simulation time t and Figs. 5 and 62 the structure and its structure compliance C = uˆ · fˆ for the marked points 1, 2, 3a, 4a, and 5a in Fig. 4. Then, we pick some intermediate results that we want to improve by holding the structure volume constant for additional 150 iteration steps. We restart the simulation at the picked points 3a, 4a, and 5a by simply using the resepctive results for the density field as initial conditions and apply a growth function according to Eq. (50) with ϱ(t (i+1) ) = const. (Fig. 6 3b, 4b, and 5b). Due to the fact that the model only contains the growth rate ϱ˙ as direct parameter, rounding errors would cause a deviation from the target volume by only defining ϱ˙ = 0. In contrast, if Eq. (50) is used, the growth rate will be adapted to the current structure volume to reduce this deviation. The average calculation time for the mesh with 4900 nodes on 2304 elements (60 elements alongside the longest edge) was 0.75 s per iteration step. 1 Our model is implemented in 3D. To simulate two-dimensional problems, we apply meshes with thickness of exactly one finite element. Thus, the number of degree of freedoms is higher than for a true two-dimensional problem (twice as much nodes and three instead of two unknowns per node) accompanied by higher calculation effort. 2 No filter techniques are used for visualization. The values for the density field ρ = 1/ f (χ ) given in the Gauß-points are extrapolated to the nodes according to the shape function in Eq. (25). Thus, each nodal value of the density field is the average density of all neighboring Gauß-point values.

792

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

60 50 40 30 20 10 200

400

600

800

1000

1200

1400

Fig. 4. Evolution of the relative structure volume ϱ(t) in % over time. 1, 2, 3a, 4a, 5a: for free growth given by the Lagrange shift approach; 3b, 4b, 5b: growth function according to Eq. (50) with ϱ(t (i+1) ) = const.

Fig. 5. Evolution of the structure in early simulation steps. Structure compliance: C = uˆ · fˆ .

In each case, the structure compliance C is decreasing while the structure volume is held constant. Although the changes of the structure shape is minimal, intermediate (gray) regions are reduced and the truss-like structures are more elaborated. Nevertheless, the relative improvement of the structure compliance is significant: 27% for ϱ = 28%, 15.6% for ϱ = 36% and 7% ϱ = 47%. Fig. 7 shows the direct comparison of the results from a SIMP approach given in [24] with our final result (Fig. 6 5b) and the result obtained with a finer mesh (19 012 nodes on 9216 elements) to show the mesh independency of our model. The topology (meaning the number of holes) of our results are identical (except for a little hole at the tip in the finer mesh) to the result from the SIMP approach given in [24], but the shapes are different. The SIMP approach yields a more truss like structure, whereas our model seems to provide slightly more organic shapes. The most obvious difference is the bulge at the inner corner and the thinner vertical trusses connected to the support. Thus, our model seems to prioritize regions with high stresses, like the singularity at the inner corner. In our next simulation, we will prescribe only one target structure volume and have a detailed look at the change of structure compliance while holding the structure volume constant as well as comparing the calculation time with our previous models in [22,23]. We evaluate the quasi-2D simply supported beam shown in Fig. 8 with the target volume ϱtarget = 45.65% as done in [22,23]. An overview of the calculation time and a comparison with our previous models for a mesh with 96 × 16 × 1 = 1536 elements and 3298 nodes is given in Table 1. The present approach needs about 14% more iteration steps to reach the target volume than [23] but the calculation time per iteration step is reduced by about 20%. By using an optimized mesh (better ordering of the unknowns for

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

793

Fig. 6. Evolution of the structure in later simulation steps. 3a, 4a, 5a: for free growth given by the Lagrange shift approach; 3b, 4b, 5b: structural changes in each case after holding the structure volume constant for 150 iteration steps. Structure compliance: C = uˆ · fˆ .

794

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

Fig. 7. Comparison of the results for the L-shaped cantilever (ϱ = 47%): SIMP approach given in [24] (left), Lagrange shift with same mesh (middle) and Lagrange shift with a twice as fine mesh (120 elements alongside the longest edge) (right).

Fig. 8. Boundary conditions for a simply supported beam. Table 1 Overview of calculation times for the simply supported beam with the present approach and the approaches from [22,23]. Iteration steps

Total CPU time [min]

CPU time per it. step [s]

Relative CPU time

Junker, Hackl 2015 [22] Junker, Hackl 2016 [23]

1723 626

164 27

5.7 2.6

100% 16%

Present approach

713 1000

24 33

2.0 2.0

14% 20%

Present approach (optimized mesh)

713 1000

4.6 6.4

0.39 0.39

2.8% 3.9%

the solver), we could reduce the time for each iteration step by about 80%. The main calculation effort lies within the solution of the (linear elastic) mechanical problem for each new density field. To determine the evolution of the density field, multiple integrals need to be computed but no additional non-linear equations or equation systems must be solved. Thus, each improvement to the linear solver has a huge impact on the overall calculation time. In contrast, more iteration steps are needed as e.g. compared to SIMP approaches. However, it is worth to mention that our approach operates only using a finite element code — there is no need for any (external) optimization algorithm. Fig. 9 shows the evolution of the structure volume ϱ(t) and the structure compliance C, including a closer look at the steps t ≥ 713 in which the volume structure is held constant. The results using an isovolume filter3 for the 3 Employed software: ParaView. The isovolume filter displays every point of a given scalar as solid if its value is within a given interval. To display the density field in ParaView, we extrapolated the density values ρ = 1/ f (χ ) within the Gauß points to the nodal points. Thus, the density value within a node between two neighboring Gauß point densities ρ = {0, 1} will become 0.5, which matches the chosen isovolume filter.

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

795

Fig. 9. Structure volume ϱ(t) and the structure compliance C = uˆ · fˆ for the quasi-2D problem with the prescribed structure volume ϱtarget = 45.65%.

evolution of the structure over time are given in Fig. 10 for the growth process and in Fig. 11 for the structure at the time step t = 713, at which the target volume is exceeded for the first time, and at the end of the simulation t = 1000. The target structure volume ϱ = 45.65% is reached in step t = 713. The maximal absolute deviation of the structure volume was ∆ϱ ≈ 0.04%. While holding the structure volume constant, the structure is still changing and its structure compliance is reduced from 6.692 × 10−4 to 6.191 × 10−4 . This is a relative improvement of about 7%. The structure seems also “more optimal”: the thickness of single beam-like connection parts is more homogeneous and the lower arch seems more circular. We will now compare the Lagrange shift approach with an exponential growth function and also investigate different values of the parameter λS for the Lagrange shift approach. We re-introduce the bending problem from [22,23] shown in Fig. 12. We compare an exponential growth according to Eq. (51) with η = 0.01 and ∆ϱmax = 0.001 and the Lagrange shift approach applying the values λS = {0.1, 0.5, 0.9}. For both growth functions, we do not define a prescribed structure volume so that all simulations will continuously add mass to the system and will (slowly) tend to the trivial solution ϱ = 1.0. An overview over the calculation times and comparison with our previous models for a mesh with 32 × 20 × 6 = 3840 elements and 4851 nodes and ϱ = 33.6% is given in Table 2. For this example, the calculation time per iteration step of the present approach is reduced by about 6% in comparison to [23] without mesh optimization. By applying an optimized mesh, the calculation time can be reduced by about 40%. Of course, the number of required iteration steps depends strongly on the chosen growth function (and the desired target structure volume). To reach the structure volume ϱ = 33.6%, the exponential growth and Lagrange shift with λS = 0.1 require about 9% less iteration steps compared to [23]. In turn, this results in an overall reduced calculation time of about 10%, whereas Lagrange shift with λS = 0.5 requires 50% more steps and with λS = 0.9 even 87% more steps.

796

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

Fig. 10. Structural changes during the growth process ϱ < 45.65% (isovolume filter: ϱ > 0.5).

Fig. 11. Structural change for constant ϱ = 45.65% after the growth process (isovolume filter: ϱ > 0.5). Table 2 Overview of calculation times for the bending problem (ϱ = 33.6%) with the present approach and the approaches from [22,23]. Growth type

Iteration steps

Total CPU time [min]

CPU time per it. step [s]

Relative CPU time

Junker, Hackl 2015 [22] Junker, Hackl 2016 [23]

– –

2566 616

595 56

13.9 5.5

100% 9.4%

Present approach

exp. growth λS = 0.1 λS = 0.5 λS = 0.9

559 574 937 1151

48 50 82 100

5.2 5.2 5.2 5.2

8.1% 8.4% 13.8% 16.8%

Present approach (optimized mesh)

exp. growth λS = 0.1 λS = 0.5 λS = 0.9

559 574 937 1151

28 29 47 58

3.0 3.0 3.0 3.0

4.7% 4.9% 7.9% 9.7%

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

797

Fig. 12. Boundary conditions for the bending problem.

The evolution of the structure volume ϱ(t) and the structure compliances C with respect to the individual structure volume are displayed in Fig. 13. The corresponding structures for an isovolume filter with ϱ(t) ≥ 0.5 without any further filtering or smoothing are given in Fig. 14. The structure volume ϱ(t) for the exponential growth follows an exponential function until the maximum volume increment is reached at t ≈ 350. The following constant volume increase is intended, whereas flattening of the curve for t & 800 results from the very small driving forces. For structures with a high volume, additional material only improves the structure insignificantly, resulting in minor driving forces and density change. This causes the actual growth to be slower than that prescribed and finally leads to unstable behavior at t & 1200, although this is not critical since the trivial solution is nearly reached. The Lagrange shift approach also shows a comparable evolution, although flattening of the curve is significantly higher. The stronger flattening is reasonable since the Lagrange shift approach is based on a growth proportional to the actual driving forces. For greater values of λS , the growth is slower due to the fact that less volume is allowed to be added per iteration step. The structure compliance decreases with increasing λS , because the model has more time to optimize structural details. Overall, all growth functions tend to similar structures and structure compliances. With greater λS , the structure compliance improves for early iteration steps and therefore for smaller structure volumes, but this greatly increases the number of required iteration steps. The results seem reasonable: the truss-like structure for ϱ < 50% is similar to the results from [22], whereas the structure for ϱ ≈ 50% tends to the results for the variable thickness problem given by [29]. The calculation effort is comparable to our previous model in [23]. Thus, the computation time is strongly reduced compared to [22] (about 90%) and much coarser meshes are needed compared to [29] (3840 vs. 245,760 elements). We refer to [23] for further details. Although the exponential growth converges at least as quickly as the fastest Lagrange shift approach (λS = 0.1) and the results are at least as good as the best Lagrange shift approach (λS = 0.9), the effort to find appropriate parameters η and ∆ϱmax is quite high compared to the Lagrange shift approach. For the Lagrange shift, approach only one, clearly bounded parameter needs to be specified 0 < λS < 1. Furthermore, for all its values, the simulations are stable and do not contain any trivial solution as is possible for an arbitrarily defined growth function (see Section 3.4). For either value of λS , the resulting structures can become relatively fine and complex, particularly for three dimensional problems with spacious model spaces. Fig. 15 shows the results for the boundary conditions given by the bending problem in Fig. 12 but with a doubled model space width for λS = 0.1. In this case, more substructures occur, and the overall structure becomes more complex with multiple branches. Nevertheless, the structure is symmetric so that the result can be considered to be non-random and numerically stable.

798

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

Fig. 13. Structure volume ϱ(t) and the structure compliance C = uˆ · fˆ for different growth functions applied to the bending problem.

5. Conclusions and outlook In this paper, we derived a variational growth approach for topology optimization that includes direct control of global mass generation. The governing equations for the displacement field and the internal variable (compliance parameter) were derived using the Hamilton principle without any further assumptions. By introducing a Lagrange multiplier, we were able to control the structure volume within each iteration step by a given growth function. The Lagrange shift approach for a growth function led to easy-to-handle and numerically stable growth. We also introduced a finite element implementation of the model containing a discontinuous Galerkin approach for regularization from [23]. With the discontinuous Galerkin approach, we were able to provide solutions that were free of (intra-element) checkerboarding by penalizing the gradient of the internal variable describing the density on the element level. With this approach, it was also possible to find an analytical formulation for the Lagrange multiplier and a localized evolution equation so that neither the Lagrange multiplier nor the density field were additional field unknowns in the finite element scheme. We also introduced adaptive numerical parameters to stabilize the numerical behavior and increase user friendliness. The numerical results showed overall fine, smooth and, for even relatively simple boundary problems complex structures. The resulting structures seem reasonable, and in some cases they are comparable to natural and bone-like structures. Relatively coarse meshes were used, and in combination with the efficient regularization, the calculation times were only several seconds per iteration step using a single CPU thread. The number of required iterations steps to reach a reasonable structure volume depends on the growth function, although faster growth leads to “less optimal” solutions. Therefore, we introduced a method to find “more optimal” solutions for prescribed volume structures. In

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

799

Fig. 14. Snapshots for similar structure volumes from different growth functions during the simulation in which the structure volume increases continuously over time.

the given example, it was possible to make the resulting structure about 7% stiffer whereas the divergence between the actual and the prescribed structure volume were less than 0.1%. In future studies, the prescribed volume structure could be held constant over several iteration steps for multiple values. Thus, “more optimal” structures for different prescribed structure volumes could be determined within one simulation. Because we have full control of the volume structure over time, the model could be used to find an optimized structure that fulfills a failure criterion with minimum usage of material: structure growths as usual until the failure criterion is fulfilled the first time. Subsequently, the volume structure is held constant or reduced if the failure criterion is over-fulfilled by a given tolerance. These are minor additions to the presented model; however, in the long term different material models could be added to our model. For example, hyperelasticity or anisotropic materials with material orientation as an additional design variable. FE2 or equivalent multiscale approaches are also

800

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

Fig. 15. Result for the bending problem with doubled width after 626 iteration steps. (Structure volume: ϱ = 24.84%, structure compliance C = 0.473, average calculation time per step is 13.4 s).

imaginable. Since these approaches require high calculation effort, the numerical efficiency of the presented approach is promising due to localized regularization and embedding into a single finite element environment. In this work, the model was implemented in a single thread system; however, future work will consider parallelization techniques to shorten the calculation times further. References [1] M.P. Bendsøe, Optimal shape design as a material distribution problem, Struct. Optim. 1 (4) (1989) 193–202. [2] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (9–10) (1999) 635–654. [3] G. Rozvany, M. Zhou, The coc algorithm, part i: cross-section optimization or sizing, Comput. Methods Appl. Mech. Engrg. 89 (1) (1991) 281–308. [4] T.P. Harrigan, J.J. Hamilton, Bone remodeling and structural optimization, J. Biomech. 27 (3) (1994) 323–328. [5] E. Kuhl, A. Menzel, P. Steinmann, Computational modeling of growth, Comput. Mech. 32 (1–2) (2003) 71–88. [6] T. Waffenschmidt, A. Menzel, Application of an anisotropic growth and remodelling formulation to computational structural design, Mech. Res. Commun. 42 (2012) 77–86. [7] A. Klarbring, B. Torstenfelt, Dynamical systems and topology optimization, Struct. Multidiscip. Optim. 42 (2) (2010) 179–192. [8] A. Klarbring, B. Torstenfelt, Dynamical systems, simp, bone remodeling and time dependent loads, Struct. Multidiscip. Optim. 45 (3) (2012) 359–366. [9] A. Klarbring, B. Torstenfelt, Lazy zone bone remodeling theory and its relation to topology optimization, Ann. Solid Struct. Mech. 4 (1–2) (2012) 25–32. [10] J.H. Zhu, W.H. Zhang, K.P. Qiu, Bi-directional evolutionary topology optimization using element replaceable method, Comput. Mech. 40 (1) (2006) 97–109. http://dx.doi.org/10.1007/s00466-006-0087-0. [11] X. Huang, Y.M. Xie, Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials, Comput. Mech. 43 (3) (2008) 393–401. http://dx.doi.org/10.1007/s00466-008-0312-0. [12] S.D. Rajan, Sizing, shape, and topology design optimization of trusses using genetic algorithm, J. Struct. Eng. 121 (10) (1995) 1480–1487. http://dx.doi.org/10.1061/(ASCE)0733-9445(1995)121:10(1480). [13] P. Hajela, E. Lee, C.-Y. Lin, Genetic algorithms in structural topology optimization, in: Topology Design of Structures, Springer Netherlands, Dordrecht, 1993, pp. 117–133. http://dx.doi.org/10.1007/978-94-011-1804-0-10. [14] G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (1) (2004) 363–393. http://dx.doi.org/10.1016/j.jcp.2003.09.032. URL http://www.sciencedirect.com/science/article/pii/S002199910300487X. [15] Q. Xia, M.Y. Wang, Topology optimization of thermoelastic structures using level set method, Comput. Mech. 42 (6) (2008) 837–857. http://dx.doi.org/10.1007/s00466-008-0287-x. [16] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Engrg. 192 (1–2) (2003) 227–246. http://dx.doi.org/10.1016/S0045-7825(02)00559-5. URL http://www.sciencedirect.com/science/article/pii/S0045782502005595. [17] B. Bourdin, A. Chambolle, Design-dependent loads in topology optimization, ESAIM Control Optim. Calc. Var. 9 (2003) 19–48. http://dx.doi.org/10.1051/cocv:2002070. [18] L. Blank, H. Garcke, L. Sarbu, T. Srisupattarawanit, V. Styles, A. Voigt, Phase-field approaches to structural topology optimization, in: Constrained Optimization and Optimal Control for Partial Differential Equations, Springer Basel, Basel, 2012, pp. 245–256. http://dx. doi.org/10.1007/978-3-0348-0133-1-13. [19] D.J. Munk, G.A. Vio, G.P. Steven, Topology and shape optimization methods using evolutionary algorithms: a review, Struct. Multidiscip. Optim. 52 (3) (2015) 613–631. http://dx.doi.org/10.1007/s00158-015-1261-9. [20] J.D. Deaton, R.V. Grandhi, A survey of structural and multidisciplinary continuum topology optimization: post 2000, Struct. Multidiscip. Optim. 49 (1) (2014) 1–38. [21] G. Rozvany, A critical review of established methods of structural topology optimization, Struct. Multidiscip. Optim. 37 (3) (2009) 217–237.

D.R. Jantos et al. / Comput. Methods Appl. Mech. Engrg. 310 (2016) 780–801

801

[22] P. Junker, K. Hackl, A variational growth approach to topology optimization, Struct. Multidiscip. Optim. (2015) 1–12. [23] P. Junker, K. Hackl, A discontinuous phase field approach to variational growth-based topology optimization, Struct. Multidiscip. Optim. (2016) 1–14. [24] M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods and Applications, Springer, 2003. [25] P. Junker, J. Makowski, K. Hackl, The principle of the minimum of the dissipation potential for non-isothermal processes, Contin. Mech. Thermodyn. 26 (3) (2014) 259–268. [26] F. van Keulen, R. Haftka, N. Kim, Review of options for structural design sensitivity analysis. Part 1: Linear systems, Comput. Methods Appl. Mech. Engrg. 194 (30–33) (2005) 3213–3243. Struct. Des. Optim. http://dx.doi.org/10.1016/j.cma.2005.02.002. URL http://www.sciencedirect.com/science/article/pii/S0045782505000459. [27] F.-J. Barthold, A short guide to variational design sensitivity analysis, in: 6th WCSMO–World Congress on Structural and Multidisciplinary Optimization, Citeseer, 2005. [28] K. Schittkowski, C. Zillober, R. Zotemantel, Numerical comparison of nonlinear programming algorithms for structural optimization, Struct. Optim. 7 (1–2) (1994) 1–19. [29] T. Borrvall, J. Petersson, Large-scale topology optimization in 3d using parallel computing, Comput. Methods Appl. Mech. Engrg. 190 (46–47) (2001) 6201–6229. http://dx.doi.org/10.1016/S0045-7825(01)00216-X. URL http://www.sciencedirect.com/science/article/pii/S004578250100216X.