1 Introduction
In 1988 Bendsøe and Kikuchi published their seminal work on topology optimization Bendsoe1988 and seven years later, Sigmund applied topology optimization to design materials with prescribed elastic properties Sigmund1994 ; Sigmund1995 . Since then tailoring material properties with topology optimization has evolved into a large research area (see Ref. Osanov2016 for a review). One ongoing research topic is to reduce the high computational cost of topology optimization problems (e.g. Refs. Sigmund2013 ; Peetz2021 ). A large amount of this cost comes from the solution of the physical equilibrium, usually carried out with FiniteElementAnalysis. As an alternative, Fourierbased solvers show promise in becoming an efficient tool for the solution of the static mechanical homogenization problems Roters2012 . Fourierbased solvers have not yet been employed in the context of topology optimization, and we demonstrate their efficiency and versatility in the present article.
Fourierbased solvers for mechanical equilibrium trace back to the works of Moulinec and Suquet in the nineties Moulinec1994 ; Moulinec1998 . Since then, Fourierbased solvers have greatly evolved (see Ref. Schneider2021 for a current review) and have been successfully used for complex homogenization problems, e.g. Refs. deGeus2016 ; Lucarini2019
. Fourierbased solvers rely on efficient FastFourierTransformation implementations (like FFTW
Frigo2005 ) for computational efficiency. A variant of the method using compatibility projection was developed by Vondrejc, Zeman, de Geus et al. during the last decade Vondrejc2014 ; Zeman2017 ; deGeus2017 . These solvers work with some strain measure (e.g. the deformation gradient) rather than the displacements, as the unknown variable. The compatibility of this strain measure is enforced by projection, so that no additional compatibility equation must be solved. Here we will use the finiteelement projection recently proposed by us Leute2021 to avoid ringing problems, a phenomenon of artificial oscillations around discontinuities Gottlieb1997 persisting in many Fourierbased solvers (e.g. Schneider2016 ; Ma2021 ) and to enable modeling of internal free surfaces. For the topology optimization context, this means that void can be modelled with zero stiffness and does not need to be approximated by a very weak material.In this article we show how to use a Fourierbased solver with compatibility projection in gradientbased topology optimization. There are different approaches to formulate the topology optimization problem as a wellposed, differentiable optimization problem (e.g. Refs. Bendsoe2003 ; Sigmund2013 ). We will use the phasefield method proposed by Bourdin et al. in 2003 Bourdin2003 and since then successfully used in different optimization problems Wang2004 ; Burger2006 ; Wallin2012 ; Dondl2019 . To calculate the sensitivity, i.e. the gradient with respect to the design parameters, we will work with the discrete adjoint method. The adjoint method is the most efficient way to perform a sensitivity analysis of problems with many design variables (e.g. Ref. Tortorelli1994 ). It has therefore become a standard for FiniteElementMethod based topology optimization problems (e.g. Refs. Bendsoe2003 ; Osanov2016 ). Nevertheless, to the best of our knowledge, the adjoint method has not yet been formulated for a Fourierbased solver.
As exemplary topology optimization problem we have chosen two relatively simple, smallstrain 2D problems: First we optimize a composite material for a specific shear modulus; second we optimize an auxetic metamaterial with negative Poisson’s ratio.
2 Equilibrium solver with compatibility projection
In this section, we will briefly recapitulate the theory of a Fourierbased solver with compatibility projection.
We are working in a periodic unit cell of (undeformed) volume that contains our representative volume element. The aim is to solve the static mechanical equilibrium given by
(1) 
where is the nabla operator with respect to the undeformed positions , the operator is the dot product where here and in what follows we assume implicit summation over repeated indices (Einstein summation convention). is the first PiolaKirchhoff stress and the superscript
designates the tensor transpose. Here and in the following bold symbols represent second order tensors and arrows represent a vector (see appendix
B for our notation conventions). depends on the local deformation gradient through an appropriate constitutive expression. is defined by where is the deformed position.We now formulate the weak form of the equilibrium equation following the procedure outlined in Ref. deGeus2017 . Given a periodic test function , we rewrite Eq. (1) in the weak form,
(2) 
where we have applied the divergence theorem and used the periodic boundary conditions to eliminate the surface term. The operator designates the double contraction, . The operator denotes the outer product, i.e. . The tensor can be interpreted as a virtual deformation gradient , i.e. a periodic and compatible test function. By compatible we mean, that it is given by the gradient of a potential function.
Compatibility can be enforced by a convolution with the selfadjoint compatibility projection operator (see e.g. Refs. Vondrejc2014 ; Leute2021 ), so that we can write:
(3)  
where is the application of the projection operator to the right handside and a periodic test function, which no longer has to be compatible. The operation is a convolution in realspace but becomes a multiplication in Fourier space. The numerical methods behind this compatibility projection therefore employ fast Fourier transforms to accelerate the computation of (see Refs. deGeus2017 ; Zeman2017 ; Leute2021 ).
Eq. (3) is discretized on a regular grid using a Galerkin scheme, i.e. the test functions and the unknown are expressed with the same set of basis functions. This leads to the discretized equilibrium equation:
(4) 
where small greek letters refer to the element of the discretized domain. Note that in the discretized form, is a fourth order tensor and the double contraction yields , where we have omitted the greek element indices for brevity.
The projection operator projects every tensor field onto its compatible part. Suitable definitions employing Fourierderivatives for are given for finitestrain in Ref. deGeus2017 and for smallstrain in Ref. Zeman2017 . In Ref. Leute2021 , the formulations for both cases are extended for arbitrary discretization of the gradient operator. In all cases is blockdiagonal in Fourierspace, so that Eq. (4) is always evaluated in Fourierspace for the sake of numerical efficiency. Please note that whereas the specific form of varies, the form of Eq. (4) remains the same and the derivation in section 3 is independent of the specific form of that is used. The only properties of that are required are its selfadjointness and that it projects every tensor field unto its compatible part.
Equation (4) is usually nonlinear. It must therefore be solved iteratively, e.g. by Newton iteration (see Ref. deGeus2017 ):
(5)  
Here is the assembly of the lefttransposed, i.e. , of the local tangent stiffness tensor , defined by the linearization of the stresses around , , and denotes the step of the Newton iteration. The tangent stiffness tensor is blockdiagonal, so that the double dot product in Eq. (5) can be efficiently implemented. As shown in Refs. Vondrejc2014 ; Mishra2016 ; Zeman2017 , the increments are themselves compatible, which ensures the compatibility of the deformation gradient obtained by this iterative procedure. It is therefore not necessary to solve an equation for satisfying compatibility of the deformation gradient.
3 Discrete adjoint method
In this section, we will formulate the discrete adjoint method for topology optimization problems using a Fourierbased solver with compatibility projection.
The aim of the discrete adjoint method is to calculate the sensitivity (the gradient with respect to the design variable) of the optimization problem in the discretized form:
(6)  
subject to 
with as the aim function and as the number of equilibrium constraints to which the optimization problem is subjected, typically differing in boundary condition only. Here and in the following, a capital Greek letter index refers to the number of equilibrium constraints and a comma in the index serves to separate different kind of indices for readability. is the design variable, in our case the material density. Note that the problem given in Eq. (6) is in the discretized form.
3.1 Derivation of the adjoint method
For the derivation of the adjoint method, we follow Refs. Tortorelli1994 ; Giles2000 . We introduce a set of second order tensors as Lagrangian multipliers to fulfill the constraints. The optimization problem Eq. (6) is then equivalent to
(7) 
as long as the equilibrium constraints are satisfied separately. The superscript designates the complex conjugate of a variable. Since the design parameters are independent of each other, the definition of the sensitivity at element leads to:
(8)  
The most complicated terms in this equation are the total derivatives of the strains with respect to the material density. They are defined implicitly by the equilibrium constraints, i.e. they can not usually be calculated analytically. To avoid a costly numerical computation of these terms, we choose so that the expression in the parenthesis in front of these terms vanishes.
In other words, must fulfill the adjoint equations:
(9)  
For , the partial derivative of the stress with respect to the strain is 0. In the other cases, corresponds to the tangent stiffness tensor which is used in Eq. (5),
(10) 
where is the double contraction of two forth order tensors. We confine ourselves to realvalued aim functions. The complex conjugate of Eq. (10) is then:
(11) 
since the projection operator and the tangent stiffness tensor are selfadjoint.
projects to a compatible solution space. In consequence, we can set the noncompatible part of to zero without loss of generality, so that . By applying to the lefthand side of Eq. (11) we get the final adjoint equations:
(12) 
Once the Lagrangian multipliers are determined by Eqs. (12), the sensitivity can be calculated with the remaining terms of Eq. (8):
(13)  
where we have used the fact that the stress of an element does not directly depend on the material of an element . Note that the partial derivatives can easily be calculated analytically. With the selfadjointness of the projection operator, Eq. (13) can be reformulated as:
(14)  
In the second equality we have made use of the fact that are the solution of the adjoint Eqs. (12) and therefore compatible and realvalued. Note that the structur of Eq. (14) is similar to the corresponding sensitivity equations in FiniteElementAnalysis based topology optimization problems, see e.g. Bendsoe2003 .
3.2 Efficiency of the adjoint method
The most computationally costly operation in the sensitivity analysis with the discrete adjoint method is the solution of the adjoint Eqs. (12). The adjoint equation is structurally identical to the Newton iteration of the equilibrium solver, Eq. (5). The cost of solving one adjoint equation is therefore identical to one Newton iteration. This means the cost of the sensitivity analysis is typically smaller than the cost of solving a (nonlinear) equilibrium problem. An additional advantage of this form of the adjoint equations is that we can use the solver implemented for the Newtoniterations to compute the Lagrange multipliers.
4 Validation and application
In this section, we test our adjoint formulation on a 2D smallstrain example. We have implemented the method in the opensource Fourieraccelerated micromechanics solver
Spectre muspectre .4.1 The test cases
As test case, we optimize a (periodic) 2D unit cell for a target effective shear modulus in a smallstrain situation. As aim function we choose the least square difference between the average stress in the unit cell and the target stress , which is the Cauchystress in the hypothetical homogeneous unit cell with the target material. To ensure isotropy, we consider this error for three linearly independent load cases,
(15) 
where is the imposed average small strain and .
The topology optimization problem needs some form of regularization that picks specific geometries from the variety of geometries that minimize our aim function. We use the phasefield approach of Refs. Bourdin2003 ; Wallin2012 , that selects solutions with minimal interface area. The idea of the phasefield approach is to represent the material distribution by some continuous phasefield function that varies between and to represent two materials. The aim function penalizes interfaces through the expression
where () are the length of the unit cell in the first (second) dimension and is a weighting parameter controlling the width of the interface: A smaller penalizes values of between and , while allowing steeper gradients of , so that the width of the diffuse interface decreases as decreases. It can be shown that in the limit , this functional becomes proportional to the total area (or in twodimensions, length) of the interface Modica1977 ; Modica1987 . In our case, of the length of the unit cell in xdirection results in a good interface width.
In the topology optimization context, the doublewell potential serves to penalize intermediate terms of , i.e. the solution of the continuous optimization problem will converge towards the solution of the discrete optimization problem. The penalization of the gradient avoids the meshdependency of the typical topology optimization problem (e.g. Bendsoe2003 chapter 1.3). Note that, contrary to the ’Solid isotropic material with penalization’approach often used in topology optimization (e.g. Bendsoe2003 ), the phasefield approach does not require an active volume constraint to converge to discrete solutions. To keep the test cases as simple as possible, we have therefore refrained from adding a volume constraint.
The complete aim function is now:
(16)  
with a second weighting parameter that determines the effective penalty of maintaining an interface. To get sensible results, we have fixed with
the Youngs modulus of the solid phase (see below). Because of the relative smallness and simplicity of the optimization problem, we have refrained from a hyperparameter optimization (e.g.
Lynch2019 ) but simply used educated guesses.We use two different discretizations: A square unit cell discretized with a regular square grid or a hexagonal unit cell discretized with a regular hexagonal grid. In both cases, we use 31x31 grid points with two linear elements per pixel and linear shape functions. The discretized gradients are derived in detail in Leute2021 for the square grid and in A for the hexagonal grid. The material density is defined per grid point, so that both elements at one grid point always have the same material density.
The system consists of a hypothetical elastic, isotropic material (Youngs modulus and Poissons ratio ) and of void (Youngs modulus and Poissons ratio
) with a linear interpolation of the material properties in the interface, e.g.
. We want to emphasize that the proposed equilibrium solver can handle elements with zero stiffness, so that we do not need to represent the void by a very weak material. Both materials are modelled with the standard 2dimensional Hooke’s law. The partial derivatives and can easily be calculated analytically.We solve the equilibrium equation with the NewtonCG solver of Spectre. We optimize the aim function with a standard LBFGSB optimizer Nocedal2006 as implemented in scipy with 0 and 1 as bounds on the material density.
4.2 Validation of the sensitivity analysis
First, we validated the implemented adjoint method by a comparison with a finite difference calculation of the sensitivity. Figure 1 shows the norm of the difference between the finite difference calculation of the sensitivity and the sensitivity calculated with the discrete adjoint method for a square grid discretization. The difference is given for two different initial phases and two different target shear moduli. In each case, the norm decreases linearly with the finite difference , until very small values of are reached. The results for a hexagonal grid (not shown) look similar. We conclude that the finite difference calculation converges linearly towards the adjoint calculation until numerical errors get dominant, confirming our expressions and their implementation.
4.3 Optimization results
The results of the optimizations with the square grid are shown in Fig. 2, the results with the hexagonal grid in Fig. 3. Note that for comparability we have shifted the results of the different optimizations so that the hole is always in the center of the unit cell. We observe that the eight optimizations all lead to a single hole in the unit cell. For the sinus initial phase the hole is roughly square shaped in the case of the square grid (Fig. 2 b and c) and circular in the case of the hexagonal grid (Fig. 3 b and c). For both grids, the hole is smaller for a larger target shear modulus. For the random initial phase and a target shear modulus of 7/20, we get the same results as for the sinus initial phase for both grids (Fig. 2 f and Fig. 3 f). However, for a random initial phase and target shear modulus of 3/10 the hole is a quadrangle (Fig. 2 e) in the case of the square grid and triangular (Fig. 3 e) in the case of the hexagonal grid.
The results of the optimization for a sine wave as initial phase agree with our expectations: Larger target shear modulus result in smaller holes and circular holes arranged in a triagonal hexagonal grid are isotropic with respect to the elasticity tensor Fil1964 ; Hu2000 . Since the target homogenized elastic constants were isotropic, the rectangular form of the holes in the case of the square grid compensates the anisotropy introduced by the square simulation cell.
For , we get the same results for a random initial phase as for the sinus initial phase. For
, the optimization finds other minimums with a random initial phase. Those are probably local minima.
4.4 Auxetic metamaterials
As a second example we want to optimize for an auxetic structure, i.e. a structure with a negative Poisson’s ratio. The approach is the same as in section 4.1 with two changes: Firstly, we prescribe a Poisson’s ratio of in addition to the target shear modulus . Secondly, we optimize only for the load case from section 4.1. With educated guessing we fixed the weighting parameters as , all other parameters are the same as in section 4.1.
Fig. 4 represents the result of the topology optimization for a random initial phase and a square grid discretization. It consists of a pattern of alternating horizontal and vertical slits. A similar pattern has been proposed by Taylor et al. to achieve auxetic structures with low porosity Taylor2014 . We conclude that the proposed method can be used for the design of simple auxetic structures.
5 Conclusion
In this paper, we have developed the theory for the discrete adjoint method in the context of Fourierbased micromechanics with compatibility projection. We implemented and validated this method within Spectre muspectre . Our optimization examples with inclusions of zero stiffness demonstrate that this method can be used for the design of simple metamaterials. Furthermore, the concept of compatibility projection can be extended to construct preconditioners for standard finite element formulations. Our derived expressions are extremely simple and efficient to calculate within existing solution frameworks, fostering the way for the rational design of complex composites or metamaterials.
Acknowledgements
We thank W. Beck Andrews, Andrea Codrignani, Martin Ladecký, Ivana Pultarová, Antoine Sanner and Jan Zeman for enlightening discussion. We acknowledge funding by the Carl Zeiss Foundation (Research cluster “Interactive and Programmable Materials  IPROM”), the European Research Council (StG757343), the Deutsche Forschungsgemeinschaft (EXC 2193/1  390951807) and the Swiss National Science Foundation (Ambizione grant 174105).
Appendix A Discrete gradient for a hexagonal grid
In this appendix we derive the gradient for a hexagonal grid and linear shape functions. The unit cell of this grid is represented in Fig. 5
. The Cartesian coordinate system
and the element system are connected by:We use the linear shape functions:
where the superscript designates the number of the triangular element (see Fig. 5). The partial derivatives of a function follow directly as:
with the value of at the grid point , and the triangular element .
Appendix B Notation conventions
Vector  
Second order tensor  
Forth order tensor  
A at element  Small greek letter index  
A in equilibrium constraint  Capital greek letter index  
Tensor transpose  
Tensor transpose  
Lefthand transpose  
Outer product  
Dot product  
Double contraction  
Double contraction  
Double contraction  
Complex conjugate 
References
 (1) M. P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Computer Methods in Applied Mechanics and Engineering 71 (2) (1988) 197–224.
 (2) O. Sigmund, Materials with prescribed constitutive parameters: an inverse homogenization problem, International Journal of Solids and Structures 31 (17) (1994) 2313–2329.
 (3) O. Sigmund, Tailoring materials with prescribed elastic properties, Mechanics of Materials 20 (4) (1995) 351–368.
 (4) M. Osanov, J. K. Guest, Topology optimization for architected materials design, Annual Review of Materials Research 46 (2016) 211–233.
 (5) O. Sigmund, K. Maute, Topology optimization approaches, Structural and Multidisciplinary Optimization 48 (6) (2013) 1031–1055.
 (6) D. Peetz, A. Elbanna, On the use of multigrid preconditioners for topology optimization, Structural and Multidisciplinary Optimization 63 (2) (2021) 835–853.
 (7) F. Roters, P. Eisenlohr, C. Kords, D. Tjahjanto, M. Diehl, D. Raabe, DAMASK: the Düsseldorf Advanced MAterial Simulation Kit for studying crystal plasticity using an FE based or a spectral numerical solver, Procedia IUTAM 3 (2012) 3–10.
 (8) H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, Comptes Rendus de l’Académie des Sciences. Série II. Mécanique, Physique, Chimie, Astronomie. 318 (1994) 1417–1423.
 (9) H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Computer Methods in Applied Mechanics and Engineering 157 (12) (1998) 69–94.
 (10) M. Schneider, A review of nonlinear FFTbased computational homogenization methods, Acta Mechanica 323 (2021) 2051–2100.
 (11) T. W. J. de Geus, J. E. P. van Duuren, R. H. J. Peerlings, M. G. D. Geers, Fracture initiation in multiphase materials: A statistical characterization of microstructural damage sites, Materials Science and Engineering: A 673 (2016) 551–556.
 (12) S. Lucarini, J. Segurado, On the accuracy of spectral solvers for micromechanics based fatigue modeling, Computational Mechanics 63 (2) (2019) 365–382.
 (13) M. Frigo, S. G. Johnson, The design and implementation of FFTW3, P. IEEE 93 (2) (2005) 216–231.
 (14) J. Vondřejc, J. Zeman, I. Marek, An FFTbased Galerkin method for homogenization of periodic media, Computers & Mathematics with Applications 68 (3) (2014) 156–173.
 (15) J. Zeman, T. W. J. de Geus, J. Vondřejc, R. H. J. Peerlings, M. G. D. Geers, A finite element perspective on nonlinear FFTbased micromechanical simulations, International Journal for Numerical Methods in Engineering 111 (10) (2017) 903–926.
 (16) T. W. J. de Geus, J. Vondřejc, J. Zeman, R. H. J. Peerlings, M. G. D. Geers, Finite strain FFTbased nonlinear solvers made simple, Computer Methods in Applied Mechanics and Engineering 318 (2017) 412–430.
 (17) R. J. Leute, M. Ladecký, A. Falsafi, I. Jödicke, I. Pultarová, J. Zeman, T. Junge, L. Pastewka, Elimination of ringing artifacts by finiteelement projection in FFTbased homogenization, arXiv:2105.03297 (2021).
 (18) D. Gottlieb, C.W. Shu, On the Gibbs phenomenon and its resolution, SIAM Review 39 (4) (1997) 644–668.
 (19) M. Schneider, F. Ospald, M. Kabel, Computational homogenization of elasticity on a staggered grid, International Journal for Numerical Methods in Engineering 105 (9) (2016) 693–720.
 (20) X. Ma, M. Shakoor, D. Vasiukov, S. V. Lomov, C. H. Park, Numerical artifacts of fast fourier transform solvers for elastic problems of multiphase materials: their causes and reduction methods, Computational Mechanics 67 (6) (2021) 1661–1683.
 (21) M. P. Bendsoe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer Science & Business Media, 2003.
 (22) B. Bourdin, A. Chambolle, Designdependent loads in topology optimization, ESAIM: Control, Optimisation and Calculus of Variations 9 (2003) 19–48.
 (23) M. Y. Wang, S. Zhou, Phase field: a variational method for structural topology optimization, CMESComputer Modeling in Engineering and Sciences 6 (6) (2004) 547.
 (24) M. Burger, R. Stainko, Phasefield relaxation of topology optimization with local stress constraints, SIAM Journal on Control and Optimization 45 (4) (2006) 1447–1466.
 (25) M. Wallin, M. Ristinmaa, H. Askfelt, Optimal topologies derived from a phasefield method, Structural and Multidisciplinary Optimization 45 (2) (2012) 171–183.
 (26) P. Dondl, P. S. Poh, M. Rumpf, S. Simon, Simultaneous elastic shape optimization for a domain splitting in bone tissue engineering, Proceedings of the Royal Society A 475 (2227) (2019) 20180718.
 (27) D. A. Tortorelli, P. Michaleris, Design sensitivity analysis: overview and review, Inverse Problems in Engineering 1 (1) (1994) 71–105.
 (28) N. Mishra, J. Vondřejc, J. Zeman, A comparative study on lowmemory iterative solvers for FFTbased homogenization of periodic media, Journal of Computational Physics 321 (2016) 151–168.
 (29) M. B. Giles, N. A. Pierce, An introduction to the adjoint approach to design, Flow, Turbulence and Combustion 65 (34) (2000) 393–415.
 (30) https://gitlab.com/muspectre/muspectre.
 (31) L. Modica, S. Mortola, Il limite nella convergenza di una famiglia di funzionali ellittici, Bollettino dell’Unione Matematica Italiana A 14 (1977) 526–529.

(32)
L. Modica, The gradient theory of phase transitions and the minimal interface criterion, Archive for Rational Mechanics and Analysis 98 (2) (1987) 123–142.

(33)
M. E. Lynch, S. Sarkar, K. Maute, Machine learning to aid tuning of numerical parameters in topology optimization, Journal of Mechanical Design 141 (11) (2019).
 (34) J. Nocedal, S. J. Wright, Numerical optimization, 2nd Edition, Springer series in operations research, Springer, New York, 2006.
 (35) L. Fil’Shtinskii, Stresses and displacements in an elastic sheet weakened by a doublyperiodic set of equal circular holes, Journal of Applied Mathematics and Mechanics 28 (3) (1964) 530–543.
 (36) N. Hu, B. Wang, G. Tan, Z. Yao, W. Yuan, Effective elastic properties of 2d solids with circular holes: numerical simulations, Composites Science and Technology 60 (9) (2000) 1811–1823.
 (37) M. Taylor, L. Francesconi, M. Gerendás, A. Shanian, C. Carson, K. Bertoldi, Low porosity metallic periodic structures with negative poisson’s ratio, Advanced Materials 26 (15) (2014) 2365–2370.
Comments
There are no comments yet.