A sequential global programming approach for two-scale optimization of homogenized multiphysics problems with application to Biot porous media

We present a new approach and an algorithm for solving two-scale material optimization problems to optimize the behaviour of a fluid-saturated porous medium in a given domain. While the state problem is governed by the Biot model describing the fluid–structure interaction in homogenized poroelastic structures, the approach is widely applicable to multiphysics problems involving several macroscopic fields in which homogenization provides the relationship between the microconfigurations and the macroscopic mathematical model. The optimization variables describe the local microstructure design by virtue of the pore shape which determines the effective medium properties, namely the material coefficients, computed by the homogenization method. The numerical optimization strategy involves (a) precomputing a database of the material coefficients associated with the geometric parameters and (b) applying the sequential global programming (SGP) method for solving the problem of macroscopically optimized distribution of material coefficients. Although there are similarities to the free material optimization (FMO) approach, only effective material coefficients are considered admissible, for which a well-defined set of corresponding configurable microstructures exists. Due to the flexibility of the SGP approach, different types of microstructures with fully independent parametrizations can easily be handled. The efficiency of the concept is demonstrated by a series of numerical experiments that show that the SGP method can simultaneously handle multiple types of microstructures with nontrivial parametrizations using a considerably low and stable number of state problems to be solved.


Introduction
The design of fluid-saturated poroelastic media (FSPM) present a gradually increasing topic of research interest due to its mathematical complexity and a great application potential.Although the theory of FSPM has been developed in the context of geomechanics and civil engineering, nowadays theses types of materials are abundant in many engineering applications.A convenient design of microstructures can provide a metamaterial property related to controllable fluid transport, or elasticity.In particular, soft robots can be designed as inflatable porous structures generating a motion and force due to variable fluid content, e.g., [1].To this aim, the behaviour of the fluid-saturated porous materials is described by the Biot model [2], within the small strain theory, which was postulated using a phenomenological approach.The homogenization method enabled the derivation of the quasistatic Biot's equations [3].Since then, a number of works extended the results for the dynamic case, which is important for treating wave propagations, see e.g., [4].As an extension beyond the linear theory, a modified Biot model with strain-dependent poroelastic and permeability coefficients was proposed in [5].
Topology optimization of microstructures constituting the FSPM was treated in [6] and [1].Therein, the fluid-structure interaction problem was handled in the homogenization framework and an approximation towards computational simplification was proposed.
In this paper, we aim at a two-scale approach optimization allowing for a spatial grading of the microstructure design.Two-scale optimization problems have been already extensively discussed in literature before.The whole idea started with the seminal paper of Bendsøe and Kikuchi [7], in which the following concept was suggested: for a given parametrization of the unit cell, carry out the homogenization procedure on a fixed parameter grid in a preprocessing step.Then, in every step of the optimization, first retrieve, for each design element, (approximate) effective material coefficients by interpolation.Next, plug these coefficients into the state equation, solve the latter and evaluate the cost.The other way round, sensitivities are computed by the chain rule, i.e. first differentiate the quantity of interest with respect to the material coefficients and then differentiate the material coefficients with respect to to the design parametrization.This procedure opens the way for the application of any suitable gradient based optimization solver, like, e.g., OCM [8], MMA [9] or SnOpt [10], to name only those, which are most prominently used in structural topology and material optimization.
While this concept essentially carries over to other classes of problems, as it is done by [11][12][13] for thermomechanical settings, we opted to follow a slightly different avenue in this paper.There are several reasons: First, the concept depends, by its nature, to a large extent on the chosen parametrization.If the parameters enter the homogenized properties in a substantially nonconvex way (as it is the case, if, e.g., rotations of the base cells are allowed), many local minima might be introduced and additional measures must be taken to avoid getting trapped in one of them.Second, it is not easy to extend the original concept with respect to the use of completely independent types of unit cells, either characterized by different geometries or material configurations.In this case, specifying a smooth parametrization is non-trivial.The typical idea would be to first introduce an independent parametrization for either cell types (for example using sizing variables) and then add on top a smooth interpolation scheme for the effective tensors as used, for instance, in multi-material optimization (see [14]).The problems with that is however, that the second level of interpolation introduces material coefficients, for which typically no interpretation in terms of a microstructure exists.Thus, an additional penalization strategy is required, which ensures that those unphysical choices do not remain in the optimal solution.Such an approach was successfully demonstrated in the recent work [15].In another recent article, [16] chose two unit cell types, described via level-set functions, such that the mixture of their geometric parameters can be directly interpreted as a third unit cell type.[17] also opted for level-set functions to describe the geometry of the microstructures.But, with respect to the handling of multiple material classes, the authors defined floating patches, where each patch is a subdomain of the design domain and only occupied by one microstructure type.Then, the layout of these patches are optimized on the macroscopic level and their overlaps are combined via a differentiable maximum operator.
In our paper, we describe, how these disadvantages can be circumvented using the SGP concept.The basic idea has been already introduced in [18] and is now generalized to a multiphysics, two-scale setting.This involves an extension of an MMAtype block-separable model function (see [19]) to the poroelastic setting, a split of the computations into an offline and an online phase, which is particularly suited for homogenization based problems, and a numerical solution scheme for the nearly global optimization of block-separable subproblems.We would like to note here that the term block-separable implies that the minimization can be carried out separately for each design element, however a design element itself can be described by multiple design degrees of freedom.For a further motivation of the SGP method, we refer to the first paragraph in section 3. Here, we just like to add that, in the whole optimization process, two different types of sensitivities are relevant.First, there are the sensitivities of constraint or cost functions with respect to the effective material coefficients.These constitute a substantial ingredient of the block-separable model used in the heart of the SGP method.Second, there are the sensitivities of the material coefficients with respect to the chosen parametrization.In the context of the suggested two-scale SGP framework, the latter ones are not strictly required, but can help to come up with an improved interpolation model used in the offline phase.In any case, the derivation of sensitivities presented in this paper, for the particular context of fluid saturated porous media, relies on derivations in [20], where also the sensitivity of the homogenized coefficients were reported, see also [5].
Finally, we would like to comment on the generality of the presented approach.Although the SGP concept outlined in our paper can be applied to a large range of multiphysics two-scale material optimization problems, the Biot model of fluid saturated porous media provides an ideal test bed for the method.This is for several reasons: first, the physical coupling is non-trivial.Second, it is very natural to set up competing objective functions, such as the structural compliance on the one hand and the enhanced fluid flow through an outflow boundary, on the other hand.And third, configurable types of microstructures supporting either the first or the second goal can be deduced in a straightforward manner.
The structure of the remainder of this paper is as follows: In section 2 all ingredients of the two-scale problem are described.To these belong a brief repetition of the constitutive laws for the Biot model (section 2.1), the poroelastic state problem in variational form (section 2.2), a generic sketch of the two-scale problem constrained by the poroelasticity equations (section 2.3) and an adjoint analysis providing sensitivities with respect to effective material coefficients, as used later by the SGP method (section 2.4).Finally, two types microstructures are suggested in form of configurable unit cells (section 2.5).In section 3 the SGP concept for the solution of two-scale optimization problems is introduced in greater detail.For this, the two-scale problem is discretized and extended for the use of multiple types of unit cells (section 3.1).Then, a separable sequential approximation concept is suggested (section 3.2) and last the SGP method is presented in an algorithmic form (section 3.3).In section 4, the advantages of the SGP algorithm will be discussed using various types of two-scale problems.

