Topology optimization of freely floating elastic continua using the inertia relief method

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. xxx (xxxx) xxx www.elsevier.com/locate/cma

Topology optimization of freely floating elastic continua using the inertia relief method Carl-Johan Thore Division of Solid Mechanics, Linköping University, SE 581 83 Linköping, Sweden Received 18 December 2018; received in revised form 22 September 2019; accepted 5 November 2019 Available online xxxx

Abstract In many applications, it is of interest to perform static finite element analyses on freely floating bodies that are not in quasistatic equilibrium; airplanes and helicopters maneuvering in flight for example. This is particularly so if topology optimization (TO) is to be used, since TO with dynamic analyses can be very computationally expensive. The so-called inertia relief method, which essentially entails computing the rigid body inertia and subtracting it from the given loads to make the system of loads self-equilibrating, can sometimes be used to replace a dynamic analysis with a static, thus enabling the use of high-resolution TO. We derive the inertia relief method for elastic continua and obtain a static variational problem which require that we can suppress (linearized) rigid body motions without affecting the deformation. Three methods for doing this are investigated. Based on the static variational problem we consider maximizing stiffness using TO. Numerical examples show that all three methods for suppressing rigid body motion work, and indicate that optimal designs for freely floating structures undergoing rigid acceleration can differ significantly from designs optimized under static conditions. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Topology optimization; Elastic continua; Freely floating; Inertia relief; Singular systems

1. Introduction Topology optimization (TO) [6,16] is a computer-based method for automatic design of maximally efficient structures. Development of TO for elastic continua dates back to the late 1980s [5]. Nowadays, this field of research attracts wide interest both from the academy (recent reviews include [22,48]) and from the industry, as evidenced by several papers on industrial applications of TO (see e.g. [58]) and the existence of commercial software such as TOSCA (Dassault Systems) and Altair OptiStruct. In this work, we consider TO of linearly elastic continua that are freely floating in the sense that there are no supports anchoring the body to a fixed frame. Such models arise when considering for instance airplanes, helicopters or rockets. Designing for scenarios such as landing or lift-off of an airplane, under which the plane is not in (quasi-) static equilibrium seems to necessitate the use of dynamics analyses. However, TO typically requires hundreds or even thousands of finite element (FE) analyses and therefore computationally costly dynamic analyses is presently not a viable option, at least not for high-resolution models. A potential remedy is offered by the so-called inertia relief method, which may enable approximation of dynamic problems with static problems. E-mail address: [email protected]. https://doi.org/10.1016/j.cma.2019.112733 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

2

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

The basic idea behind inertia relief is to compute (an approximation of) the rigid body inertia and subtract it from the given loads to achieve a self-equilibrating system of loads. In standard FE-notation (cf. Section 5) the equilibrium equation with inertia relief is K u = f − M ar , where ar contains the approximate rigid body, translational and rotational, accelerations caused by the application of the given load f . A graphical illustration is given in the bottom panel of Fig. 1 in Section 6 below. Inertia relief is often available as an option in FE software, and its use for FE analysis has been treated in the engineering literature [38,43]. It may be mentioned that inertia relief is sometimes also used by engineers for bodies which are not unconstrained but where only the forces acting on the body are available, e.g. from measurements or from multibody dynamics simulations [46]. As for TO in combination with inertia relief there is not much published [17,41,50], and nothing particularly detailed. To the best of our knowledge, inertia relief has not been treated in the mathematical literature, but mathematical modeling of freely floating continua leads to so-called singular Neumann problems, and such problems have attracted much interest [7,8,33–35,42,47,57]. As shown herein, inertia relief can from a mathematical point of view be seen as a way to ensure existence of solutions to the singular (elastic) Neumann problem. In the mathematical literature, conditions for existence of solutions are assumed to be satisfied, and the focus is instead on how to handle the fact that solutions are not unique. Non-uniqueness is a central issue also herein, and is further discussed in Section 2.3 below. The present paper provides a detailed derivation of the inertia relief method for elastic continua. The resulting static variational problem (20) is well-posed, but require that we can suppress rigid body motions without affecting the deformation. To handle this, we consider a relaxed version of the problem. The use of inertia relief implies that the relaxed problem is solvable, but not uniquely so since an arbitrary (linearized) rigid body motion can be superimposed to get a new solution. To remedy the situation we investigate three different methods to obtain a unique solution, each method suitable for use with iterative linear systems solvers necessary for handling large-scale TO problems. The considered TO problem is the common problem of maximizing stiffness with an upper bound on the total volume [6,16], but here with a freely floating design domain undergoing rigid acceleration. Finite element-discretized versions of the TO problems are then solved in Section 6 using a multi-grid pre-conditioned conjugate gradient (CG) method to solve the equilibrium equations. The examples indicate the relative performance of the above-mentioned three methods to handle non-uniqueness (manifested in a singular stiffness matrix) and show optimized designs for freely floating design domains undergoing rigid acceleration. 1.1. Notation and basic facts We shall only be concerned with standard L p and Sobolev spaces [1]. Spaces of vector- or matrix-valued functions in which each component is in a function space X is denoted by bold-face X. Specifically, we use, for measurable sets Ω ⊂ Rn , the space L ∞ (Ω ) of essentially bounded functions; the space L 2 (Ω ) with inner product and norm ∫ √ (v, w)Ω = v : w, ∥u∥ L 2 (Ω) = (u, u)Ω , Ω

where : is the Frobenius inner product; and the Sobolev space H 1 (Ω ). Throughout the article, we do not explicitly specify the measure symbol in integrals unless there is a risk of confusion, and we often write (v, w) instead of (v, w)Ω when the domain is obvious. The Euclidean norm on Rn is denoted by ∥ · ∥. For a set Ω ⊂ Rn , |Ω | denotes its Lebesgue measure, and Ω its closure. If A ∈ Rd×d , d = 2, 3, is skew-symmetric, i.e. A = − AT , then there exists a ∈ Rd such that Ax = a × x, where × is the cross product. For d = 2, a × x should be interpreted as the scalar function a1 x2 − a2 x1 . In the appendices, c denotes a generic positive constant whose value may change from line to line. 2. The state problem 2.1. Design parametrization We consider here the material distribution method for TO [6]. The design domain is a connected, bounded and open set Ω ⊂ Rd , where d = 2 or 3, with Lipschitz boundary Γ , and the design is parameterized by a scalar field Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

3

ρ ∈ L ∞ (Ω ) indicating where to put material (ρ = 1) and where not to (ρ = 0). This field may be interpreted as relative density. The set of admissible designs is defined as ∫ { } H = ρ ∈ L ∞ (Ω ) | 0 ≤ ρ ≤ 1 a.e. in Ω , S(ρ) ≤ γ |Ω | , (1) Ω

where a.e. is short for almost everywhere, γ ∈ (0, 1), and |Ω | is the area (2D) or volume (3D) of the design domain (for the 2D problems we assume, as usual, a constant thickness). The structural stiffness is not taken to depend directly on ρ but rather on the “physical design” ρ˜ : L ∞ (Ω ) → ∞ L (Ω ) defined as ρ(ρ) ˜ = ρ + g(S(ρ))(1 − ρ)

∀x ∈ Ω

(2)

where the positive constant ρ ≪ 1 is the relative density of “void”. To enforce “black-and-white” solutions, the so-called Rational Approximation of Material Properties (RAMP) scheme [51] is used, so ρ , (3) g(ρ) = 1 + q(1 − ρ) where q > 0. To avoid problems such as non-existence of solutions to the infinite-dimensional problem and checkerboards for certain discretizations [49] we use a linear design filter [10] defined by ∫ S(ρ) = ρ( y)Ψ (x, y)d y ∀x ∈ Ω , (4) Ω

where the kernel is chosen as max (0, R − ∥x − y∥) , Ψ (x, y) = ∫ Ω max (0, R − ∥x − z∥) d z

(5)

in which R > 0 is the filter radius. This choice of kernel implies that 0 ≤ S(ρ) ≤ 1 for all x ∈ Ω and ρ ∈ H [9]. The density of the material is given by ρ M (ρ) = ρsolid S(ρ),

(6)

where the constant ρsolid > 0 is the absolute density of the solid material. The reason for using linear interpolation in this expression, and not the RAMP scheme, is to reduce problems with “gray” regions of intermediate design values in optimized designs; cf. [37, Fig. 1] and [14, Fig. 11]. (These works used SIMP, but we have found numerically that RAMP works better.) The appearance of such regions has been attributed to large ratios between stiffness and mass in low relative density regions, leading to very large deformations [roughly proportional to ρ(ρ)/ρ ˜ M (ρ)] [14]. To minimize such deformations, the optimization solver adds small amounts of material, generating “gray” regions in the design domain. 2.2. The state problem Let I = (0, T ) be a time interval, and let two superposed dots denote a second derivative with respect to time. Neglecting damping and friction, the weak form of the equation of motion is ¨ v) + a(ρ(ρ), (ρ M (ρ)u, ˜ u, v) = ℓs (v) ∀(v, t) ∈ H 1 (Ω ) × I.

(7)

Using Hooke’s law, the internal elastic work can be written as ∫ ∫ a(ρ(ρ), ˜ u, v) = ρ(ρ)σ ˜ (u) : ε(v) = ρ(ρ)Eε(u) ˜ : ε(v),

(8)

Ω

Ω

where σ (u) = Eε(u) is the stress tensor, and ε(u) = 12 (∇uT + ∇u) is the linearized strain tensor. The components of the fourth-order elasticity tensor E are assumed to lie in L ∞ (Ω ) and satisfy the symmetry conditions E i jkl = E jikl = E kli j . It is also assumed that there exists β > 0 such that E(x)ε : ε ≥ βε : ε, for all symmetric ε ∈ R [E(x)ε]i j =

(9)

d×d

d ∑ d ∑

and almost all x ∈ Ω . The components of E(x)ε are given by

E i jkl (x)εkl .

k=1 l=1 Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

4

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

The load functional in (7) finally, is assumed to be constant with respect to time, linear and continuous over H 1 (Ω ); it is here taken as ℓs (v) = (t, v)Γ ,

(10) 2

where the surface load t ∈ L (Γ ). Since the right-hand side in (7) is constant in time (or can be assumed constant over I ), we can under some simplifying assumptions reduce the equation to a static one. To this end, first note that since (9) holds and ρ > 0 a.e., the internal work is zero if and only if ε(u) = 0. To see this, take u such that a(ρ(ρ), ˜ u, v) vanishes for all test functions. Then, ∫ ∫ 0 = a(ρ(ρ), ˜ u, u) = ρ(ρ)Eε(u) ˜ : ε(u) ≥ ρβ (11) ε(u) : ε(u) = ρβ∥ε(u)∥2L 2 (Ω) ≥ 0, Ω

Ω

implying ε(u) = 0 a.e. in Ω . The functions satisfying this conditions are referred to as linearized1 (or sometimes infinitesimal) rigid body motions, and form the space [19, Theorem 6.15-2] R M = {u ∈ H 1 (Ω ) | ε(u) = 0 in Ω } = {u ∈ H 1 (Ω ) | u(x) = Ax + b, A = − AT ∈ Rd×d and b ∈ Rd }.

(12)

This is a finite-dimensional subspace of H 1 (Ω ), with dimension p = 21 d(d + 1). Being finite-dimensional, R M is a closed subspace of H 1 (Ω ), and thus also of L 2 (Ω ). It follows that any displacement in H 1 (Ω ) can be uniquely decomposed as [19, Theorem 4.5-2] u = ur + u d

(13)

∀t ∈ I,

where ur ∈ R M is the rigid body motion and ud ∈ R M Theorem 4.5-2]

⊥