Formulation of the two-scale optimization problem
In this section, we explain our optimization strategy.Although it can be applied to similar problems involving several physical fields or multiphysics problems, in this paper, we consider the fluid saturated porous media represented by the Biot model which can be derived using the homogenization of the fluid-structure interaction problem restricted to small deformation kinematics, see e.g., [3,21,22].In the next section we report the homogenization result presented

Notation
We employ the following notation.Since we deal with a two-scale problem, we distinguish the "macroscopic" and "microscopic" coordinates, x and y, respectively.We use ∇ x = (∂ x i ) and ∇ y = (∂ y i ) when differentiation with respect to coordinate x and y is used, respectively, whereby ∇ ≡ ∇ x .By e(u) = 1/2[(∇u) T + ∇u], we denote the strain of a vectorial function u, where the transpose operator is indicated by the superscript T .The Lebesgue spaces of 2nd-power integrable functions on an open bounded domain D ⊂ R 3 is denoted by L 2 (D), the Sobolev space W 1,2 (D) of the square integrable vector-valued functions on D including the first order generalized derivative, is abbreviated by H 1 (D).Further, H 1 # (Y m ) is the Sobolev space of vector-valued Y-periodic functions (the subscript #).

The homogenized Biot -Darcy model
We report the homogenization result presented e.g., in [5], cf.[20], where the problem of locally optimized microstructures has been described.The homogenized model of the porous elastic medium incorporates local problems for characteristic responses which are employed to compute the effective material coefficients of the Biot model.The local problems specified below, related to the homogenized model, are defined at the microscopic representative unit cell . which splits into the solid part occupying domain Y m and the complementary channel part Y c .Thus, where where ID = (D ijkl ) is the elasticity tensor satisfying the usual symmetries, ) is the linear strain tensor associated with the displacement field v.
In what follows, by the microstructure Y(x), we mean the decomposition eq. ( 1) of the representative cell Y and the material properties, as represented by the elasticity ID only in our case.If the structure is perfectly periodic, microstructures Y ≡ Y(x) are independent of the macroscopic position x ∈ Ω.Otherwise, the local problems must be considered at any macroscopic position, i.e. for almost any x ∈ Ω, see e.g., [21] in the context of slowly varying "quasi-periodic" microstructures.It should be pointed out, that this issue is of a special importance when dealing with homogenization-based material design optimization; as will be explained below, a regularization is required to control the design variation within Ω.
The local microstructural response is obtained by solving the following decoupled problems: where • ∀v ∈ H 1 # (Y c ) and ∀q ∈ L 2 (Y c ). Effective material properties of the homogenized deformable fluid-saturated porous medium are described in terms of homogenized poroelastic coefficients: the drained elasticity A A, the stress coupling C and the compressibility N , all being related to the solid skeleton.All these coefficients including the intrinsinc hydraulic permeability K are computed using the characteristic microscopic responses eqs.(3) to (5) substituted in following expressions: Obviously, the tensors A A = (A ijkl ), C = (C ij ) and K = (K ij ) are symmetric, A A adheres all the symmetries of ID; moreover A A is positive definite and N > 0. The hydraulic permeability K is, in general, positive semi-definite.It is positive definite whenever the channels constitute a simply connected domain generated as the periodic lattice by Y c ; for this, denoting by Γ k Y ⊂ ∂Y , k = 1, . . ., 6 the faces of Y , it must hold that Γ k Y ∩ ∂Y c = ∅ for all k = 1, . . ., 6.

Coupled flow deformation problem
The Biot-Darcy model of poroelastic media for quasi-static, evolutionary problems imposed in Ω is constituted by the following equations involving stress σ, displacement u, strain e(u), fluid pressure p and the seepage velocity w: where the homogenized coefficients are given by eq. ( 6) and Above, η is the relative fluid viscosity, γ is the fluid compressibility and φ = |Y c |/|Y | is the porosity (volume fraction of the fluid-filled channels).The effective volume forces in eq. ( 7), acting in the solid and fluid phases, are denoted by f s and f f , respectively.It is important to note that η = η phys /ε 2 0 is defined for a given fluid (η phys ) and microstructures scale: ε 0 = 0 /L where L is a characteristic macroscopic length, and 0 is the characteristic microstructure size, typically given by the "pore diameter".Thus, for a given fluid, the effective permeability K/η is proportional to ε 2 0 , i.e. reflecting the microstructure size.In contrast, all other coefficients are scale-independent (when the scale separation holds, i.e. ε 0 being small enough).
Remark 1.In this paper, we only consider steady state problems for the Biot medium, such that all time derivatives in eq.(7) vanish.Consequently, the Biot compressibility M is not involved, as far as the porous phase, generated as a periodic lattice by channels Y c , is connected.For any microstructure with disconnected pores, such that Y c ⊂ Y , thus, Y c constitute one, or more inclusions with one cell Y , see [22], the permeability vanishes.Then, the time integration in eq. ( 7) leads to the mass conservation equation in the form B : e(u) + M p = 0, assuming an undeformed initial configuration with the zero pressure in the inclusions.In the optimization problem, besides microstructures with nondegenerate permeabilities, we shall consider also microstructures with spherical, thus, disconnected pores, constituting impermeable material.For this case, one can choose either fluid filled pores, or empty pores; the only difference is the use of the so-called undrained material elasticity, A A U = A A + M −1 B ⊗ B, or the elasticity A A describing effective elasticity of the "drained" skeleton, with empty pores.

State problem formulation
Let Ω ⊂ R 3 be an open bounded domain.Its boundary ∂Ω splits, as follows: p , and Γ 1 p ∩ Γ 2 p = ∅.We consider the steady state problems for the linear Biot continuum occupying domain Ω.The poroelastic material parameters and the hydraulic permeability referred to as the homogenized coefficients, in general, are given by the locally defined microstructures Y(x) which can vary with x ∈ Ω.The two-scale optimization approach proposed in this paper enables to combine microstructures characterized by connected and disconnected pores, the latter characterized by a vanishing permeability.To this aim, the domain Ω = Ω 0 ∪ Ω + is decomposed into in two parts: the permeable Ω + and the impermeable Ω 0 , which may not constitute connected domains, being split into more disconnected subparts.Consequently, the interface Γ + = ∂Ω + ∩ ∂Ω 0 is impermeable.Regarding the boundary decomposition, we assume that Γ k p+ := Γ k p ∩ ∂Ω + = ∅, for k = 1, 2, so that the porous structure permits the fluid transport through domain Ω + , if this one connects Γ 1 p+ and Γ 2 p+ .We consider the following macroscopic problem: Given the traction surface forces g, and pressures pk on boundaries Γ k p , find displacements u and the hydraulic pressure P which satisfy where P = 0 in Ω 0 .Whereas, in Ω + , P satisfies For the steady state problem the set of equations eq. ( 7) yields the two problems eq.( 9) and eq.( 10) as a decoupled system: first, eq. ( 10) can be solved for P , then eq. ( 9) is solved for u.Moreover, for the considered type of the boundary conditions and since volume forces are not involved, the solutions are independent of the viscosity η, see eq. (7).
Further, we consider an extension of pk from boundary Γ k p to the whole domain Ω, such that pk = 0 on Γ l p (in the sense of traces) for l = k.Then P = p + k pk in Ω + , such that p = 0 on Γ p+ .Note that p can be simply extended by 0 in Ω 0 .For the sake of notational simplicity, we introduce p = k pk .By virtue of the Dirichlet boundary conditions for u and p, we introduce the following spaces: We employ the bilinear forms and the linear functional g, In order to define the state problem in the context of two-scale optimization, we employ the weak formulation which reads, as follows: Find u ∈ V 0 and p ∈ Q 0 , such that, for all v ∈ V 0 and q ∈ Q 0 , To define p uniquely in Ω, p ≡ 0 in Ω 0 = Ω \ Ω + .
Since the two fields are decoupled, first p is solved from eq. ( 13) 2 , then u is solved from eq. ( 13) 1 , where p is already known.
Remark 2. In the context of the undrained porosity defined by fluid-filled closed pores Y c ⊂ Y , see Remark 1, formulation eq. ( 13) is consistent also with this microstructure class type Y 0 with A A U replacing A A in the elasticity bilinear form eq. ( 12) 1 .Pressure is then defined pointwise in Ω 0 by P := −B : e(u)/M .
By α(x) we denote an abstract optimization variable which determines the homogenized coefficients for any position x ∈ Ω. Below we consider α representing several geometrical parameters characterizing microstructures Y(x) of a given type.Although, in this section, we disregard some particular details related to the treatment of multiple types of Y, we bear in mind the existence of two microstructure classes, Y + and Y 0 , associated with the pore connectivity type, as discussed above.The "permeable" domain Ω + is occupied by the material given pointwise by Y(x) ∈ Y + for all x ∈ Ω + .Hence, both the subdomains of Ω are defined implicitly by the microstructure type: Ω i is the set of x ∈ Ω, such that Y(x) ∈ Y i , where i = +, 0.
In the next section, we shall consider a twoscale optimization problem which is characterized by the following features: • Geometrical restrictions are stated in respective definitions of the admissibility designs sets for a chosen type of microstructure.For the sake of brevity, let A be the set of admissible designs, further we consider α(x) ∈ A for any x ∈ Ω. • We consider multiple optimization criteria which perform as the objective functions, or equality constraints.Without loss of generality, we confine ourselves to the two criteria Φ α (u) and Ψ α (p) that are defined, as follows: While Φ α (u) expresses the structural compliance, criterion function Ψ α (p) expresses the amount of the fluid flow through surface Γ 2 p due to the pressure difference p1 − p2 , see the boundary condition eq. ( 10) 2 .These two criteria are antagonist: the pore volume reduction leads naturally to stiffening the structure, but reduces the permeability.Hence, for the objective function Φ α , function Ψ α serves as a constraint and vice versa.

Two-scale optimization problem
Here, for the ease of notation, we restrict to one microstructure type only, namely Y(x) ∈ Y + , so that we may consider Ω ≡ Ω + .Hence, all the bilinear forms in eq. ( 12) are defined by integration in Ω.Later, in section 3, we will consider microstructures characterized by different unit cell types of classes Y + and Y 0 , however, the formulations introduced below can be adapted easily.
We first define the direct optimization problem to find design α(Ω) that minimizes a cost functional based on the criteria defined in eq. ( 14).Further, we introduce the set T = S 6 × S 3 × S 3 × R × R and denote by IH = (A A, B, K, ρ m , R) ∈ T the (local) material parameters involing the effective (homogenized) material coefficients, the solid part volume ρ m = 1−φ = |Y m |/|Y |, and a regularization parameter R, which typically depends only on the design.We note that the dimension of the regularization label R is, for ease of notation, chosen as 1 for now, although later in section 4.3 more general regularization labels are used.Obviously, IH is given uniquely by the local admissible design α(x) ∈ A, x ∈ Ω, whereby for a suitably chosen parametrization, the admissibility set is given simply by Examples for such parametrizations along with a description of the lower and upper bounds a, a ∈ R n are presented in section 2.5.For a given admissible design α(Ω), the state z = (u, p) is the solution of eq. ( 13), where the homogenized coefficients IH(α) are given in eq. ( 6) using the characteristic responses W (α) := (ω ij , ω P , ψ k , π k ).W (α) are the solutions of eqs.(3) to (5), which depend on α(x) in terms of the microconfigurations Y(x).In this way, mapping S : α(Ω) → z(Ω) introduces the admissible state.
It can be defined by a composition map, S = Z • E • W, where W represents the resolvents of the characteristic problems imposed on the local microconfigurations, E provides the homogenized material, and Z is the resolvent of the macroscopic state problem, so that Further, we employ the mapping The macroscopic state problem is the implicit form of the mapping where S 0 is the space of admissible state problem solutions.For the Biot medium problem, eq. ( 16) is identified with eq. ( 13).