the deformation. The latter space is given by [19,

R M ⊥ = {u ∈ H 1 (Ω ) | (u, v) = 0 ∀v ∈ R M},

(14)

2

where we recall that (·, ·) is the L (Ω ) inner product. Note that since R M only approximates the space of “true” rigid body motions the terms “rigid body motion” and “deformation” are not entirely accurate. We shall nevertheless use them in the following. Now substituting (13) into (7) gives (temporarily suppressing the dependence on ρ) (ρ M u¨ d , v) + (ρ M u¨ r , v) + a(ud , v) + a(ur , v) = ℓs (v) ∀(v, t) ∈ H 1 (Ω ) × I. Since the right-hand side in this equation is constant we assume that any variations of deformation with time is eventually damped out, so that u¨ d = 0. Then noting that a(ur , v) = 0, the equation reduces to (ρ M (ρ)u¨ r , v) + a(ρ(ρ), ˜ ud , v) = ℓs (v) ∀(v, t) ∈ H 1 (Ω ) × I

(15)

from which u¨ r must somehow be eliminated to get a static problem. To this end we note that for all test functions v ∈ R M, (15) reduces to (ρ M (ρ)u¨ r , v) = ℓs (v).

(16)

Seeking to express u¨ r in terms of known quantities we first assume a basis for R M is given. In the case d = 3 for instance in the form of the functions e1 , e2 , e3 , x × e1 , x × e2 , x × e3 ,

(17)

where ei denotes the standard basis vectors in R3 (Appendix A). The basis functions are denoted by z i = z i (x), i = 1, . . . , p. Now substituting v(x) =

p ∑ j=1

1

z j (x)d j

and

u¨ r (x, t) =

p ∑

z i (x)αi (t)

(18)

i=1

They approximate the “true” rigid body motions (see e.g. [31, section 12.3]) for small rotations.

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

5

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

in (16) yields, after some rearrangements, p p p ∑ ∑ ∑ (ρ M z i , z j )αi d j = ℓs (z j )d j . i=1 j=1

j=1

Since the coefficients d j are arbitrary, this is equivalent to requiring that p ∑ (ρ M (ρ)z i , z j )αi = ℓs (z j ),

j = 1, . . . , p.

(19)

i=1

This system of p equations for the rigid body accelerations αi is uniquely solvable for every ρ ∈ H such that ρ > 0 in some open set ω ⊂ Ω (Appendix B), i.e. whenever there is at least some mass in the design domain. With the basis (17), α1 , α2 and α3 are translational accelerations; and α4 , α5 and α6 are angular accelerations. It is clear from (19) that α = (α1 . . . α p )T is constant in time, so substituting the second of (18) into (15), and dropping the subscripted d on ud , leads to the static variational problem find u ∈ R M ⊥ :

a(ρ(ρ), ˜ u, v) = ℓ(ρ M (ρ), v) ∀v ∈ R M ⊥ ,

(20)

in which p ∑ ℓ(ρ M (ρ), v) = ℓs (v) − (ρ M (ρ)z i , v)αi (ρ M (ρ)),

(21)

i=1

where αi (ρ M (ρ)) solves (19). Existence and uniqueness of solutions to (20) is proved in Appendix C. The last term in (21) can be described as the inertia relief correction, and is the (virtual) work done by the approximate inertial load (density times acceleration) ρ M (ρ)

p ∑

z i αi (ρ M (ρ)),

i=1

acting in each point in the design domain. This load constitutes a body load that counteracts the given loads such that the body is in equilibrium under the total load; cf. e.g. the bottom plot in Fig. 1. Remark 1. Since ρ > 0, the only rigid body motions that can occur are for the design domain Ω as whole; i.e. we do not get parts inside Ω freely floating relative to each other; see e.g. [13,54] for this case. 2.3. Suppressing rigid body motions The fact that the solution and the test functions in (20) belong to R M ⊥ means that we must somehow suppress rigid body motions without affecting the deformations. For sake of discussion we consider first the relaxed problem find u ∈ H 1 (Ω ) :

a(ρ(ρ), ˜ u, v) = ℓ(ρ M (ρ), v) ∀v ∈ H 1 (Ω ).

(22)

It follows from the existence of solutions to (20) that this problem is solvable (for every ρ ∈ H) if and only if ℓ(ρ M (ρ), v) = 0

∀v ∈ R M,

(23)

a condition which holds here because of the inertia relief correction in (21). This condition is equivalent to a linearized version of Euler’s laws of motion written for the entire design domain (Appendix D), so the physical interpretation is that (22) is solvable because the body is in static equilibrium. Although problem (22) is solvable and the deformation for a given load is unique, the total solution is not unique since one can always add a function from R M and get a new solution. The practical implication is that the stiffness matrix arising from a standard Galerkin FE discretization is singular. There exist several different methods for handling the non-uniqueness [34] and the implied singularity of the stiffness matrix. We shall here consider three such methods, all of which can be used with iterative linear solvers, a necessity for large-scale TO problems. Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

6

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

2.3.1. The CG method applied to the singular problem The first method is to apply the iterative linear solver, here the CG method, directly to the discrete, singular equilibrium equation (see (37) below). Convergence of the method can be shown [30] provided that the load vector is orthogonal to the nullspace of the stiffness matrix (this is ensured by using the inertia relief method), the initial guess has no component in the kernel of the stiffness matrix, and a positive definite preconditioner is used. Iterative linear solvers applied to singular systems in TO has previously been studied by Washizawa et al. [54]. 2.3.2. Pin-pointing Another method, sometimes referred to as pin-pointing, is to prescribe p nodal degrees of freedom in the FE mesh so as to prevent rigid body motions while leaving the deformation unaffected. Let u denote the solution to (20) and uh the solution to the FE-discretized state problem (29) given below. Then combining Cea’s lemma, i.e. ∥u − uh ∥ H 1 (Ω) ≤ c

inf

v h ∈V h ∩R M ⊥

∥u − v h ∥ H 1 (Ω) ,

with Lemma 19 of Appendix G shows that a convergent FE method (c.f. also [23, Proposition 3.84]) can be obtained even if u is only in H 1 (Ω ), hence is not necessarily continuous. The only potential problem is ill-conditioning of the stiffness matrix resulting from prescribing point values of the FE displacement [7, Theorem 5.2]. This problem may be exacerbated in TO, where supports can end up in parts of Ω with low relative density (ρ ≈ ρ). In applications however, one would either use a direct solver or apply a preconditioner to an iterative method, so it is not clear from the outset that ill-conditioning is a practical problem. The fact that pin-pointing can be used with standard direct solvers is of course an advantage compared to the other two methods that we discuss. 2.3.3. Penalty Uniqueness can also be obtained by adding a term to the internal work which penalizes rigid motions without affecting the deformation [8,34]. For a specific choice of penalty term the modified internal work then reads a p (ρ(ρ), ˜ u, v) = a(ρ(ρ), ˜ u, v) +

p ∑ (u, z i )ω (v, z i )ω .

(24)

i=1

where ω is taken as Ω following [34]. The associated state problem is to find u ∈ H 1 (Ω ) :

a p (ρ(ρ), ˜ u, v) = ℓ(ρ M (ρ), v)

∀v ∈ H 1 (Ω ).

(25)

As shown in Appendix E, this problem is well-posed and provides a solution to (20). A drawback of this penalty method is that it leads to a dense stiffness matrix (see Section 5). This precludes the use of direct linear solvers, but matrix–vector products can still be computed efficiently so that iterative solvers can be applied. A positive note is that there is no additional free parameter introduced with the penalty term in (24). Remark 2. While solving the singular system using the CG method or solving the penalized problem (25) both give a solution in R M ⊥ , i.e. a deformation with no superimposed rigid motion, the pin-pointing method yields the same solution but with a superimposed non-zero rigid body motion. However, if desired it is straightforward to compute and subtract the part of the solution in R M, so therefore we may in practice treat the solution obtained with the pin-pointing method as being in R M ⊥ . Note that this implies that the choice of the p constraints for the pin-pointing method is arbitrary; as will be seen in Section 6, even constraining nodes that happen to end up in “voids” (the location of which are generally not known before-hand) is not an issue. The numerical performance of the three methods for suppressing rigid body motions will be assessed in Section 6. 3. The design problem In order to maximize stiffness we minimize the so-called compliance, a commonly used inverse measure of average stiffness [6,16]. The design problem thus reads inf ℓ(ρ M (ρ), u(ρ)),

ρ∈H

(26)

where u(ρ) solves (20), and the explicit design-dependency in the objective comes from the inertia relief correction in (21). Existence of a globally optimal solution to problem (26) is proved in Appendix F. Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

7

4. Finite-dimensional approximations Finite-dimensional approximations of the design problem (26) are obtained using a standard Galerkin FE method with an element-wise constant approximation for the design. In the numerical examples in Section 6 we use rectangular design domains meshed with equal-sized squares or cubes, but the theory in Appendix G is not complicated by generalizing to polyhedral domains and standard, affine-equivalent Lagrange C 0 -elements forming m quasi-uniform [24, { Def. 1.140] meshes} {Ωe }e=1 with convex elements. Here m = m(h), where the mesh parameter h = maxe∈{1,...,m} max x, y∈Ω e ∥x − y∥ . The approximation spaces for the displacement are based on the space V h = {v ∈ H 1 (Ω ) | vi ∈ C 0 (Ω ) : vi |Ωe ∈ P(Ωe ), e = 1, . . . , m}, where P(Ωe ) denotes a space of polynomials on Ωe . For the pin-pointing method of Section 2.3.2 the approximation space is of the form {v ∈ V h | vi (x j ) = 0 for some i, j}, where displacement components i and nodes j are chosen to prevent rigid motion while leaving the deformation unaffected (cf. [23, p. 158]). The total number of point-wise constraints on displacement components must of course be exactly p since this is the dimension of R M. The design is approximated as element-wise constant, i.e. we use the approximation space L h = {ρ ∈ L ∞ (Ω ) | ρ|Ωe = ρe , e = 1, . . . , m}. The constant values ρe , e = 1, . . . , m, are collected in a vector ρ ∈ Rm , which then uniquely determines a function ρh ∈ L h . The filtered design S(ρh ) is also approximated as an element-wise constant function defined by ∑ f ∈N ρ f |Ω f |Ψ (x e , x f ) , e = 1, . . . , m. (27) Sh (ρh )|Ωe = ∑ e f ∈Ne |Ω f |Ψ (x e , x f ) Here x e is the centroid of element e and Ne is the set of elements whose centroid is at a distance from x e less than R. It is readily verified that 0 ≤ Sh (ρh ) ≤ 1 a.e. provided ρ f ∈ [0, 1], f = 1, . . . , m. The finite-dimensional approximation of (26) at discretization level h now reads inf

ρh ∈Hh

ℓ(ρ Mh (ρh ), uh (ρh )),

(28)

where ∫ } { Sh (ρ) ≤ γ |Ω | , Hh = ρ ∈ L h | 0 ≤ ρ ≤ 1 a.e. in Ω , Ω

and uh (ρh ) solves the discretized state problem a(ρ˜h (ρh ), uh , v h ) = ℓ(ρ Mh (ρh ), v h ) ∀v h ∈ V h ∩ R M ⊥ ,

(29)

where ρ˜h (ρh ) = ρ + g(Sh (ρh ))(1 − ρ),

(30)

ρ Mh (ρh ) = ρsolid Sh (ρh ).

(31)

and

Existence of solutions to problem (28) follows from the equivalence of this problem with the matrix–vector problem (39) below which admits a globally optimal solution by Weierstrass’ theorem [4]. Regarding the relation between problems (28) and (26) it is shown in Appendix G that any sequence {ρh∗ } of globally optimal solutions to (28) has a subsequence converging to a global solution to (26). Since the sequence {ρh∗ } is arbitrary, it follows that every sequence converges to a solution to (26). Analogous results are expected to hold for sequences of stationary points [25,26]. Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

8

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

5. The matrix–vector problem The dimension of the approximation space V h is n = N d, where N is the number of nodes in the FE mesh. n Given a (vector) basis {ϕ i }i=1 and a function v h ∈ V h we denote by vˆ = (vˆ1 , . . . , vˆn )T ∈ Rn the coordinate vector relative to this basis, i.e. v h (x) =

n ∑

ϕ i (x)vˆi =

N ∑

i=1

ϕ j (x)ˆv j ,

j=1

where the vector basis is related to the scalar basis via ϕ d(i−1)+k = ϕi ek , i = 1, . . . , N , k = 1, . . . , d. The space R M defined in (12) consists of linear functions only. Therefore R M ⊂ V h , and the basis functions can be represented exactly as N ∑

z j (x) =

ϕi (x)ˆz ji ,

i=1

where zˆ ji = z j (x i ) ∈ Rd . With the basis (17) specialized to d = 2 the coordinate vectors become zˆ 1 = (1 0 1 0 . . . 1 0 )T zˆ 2 = (0 1 0 1 . . . 0 1 )T

(32)

zˆ 3 = (−y1 x1 − y2 x2 . . . − y N x N ) , T

where x1 , y1 , . . . , x N , y N are coordinates of the nodes in the FE mesh. Introducing Z = (ˆz 1 . . . zˆ p ) ∈ Rn× p , substitution of the FE functions into (19) and some algebra yields Z T M(ρ)Zα = Z T f ,

(33)

where the mass matrix is defined as [M(ρ)]i j = (ρ Mh (ρ)ϕ i , ϕ j ) for all i, j = 1, . . . , n, in which ρ Mh (ρ) = ρ Mh (ρh (ρ)); and )T (∫ ∫ T T t ϕn . t ϕ1 . . . f = Γ

Γ

Starting from the discrete version of (25), i.e. p p ∑ ∑ a(ρ˜h (ρ), uh , v h ) + (uh , z i )(v h , z i ) = ℓs (ρ Mh (ρ), v h ) − (ρ Mh (ρ)z i , v h )αi i=1

∀v h ∈ V h ,

i=1

where ρ˜h (ρ) = ρ˜h (ρh (ρ)), we derive the equation [K (ρ) +

p ∑

M 0 zˆ i (M 0 zˆ i )T ]uˆ = f − M(ρ)Zα,

(34)

i=1

in which the symmetric and positive semi-definite stiffness matrix K (ρ) ∈ Rn×n is defined by [K (ρ)]i j = a(ρ˜h (ρ), ϕ i , ϕ j ) ∀i, j = 1, . . . , n,

(35)

1 M(1), in which 1 is an m × 1 vector of ones. and M 0 = ρsolid Now define the inertia relief-corrected right-hand side as

f α (ρ) = f − M(ρ)Zα(ρ),

(36)

where α(ρ) solves (33). Omitting the penalty term, the equilibrium equation (34) then becomes K (ρ)uˆ = f α (ρ),

(37)

where the singular matrix K (ρ) is replaced by K pen (ρ) = K (ρ) +

p ∑

M 0 zˆ i (M 0 zˆ i )T

(38)

i=1 Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

9

for the penalty approach of Section 2.3.3 and K p− p (ρ) when using pin-pointing. Here the latter matrix is obtained from K (ρ) by zeroing all entries in the p rows and p columns associated with the prescribed displacement components and setting the corresponding diagonal entries to one. The right-hand side f α (ρ) is also modified by setting the entries corresponding to prescribed displacement components to zero. An alternative is to reduce the size of the system, but this results in a more complex code for the multi-grid preconditioner described below. The matrix–vector version of the design problem (26) reads ˆ min f α (ρ)T u(ρ),

(39)

ρ∈Hd

ˆ where u(ρ) solves (37) and m { } ∑ Hd = ρ ∈ Rm | ρe ∈ [0, 1], e = 1, . . . , m, ρ˜e (ρ)/|Ω | ≤ γ ,

(40)

e=1

where ρ˜e (ρ) =

∑

He f = ∑

f ∈Ne

He f ρ f , with

|Ω f |Ψ (x e , x f ) . f ∈Ne |Ω f |Ψ (x e , x f )

Division by |Ω | in the volume constraint in (40) is done to improve numerical performance (the code for solving the subproblems in the Method of Moving Asymptotes [52] used here is sensitive to scaling of the constraints). For ˆ 0 ))−1 , where ρ 0 is the starting guess the same reason we also scale the objective function in (39) with ( f α (ρ 0 )T u(ρ for the design, in the computer implementation. Remark 3. It is relatively straightforward to show (see Appendix H) that the vectors zˆ j , j = 1, . . . , p, recall (32), constitute a basis for the null space N (K ) of K . Using (36) and (33) we get Z T f d (ρ) = Z T ( f − M(ρ)Zα(ρ)) = Z T f − Z T f = 0, i.e. the inertia relief corrected load vector is orthogonal to N (K ). This guarantees that (37) is solvable (though not uniquely so) for every admissible design. Introducing P = I − M Z(Z T M Z)−1 Z T , where I is an identity matrix, we can write f α (ρ) = P f . The action of (the non-orthogonal projection) P is to project vectors in Rn onto the range space of K . If the vector is interpreted as a load vector, then the physical interpretation of this is that P modifies the vector to ensure static equilibrium of the structure. 5.1. Solving the equilibrium equation For large-scale TO problems, the stiffness matrix in (37) is very large. Therefore, although, except in its penalized version (38), it is sparse, it is necessary to use iterative linear solvers. Since all of the stiffness matrices used here are symmetric and positive semi-definite we use a pre-conditioned CG method to solve (37). Following recent work [2,32] on iterative solvers in TO we use one step of a V-cycle multigrid method as a preconditioner. Given a function mg = mg( A, b) representing the multi-grid preconditioner, the CG method for a given design ρ can now be described as follows [40] (here Kˆ is one of the three versions of the stiffness matrix considered): ˆ Let r = Kˆ uˆ − f α and compute − p = y = mg( A, r). Then Solve (33) and define an initial guess u. while true κ1 = r T y/ pT Kˆ p uˆ = uˆ + κ1 p r˜ = r + κ1 Kˆ p ˜y = mg( A, r˜ ) κ2 = r˜ T ˜y/r T y p = − ˜y + κ2 p check stopping criteria r = r˜ , y = ˜y end Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

10

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

The matrix A supplied to the preconditioner is for the pin-pointing method taken as K p− p , whereas for the other two methods it is taken as K + M, a choice suggested in [34] and which we found to work better than K p− p . For the penalty method, the matrix products in the listed algorithm are computed as Kˆ p = K p +

p ∑ (M 0 zˆ i )((M 0 zˆ i )T p),

(41)

i=1

where the vectors M 0 zˆ i , i = 1, . . . , p, are computed and stored once at the beginning of the optimization process. The initial guess uˆ for the CG method is taken as the equilibrium solution from the previous optimization iterate (and 0 at the beginning of the optimization). As stopping criteria we use ∥˜r ∥ < CGtol∥ f α ∥ and a maximum number of iterations. We remark that early termination of CG methods in TO has been subject to some research [2,3,39], but such investigations are beyond the scope of the present paper. When applying the CG method to the singular system, as discussed in Section 2.3.1, it occasionally happened that the method failed, despite the linear system being well-scaled, to converge to the desired tolerance within a reasonable amount of iterations (<1000 say). The problem was found to be existence of a fluctuating, nonvanishing null-space component of the residual. As a remedy, we remove any null-space components by setting r˜ = r˜ − Z(Z T Z)−1 Z T r˜ just prior to the application of the preconditioner in the algorithm listed above. By doing this the first and then every 50:th CG iteration, the CG method never failed in any of our numerical experiments. 6. Numerical examples The code for solving problem (39) is implemented in Matlab (R2018b). As optimization solver we use the Method of Moving Asymptotes [52], version 2007 [53] with default settings. The implementation of the CG method and the multi-grid preconditioner for solving the equilibrium equation (37) is based on the Matlab code by Amir et al. [2]. The derivative of the objective, f = f (ρ), in (39) with respect to optimization variable ρi is given by ( ) m ∑ ∂M ∂f ∂K ∂ f ∂ ρ˜e ∂ f ∂M ˆ u, = , = −2 uˆ T Zα − λT Z T Zα − uˆ T ∂ρi ∂ ρ ˜ ∂ρ ∂ ρ ˜ ∂ ρ ˜ ∂ ρ ˜ ∂ ρ˜e e i e e e e=1 ˆ and ∂ ρ˜e /∂ρi = Hei . where the adjoint vector λ solves the (very small) p × p system Z T M Zλ = Z T M u, In all the examples, the volume fraction is set to γ = 0.3, and the relative density of void to ρ = 10−9 . The material is assumed to be isotropic, with Young’s modulus E = 1 and Poisson’s ratio ν = 0.3. A state of plane stress is assumed for the 2D examples, for which the thickness is set to 0.1. The absolute density of the solid material is set to ρsolid = 1. The penalty parameter q in (3) is set to 4. The multi-grid preconditioner uses three grids and a damping factor for the Jacobi smoother set to 0.8 for the 2D examples and 0.6 for the 3D example [2]. The stopping tolerance for the CG method is CGtol = 10−10 and the maximum number of CG iterations 1000. The numerical calculations are carried out by an Intel Core i7-4712MQ CPU. The examples below are solved using the three methods for suppressing rigid body motions introduced in Section 2.3: CG on the singular system, pin-pointing, and penalty. For the pin-pointing we lock all displacement components for the FE node closest to the geometric center of Ω and one component in either one (2D) or three (3D) neighboring nodes to prevent rotation. 6.1. MMB-beam Our first example is based on the so-called MBB-beam shown in Fig. 1. The rectangular design domain Ω is of size 8 times 2 and the FE mesh consists of 416 × 104 equal-sized, four-noded, bi-linear elements. The filter radius R is set to 0.6% of the domain width. As a reference, the top plot shows a “standard”, static minimum compliance design where the left corner is fixed and the right is prevented from moving vertically. The middle plot shows a design obtained with the supports in the lower corners removed so that the design domain is freely floating. The reaction forces in the supports are replaced by prescribed loads with magnitude κt, where κ > 0 is chosen such that the design domain accelerates upwards. The bottom plot in Fig. 1 shows the total load acting on the design: the prescribed loads in black and the inertia relief load −M(ρ)Zα(ρ) acting as a body load in red. The presence of the inertia relief load (mass times acceleration) means that the total system of loads is self-equilibrating, thus enabling static FE-analyses to be performed. Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

11

Fig. 1. MBB-beam-like examples. Top: Reference design optimized with fixed displacements at the lower corners. Middle: The supports are removed and replaced by prescribed loads causing an upwards acceleration. The picture shows an optimized design for this setup. Bottom: Total load f α (ρ) on the freely floating design consisting of prescribed loads f in black and computed inertia relief load −M(ρ)Zα(ρ) in red. The magnitude of the components in the inertia relief load is at least 1000 times smaller than for the components of f ; hence the black and red arrows are not drawn to scale. The small red dots are inertia relief loads in the “void” where ρ ≈ ρ. To avoid cluttering, the inertia relief load is only plotted for a few of the 43 785 nodes in the domain. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 shows further examples based on the same setup but with increasingly large upwards acceleration. Compared to the middle plot of Fig. 1 which shows a design quite similar to the static reference design, the designs in Fig. 2 are markedly different. There is a clear tendency, in particular for the lower design, to move material away from the middle parts towards the left and right sides of the design domain. This can be explained by noting that this leads to large inertial loads (−M(ρ)Zα(ρ)) pointing downwards to counteract the upwards-acting prescribed loads, thus reducing the compliance in these regions. The downwards-acting prescribed load, applied in the upper middle part of the domain, is now comparatively low, making it unfavorable to place material in the inner parts of the design domain. Another notable difference between the designs in Fig. 2 and the static reference design is the curvature of some of the structural members. This effect is due to the inertia relief load acting as a body load; a similar effect is for example seen in works on TO with self-weight. We mention in particular Fig. 17 in [14] showing results for an MBB-setup similar to the one used here. There the same effect of material being removed for the middle parts and some of the members becoming more curved as the relative importance of the self-weight load increases can be seen. A difference compared to the designs in Fig. 2 is that in [14] there is no large build-up Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

12

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 2. Additional MBB-beam-like examples with increasing upwards acceleration going from top to bottom. The plots show designs and corresponding inertia relief loads for two different upwards loads.

of “pillars” at the sides of the design domain. This difference is explained by the fact that in [14], the upwards acting forces are reaction forces developed at the supports and adjusting to the loading and material distribution, whereas in our case, where the structure is freely floating without supports, the upwards acting forces are prescribed and will not adjust to other loads, resulting instead in an adjustment of the material distribution. The designs in the plots of Fig. 2 look unusual when compared with “standard” static designs such as the one in Fig. 1 (top), possibly raising some concerns about the optimality of the results. To assess this, we take the design, Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

13

Fig. 3. MBB-beam-like examples with prescribed loading inducing rotational acceleration. Note the slight curvature of the inertia relief load (red arrows) in the bottom plot. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 1 Data for three different starting points for an MBB-type setup. Columns show (i) minimum objective function value; (ii) maximum objective function value; (iii) average number of CG iterations; and (iv) average time per CG iteration. Method

min. obj.

max. obj.

CG iter.

CG time [s]

CG on sing. system Pin-pointing Penalty

0.99 0.99 0.99

1.00 1.00 1.00

69 27 121

1.09 0.46 1.93

call it ρ s , obtained in the static case in Fig. 1 (top) and subject it to the prescribed loading illustrated in Fig. 2 (bottom). The total load on the static design thus becomes f α (ρ s ) = f − M(ρ s )Zα(ρ s ). The compliance for this design is then almost two times higher than for the dynamic design in Fig. 2 (bottom), thus indicating the optimality of the latter design. Fig. 3 shows a case where lower left load is larger (κ1 = 1.3κ) than the lower right (κ2 = 0.7κ), thus inducing (clockwise) rotational acceleration. A large fraction of material is now placed at the left side to counteract the effect of the prescribed load acting there. Compared to Figs. 1 and 2, the inertia relief load in Fig. 3 (bottom) now has a slight curvature due to the rotational acceleration. Table 1 shows a comparison of the performance of the three methods for suppressing rigid body motions introduced in Section 2.3 when used in TO. The design domain and prescribed loading is the one used to obtain the bottom plot of Fig. 2. For each method, the TO problem is solved using three initial designs as starting guess: one uniform, feasible design and two designs obtained from a uniform distribution over (0, 1)m . One of the random starting points is re-scaled to satisfy the volume constraint in (40), while the other one does not satisfy this constraint. The optimization runs are carried out with a maximum computational time of 0.5 h and 1000 MMA iterations as the only stopping criteria. Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

14

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 4. Objective function value versus iteration number for the three different starting points. Table 2 Data for three different starting points for the 2D cantilever example. Columns show (i) minimum objective function value; (ii) maximum objective function value; (iii) average number of CG iterations; and (iv) average time per CG iteration. Method

min. obj.

max. obj.

CG iter.

CG time [s]

CG on sing. system Pin-pointing Penalty

0.98 0.98 0.98

0.99 0.99 1.00

20 21 32

2.32 2.24 3.80

The first thing that may be noted from Table 1 is that for each starting point, the final design is the same for all three methods. This is expected since the compliance depends only on the deformation which is unique, hence independent on how rigid motions are treated. As for the average number of CG iterations and time, the pin-pointing method is the most efficient despite the two nodes with prescribed displacements ending up in the void part in the middle. CG on the singular system is only about half as fast. The penalty method required most time owing the relatively large number of CG iterations together with the additional matrix multiplications in (41). It should of course be noted that the Matlab implementation is certainly not optimal in terms of computational performance, so it is mainly the relative timings and the number of CG iterations that are of interest. Iteration histories for each method using all three starting points are shown in Fig. 4. For a given start point, the graphs coincide almost completely for the three methods, hence cannot be distinguished visually. As is quite common in TO, the most notable changes of the objective function value occur during the first 50–100 iterations or so, whereafter the change is much slower. 6.2. Cantilever The setup for the second example is shown in the top plot of Fig. 5, including an optimized design. The loads are chosen to induce rotational acceleration. The design domain Ω is a rectangle of size 2 times 1 and the FE mesh consists of 720 × 360 equal-sized, four-noded bi-linear elements. The filter radius R is set to 0.8% of the domain width. The lengths L are set to 0.2 and the magnitudes of the surface load to t1 = t = 0.1. The optimized design in the top plot of Fig. 5 resembles classical designs [6], but there is slight curvature of the structural members caused by the inertia relief load. Fig. 6 shows a design obtained by increasing the magnitude t1 by a factor 1.2. The resulting design differs notably from the one in Fig. 5. Table 2 shows a comparison of the performance when using one of the three methods for suppressing rigid body motions and, for each method, three different starting points chosen as for Table 1. It can be seen that all three methods perform similarly, for each starting point finding essentially the same design. For this problem, direct application of the CG method to the singular system is the fastest method. This notable improvement of performance compared to the MBB-example (Table 1) could be due to the much lower aspect ratio (2:1 compared to 4:1) of the Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

15

Fig. 5. Top: Setup and an optimized design. Bottom: Total load f α (ρ) consisting of prescribed loads f in black and computed inertia relief load −M(ρ)Zα(ρ) in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

design domain. The pin-pointing method again worked well despite supports ending up in a void area. The penalty method required most time owing to the additional matrix multiplications in (41) and a slightly larger number of CG iterations. Fig. 7 shows a 3D-version of the cantilever example. The design domain is a rectangular box with dimensions 2 × 1 × 1, meshed using 96 × 48 × 48 equal-sized eight-noded, tri-linear brick elements. Uniform loads are applied over square regions of Γ , with magnitudes t1 /1.8 = t = 0.1. (Exploiting symmetry around the plane y = 0.5 could be an option but was not implemented here.). The filter radius R is set to 0.08. The so-called boundary extension method [20] is implemented by prescribing the design to zero in a thin layer at the boundary of Ω except near and on the loaded parts; the effect of this is to avoid an unnatural sticking to the boundary which can be notable in 3D. A uniform initial guess is used with a maximum CPU-time of 3 h as the only stopping criteria for the optimization process. An optimized design is shown in Fig. 7 (top right). The prescribed loads have been chosen to induce rotational acceleration, as indicated by the curvature of the inertia relief force field shown in Fig. 7 (bottom). Performance data from the optimization runs is shown in Table 3. Again all three methods perform well, but CG on the singular system is the fastest (this explains the slightly lower objective function value, 0.99 versus 1.0 – there was simply a Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

16

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 6. Top: Optimized design with a larger magnitude for the bottom right load. Bottom: Total load.

Table 3 Data for the 3D example. Columns show (i) objective function value; (ii) average number of CG iterations; and (iii) average time per CG iteration. Method

obj.

CG iter.

CG time [s]

CG on sing. system Pin-pointing Penalty

0.99 1.0 1.0

19 29 32

10.72 12.93 15.20

few more optimization iterations taken within the time limit). The total number of optimization iterations is roughly between 500 and 600, with rapid changes during the first 50–100 iterations and thereafter slower; c.f. Fig. 4. As a final example we consider a version of the 3D cantilever with a slightly longer design domain having dimensions 3 × 1 × 1, meshed using 120 × 40 × 40 FEs of the same type as before. The magnitude of the loads are now set to t1 /3 = t = 0.1. An optimized design can be seen in Fig. 8. Similarly to the MBB-beam example (see e.g., Fig. 2) there is a large build-up of material around the region of application of one of the loads (the one with magnitude t1 ), creating a structure loosely resembling a hammer. This leads to a large upwards-acting load (right plot in Fig. 8) which counteracts the bending effect of the prescribed loads in order to minimize compliance. Performance data for the three different methods, again with a uniform initial guess and a maximum CPU-time of 3 h, is seen in Table 4. In Table 3 the methods performed very similarly, but here CG on the singular system performs markedly better, in particular in comparison with the penalty method. Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

17

Fig. 7. Top left: Design domain and loading for the 3D examples. Top right: An optimized design. Elements whose density is above 0.7 are shown. Bottom: Inertia relief load for the top right design.

Fig. 8. Left: An optimized design with longer design domain. Elements whose density is above 0.7 are shown. Right: Inertia relief load.

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

18

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx Table 4 Data for 3D example with longer design domain. Columns show (i) objective function value; (ii) average number of CG iterations; and (iii) average time per CG iteration. Method

obj.

CG iter.

CG time [s]

CG on sing. system Pin-pointing Penalty

0.98 0.99 1.0

27 34 61

10.57 15.92 22.91

7. Concluding remarks We have given a detailed theoretical account of TO with the inertia relief method for design of freely floating elastic continua undergoing rigid acceleration. The use of inertia relief leads to a well-posed static variational problem which require that rigid body motions can be suppressed without affecting the deformation. Three methods to accomplish this has been tested numerically. The numerical examples in Section 6 indicate that the commonly used pin-pointing method for suppressing rigid motion is applicable also for TO despite the possibility of supports ending up in low relative density regions. Direct solution of the singular equilibrium equation was slower on the MBB-example but faster for the cantilever-examples. A drawback compared to the pin-pointing method is the need to modify the CG-code, as discussed at the end of Section 5.1, to ensure a reasonable rate of convergence in all cases (pin-pointing requires modification of the stiffness matrix, but it can be argued that this is a more standard task). The penalty method was the slowest in all the numerical examples. These observations support the conclusion that all methods work, and the recommendation that pin-pointing be used as the first choice. As for the optimized designs, Figs. 2, 6 and 8 in particular indicate that, in cases where the relative importance of inertial loads compared to other loads is high, optimal designs for freely floating structure undergoing rigid acceleration differ significantly from “traditional” designs optimized for static conditions; cf. e.g. Figs. 1 and 2 or 3. As a final remark we mention that although the numerical examples are based on maximizing stiffness only, hence do not take into account important aspects of mechanical design such as stability, uncertainty, manufacturing constraints, stress levels and fatigue, we see no principle difficulty in combining inertia relief with such things in TO. Acknowledgments This research was supported in part by the Swedish Foundation for Strategic Research Grant No. AM13-0029, and in part by Bernt J¨armarks stiftelse f¨or vetenskaplig forskning, Sweden. Appendix A. A basis for R M Theorem 4.

The functions

e1 , e2 , e3 , x × e1 , x × e2 , x × e3 ,

(A.1)

constitute a basis for R M when d = 3. Proof. Let the functions in (A.1) be denoted z i = z i (x), i = 1, . . . , 6. To see that these functions span R M, note that for every v ∈ R M, there exists vectors a and b such that v(x) = a × x + b = (a1 e1 + a2 e2 + a3 e3 ) × x + (b1 e1 + b2 e2 + b3 e3 ) =

6 ∑

z i (x)vi ,

i=1

where (v1 , . . . , v6 )T = (b1 , b2 , b3 , −a1 , −a2 , −a3 )T . To see that the functions in (A.1) are linearly independent, we consider a function v ∈ R M such that v(x) =

6 ∑

z i (x)vi = 0

∀x ∈ Ω

⇔

i=1

v1 e1 + v2 e2 + v3 e3 + v4 (x × e1 ) + v5 (x × e2 ) + v6 (x × e3 ) = 0

∀x ∈ Ω .

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

19

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Taking derivatives with respect to x and y gives v4 = v5 = v6 = 0. It then follows immediately that v1 = v2 = v3 = 0, thus showing that the z i ’s are linearly independent. □ Appendix B. Existence and uniqueness of rigid accelerations Theorem 5. The system (19) has a unique solution for every ρ ∈ H such that ρ > 0 in some open set ω ⊂ Ω , and this solution satisfies |αi | ≤ c∥t∥ L 2 (Γ ) , i = 1, . . . , p. ˜ with components Mi j = (ρ M z i , z j ), and let c˜ = Proof. Introduce the p × p matrix M vectors c ∈ R p , ∫ p p ∑ ∑ ˜ c= cT M ci c j (ρ M z i , z j ) = (ρ M c˜ , c˜ ) = ρ M c˜ T c˜ ≥ 0,

∑p

i=1 ci z i .

Then for arbitrary

Ω

i=1 j=1

since ρ M ≥ 0 a.e. Now consider c such that ˜ c = (ρ M c˜ , c˜ ) = (ρ M c˜ , c˜ )ω . 0 = cT M Since ρ M > 0 in ω, this implies ∥˜c∥ L 2 (ω) = 0, hence c˜ = 0 a.e. and thus c = 0 since the z i ’s are linearly independent ˜ is positive definite, hence non-singular. Therefore on ω (cf. the proof of Theorem 4). This shows that the matrix M (19) has the unique solution αi =

p ∑ (M −1 )i j ℓs (z j ),

i = 1, . . . , p.

j=1

Recalling (10) we get |αi | ≤ p max |(M −1 )i j ||ℓs (z j )| ≤ c max ∥z j ∥ L 2 (Γ ) ∥t∥ L 2 (Γ ) j=1,..., p

j=1,..., p

using a trace theorem [56, Theorem 8.7] and Cauchy–Schwarz inequality. □ Appendix C. Existence of solutions to problem (20) To prove existence of solutions to the state problem (20) we make use of the following Korn’s inequality: Lemma 6.

There exists a constant c > 0 such that for all u ∈ R M ⊥ ,

∥ε(u)∥ L 2 (Ω) ≥ c∥u∥ H 1 (Ω) .

(C.1)

Proof (Cf. Section 3.3 in [11].). To get a contradiction, assume that there is no c > 0 for which (C.1) holds. Then there exists a sequence {v n } ⊂ R M ⊥ such that 1 and ∥v n ∥ H 1 (Ω) = 1 (C.2) ∥ε(v n )∥ L 2 (Ω) < n for all n. This sequence is bounded in H 1 (Ω ) and therefore admits a weakly convergent subsequence {v n k }. It follows from the Rellich–Kondrachov theorem [12, Theorem 9.16] and uniqueness of limits that v n k converges strongly in L 2 (Ω ). Korn’s (first) inequality [19, Theorem 6.15-1] gives c∥v n k − v m k ∥2H 1 (Ω) ≤ ∥ε(v n k ) − ε(v m k )∥2L 2 (Ω) + ∥v n k − v m k ∥2L 2 (Ω) ≤ 2∥ε(v n k )∥ L 2 (Ω) + 2∥ε(v m k )∥ L 2 (Ω) + ∥v n k − v m k ∥2L 2 (Ω) for some c > 0 (here we have used Young’s inequality in the last step). Invoking the first of (C.2) gives ( ) 1 2 2 2 2 ∥v n k − v m k ∥ H 1 (Ω) < + + ∥v n k − v m k ∥ L 2 (Ω) . c nk mk Since the right-hand side can be made arbitrarily small it follows that {v n k } is a Cauchy sequence in H 1 (Ω ), hence converges strongly in this space to some v ∈ H 1 (Ω ). It follows that ∥ε(v n k )∥ L 2 (Ω) → ∥ε(v)∥ L 2 (Ω) , Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

20

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

and then the first of (C.2) implies that ε(v) = 0 by uniqueness of limits. Since R M ⊥ is a closed subset, the limit v ∈ R M ⊥ . This implies v = 0, which contradicts the second of (C.2). It follows that there exists a c > 0 such that (C.1) holds. □ With this we can now prove Theorem 7.

Problem (20) has a unique solution for every ρ ∈ H.

Proof. We verify that a, ℓ and R M ⊥ satisfy the conditions of the Lax–Milgram lemma [19, Theorem 6.2-1], from which existence and uniqueness then follows. R M ⊥ being linear and closed (the latter follows from the continuity of the L 2 (Ω ) inner product) in the Hilbert space H 1 (Ω ) is itself a Hilbert space. To see that the linear form ℓ is bounded, hence continuous, we first note that |ℓs (v)| ≤ c∥t∥ L 2 (Γ ) ∥v∥ H 1 (Ω) from Cauchy–Schwarz inequality and a trace theorem [56, Theorem 8.7]. Then, using the H¨older inequality and the boundedness of αi (see Appendix B), |ℓ(v)| = |ℓs (v)| +

p ∑

|(ρ M z i , v)αi | ≤

i=1

c∥v∥ H 1 (Ω) +

p ∑

∥ρ M ∥ L ∞ (Ω) ∥z i ∥ H 1 (Ω) ∥v∥ H 1 (Ω) |αi | ≤ c∥v∥ H 1 (Ω) .

(C.3)

i=1

Symmetry and continuity of a is almost immediate. Coercivity then follows using (9), (2) and Lemma 6 to get a(ρ(ρ), ˜ u, u) ≥ cρβ∥u∥2H 1 (Ω) . □ Appendix D. Euler’s laws of motion The compatibility condition (23) is equivalent to linearized versions of Euler’s laws of motion (Euler I and II) holding for the entire design domain. To see this, first choose v in (23) as, in turn, e1 , e2 , e3 , e1 × x, e2 × x, and e3 × x. This gives Euler I: ( p ) ∫ ∫ ∑ t. (D.1) ρM z i αi = Ω

Γ

i=1

Then, using aT (b × c) = bT (c × a), we get Euler II: ( p ) ∫ ∫ ∑ x × t. ρM x × z i αi = Ω

(D.2)

Γ

i=1

As for the other direction, if v ∈ R M then, recalling (12), ∫ ∫ p p ∑ ∑ ℓ(v) = t Tv − (ρ M z i αi , v) = t T (a × x + b) − (ρ M z i αi , a × x + b). Γ

Γ

i=1

i=1

Continuing on we find ∫ ∫ p p ∑ ∑ T T t (a × x) + t b− (ρ M z i αi , a × x) − (ρ M z i αi , b) = Γ

a

Γ T

(∫

i=1

∫

ρM x ×

x×t− Γ

Ω

i=1

( p ∑ i=1

)) z i αi

T

(∫

+b

t− Γ

p ∫ ∑ i=1

) ρ M z i αi

= 0,

Ω

where the last equality follows since (D.1) and (D.2) are assumed to hold. Appendix E. Existence of a solution to the state problem with penalty Theorem 8. Problem (25) has a unique solution u ∈ H 1 (Ω ) for every ρ ∈ H, and since ℓ satisfies (23), u ∈ R M ⊥ and therefore solves (20). Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

21

Proof. We know from the proof of Theorem 7 that ℓ and R M ⊥ satisfy the conditions of the Lax–Milgram lemma [19, Theorem 6.2-1], so if we can show that a p is symmetric, continuous and coercive on H 1 (Ω ) we are done with the existence part. Symmetry of a p is immediate, and continuity follows from the continuity of a and the Cauchy–Schwarz inequality: |a p (u, v)| ≤ |a(u, v)| +

p ∑

|(u, z i )| |(v, z i )| ≤ c∥u∥ H 1 (Ω) ∥v∥ H 1 (Ω) .

i=1

For coercivity we recall the decomposition (13), by which u = ud + ur , with ud ∈ R M ⊥ and ur ∈ R M. Now write a p (u, u) = a p (ud , ud ) + 2a p (ud , ur ) + a p (ur , ur ) = p p ∑ ∑ (ur , z i )2 , a(ud , ud ) + (ur , z i )2 ≥ ρβc∥ud ∥2H 1 (Ω) +

(E.1)

i=1

i=1

where the inequality follows from Lemma 6. To handle the second term in the last of (E.1) we note that R M is finite-dimensional, so all norms on R M are equivalent [19, Theorem 2.7], in particular to the H 1 (Ω )-norm. Let   p ∑ f (u) = √ (u, z i )2 . (E.2) i=1

This function is non-negative and absolutely homogeneous. Using Minkowski’s inequality, we see that it satisfies the triangle inequality, i.e.    p  p ∑ ∑ √ 2 (u1 + u2 , z i ) = √ |(u1 , z i ) + (u2 , z i )|2 ≤ f (u1 + u2 ) = i=1

i=1

   p  p ∑ ∑ 2 √ |(u1 , z i )| + √ |(u2 , z i )|2 = f (u1 ) + f (u2 ). i=1

i=1

Finally, f (u) = 0 ⇔ (u, z i ) = 0, i = 1, . . . , p. The latter implies u ∈ R M ⊥ . But if also u ∈ R M, it follows that u = 0. Having thus verified that f defined in (E.2) is a norm on R M, it follows by equivalence of norms that f (ur )2 ≥ c∥ur ∥2H 1 (Ω) . Substitution in (E.1) yields a p (u, u) ≥ c1 ∥ud ∥2H 1 (Ω) + c2 ∥ur ∥2H 1 (Ω) ≥ c(∥ud ∥2H 1 (Ω) + ∥ur ∥2H 1 (Ω) ),

(E.3)

where c = min{c1 , c2 }. Now, ∥u∥2H 1 (Ω) = ∥ud + ur ∥2H 1 (Ω) ≤ (∥ud ∥ H 1 (Ω) + ∥ur ∥ H 1 (Ω) )2 ≤ 2(∥ud ∥2H 1 (Ω) + ∥ur ∥2H 1 (Ω) ), using Young’s inequality in the last step. Finally, substitution in (E.3) gives a p (u, u) ≥ c∥u∥2H 1 (Ω)

∀u ∈ H 1 (Ω ).

We now show that the solution u of (25) of belongs in R M ⊥ . Since the test function is arbitrary we may choose v = z j , j = 1, . . . , p. Then a p (u, z j ) = a(u, z j ) +

p ∑ (u, z i )(z j , z i ) = ℓ(z j ) = 0, i=1

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

22

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

where the last equality follows from (23). This gives p ∑ (u, z i )(z j , z i ) = 0,

j = 1, . . . , p.

i=1

Multiplying by (u, z j ) and summing over all j:s lead to ( p ) ∫ p p p p ∑ ∑ ∑ ∑ ∑ (u, z j ) (u, z i )(z j , z i ) = ( (u, z j )z j , (u, z i )z i ) = ∥ (u, z i )z i ∥2 = 0, j=1

i=1

j=1

i=1

Ω

i=1

∑p

which implies i=1 (u, z i )z i = 0 in Ω . Since the z i :s constitute a basis we conclude that (u, z i ) = 0, i = 1, . . . , p; i.e., u ∈ RM⊥ . □ Appendix F. Existence of solutions to problem (26) Existence of solutions for related problems has been shown with similar techniques in e.g. refs. [9,10,15,27,44, 45,55]. Anticipating their use for the FE convergence proof in Appendix G, some of the results below are stated in a little more generality than what is needed here. In this appendix and the next we will make use of the following relation between two types of convergence: Lemma 9. Let { f n } be a sequence of functions in L ∞ (Ω ) that are bounded a.e. and converges to zero a.e. in Ω . ∗ Then f n ⇀ 0 in L ∞ (Ω ). Proof. Let ϕ ∈ L 1 (Ω ) be arbitrary. Then since f n is bounded (by M), | f n ϕ| ≤ M|ϕ| whence ∫ |

f n ϕ| ≤

Ω

∫

a.e. in Ω , | f n ϕ| → 0

Ω

by Lebesgue’s dominated convergence theorem. The conclusion follows from the definition of weak∗ -convergence in L ∞ (Ω ). □ To prove Theorem 14 below we follow essentially the direct method of calculus of variations [21], beginning with a proof of compactness of the feasible set and then proceeding to show continuity (hence lower semi-continuity) in a suitable sense of the objective function over this set. Lemma 10.

The feasible set H defined in (1) is weakly* sequentially compact.

Proof. The pointwise bounds on ρ imply that H is bounded in L ∞ (Ω ) so it follows from the sequential Banach– ∗ Alaoglu theorem [36, Theorem 12, p. 107] that there exists a subsequence {ρn k } ⊂ H such that ρn k ⇀ ρ in L ∞ (Ω ). We now verify that the limit ρ ∈ H. That ρ ∈ [0, 1] a.e. follows from [9, Lemma 3.1]. To see that the volume constraint is satisfied we first note that the filtered design converges pointwise, i.e. ∫ [ ] S(ρn k ) − S(ρ) = ρn k ( y) − ρ( y) Ψ (x, y)d y → 0 ∀x ∈ Ω , (F.1) Ω

since Ψ (x, ·) is continuous with compact support, hence belong in L 1 (Ω ), whence convergence follows by definition of the weak∗ convergence. Since ρ ≥ 0 a.e., S(ρ) ≥ 0 in Ω , so Fatou’s lemma gives ∫ ∫ S(ρ) ≤ lim inf S(ρn k ) ≤ γ |Ω |, n k →∞

Ω

Ω

and thus ρ ∈ H. □ ∗

Let ρn ⇀ ρ in L ∞ (Ω ) and let, for each n = 0, 1, . . ., ρ Mn and ρ˜n be functions from L ∞ (Ω ) to L ∞ (Ω ) that take on positive values a.e. in Ω , are bounded in L ∞ (Ω ), and satisfy ∗

ρ Mn (ρn ) ⇀ ρ M (ρ) in L ∞ (Ω )

(F.2)

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

23

and ρ˜n (ρn ) → ρ(ρ) ˜ a.e. in Ω ,

(F.3)

respectively. Since the penalization function g in (2) is continuous it follows immediately from the pointwise convergence (F.1) that ρ(ρ ˜ n ) → ρ(ρ) ˜ in Ω . Recalling (6) it also follows that ρ M (ρn ) → ρ M (ρ) in Ω . Since |ρ M (ρn ) − ρ M (ρ)| ≤ ρsolid this implies, by Lemma 9, ∗

ρ M (ρn ) ⇀ ρ M (ρ) in L ∞ (Ω ). For the purposes of this appendix we could thus take ρ Mn = ρ M and ρ˜n = ρ˜ for all n, but in the FE convergence proof in Appendix G the functions themselves (and not only the design) will actually vary with n (corresponding to the mesh parameter h). We now continue to show convergence of the rigid accelerations: Lemma 11.

∗

Let ρn ⇀ ρ in L ∞ (Ω ). Consider a sequence {α n } of solutions to

p ∑ (ρ Mn (ρn )z i , z j )αin = ℓs (z j ),

j = 1, . . . , p.

(F.4)

i=1

This sequence converges strongly to the solution α = α(ρ) of (19). Proof. For each ρ Mn (ρn ), Theorem 5 in Appendix B ensures the existence of a unique α n satisfying (F.4). The sequence {α n } ⊂ R p is thus bounded and therefore admits a subsequence converging weakly, hence strongly since p < ∞, to some α. Now consider an arbitrary convergent subsequence, again denoted {α¨ n }, and use (F.4) to write, for j = 1, . . . , p, p p p ∑ ∑ ∑ (ρ M (ρ)z i , z j )αi − ℓs (z j ) = (ρ M (ρ)z i , z j )αi − (ρ Mn (ρn )z i , z j )αin = i=1

i=1

i=1

p p ∑ ∑ ([ρ M (ρ) − ρ Mn (ρn )]z i , z j )αi + (ρ Mn (ρ)z i , z j )(αi − αin ). i=1

(F.5)

i=1 ∗

The first term on the second line side tends to zero since ρ Mn (ρn ) ⇀ ρ M (ρ) in L ∞ (Ω ) and z iT z j ∈ L 1 (Ω ). Using Cauchy–Schwarz, |(ρ Mn (ρ)z i , z j )| ≤ ∥ρ Mn (ρ)∥ L ∞ (Ω) ∥z i ∥ L 2 (Ω) ∥z j ∥ L 2 (Ω) , so (ρ Mn (ρ)z i , z j ) is bounded for all n. Therefore αin → αi , i = 1, . . . , p, implies that the second term in the last line of (F.5) tends to zero. This shows that the limit α solves (19). Since, by Theorem 5, this equation has a unique solution and the chosen subsequence was arbitrary it follows that the entire sequence converges to α = α(ρ). □ Next we establish continuity of the compliance in a certain sense. Lemma 12.

∗

Let ρ Mn (ρn ) ⇀ ρ M (ρ) in L ∞ (Ω ) and v n ⇀ v in H 1 (Ω ). Then

ℓ(ρ Mn (ρn ), v n ) → ℓ(ρ M (ρ), v).

(F.6)

Proof. The left-hand side in (F.6) is given by p ∑ ℓ(ρ Mn (ρn ), v n ) = ℓs (v n ) − (ρ Mn (ρn )z i , v n )αin ,

(F.7)

i=1

where αin solves (F.4). Here the first term on the right-hand side tends to ℓs (v) since ℓs is a continuous linear functional over H 1 (Ω ). As for the second term we write, for i = 1, . . . , p, (ρ Mn (ρn )z i , v n ) = (ρ Mn (ρn )z i , v n − v) + (ρ Mn (ρn )z i , v).

(F.8)

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

24

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

It follows from the Rellich–Kondrachov theorem [12, Theorem 9.16] and uniqueness of limits that v n → v in L 2 (Ω ), so that by Cauchy–Schwarz, |(ρ Mn (ρh )z i , v n − v)| ≤ ∥ρ Mn (ρn )z i ∥ L 2 (Ω) ∥v n − v∥ L 2 (Ω) → 0. Combining (F.8), the weak∗ convergence (F.2), and Lemma 11 we get (ρ Mn (ρn )z i , v n )αin → (ρ M (ρ)z i , v)αi , Now (F.6) follows from (F.7).

i = 1, . . . , p.

□

Now introduce a sequence of spaces {V n }, V n ⊂ R M ⊥ , that are dense in R M ⊥ in the sense that for every v ∈ R M ⊥, lim inf ∥v − v n ∥ H 1 (Ω) = 0.

n→∞ v n ∈V n

(F.9)

In other words, for every v ∈ R M ⊥ there exists a sequence {v n }, v n ∈ V n , which converges strongly to v in H 1 (Ω ). To prove Theorem 14 we only use V n = R M ⊥ for all n, but for the FE convergence in Appendix G the space will of course vary depending on the mesh parameter. Lemma 13.

∗

Let ρn ⇀ ρ in L ∞ (Ω ). Then the sequence of displacements {un }, un ∈ V n , solving

a(ρ˜n (ρn ), un , v n ) = ℓ(ρ Mn (ρn ), v n ) ∀v n ∈ V n

(F.10)

converges weakly in H 1 (Ω ) to the solution u = u(ρ) of (20). Proof. Since V n ⊂ R M ⊥ and a(ρ(ρ), ˜ ·, ·) is coercive on R M ⊥ it follows that a(ρ˜n (ρn ), ·, ·) is coercive on V n . This gives (recalling (11) and Lemma 6) βρc∥un ∥2H 1 (Ω) ≤ a(ρ˜n (ρn ), un , un ) ≤ |ℓ(ρ Mn (ρn ), un )| ≤ c2 ∥un ∥ H 1 (Ω) , where the last inequality follows from (C.3), so that ∥un ∥ H 1 (Ω) ≤ c2 /(βρc). The sequence {un } is thus bounded and therefore admits a weakly convergent subsequence, again denoted {un }. Let v be an arbitrary function in R M ⊥ . Since, by assumption, V n is dense in R M ⊥ there exists a sequence {v n } converging strongly in H 1 (Ω ) to v. Now write the left-hand side in (F.10) as a(ρ˜n (ρn ), un , v n − v) + a(ρ˜n (ρn ) − ρ(ρ), ˜ un , v) + a(ρ(ρ), ˜ un , v).

(F.11)

Due to the continuity of a we have |a(ρ˜n (ρn ), un , v n − v)| ≤ c∥ρ˜n (ρn )∥ L ∞ (Ω) ∥un ∥ H 1 (Ω) ∥v n − v∥ H 1 (Ω) , implying that a(ρ˜n , un , v n − v) → 0 due to the strong convergence of v n to v in H 1 (Ω ) and the boundedness of un and ρ˜n (ρn ) in H 1 (Ω ) and L ∞ (Ω ), respectively. The symmetry of a together with Cauchy–Schwarz inequality give ∫ |a(ρ˜n (ρn ) − ρ(ρ), ˜ un , v)| ≤ |(ρ˜n (ρn ) − ρ(ρ))σ ˜ (v) : ε(un )| ≤ Ω d ∑ d ∑

∥(ρ˜n (ρn ) − ρ(ρ))σ ˜ i j (v)∥ L 2 (Ω) ∥εi j (un )∥ L 2 (Ω) .

i=1 j=1

Since σi j (v) ∈ L 2 (Ω ) and thus are finite almost everywhere it follows from the assumed pointwise convergence (F.3) that (ρ˜n (ρn ) − ρ(ρ))σ ˜ i j (v) → 0 a.e. in Ω . In addition, |(ρ˜n (ρn ) − ρ(ρ))σ ˜ i j (v)| ≤ |σi j (v)| a.e. in Ω , Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

25

so by Lebesgue’s dominated convergence theorem and the boundedness of ∥εi j (un )∥ L 2 (Ω) we get |a(ρ˜n (ρn ) − ρ(ρ), ˜ un , v)| → 0. Finally the last term in (F.11) tends to a(ρ(ρ), ˜ u, v) due to the weak convergence of un to u. The proof is concluded by applying Lemma 12 to the right-hand side of (F.10) and noting that the entire sequence {un } converges to u since the state problem (20) has a unique solution. □ We are finally ready to prove the main theorem: Theorem 14.

Problem (26) has at least one globally optimal solution.

Proof. It follows from (11) that ℓ(ρ M (ρ), u, u) = a(ρ(ρ), ˜ u, u) ≥ 0

∀(ρ, u) ∈ H × H 1 (Ω );

i.e. the compliance is bounded below. Now let {ρn } ⊂ H be a minimizing sequence; i.e., lim ℓ(ρ M (ρn ), u(ρn )) = inf ℓ(ρ M (ρ), u(ρ)). ρ∈H

n→∞

By Lemma 10, this sequence has a subsequence {ρn k } converging weakly∗ in L ∞ (Ω ) to some ρ∗ ∈ H. According to Lemma 13 the corresponding sequence of displacements {u(ρn k )} converges weakly in H 1 (Ω ) to some u = u(ρ∗ ). Using Lemma 12 we therefore obtain ℓ(ρ M (ρn k ), u(ρn k )) → ℓ(ρ M (ρ∗ ), u(ρ∗ )). Since {ρn k } is a minimizing sequence it follows, by uniqueness of limits, that ℓ(ρ M (ρ∗ ), u(ρ∗ )) = inf ℓ(ρ M (ρ), u(ρ)), ρ∈H

i.e. ρ∗ is a globally optimal solution to (26). □ Appendix G. FE convergence We show here that any sequence of (globally) optimal solutions to problem (28) admits a convergent subsequence, the limit of which is a globally optimal solution to (26). Similar proofs can be found in e.g. refs. [9,10,28,44,45]. The discretized design and the discretized filtered design are both element-wise constant, thus raising some concerns regarding what happens at the boundaries between elements. However, the boundary points of the mesh constitute a set of measure zero and therefore what happens there is of no consequence here. More precisely: Let ∪m e=1 ∂Ωe denote the set of points in Ω belonging to an element side. Since there are, for every h ≥ 0, countably many elements, each with a boundary ∂Ωe of measure zero, it follows that ∪m e=1 ∂Ωe also has measure zero, hence that almost every point in Ω belongs in ∪m Ω , the set of points that are in the interior of some element. In the e e=1 following, pointwise convergence a.e. suffices, hence we can, for the proofs that follow, confine ourselves to points in ∪m e=1 Ωe . We start by showing the pointwise a.e. convergence of the discretized filtered design. Lemma 15.

∗

With Sh defined in (27), let ρh ⇀ ρ, ρh ∈ Hh , in L ∞ (Ω ). Then

Sh (ρh ) → S(ρ) a.e. in Ω .

(G.1)

(Here and in the following, convergence refers to sequences obtained as the mesh parameter h → 0 from above.) Proof. Let x be an arbitrary point in ∪m e=1 Ωe . At this point we have |Sh (ρh ) − S(ρ)| ≤ |S(ρh ) − S(ρ)| + |Sh (ρh ) − S(ρh )|.

(G.2)

The first term to the right tends to zero according to (F.1). To see that the second term tends to zero we first note that, since the elements are disjoint, x ∈ Ωe for some e. Let ρh ∈ [0, 1] a.e. be arbitrary. Then, making use of the Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

26

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

integral mean value theorem in the last line, we get ∫ ⏐ ∑ ⏐ ⏐ ⏐ ρ f |Ω f |Ψ (x e , x f ) − ρh ( y)Ψ (x, y)d y⏐ = ⏐ Ω

f ∈Ne

m ∫ ⏐ ∑ ∑ ⏐ ρ f |Ω f |Ψ (x e , x f ) − ⏐ f ∈Ne

⏐ ∑ ⏐ ρ f |Ω f |Ψ (x e , x f ) − ⏐

f =1 Ω f ∫ m ∑

ρf

f =1

f ∈Ne

∑

⏐ ⏐ ρh ( y)Ψ (x, y)d y⏐ =

Ωf

⏐ ⏐ Ψ (x, y)d y⏐ ≤ ∑

ρ f |Ω f ∥Ψ (x e , x f ) − Ψ (x, ξ f )| +

f ∈Ne

ρ f |Ω f ∥Ψ (x, ξ f )|,

(G.3)

f ∈{1,...,m}\Ne

for some ξ f ∈ Ω f . The first term in the last line comes from those elements whose centroid lies within the support of the filter, and the other from those whose centroid lies outside. Since the elements are convex, x → x e and ξ f → x f as h tends to zero. Given ε > 0, for small enough h the continuity of Ψ gives |Ψ (x, ξ f ) − Ψ (x e , x f )| < ε/2|Ne | |Ω f | for f ∈ Ne and |Ψ (x, ξ f )| < ε/2(m − |Ne |)|Ω f | for f ∈ / Ne . Here |Ne | denotes the number of elements in Ne . Substitution in (G.3) together with ρ f ≤ 1 leads to ∫ ⏐ ⏐ ∑ ⏐ ⏐ ρ f |Ω f |Ψ (x e , x f ) − ρh ( y)Ψ (x, y)d y⏐ < ⏐ Ω

f ∈Ne

∑ f ∈Ne

ε|Ω f | + 2|Ne | |Ω f |

∑ f ∈{1,...,m}\Ne

ε|Ω f | < ε. 2(m − |Ne |)|Ω f |

From this and the fact that every x ∈ ∪m e=1 Ωe belongs to some Ωe we conclude that ∫ ∑ Ω ρh ( y)Ψ (x, y)d y f ∈Ne ρ f |Ω f |Ψ (x e , x f ) Sh (ρh )(x) = ∑ → S = S(ρh )(x) a.e. in Ω , 1 f ∈Ne |Ω f |Ψ (x e , x f ) using the fact that ∫ Ψ (x, y)d y = 1

∀x ∈ Ω

Ω

with the kernel chosen as in (5). Now (G.1) follows from (G.2).

□

From the pointwise convergence of the filtered design the following result is almost immediate: Lemma 16.

∗

Let ρh ⇀ ρ in L ∞ (Ω ), ρh ∈ Hh . Then ∗

ρ˜h (ρh ) → ρ(ρ) ˜ a.e. in Ω and ρ Mh (ρh ) ⇀ ρ M (ρ) in L ∞ (Ω ), where ρ˜h and ρ Mh are defined in (30) and (31), respectively. Proof. Since the penalization function g in (2) is continuous it follows immediately from the pointwise convergence (G.1) that ρ˜h (ρh ) → ρ(ρ) ˜ a.e. in Ω . Recalling (6) it follows that ρ Mh (ρh ) → ρ M (ρ) a.e. in Ω . Since |ρ Mh (ρh ) − ρ M (ρ)| ≤ ρsolid this implies, by ∗ Lemma 9, that ρ Mh (ρh ) ⇀ ρ M (ρ) in L ∞ (Ω ). □ In order to prove the FE convergence in Theorem 20 we need density results for the design and the displacements; i.e., that functions in H and R M ⊥ can be approximated to arbitrary precision by functions from the respective FE spaces. For the design we have the following result: Lemma 17.

∗

For every ρ ∈ H there exists a sequence {ρh }, ρh ∈ Hh , such that ρh ⇀ ρ in L ∞ (Ω ).

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

27

Proof. Given ρ ∈ H, introduce the interpolation operator Πh defined by ∫ 1 Πh ρ|Ωe = ρ, e = 1, . . . , m. |Ωe | Ωe We start by showing that Πh ρ → ρ a.e. in Ω

(G.4)

and then modify Πh to obtain an interpolation operator mapping ρ into Hh for every h. For each x ∈ ∪m e=1 Ωe there holds that ∫ 1 Πh ρ(x) = ρ, |Ωe | Ωe for some e since the elements do not overlap. Let B(x, h) be the open (Euclidean) ball of radius h centered at x. Then ∥x − y∥ ≤ sup x, y∈Ωe ∥x − y∥ = h, so that y ∈ B(x, h), and thus Ωe ⊂ B(x, h). Since the meshes are quasi-uniform, there holds that h ≤ cre , where re is the radius of the largest ball inscribed in Ωe . It follows that ( ) 1 3 max(0,d−2) 1 d d |Ωe | ≥ πre ≥ d π h = d |B(x, h)|, c c 4 so now a variant of Lebesgue’s differentiation theorem [29, p. 460] and the fact that almost every point in ∪m e=1 Ωe is a Lebesgue point [29, p. 458] for Πh ρ gives (G.4). Straightforward calculations show that Πh ρ ∈ [0, 1] a.e. in Ω . It therefore only remains to ensure satisfaction of the volume constraint in Hh . To this end we multiply Πh by { γ |Ω | } . (G.5) ch = min 1, ∫ Ω Sh (Πh ρ) Clearly ch ∈ [0, 1], so ch Πh ρ ∈ [0, 1] a.e. in Ω . We also have ∫ Sh (ch Πh ρ) ≤ γ |Ω |.

(G.6)

Ω

The modified interpolation operator thus satisfies ch Πh ρ ∈ [0, 1] and the volume constraint in Hh . It remains to ∗ show that ch Πh ρ ⇀ ρ in L ∞ (Ω ). ∗ The pointwise convergence of Πh ρ to ρ together with |Πh ρ − ρ| ≤ 1 a.e. imply (Lemma 9) that Πh ρ ⇀ ρ in L ∞∫(Ω ), hence by Lemma 15 that Sh (Πh ρ) → S(ρ) a.e. Therefore Lebesgue’s dominated convergence theorem ∫ gives Ω Sh (Πh ρ) → Ω S(ρ). This in turn implies that { { γ |Ω | } γ |Ω | } ch = min 1, ∫ → min 1, ∫ = 1, (G.7) Ω Sh (Πh ρ) Ω S(ρ) ∫ where the last equality follows from Ω S(ρ) ≤ γ |Ω |. Using (G.4) and (G.7) we therefore get ch Πh ρ → 1 · ρ = ρ a.e. in Ω . ∗

Applying Lemma 9 gives ch Πh ρ ⇀ ρ in L ∞ (Ω ); i.e, {ch Πh ρ} is the sequence whose existence was asserted in the statement of the lemma. □ For the displacement we now prove density in the sense of (F.9) for the sequence of FE spaces {V h ∩ R M ⊥ } in R M ⊥ . To this end we first show that every function in R M ⊥ can be approximated to arbitrary precision by a function in R M ⊥ ∩ C ∞ (Ω ), then use this result to prove the density of {V h ∩ R M ⊥ } in R M ⊥ in Lemma 19. Lemma 18. For every v ∈ R M ⊥ and ε > 0 there exists a function v ε ∈ R M ⊥ ∩ C ∞ (Ω ) such that ∥v − v ε ∥ H 1 (Ω) < ε. Proof. The space C ∞ (Ω ) is dense in H 1 (Ω ) [1, Theorem 3.18], so there exists for every ε > 0, a function v˜ ε ∈ C ∞ (Ω ) such that ∥v − v˜ ε ∥ H 1 (Ω) < ε.

(G.8)

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

28

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

However, it is not necessarily the case that v˜ ε ∈ R M ⊥ . To obtain such a function we note that (13) gives the unique decomposition v˜ ε = v εd + v εr , with v εd ∈ R M ⊥ and v εr ∈ R M. Since v εr can be written as v εr = Ax + b, see (12), it is clear that v εr ∈ C ∞ (Ω ). Then the function v ε = v˜ ε − v εr belongs in R M ⊥ ∩ C ∞ (Ω ), and we now proceed to show that this function is the one whose existence is asserted by the lemma. It is straightforward to show (cf. Appendix D) the alternative characterization ∫ ∫ ⊥ 1 R M = {u ∈ H (Ω ) | u = 0, x × u = 0}. (G.9) Ω

Ω

Since v ε ∈ R M , (G.9) gives ∫ ∫ ∫ vε = [˜v ε − v εr ] = 0 and x × vε ⊥

Ω

Ω

∫

x × [˜v ε − v εr ] = 0.

=

Ω

(G.10)

Ω

p

Given a basis {z i }i=1 for R M we can write p ∑

v εr =

ai z i .

(G.11)

i=1

Substitution into ∫ ∫ (G.10) shows that the ai ’s can be expressed (uniquely) as linear combinations of the components of Ω v˜ ε and Ω x × v˜ ε ; i.e., ai =

d ∑

∫

v˜εj +

ai j Ω

j=1

d ∑

∫

[x × v˜ ε ] j

bi j

(G.12)

Ω

j=1

for real numbers ai j and bi j . Here we define [x × v˜ ε ]2 = 0 if d = 2. Using the triangle inequality and (G.11) we get ∥v − v ε ∥ H 1 (Ω) = ∥v − v˜ ε + v εr ∥ H 1 (Ω) ≤ ∥v − v˜ ε ∥ H 1 (Ω) + ∥v εr ∥ H 1 (Ω) ≤ ∥v − v˜ ε ∥ H 1 (Ω) +

p ∑

(G.13)

|ai | ∥z i ∥ H 1 (Ω)

i=1

We now try to bound |ai | from above in terms of ∥v − v˜ ε ∥ H 1 (Ω) . To this end we use (G.12) and the inequality (a + b)2 ≤ 2a 2 + 2b2 to write ∫ ∫ ∫ ∫ d d d d ∑ ∑ ∑ ∑ ai j bi j ai j bi j |ai |2 = | v˜εj + [x × v˜ ε ] j |2 ≤ 2| v˜εj |2 + 2| [x × v˜ ε ] j |2 . (G.14) Ω

j=1

j=1

Ω

j=1

Ω

j=1

Ω

Now consider |

d ∑

⎛ ⎞2 ∫ d ∑ v˜εj |2 ≤ ⎝ |ai j | | v˜εj |⎠ .

∫ ai j Ω

j=1

j=1

Ω

∫ Since v ∈ R M ⊥ , and thus Ω v j = 0, we can write ∫ ∫ | v˜εj | = | [v j − v˜εj ]| ≤ |Ω |1/2 ∥v j − v˜εj ∥ L 2 (Ω) ≤ |Ω |1/2 ∥v j − v˜εj ∥ H 1 (Ω) , Ω

Ω

using Cauchy–Schwarz for the first inequality. Then ∫ d d ∑ ∑ |ai j | | v˜εj | ≤ |ai j | |Ω |1/2 ∥v j − v˜εj ∥ H 1 (Ω) ≤ Ω

j=1 1/2

|Ω |

max |ai j |

j=1,...,d

j=1 d ∑ j=1

∥v j − v˜εj ∥ H 1 (Ω) = ci

d ∑

∥v j − v˜εj ∥ H 1 (Ω) ,

j=1

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

∑d

29

∑d

which, using the inequality ( i=1 xi )2 ≤ 4 i=1 xi2 valid for d = 2 or 3, gives ⎛ ⎞2 ⎛ ⎞2 ∫ ∫ d d d ∑ ∑ ∑ | ai j v˜εj |2 ≤ ⎝ |ai j | | v˜εj |⎠ ≤ ci2 ⎝ ∥v j − v˜εj ∥ H 1 (Ω) ⎠ ≤ Ω

j=1

Ω

j=1

4ci2

d ∑

j=1

∥v j − v˜εj ∥2H 1 (Ω) = Ci ∥v − v˜ ε ∥2H 1 (Ω) .

j=1

Substitution in (G.14) gives d ∑

|ai |2 ≤ 2|

∫

v˜εj |2 + 2|

ai j

d ∑

Ω

j=1

∫ bi j Ω

j=1

[x × v˜ ε ] j |2 ≤ ci ∥v − v˜ ε ∥2H 1 (Ω) + 2|

d ∑

∫ bi j

j=1

[x × v˜ ε ] j |2 .

(G.15)

Ω

To bound the second term in (G.15) we note first that ⎞2 ⎛ ∫ ∫ d d ∑ ∑ bi j [x × v˜ ε ] j |2 ≤ ⎝ |bi j | | [x × v˜ ε ] j |⎠ . | Ω

j=1

Since v ∈ R M ⊥ , it holds that d ∑

∫

∫

Ω [x

d ∑

[x × v˜ ε ] j | =

|bi j ∥ Ω

j=1

Ω

j=1

× v] j = 0, so that ∫

[x × (˜v ε − v)] j |.

|bi j ∥ Ω

j=1

Now writing ∫

[x × (˜v ε − v)] j | = |

|

∫ ∑ d

Ω

X jk (v˜εk − vk )|,

Ω k=1

where the matrix-valued function X is such that Xv = x × v, we get ∫

[x × (˜v ε − v)] j | ≤

| Ω

d ∫ ∑ k=1

max ∥X jk ∥ L 2 (Ω)

k=1,...,d

d ∑

|X jk (v˜εk − vk )| ≤

Ω

d ∑

∥X jk ∥ L 2 (Ω) ∥v˜εk − vk ∥ H 1 (Ω) ≤

k=1

∥v˜εk − vk ∥ H 1 (Ω) = c j

k=1

d ∑

∥v˜εk − vk ∥ H 1 (Ω) ,

k=1

so that d ∑

∫

[x × v˜ ε ] j | ≤

|bi j | | Ω

j=1

d ∑

( |bi j | c j

j=1

d ∑

) ∥v˜εk − vk ∥ H 1 (Ω)

≤ ci

k=1

d ∑

∥v˜εk − vk ∥ H 1 (Ω) .

k=1

This gives |

d ∑ j=1

∫

⎞2 ⎛ ∫ d d ∑ ∑ 2 ⎠ ⎝ |bi j | | [x × v˜ ε ] j | ≤ 4ci ∥v˜εk − vk ∥2H 1 (Ω) = ci ∥v − v˜ ε ∥2H 1 (Ω) , [x × v˜ ε ] j | ≤ 2

bi j Ω

Ω

j=1

k=1

which when inserted in (G.15) yields |ai |2 ≤ 2|

d ∑

∫ ai j

j=1

ci ∥v −

vεj |2 + 2|

Ω

v˜ ε ∥2H 1 (Ω)

d ∑ j=1

∫

[x × v˜ ε ] j |2 ≤

bi j Ω

+ 2di ∥v − v˜ ε ∥2H 1 (Ω) = ci ∥v − v˜ ε ∥2H 1 (Ω) .

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

30

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Substituting this result in (G.13) gives ∥v − v ε ∥ H 1 (Ω) ≤ ∥v − v˜ ε ∥ H 1 (Ω) +

p ∑

|ai | ∥z i ∥ H 1 (Ω) ≤

i=1 p

∥v − v˜ ε ∥ H 1 (Ω) +

∑

1/2

ci ∥v − v˜ ε ∥ H 1 (Ω) ∥z i ∥ H 1 (Ω) ≤ (1 + c)∥v − v˜ ε ∥ H 1 (Ω) .

i=1

Letting ε > 0 be arbitrary and choosing v˜ ε such that ∥v − v˜ ε ∥ H 1 (Ω) < ε/(1 + c) we find for v ε = v˜ ε − v εr that ∥v − v ε ∥ H 1 (Ω) < ε for all ε > 0. □ Given Lemma 18 we can now prove Lemma 19. lim

For every v ∈ R M ⊥ , there holds that inf

h→0+ v h ∈V h ∩R M ⊥

∥v − v h ∥ H 1 (Ω) = 0.

(G.16)

Proof. Let π h v ∈ V h denote the L 2 -projection of v ∈ R M ⊥ onto V h , i.e. (π h v, v h ) = (v, v h ) ∀v h ∈ V h . Since R M ⊂ V h it follows that for every v ∈ R M ⊥ , (π h v, v h ) = (v, v h ) = 0

∀v h ∈ R M,

meaning that π h v ∈ R M ⊥ . Since the FE meshes are quasi-uniform, and since the components of the L 2 -projection are the L 2 -projections of the components, we get the estimate [23, Prop. 1.134] ∥v − π h v∥ H 1 (Ω) ≤ ch|v| H 2 (Ω)

∀v ∈ C ∞ (Ω ) ∩ H 1 (Ω ).

By Lemma 18 we can for every ε > 0 find a function v ε ∈ R M ⊥ ∩ C ∞ (Ω ) such that ∥v − v ε ∥ H 1 (Ω) < ε/2. Then π h v ε ∈ V h ∩ R M ⊥ and we get inf

v h ∈V h ∩R M ⊥

∥v − v h ∥ H 1 (Ω) ≤ ∥v − π h v ε ∥ H 1 (Ω) ≤

∥v − v ε ∥ H 1 (Ω) + ∥v ε − π h v ε ∥ H 1 (Ω) < ε/2 + ch|v ε | H 2 (Ω) < ε for every h < ε/(2c|v ε | H 2 (Ω) ). Now (G.16) follows since ε > 0 is arbitrary. □ The proof of Lemma 19 is quite similar to e.g. that of Theorem 3.2.3 in [18], but here we used the L 2 -projection instead of the Lagrange interpolant to ensure π h v ∈ R M ⊥ . We are now ready to prove the main result. Theorem 20.

Any sequence {ρh∗ } of solutions to (28) admits a convergent subsequence whose limit solves (26).

Proof. Consider an arbitrary sequence {ρh }, ρh ∈ Hh . Since Hh is bounded in L ∞ (Ω ) for every h > 0, the sequential Banach–Alaoglu theorem [36, Theorem 12, p. 107] asserts the existence of a subsequence, again denoted {ρh }, converging weakly∗ to some ρ ∈ L ∞ (Ω ). That ρ ∈ H then follows as in Lemma 10, using the pointwise convergence of Sh (ρh ) to S(ρ) in (G.1). By combining Lemmas 16 and 19 with Lemma 13 one infers that the corresponding sequence of displacements {uh (ρh )} converges weakly to the solution u(ρ) of (20). The preceding paragraph shows that any sequence {ρh } feasible in (28) admits a convergent subsequence whose limit is feasible in (26). This holds in particular for the limit ρ ∗ of any such sequence {ρh∗ } of solutions to the FE discretized problem (28). It therefore only remains to show that such a ρ ∗ is also optimal in (26). Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

31

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx ∗

Let ρ ∈ H be arbitrary. By Lemma 17 there exists a sequence {ρh }, ρh ∈ Hh , such that ρh ⇀ ρ in L ∞ (Ω ). Since ρh∗ solves (28) we have ℓ(ρ Mh (ρh∗ ), uh (ρh∗ )) ≤ ℓ(ρ Mh (ρh ), uh (ρh ))

∀ρ ∈ H. ∗

∗

Now, since ρ Mh (ρh∗ ) ⇀ ρ M (ρ ∗ ) in L ∞ (Ω ), uh (ρh∗ ) ⇀ u(ρ ∗ ) in H 1 (Ω ), ρ Mh (ρh ) ⇀ ρ M (ρ) in L ∞ (Ω ), and uh (ρh ) ⇀ u(ρ) in H 1 (Ω ), Lemma 12 implies that ℓ(ρ M (ρ ∗ ), u(ρ ∗ )) ≤ ℓ(ρ M (ρ), u(ρ)) ∀ρ ∈ H. In other words, ρ ∗ solves problem (26). □ It may be noted that the proof of Theorem 20 of course constitutes an alternative existence proof for problem (26), under somewhat more restrictive assumptions than those of Theorem 14. If desired, one can follow for instance Petersson [44] and also prove strong convergence in H 1 (Ω ) for the displacements. Strong convergence in L p (Ω ), p ∈ [1, ∞), of the filtered design can also be shown. Appendix H. The null space of K p

Theorem 21. The vectors {ˆz i }i=1 constitute a basis for N (K ). Proof. By definition of a basis we need to show that the zˆ i :s span N (K ) and are linearly independent. We begin with the spanning property, i.e. we show that vˆ ∈ N (K ) ⇒ ∃κ ∈ R p : vˆ =

p ∑

zˆ i κi .

i=1

It is straightforward to show that the stiffness matrix defined in (35) may equivalently be characterized by ˆ a(v h , w h ) = vˆ T K w

∀v h , w h ∈ V h .

Since vˆ ∈ N (K ) we therefore get ˆ = a(v h , w h ) ∀w h ∈ V h . 0 = vˆ T K w Since w h is arbitrary it follows that v h ∈ R M, i.e. ∃κ ∈ R p : v h =

p ∑

z i κi .

i=1

∑ Substituting the representation z i = nj=1 ϕ j [ˆz i ] j gives ⎞ ⎛ ( p ) p n n ∑ ∑ ∑ ∑ ⎠ ⎝ vh = ϕ j [ˆz i ] j κi = ϕj [ˆz i ] j κi . i=1

Since

{ϕ j }nj=1 vˆ j =

j=1

j=1

i=1

is a basis for V h it follows by comparison with v h =

p ∑

[ˆz i ] j κi , j = 1, . . . , n ⇔ vˆ =

i=1

p ∑

∑n

j=1

ϕ j vˆ j that

zˆ i κi ,

i=1 p

thus showing that {ˆz i }i=1 spans N (K ). ∑p p To show that {ˆz i }i=1 are linearly independent we assume that i=1 κi zˆ i = 0. Then ⎛ ⎞ ( p ) p p n n ∑ ∑ ∑ ∑ ∑ κi z i = κi ⎝ ϕ j [ˆz i ] j ⎠ = ϕj κi [ˆz i ] j = 0. i=1

i=1

j=1

j=1

i=1

p

But now, since {z i }i=1 is a basis for R M, this implies κi = 0, i = 1, . . . , p; i.e. the zˆ i ’s are linearly independent. □ Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

32

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

References [1] R. Adams, Sobolev Spaces, Academic Press, 1975. [2] O. Amir, N. Aage, B. Lazarov, On multigrid-CG for efficient topology optimization, Struct. Multidiscip. Optim. 49 (5) (2014) 815–829. [3] O. Amir, M. Stolpe, O. Sigmund, Efficient use of iterative solvers in nested topology optimization, Struct. Multidiscip. Optim. 42 (1) (2010) 55–72. [4] N. Andréasson, A. Evgrafov, M. Patriksson, An Introduction to Continuous Optimization, Studentlitteratur, 2005. [5] M. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (2) (1988) 197–224. [6] M. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods and Applications, Springer-Verlag, 2003. [7] P. Bochev, R. Lehoucq, On the finite element solution of the pure Neumann problem, SIAM Rev. 47 (1) (2005) 50–66. [8] P. Bochev, R. Lehoucq, Energy principles and finite element methods for pure traction linear elasticity, Comput. Methods Appl. Math. 11 (2) (2011) 173–191. [9] T. Borrvall, J. Petersson, Topology optimization using regularized intermediate density control, Comput. Methods Appl. Mech. Engrg. 190 (37) (2001) 4911–4928. [10] B. Bourdin, Filters in topology optimization, Internat. J. Numer. Methods Engrg. 50 (9) (2001) 2143–2158. [11] D. Braess, Finite elements – Theory, fast solvers, and applications in solid mechanics, third ed., Cambridge University Press, 2007. [12] H. Brezis, Functional Analysis, Sobolev Spaces and Partial Differential Equations, Springer Science & Business Media, 2010. [13] T. Bruns, Zero density lower bounds in topology optimization, Comput. Methods Appl. Mech. Engrg. 196 (2006) 566–578. [14] M. Bruyneel, P. Duysinx, Note on topology optimization of continuum structures including self-weight, Struct. Multidiscip. Optim. 29 (4) (2005) 245–256. [15] J. Cea, K. Malanowski, An example of a max-min problem in partial differential equations, SIAM J. Control 8 (1970) 18. [16] P. Christensen, A. Klarbring, An Introduction to Structural Optimization, Springer, 2009. [17] C.-H. Chuang, S. Chen, R.-J. Yang, P. Vogiatzis, Topology optimization with additive manufacturing consideration for vehicle load path development, Internat. J. Numer. Methods Engrg. (2017). [18] P. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, 1978. [19] P. Ciarlet, Linear and Nonlinear Functional Analysis with Applications, SIAM, 2013. [20] A. Clausen, E. Andreassen, On filter boundary conditions in topology optimization, Struct. Multidiscip. Optim. (2017). [21] B. Dacoronga, Direct Methods in the Calculus of Variations, Springer, 2008. [22] J. Deaton, R. Grandhi, A survey of structural and multidisciplinary continuum topology optimization: post 2000, Struct. Multidiscip. Optim. 49 (1) (2014) 1–38. [23] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Springer, 2004. [24] A. Ern, J.-L. Guermond, Evaluation of the condition number in linear systems arising in finite element approximations, ESAIM: Math. Model. Numer. Anal. 40 (1) (2006) 29–48. [25] A. Evgrafov, M. Gregersen, M. Sørensen, Convergence of cell based finite volume discretizations for problems of control in the conduction coefficients, ESAIM Math. Model. Numer. Anal. 45 (6) (2011) 1059–1080. [26] A. Evgrafov, K. Marhadi, Control in the coefficients with variational crimes: Application to topology optimization of Kirchhoff plates, Comput. Methods Appl. Mech. Engrg. 237–240 (2012) 27–38. [27] L. Hägg, E. Wadbro, Nonlinear filters in topology optimization: existence of solutions and efficient implementation for minimum compliance problems, Struct. Multidiscip. Optim. 55 (3) (2017) 1017–1028. [28] J. Haslinger, M. Kocvara, G. Leugering, M. Stingl, Multidisciplinary free material optimization, SIAM J. Appl. Math. 70 (2010) 2709–2728. [29] F. Jones, Lebesgue Integration on Euclidean Space, Jones and Bartlett Publishers, 1993. [30] E. Kaasschieter, Preconditioned conjugate gradients for solving singular systems, J. Comput. Appl. Math. 24 (1) (1988) 265–275. [31] A. Klarbring, Models of Mechanics, Springer, 2006. [32] M. Kocvara, S. Mohammed, Primal-dual interior point multigrid method for topology optimization, SIAM J. Sci. Comput. 38 (5) (2016) B685–B709. [33] M. Kuchta, K.-A. Mardal, M. Mortensen, Characterization of the space of rigid motions in arbitrary domains, International Center for Numerical Methods in Engineering (CIMNE), Barcelona, Spain, 2015, pp. 259–274. [34] M. Kuchta, K.-A. Mardal, M. Mortensen, On the singular Neumann problem in linear elasticity, Numer. Linear Algebra Appl. (2018) e2212. [35] M. Kˇrížek, Z. Milka, On an unconventional variational method for solving the problem of linear elasticity with Neumann or periodic boundary conditions, Banach Center Publ. 29 (1) (1994) 65–77. [36] P. Lax, Functional Analysis, Wiley, 2002. [37] E. Lee, K. James, J. Martins, Stress-constrained topology optimization with design-dependent loading, Struct. Multidiscip. Optim. 46 (5) (2012) 647–661. [38] L. Liao, A study of inertia relief analysis, in: 52nd Structural Dynamics and Materials Conference. Denver, Colorado: AIAA, 2011, pp. 1–10. [39] A. Limkilde, A. Evgrafov, J. Gravesen, On reducing computational effort in topology optimization: we can go at least this far!, Struct. Multidiscip. Optim. (2018). [40] J. Nocedal, S. Wright, Numerical Optimization, second ed., Springer, 2006. [41] N. Pagaldipti, Y.-K. Shyy, Influence of inertia relief on optimal designs, in: Proceedings of the 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, American Institute of Aeronautics and Astronautics (AIAA), 2004, pp. 616–621.

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.

C.-J. Thore / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

33

[42] M. Papadrakakis, Y. Fragakis, An integrated geometric–algebraic method for solving semi-definite problems in structural mechanics, Comput. Methods Appl. Mech. Engrg. 190 (49) (2001) 6513–6532. [43] H. Pengqiu, Q. Sun, Modified inertia relief method based on accurate inertia loads, AIAA J. 55 (2017). [44] J. Petersson, Some convergence results in perimeter-controlled topology optimization, Comput. Methods Appl. Mech. Engrg. 171 (1–2) (1999) 123–140. [45] J. Petersson, O. Sigmund, Slope constrained topology optimization, Internat. J. Numer. Methods Engrg. 41 (8) (1998) 1417–1434. [46] R. Rajan, P. Naveenchandran, C. Thamotharan, Inertia relief analysis of automotive control arm, Int. Res. J. Eng. Technol. 3 (2016) 3017–3021. [47] E. Savenkov, H. Andrä, O. Iliev, An analysis of one regularization approach for solution of pure Neumann problem, Technical Report 137, Fraunhofer (ITWM), 2008. [48] O. Sigmund, K. Maute, Topology optimization approaches – A comparative review, Struct. Multidiscip. Optim. 48 (6) (2013) 1031–1055. [49] O. Sigmund, J. Petersson, Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima, Struct. Optim. 16 (1) (1998) 68–75. [50] W. Song, H. Geaand, R.-J. Yang, C.-H. Chuang, Application of regional strain energy in topology optimization with inertia relief analysis, in: ASME IDETC/CIE 2013, 2013. [51] M. Stolpe, K. Svanberg, An alternative interpolation scheme for minimum compliance topology optimization, Struct. Multidiscip. Optim. 22 (2001) 116–124. [52] K. Svanberg, The method of moving asymptotes – a new method for structural optimization, Internat. J. Numer. Methods Engrg. 24 (2) (1987) 359–373. [53] K. Svanberg, MMA and GCMMA, versions September 2007, 2007, http://www.math.kth.se/krille/gcmma07.pdf. [54] T. Washizawa, A. Asai, N. Yoshikawa, A new approach for solving singular systems in topology optimization using Krylov subspace methods, Struct. Multidiscip. Optim. 28 (5) (2004) 330–339. [55] N. Wiker, A. Klarbring, T. Borrvall, Topology optimization of regions of Darcy and Stokes flow, Int. J. Numer. Methods Eng. 69 (7) (2007) 1374–1404. [56] J. Wloka, Partial Differential Equations, Cambridge University Press, 1987. [57] A. Yeckel, J. Derby, On setting a pressure datum when computing incompressible flows, Internat. J. Numer. Methods Fluids 29 (1) (1999) 19–34. [58] W. Zhang, J. Zhu, L. Xia, Topology optimization in aircraft and aerospace structures design, Arch. Comput. Methods Eng. 23 (4) (2016) 595–622.

Please cite this article as: C.-J. Thore, Topology optimization of freely floating elastic continua using the inertia relief method, Computer Methods in Applied Mechanics and Engineering (2019) 112733, https://doi.org/10.1016/j.cma.2019.112733.