Direct two-scale optimization problem
For the given two functions of interest Φ and Ψ, both depending on the material distribution IH(x) and the state z(x), the two-scale abstract optimization problem reads: where the term Ξ(IH) in the objective is related to the design regularization, namely to parameter R, and Λ Ξ ∈ R + is a penalty parameter.Recall the chain mapping H : α(x) → IH(x) for any x ∈ Ω, then z = Z(Ω).Below, we abbreviate Φ α (z) =: Φ(H(α), z) and also Ψ α (z) =: Ψ(H(α), z).In eq. ( 14), specific examples relevant for the Biot medium optimization were given.Optimization problem eq. ( 17) is associated with the following inf-sup problem, with the Lagrangian function, where Λ = (Λ Φ , Λ Ψ ) ∈ R 2 are the Lagrange multipliers associated with the objective and constraint functionals Φ and Ψ, and z ∈ S 0 are Lagrange multipliers -the adjoint variables -associated with the constraints of the problem eq. ( 17).
For a while, we may consider material coefficients IH as the optimization variables (although they are parameterized by α ∈ A).Further, let us assume a given value Λ ∈ R 2 ; note that the entries of Λ can be positive or negative depending on the desired flow augmentation, or reduction.In the numerical examples, we chose Λ Φ > 0, whereas Λ Ψ < 0 indicates the constraint effect of Ψ relative to Φ. Upon denoting by Im(H) = H(A), the image space of all admissible designs, and defining for a.a.x ∈ Ω} , the optimization problem eq. ( 17) can be rephrased as the two-criteria minimization problem, where For the Biot medium optimization, where the two criterion functions Φ α and Ψ α are given in eq. ( 14), the Lagrangian function attains the form

Adjoint responses and the sensitivity analysis
In this section, we provide details concerning the sensitivity analysis employed in the preceding section.We consider α to represent a general optimization variable which is related to the effective medium parameters IH.It is worth to note that one may also consider α ≡ IH in the context of the free material optimization (FMO).
To obtain the adjoint equation, we consider the optimality condition for (u, p).Thus, from eq. ( 21) it follows that δ (u,p) L(α, (u, p), Λ, (ṽ, q)) where To avoid computation of the gradient ∇q on Γ 2 p ⊂ ∂Ω, we consider p ∈ H 1 (Ω) such that p = 0 on Γ \ Γ 2 p , while p = 1 on Γ 2 p , then it is easy to see that The optimality conditions eq. ( 22), related to the state admissibility, yield the adjoint state (ṽ, q) ∈ V 0 × Q 0 which satisfies the following identities: These equations can be rewritten using eq.( 23) and eq.( 24), as follows for all (ṽ, q) ∈ V 0 × Q 0 : To allow for the independence of the state adjoint on Λ, we define the split ṽ = Λ Φ θ , where θ and qk , k = 1, 2 satisfy for all (ṽ, q, q) We can compute the total variation of the Lagrangian with If the pair (u, p) solves the state problem and (ṽ, q) is its adjoint state, eq. ( 29) is equivalent to the following expression: Above, the shape derivatives δ α of the bilinear forms can be rewritten in terms of the sensitivity of the homogenized coefficients.Besides the obviously vanishing derivative δ α g(u) = 0, it holds that Using the "total pressure" P := p+p, the following tensors are employed to evaluate the expression in eq. ( 31): Now, using these tensors, eq. ( 29) is computed, as follows: Hence the variations of L with respect to A A, B and K are given by the following formulae As Ξ(IH) solely depends on the regularization parameter R, see eq. ( 47), we get for the regularization term in eq. ( 33).In the context of the finite element discretization introduced in section 3, the homogenized coefficients are supplied as constants in each element Ω e of the partitioned domain Ω.Accordingly, the expressions in eq. ( 32) are supplied elementwise at the Gauss integration points.

Design parametrization
The design of the cell Y , that is the decomposition into the solid skeleton Y m and the pores Y c , can be parameterized in a number of ways.In [20], we employed a so-called spline-box structure parameterized by design variables defining positions of the spline control polyhedron.This kind of parametrization is convenient due to its generality to handle quite arbitrary design, but leads to complicated formulations of design constraints which are needed to preserve essential geometrical requirements (e.g., positivity of channel crosssections).
In this paper, we employ two specific types of microstructures illustrated in fig. 1, where the channels are shaped as a 3D cross (type 1), or a sphere (type 2).Hence, the latter microstructure is featured by zero permeability and therefore, we consider dry pores (voids) in the mechanical model.Due to these specific geometries, we can use a rather simple parametrization, which is listed in table 1.For a unit cell of type 1, r x and r y refer to the radii of the cylinders pointing in xand y-direction respectively.The third parameter ϕ describes the cell rotation, about axis z.For the unit cell type 2, the spherical voids, whose radii are described by r s , provide an orthotropic material with nearly isotropic elastic properties.Therefore, rotations are not enabled for this cell type.Importantly, box constraints can be imposed on r x , r y and r s straightforwardly to guarantee geometric feasibility.
Table 1: The parametrization of the pore geometry for the two types of the microstructures: 1: the 3D cross, 2: the sphere.To illustrate a sensitivity of the material properties determined by the homogenized coefficients IH, In fig.2, for unit cell type 2, the elasticity as the only relevant material property is displayed as function of r s .In fig.3, for unit cell type 1, selected components of the poroelastic tensors and of the permeability are reported as functions of r y .

A Sequential Global Programming formulation
The basic description of the Sequential Global Programming algorithm along with convergence aspects were presented in [18], where SGP was applied to a multi-material optimization based on a two-dimensional time harmonic Helmholtz state equation.The setting and procedure described in this manuscript differs from the one in [18] in the following major points: first, in [18] a selection of finitely many fixed materials was considered as admissible set.In this paper, each admissible material is computed by homogenizing unit cell, which itself is configurable by a number of geometric parameters.Thus, the designer can choose in each point of the design domain from M different unit cell types and adjust the geometric parameters for the latter.Second, the SGP approach is extended to a multi-physics setting using a slightly different separable approximation and third, a different solution strategy is employed for the subproblems arising from this.This strategy does not impose any assumption on the parametrization.In particular, parametrizations can be non-analytical and non-differentiable.This leads to a greater design flexibility.Despite these differences, there is also an important feature, the approach presented here has in common with the one outlined in [18]: separable models are established in terms of (effective) material tensors IH rather than their parameterization α.
Then, the parametrization is directly treated at the level of sub-problems without further convexification.Thanks to the separable character of the chosen first order model the resulting generally non-convex sub-problems can -in principal -still be solved to global optimality.The advantages of this approach are twofold: first, due to the separable model functions being able to capture also non-convex features of the original cost function typically a low number of outer iterations, equivalently to the number of state problems to be solved, is required; and second, due to the good fit of the separable models with the cost function as well as the fact that non-convex sub-problems are solved to global optimality the overall algorithm is less start value dependent and less prone to be trapped in poor local minima.This is in contrast to traditional approaches, where a local model is established directly based on the sensitivity of cost functions with respect to the design parameterization α.
In the following we first derive a fullly discretized counterpart for a slightly generalized of problem eq. ( 20).Then we describe in detail how the separable first order approximations can be constructed and finally present a practical outline of the full SGP algorithm including a generic sub-solver allowing to compute near globally optimal solutions for sub-problems using a brute-force strategy.

A fully discretized 2-scale design problem
For the sake of simplicity, the definitions of sets and functions were introduced in sections 2.2 and 2.3 based on the assumption that there is only one type of unit cell such that M = 1.
Here, for a more general setting, we consider M unit cell types, each one with n i design parameters, and introduce index set I := {1, . . ., M }.For each unit cell type i ∈ I, the admissibility set is defined in terms of box constraints and other purely geometrical constraints.By choosing a suitable parameterization, we can identify these with (geometric) parameter sets with a i , a i ∈ R ni being lower and upper bound vectors constraining the corresponding parameter vector α i ∈ R ni .
Remark 3. We note that, while in this manuscript the parameters in eq. ( 35) are always used to vary the geometrical properties of the unit cell, variations in the material parameters could be described in the same way.Thus, SGP can handle both of these situations.
We further define for all i ∈ I map where H i (α) performs the homogenization procedure described in section 2.3.fig. 4 illustrates the components of H i (α i ).
We denote the union of the ranges of all H i by and with that generalize the set of admissible design functions to become for a.e.x ∈ Ω} .

Now the state problem operator
with displacement function u(IH) and hydraulic pressure function p(IH) reads exactly as before.
We finally use a slightly more general resource function than in sections 2.2 and 2.3 as follows: A concretization could be the total volume fraction of a specific material phase (see description of ρm in section 2.3).
Based on these definitions, we then formulate an FMO-type problem min where ρm ∈ R is the resource constraint value and cost functions and Φ, Ψ, Ξ and their weights Λ Ψ , Λ Φ , Λ Ξ have been already introduced in section 2.3).
Although problem eq. ( 40) is formulated directly in the tensor variable IH, a realization of the feasibility condition IH ∈ U ad would force us to evaluate the homogenization maps H i (i ∈ I).This has the consequence that for each evaluation of the cost function, a homogenization procedure, which contains a series of cell problems, has to be conducted.To alleviate this situation, we follow [7] and carry out the homogenization procedure only for discrete samples of the design parameter space.For each unit cell type i, we introduce a grid with nodes A nodes i ⊆ A i and effective material coefficients are only computed, via homogenization, at the sampled nodes of this grid.In addition, we define a piecewise cubic Hermite interpolator for these samples to realize the continuous mapping Hi : for all i ∈ I.We denominate this procedure as the offline phase of a two-scale optimization approach, as it can be performed independent from the online optimization procedure that is subject to constraints, that go beyond the box constraints on the parameter sets as in eq. ( 35).For the case M = 1, the conventional approach would be now, to perform the optimization based on the interpolated functions H1 over the full parameter set A 1 .This is not directly possible for M > 1.One way to get around this would be to introduce another interpolation between the different unit cell types similar as it is done in discrete material optimization (DMO) [14].Rather than that we introduce design grids for all unit cell types.Only elements of A grid i , i ∈ I will be considered in the optimization process later.This way, in general, only an approximate solution of the design problem can be computed.However it will turn out that this strategy combines well with the separable non-convex model introduced later in section 3.2.Moreover the resulting error can be easily controlled by the distance and number of samples in A grid i , i ∈ I.The relation of different grids and mappings for the material coefficients are visualized and elaborated in fig. 5.
As we only optimize on A grid i , i ∈ I, eq. ( 37) is approximated by We note that elements of H can be precomputed already in the offline phase.In general, this leads to a higher memory requirement, but additionally reduces online computation time.Finally, we briefly introduce a finite element approximation, with n el finite elements, and therefore introduce element index set E := {1, . . ., n el } to indicate a finite element distinctively by its index e ∈ E. We further assume that the design is constant on each element and can thus be represented by IH ∈ Hn el We remark that through the definition of H in eq. ( 43) this condition already states that only material tensors are eligible, for which a unit cell type i and a parameter vector α i in A grid i exists.Moreover, we replace physical functions Φ and Ψ, regularization function Ξ and solution operator Z by their discretized counterparts, e.g., where n dof is the dimension of the discrete state solution space.The discretized version of resource function ρ eq. ( 39) is The optimization problem, fully discretized in design and state space, then reads min with We note that we have eliminated the resource constraint by the Lagrange formalism.Later we will suggest to use a bisection strategy as introduced in [8] for the framework of the well known OCM method.We finally specialize the regularization term to become where F denotes a standard density filter function (see, e.g., [23]) with and R is the vector of regularization labels associated with all finite elements e ∈ E.

Construction of subproblems
For any sequential programming algorithm first a sequence of subproblems has to be defined.Here, in each iteration k, we construct separable first order approximations, about an expansion point IH k ∈ Hn el , for the components of cost function of the original optimization problem in eq. ( 46).The model problem is min where our model function is defined as with ID e := (A A e , B e , K e ) ∈ S 6 × S 3 × S 3 , In the following, we describe each component of J sep in more details.For this, we split J (IH, λ ρ ) as From tuple IH, only the effective material coefficients A A, B and K, are relevant for J phys .Consequently, for J phys , we define a separable approximation of type where J phys is the following generalization of the first-order MMA-like model suggested in [19] for functions defined in tensor variables: Here C phys is a constant that is chosen to establish the zeroth order correctness of the model and < •, • > {S 6 ,S 3 } denotes the Frobenius inner products for matrices from S 6 and S 3 , respectively.It is further mentioned that in contrast to the model in [19], we refrain from working with flexible generalized asymptotes L A A e ∈ S 6 , L B e , L K e ∈ S 3 , but simply choose all of them to be zero matrices.The partial derivatives of J phys with respect to the material coefficients A A, B and K can be easily extracted from the expressions in eq. ( 34).
The function J vol that describes the fraction of utilized matrix material, is separable by definition, and depends solely on ρ m .We accordingly choose The function J reg given in eq. ( 54) solely depends on the regularization label R ∈ R n el , which is a component of tuple IH ∈ Hn el .The separable approximation of J reg is thus of the form where In eq. ( 59), we further employ function in which the regularization label is varied only in the e-th entry by value R, and contributions of expansion point R k are used in the neighboring entries.Is is noted that eq. ( 59) can be reduced to a convex quadratic function of type by precomputing a e , b e , c e ∈ R, which are independent from R e .Finally, we implement a step size control for the design from one iteration to the next one by adding with a positive factor Λ g to the model cost function.Alternatively, a more general globalization strategy, similar to the regularization approach with regularization label R in eq. ( 59), could be pursued by introducing particular globalization labels.Here, we assume that evaluating the design step size based on the stiffness tensor A A e and A A k e is sufficient, and, in particular, the uniqueness of the globalization labels, such that is satisfied.

The SGP algorithm with a brute-force sub-solver
Having at hand the separable first-order approximations of the objective function and penalization terms, we are now able to formulate the iterative scheme that is described by algorithm 1.We make extensively use of the separable structure of and solve the subproblems, of each iteration k, for each finite element e ∈ E individually.This is done by evaluating J sep,e for all (finitely many) IH e ∈ H and, based on these evaluations, identifying a global minimizer IH * e .Note that, with each IH e , a unique geometric cell label α e is associated and thus, by determining IH * e , we also determine respective α * e and material class index i * .As mentioned already earlier a bisection strategy is applied to treat the resource constraint, see algorithm 2 for the details.To keep things simple, it is assumed that the resource constraint is always active at a minimizer.If no resource constraint is applied, the outer loop in algorithm 2 is simply omitted.
After each iteration, the original cost function J is evaluated with the current solution of the subproblems IH * e .If a descent in J was achieved, we continue the iterative process.If not, we employ the step width control, by increasing multiplier Λ g of globalization term eq. ( 60), and resolve the subproblems using algorithm 2.
Algorithm 1 Sequential Global Programming for parametrized multi-material optimization while J diff < 0 do k ← k + 1 13: end while

Numerical results
In this section, we demonstrate the abilities of SGP by means of numerical examples.It is build up successively by first increasing the design freedom to the two-scale optimization problem, while observing the respective optimized designs and then studying the effect of regularization.
Algorithm 2 Solve subproblems via brute force strategy 1: initialize λ ρ ∈ R for volume bisection 2: while volume constraint is not satisfied do IH * e ← evaluate Hi * (α * ) (see eq. ( 41)) 10: end for 11: ρ ← evaluate ρ h (IH * e ) (see eq. ( 45)); 12: if ρ > ρm then end if 17: end while In section 4.1, we start with the unit cell that is constructed by three intersection fluid channels, visualized in the top row of fig. 1, and study the impact of the micro-structure's local orientation on the performance of the optimized designs.It will be seen that, thanks to the strength of our model, we do neither have to use smart initial orientations, as proposed e.g., in [24,25] by aligning the anisotropic material with respect to principal directions of the stress tensor, nor we have to enforce artificially a regular design.
Then, we present a pareto front and investigate the influence of different weightings of compliance and fluid flux, in the cost function, on the resulting designs.When we proceed from one point on the Pareto front to the next one, we intentionally refrain from using the previous design as a warm start.Nevertheless and despite the non-convex character of our weighted cost function, Pareto curves are obtained, in which none of the points is dominated by another one.We trace this observation back to the ability of the SGP method to avoid poor local solutions.
In section 4.2, we proceed to demonstrate the ability of SGP to handle more than one unit cell type.We again compute a Pareto curve for this case.It will be observed that the new Pareto front is, due to the increase in the design freedom, is strictly dominating the previous one.It will be observed that the more complex parametrization does on average not lead to an increase in the number of state problems to be solved per optimization run.Note that for the settings presented in section 4.1 and section 4.2, it was not necessary to employ a globalization strategy to control design changes from one iteration to the next one.Thus, we set the globalization parameter Λ g = 0.
In the end, in section 4.3, we apply a filtering technique onto the design parameters to both control the speed of variation of local orientation, as well as the interface length between the two unit cell types.Here, we also employ the globalization term described in eq. ( 60).
The setting of the poroelastic problem is depicted in fig.6.It is a recapitulation of the macroscopic problem setting from [20], where the authors selected a finite element from the macroscopic domain and optimized the shape of the local microstructure via a spline box approach.In the present paper we provide an extension to this example by solving the two-scale optimization problem with the SGP method described in section 3. We note that we work with a rather coarse discretization of the macroscopic domain.The reason is that such a discretization is sufficient to demonstrate the capabilities of SGP as described above.On the other hand, it is readily seen in algorithm 2 that the number of macroscopic elements enters the computational complexity for SGP linearly.Thus, in principle there is no obstacle to work with finer discretizations.

Optimization with one unit cell type
In this section, we employ unit cell type 1, depicted in fig.For this grid, we chose a sample size of 28 for each of the two channel radii; again the samples are equally spaced.For the following optimization results with the weighted sum formulation of structural compliance and fluid flux, we employ an initial design guess, visualized in fig.8, that is neither particularly favorable for the mechanical nor for the fluid flow state.
Fig. 8: Homogeneous initial design with r x = r y = 0.15 and no cell rotation and physical performance Φ init = 28.9 and Ψ init = 0.135.
For the described setting, we choose Λ Ψ = −10 and obtain the optimized design shown in fig.9a.Note that the design domain is discretized by two finite element layers in z-direction.We made the experience that, for all numerical results presented in this paper, the differences of optimized designs at layer z = 0 and layer z = 1 are so small such that they cannot be visually discernible.For this reason, we will only show optimized designs for layer z = 0 in the rest of the paper.
SGP stopped after 19 iterations, because the difference between the objective values of the old and new design was found to be 0. We note that this comparably low number of iterations is related to the fineness of the design discretization.Thus, using more grid points could lead to a slightly larger number of iterations.On the other hand, in those experiments that we performed in this direction, the visualizations of the obtained result could be hardly distinguished, see fig. 10.This is why we do not report results for different choices of A grid i , i ∈ I.A second observation we can make is that the fluid channels in resulting designs are fully connected.This is due to the fact that no rotational design degrees of freedom were used.On the other hand we will see next that      Again, the initial guess is neither particularly favorable for the mechanical nor for the fluid flow state.After the first iteration, we see in fig.11a that some channels, close to the outflow region, are opened widely and cells closer to the mechanical support were adjusted to have narrower fluid channels to improve the mechanical performance of the design.In comparison to the solution in fig.9, where the orientation was fixed, this solution has a 1% smaller compliance and a fluid flux which is about 47% higher.We would like to emphasize that local orientation field looks rather smooth although we have neither applied a stress based warm start for the rotation variable, as proposed by [24,25], nor we have employed a regularization technique.We also can observe that the total number of iterations required did not increase after addition of the additional design degrees of freedom.
We conclude this subsection by presenting a Pareto front for this type of bicriterial weighted sum formulation in fig.13.All optimizations were based on the initial guess that is shown in fig.8.This implies that again, no warm starting technique was employed to proceed from one point to the next on the Pareto curve.Nevertheless a Pareto curve is obtained, in which none of the points is dominated by another one.This again is a hint that the SGP method is able to avoid poor local solutions.The number of outer iterations required to solve the problems corresponding to all points on the Pareto curve varied between 3 and 31.The rather low number of 3 iterations was obtained for the extreme case, where Λ Ψ = 0.The optimized designs for various choices of Λ Ψ are visualized in fig.14.It is observed that the with decreasing Λ Ψ the compliance minimized is design (fig.14a) is almost smoothly transformed into a fully flux based design (fig.14h).

Optimization with two unit cell types
We want to study the ability of SGP to handle more than one unit cell type.For this purpose, we add unit cell type 2 that comprises of a void sphere surrounded by matrix material (see second row of fig.1).The only design parameter is the radius r s ∈ [0.1, 0.4] of the void sphere in this case.The smaller the void sphere, the higher the volume fraction of the matrix phase and therefore the stiffer the cell.Thus, cells of type 2 are particularly favorable for the mechanical part of the objective.When only optimizing the compliance, we obtain the trivial solution shown in fig.15.For the fluid flow, cells of type 2 are futile as they are not permeable.However, for numerical we set the permeability of the latter cells to 0.001.Cells of type 1 have orthotropic mechanical properties and transversal isotropic permeability tensors, whereas cells of type 2 have isotropic mechanical properties and no permeability.Although cell types 1 and 2 are disjunct in their parameter spaces, the corresponding ranges of volume fractions, of the stiff matrix material, overlap.We have ρ (H 1 ([0.08,0.08])) = 87.9%,ρ (H 1 ([0.22,0.22])) = 71.54%and ρ (H 2 (0.4)) = 73.19%,ρ (H 2 (0.1)) = 99.6%.A nodes 2 , the basis for the interpolation of H 2 , consisted of 30 uniformly distributed samples for r s ∈ stress that we did not use enhanced initial designs for the computation of the points on the Pareto curve.The comparison of the new (blue) curve with the old (red) curve shows that consistently better designs are obtained.Points on the blue curve strictly dominate points on the red curve in the Pareto sense.This is not surprising as, with the addition of a new unit cell type, the design freedom is increased.Still it is worth to mention that the fact that we do not observe any outliers in this respect again underlines the stability of our SGP method.The numbers of required outer iterations varied between 4 and 40, which means that no significant increase in the number of iterations is observed, although a second cell type has been added.In fig.17, we can observe how the number of cells of type 2, in the optimized design, decreases with decreasing Λ Ψ .This is expected, as cell type 2 is completely useless for a flux favored design.
We note that so far all results presented have been computed without employing a resource constraint.Just to demonstrate that SGP can also easily handle problems, where a resource constraint is added, we briefly discuss a selected result in fig.18.

Optimization with both cell types and regularization of design labels and interface
We introduce a regularization of the optimization problem by applying a weighted-sum filter F (e.g., [23,26]), that is often used in the context of topology optimization, on regularization labels that are directly related to the unit cells' geometric parameters.For this we introduce mappings This choice of labeling has the following effects: Within type 1, the maximal distance from lower to upper label bound is 1.This is the same distance required to jump from the stiffest cell of type 1, with r x = r y = 0.08, to any cell of type 2. Therefore, the interface between cells of type    (R 1 ) 3 is employed to circumvent disambiguities for the angular variable.
Employing these regularization labels, J reg from eq. (59) changes to where R ∈ R n el collects the -the components of the regularization label assigned to each finite element, which is defined by formula eq. ( 63) or eq.( 64), if cell type 1 or cell type 2 is chosen for the corresponding finite element e, respectively.Next, we study the influence of regularization with the optimized result for the particular choice Λ Ψ = −3.The result displayed in fig.19

Conclusion and Outlook
We presented an Sequential Global Programming (SGP) approach to homogenization-based structural optimization which can be viewed as an free material optimization constrained by the set of admissible geometric material parameters.
By means of numerical examples, where we successively added more ingredients to the optimization problem, we demonstrated that the proposed SGP approach, with its first-order approximations, provides good and reasonable optimized designs without the necessity of particular design initialization or the employment of a regularization strategy for purposes of convergence.Furthermore, SGP is able to handle several material classes with disjunct parameter sets without additional interpolation and penalization strategies.We further observed that optimizing the local orientation of the microstructure brings along a significant improvement, up to 48%, of the fluid flux.We have not actively addressed the subject of connectivity within the microstructure, that is to ensure connectivity of the fluid saturated channels.However, the regularization approach presented in section 4.3 can be used to control the degree of variation of the local microstructure rotation and we have seen, by means of the presented numerical examples, that only a mild regularization has already a fair impact on the design.
Although the resolution of the finite element approximation, and thus the number of design elements, of the examples in section section 4 was chosen rather coarsely, it served the purpose of demonstrating the presented features of SGP.With regard to finer resolutions: the algorithm can be well parallelized with respect to the design elements due to the block-separability of the first-order approximations.
The brute-force approach in the subproblem solver, described in algorithm 2, can further be speeded up by employing a hierarchical scanning of the design grids A grid i : Start with a rather coarse number of samples and determine the minimizer among those.In the next level, consider only the current minimizer and its neighbors and perform the same search within this subset of A grid i , for all i ∈ I. Repeat this step until the maximum desired number of levels or some accuracy is achieved.Note that, with this strategy, the quality of the design depends on the number of samples on the coarsest grid level.An alternative would be to apply a Lipschitz optimization solver, see [27], to each design element and type in a black box manner.
Further research will focus on extending the SGP approach for homogenization-based optimization to transient problems and, in particular, to dynamic metamaterial design.Another challenge is to extend the proposed optimization approach for an approximate treatment of nonlinear two-scale problems with the homogenized coefficients depending on the macroscopic response by virtue of the sensitivity analysis as discussed in [5].
by Y d for d = m, c, we denote the closure of the open bounded domain Y d .By ∼ Y d = |Y | −1 Y d , with Y d ⊂ Y for d = m, c, we denote the local average (|Y | is the volume of domain Y ).Obviously, the unit volume |Y | = 1 can always be chosen.We employ the usual elasticity bilinear form, involving two vector fields w and v, that reads

Fig. 1 :
Fig. 1: Parametrization of unit cells: unit cell type 1 is parameterized by radii r x and r y , both ranging from 0.08 to 0.22, r z = 0.15 and r s = 0.25 are kept constant; unit cell type 2 is parameterized by radius r s ranging from 0.1 to 0.4.

Fig. 2 :
Fig. 2: Unit cell type 2: dependence of A 1111 on parameter r s .

Fig. 3 :
Fig. 3: Unit cell type 1: dependence of homogenized coefficients A A, B, and K on r y ; r x = 0.15 is fixed.

Fig. 4 :
Fig. 4: Collection of materials: each material, represented by a unit cell object, comes along with a collection of data such as geometric parameters, physical properties and further labels.

Fig. 5 :
Fig.5: Left: Sketch of parameter set A i and samples from its subsets A nodes i (blue dots), that serves as a construction basis of interpolated Hi , and A grid i (red squares), on which the optimization process is performed.In general, A nodes i λ ρ ; IH k e = e∈E J phys ID e ; ID k e + λ ρ J vol ((ρ m ) e ) + Λ Ξ J reg,e (R e ; R k e ) + Λ g J glob (A A e ; A A k e ) (51)

3 :
for all finite element e ∈ E do 4: for all unit cell types i ∈ I do minimizer among all α * i (i ∈ I)8:i * ← unit cell type index of α * 9:

1 . 6 :Fig. 7 :
Fig.7: Visualization of directional stiffness of unit cell with maximally opened fluid channels (r x = 0.22, r y = 0.22).This spherical plot was generated by drawing the entry A 1111 of the rotated material tensor A A ∈ S 6 for varying rotation angles (θ, φ) ∈ [0, 2π] 2 about z-and y-axes.For instance, the sketched arrow points to (π/2, 0) and its length of 1.9457 comes from first entry of the material tensor that is rotated by π/2 about the z-axis.

Fig. 9 :
Fig. 9: Optimization result for Λ Ψ = −10 and fixed local micro-structure orientation (no rotation) with Φ opt = 27.25 and Ψ opt = 0.275 for the optimized design in (a),(b).The initial guess is the design shown in fig.8.In (c) the mechanical state of the optimized design is visualized by deforming the domain by the physical displacements.The strain energy is shown in colors.In (e), the flow direction is visualized by equally scaled arrows and the colors indicate magnitude of the flow field.

Fig. 10 :
Fig. 10: Two optimized designs for different sample sizes of A grid 1 .(a) 10 samples each for r x and r y and 180 samples for ϕ.(b) 28 samples each for r x and r y and 180 samples for ϕ.Here, Λ Ψ = 1 and Λ Ψ = −10.The visual differences are barely perceptible, although (b) has a 1.5% lower compliance and a 1.7% higher flux than (a).

4. 1 . 1
Optimized local in-plane rotation of micro-structureWe introduce angle variable ϕ ∈ [0, π] to allow in-plane rotation, about the z-axis, of the microstructure.The effective material coefficients are rotated by ϕ with the following analytical expressions:A A rot (r x , r y , ϕ) = Q 6 (ϕ)A A(r x , r y )Q 6 (ϕ) T , B rot (r x , r y , ϕ) = Q 3 (ϕ)B(r x , r y )Q 3 (ϕ) T , K rot (r x , r y , ϕ) = Q 3 (ϕ)K(r x , r y )Q 3 (ϕ) T , (62)where Q 6 ∈ R 6×6 are rotation matrices for the stiffness tensor A A in Voigt notation and Q 3 ∈ R 3×3 are rotation matrices for the Biot coupling and permeability tensor.We note that no additional evaluation of the homogenization operators are required, as, instead of the micro-structure, the effective material tensors are rotated.ϕ is discretized with 180 steps for the brute force approach to solve the SGP subproblem with algorithm 2.Let us again set Λ Φ = 1 and Λ Ψ = −10, as in fig.9, and observe in figs.11a and 11b how the design evolves as both physical models counteract each other: the mechanical model strives for as much material as possible to minimize the compliance while the fluid flux is maximized when there is less material in the design domain.The convergence plot for the merit function J and

Fig.
Fig. Pareto front for varying Λ Ψ in weightedsum formulation F phys = Φ + Λ Ψ Ψ.The optimization was based on cells of type 1 and the initial design was always [0.15, 0.15, 0] n el .As we are minimizing Φ and maximizing Ψ, a point P = (P Φ , P ψ ) in the image space of Φ and Ψ is dominating a pointQ = (Q Φ , Q ψ ) if P Φ ≤ Q Φ and P Ψ ≥ Q Ψ .

Fig.
Fig. Comparison of Pareto curves for varying Λ Ψ .Blue: optimization with cells of type 1 and 2. Red: optimization with only cells of type 1.The blue curve clearly dominates the red curve.

Fig. 17 :
Fig. 17: Results of bicriterial optimization with cells from both type 1 and 2 for varying Λ Ψ .The designs visualized here corresponds to the labeled data points of the pareto curve in fig.16.

Fig. 18 :
Fig.18: Result of pure compliance minimization when allowing unit cells of type 1 and 2 with an active volume fraction constraint setting ρm = 0.8 on the stiff material phase.Comparing to fig.17a, it is observed that only now also cells of type 1 appear in the design.Moreover, the resource constraint leads to a variation of the parameter r s for cell type 2.
displays the changes in design with increasing regularization parameter p filt .The respective objective values are listed in table 2. The regularization of fluid channel radii can be observed well when comparing the designs in the right lower corner of fig.19b and fig.19c.With increasing p filt , the interface