A sequential mixed-integer programming method for concurrent optimization of core topology and face sheet thickness of a sandwich beam

A method is proposed that allows for the concurrent optimization of core topology and face sheet thickness of a sandwich beam under compliance constraints. The problem is solved using a novel mixed-linear extension of the Topology Optimization of Binary Structure (TOBS) topology optimization method aiming to minimize the total mass of the beam. The method has been demonstrated on a clamped beam example and the results have been compared to results from topology optimization of the core with a range of a priori fixed face sheet thicknesses. It is shown that the new method, starting from a fully populated core, finds a minimum mass that is lower than but in the neighbourhood of the best results from the topology optimization with fixed face sheet thicknesses. By varying the compliance constraint it is shown that the core topology approaches an ideal corrugated geometry as the compliance constraint is relaxed. The trends observed in the results are compared to analytical models for an idealized core.


Introduction
The concept of a sandwich structure has been used for more than 100 years to design lightweight components in a wide range of different engineering applications, e.g., aerospace, ground and sea transportation vehicles of all kinds, buildings, etc. Traditionally the core of a sandwich is a closed cell foam, a corrugated sheet or a honey comb, having the specific task of providing shear stiffness and to maintain the distance between the face sheets, which are providing the flexural stiffness. 1 With the rapidly growing interest in lightweight designs and the recent advances in structural topology optimization, not only the material selection and the sizing of the sandwich but also the actual topology of the core have come into focus in recent research. With a few exceptions, in most works the face sheets were fixed during the optimization of the core.
The work presented in this paper seeks to develop a method for concurrently optimizing the core topology and face sheet thickness of sandwich structures. The objective is to minimize the mass of a sandwich structure, assuming a set of given materials for face sheets and the core, without having to choose the face sheet thickness a priori.
Structural topology optimization has received extensive research interest since its introduction by Bendsoe and Kikuchi, 2 mainly thanks to its ability to generate components that fulfill specific constraints without the general layout of the component being known a priori. Following these initial works, a number of different optimization methods have been developed to solve topology optimization problems. In some of these, the design domain is discretized into finite elements and the optimization then modifies the properties of individual elements in order to improve the design. Most commonly, one wants to assign the value 0 or 1 to every element, where 0 means that the element is absent from the design and 1 means that the element is present. These methods can be divided into continuous density methods and discrete methods.
In continuous density methods, the integrality constraint of the design variables is relaxed, allowing the design variable to take intermediate values between 0 and 1. The mechanical properties of the element are then scaled according to the design variable, usually with a penalization of intermediate values in order to promote solutions with a 0-1 configuration. The most common of these methods is the Solid Isotropic Material with penalization (SIMP) 3 method, where usually the elastic modulus is penalized in order to make intermediate values of the design variable less efficient. The problem can then be solved using a non-linear gradient based optimization algorithm, such as the Method of Moving Asymptotes (MMA). 4 In discrete methods, elements are either removed or added to the design based on some criterion. The most researched of these methods are the Evolutionary Structural Optimization (ESO) 5 and the Bi-Directional Evolutionary Structural Optimization (BESO), 6 where each element is assigned a sensitivity number derived from the contribution of the element to the cost function and the constraints. The elements are then ranked according to their sensitivity number and removed or added based on their ranking.
Another branch of discrete methods are the sequential integer linear programming methods. 7 These methods employ more formal mathematical programming to determine what elements are removed or added, compared with the heuristic ranking of the evolutionary methods. The Topology Optimization of Binary Structures (TOBS) 8 is one such recent development, along with the canonical relaxation algorithm. 9 Other types of methods do not operate directly on the finite element grid. These include boundary variation methods, such as the level set method, 10 that modify the boundary between solid and void faces. Another type of method is the Method of Morphable Components (MMC), 11 which defines the geometry as a collection of structural members, such as trusses and beams, and then vary their position, length and size.
In terms of the chosen objective function, a number of different problems from multiple fields in engineering have been researched. Some examples are, the perhaps most commonly studied minimization of the mean or static compliance, 2 buckling 12 and maximum first principal stress. 13 Topology optimization in the design of sandwich structures have mostly been focused on the core itself. An early example of this may be found in Zhang and Sun,14 where material properties are distributed in a design domain to minimize the static compliance, wherafter topology optimization is used to design periodic microstructures to achieve these specific properties. Huang and Xie 15 studied a similar problem, where the core of a sandwich beam was topology optimized for minimum static compliance with different periodic constraints using the BESO method. The same method is also used in Zuo et al. 16 to design periodic sandwich cores with maximal eigenfrequencies and in Cameron et al. 17 to design a sandwich core with minimum mass and constraints on static deformation, eigenfrequency and buckling load. The SIMP method was used in Strek et al., 18 to optimize the topology of sandwich cores, composed of structural steel and an auxetic material, in order to minimize the static compliance. In Chu et al., 19 the MMC method was used to design sandwich panels with truss cores for minimum static compliance. Another example is found in the work of Sun et al., 20 where a hybrid core for a sandwich with carbon fiber face sheets was investigated. The core was composed of thin walls similar to a honeycomb core. By starting out with a large number of walls and then changing the thickness of the walls and deleting walls that became redundant the authors were able to minimize the weight of the panel subject to constraints on local buckling in the core when subjected to compressive loading, adhesive bonding area between the core and face sheets and the embedding volume of the core cells.
As mentioned in the beginning, usually only the core of the sandwich is considered during the optimization, in fact, to the best of the authors knowledge, there are only two examples of published works in which the face sheets have been included in the optimization. In Chu et al., 19 the face sheet thickness was part of the optimization in one of the numerical examples, and in Cameron et al. 17 a sequential approach was used, where first the core topology was optimized on its own with an a priori chosen face sheet and then a separate optimization of the face sheets was done.
To concurrently optimize the core topology and the face sheet thickness, a mixedinteger extension of the TOBS method is proposed in the current work. The paper is organised as follows: Methodology describes the methodology, Model problem and numerical results presents the model problem and the numerical experiments, Discussion discusses the results and Conclusions concludes the paper.

Problem formulation
The response of a structure to static loads are described by the finite element equation where K is the global stiffness matrix, f is the nodal load vector and u is the vector of nodal displacements.
The objective of the current work is to minimize the mass of a sandwich structure subject to constraints on the static compliance C, defined as The mass minimization problem for the sandwich structure may then be formulated as where x is the vector of core topology variables, t f is the face sheet thickness and C is the upper bound on the static compliance. Similarly m sw is the mass of the sandwich structure given by, where ρ f is the material density of the face sheet material, a f is the total face sheet area, m i is the mass of core element i when the element is solid and N e is the total amount of elements in the core.

Linearization
The current formulation presents an extension to the TOBS method, 8 here the continuous variable t f for the face sheet thickness has been added. This addition converts equation (3) into a nonlinear mixed-integer program that is quite difficult to solve. A linearization around the design of the k: th iteration, according to the principle proposed for the TOBS method, 8 results in the following mixed-integer linear program, where Δm sw , Δx k and Δt k f are the differences in mass, core topology variables and face sheet thickness in the k: th iteration. In addition, ΔC is the upper bound on the allowed change in the compliance.
Equation (5) is solved in every iteration and the design is then updated according to the equation

Relaxation of constraints
The constraint in equation (5) needs to be relaxed in order for the problem to have a feasible solution. If the current compliance is much lower than the upper bound, then ΔC is very large and this allows a large amount of material to be removed from the design, far outside the region where the linearization is valid. If on the other hand the compliance is much larger than the upper bound, ΔC is a large negative number, meaning that the compliance has to be significantly reduced in one iteration and there might not exist such a solution in the linearized model. The relaxation is done in analogy with the TOBS method as, where ε is some small positive number. In order to prevent too much core material from being removed in each iteration, a neighbourhood constraint needs to be included in the optimization. For the pure topology optimization problem this step is unnecessary since the relaxation naturally imposes an upper bound on the amount of material that may be removed in each iteration. However, including the face sheet thickness into the problem, may lead to designs where a slight increase in face sheet thickness lowers the compliance so much that, according to the linearized model, a large amount of core material can be removed without violating the relaxed constraint. This large removal of core material would be outside the region where the linearization is valid and would lead to a drastic shift in core topology, changing the nature of the problem. Thus, to prevent this, the following additional constraint is imposed in equation (5), where ϵ r and ϵ a are some small positive numbers and A solid and A void are defined by, and

Sensitivities
In order to formulate equation (5), the sensitivities of the mass and compliance with regards to the core topology variables and the face thickness are required. Starting out with the mass, the sensitivities are acquired by differentiation of equation (4), and The sensitivities of the compliance with regard to the topology variables result from differentiation of equations (2) and (1) and some rearranging of terms, here assuming that the load vector is independent of the topology, meaning that the corresponding derivatives vanish. For a more detailed derivation, see Bendsøe and Sigmund. 12 The sensitivity is given by, where u i is the element vector of nodal displacements and K 0,i is the element stiffness matrix, of element i. In order to evaluate equation (13), the element stiffness matrix is needed as a function of the design variables. In this paper, the interpolation model proposed by Huang 21 has been used, where E is the Young's modulus of the interpolated material, E s is the Young's modulus of the constituent solid material, x min is the value of the void elements, chosen such that the stiffness matrix remains non-singular even in the presence of void elements, and p is a penalty exponent.
Combining equations (13) and (14), the sensitivities can finally be expressed as, In order to prevent checker-boarding and mesh-dependent solutions, the computed topology sensitivity numbers are filtered according to the procedure described in Huang and Xie. 22 The filtered sensitivity numbers are also averaged with their history in order to stabilize the optimization process, according to the procedure in Huang and Xie. 6 The sensitivity of the compliance with regards to the face sheet thickness is also needed. To compute this, we start by differentiating equations (1) and (2) with regards to the thickness of a single face sheet element, where t f,j is the thickness of face sheet element j. It is assumed that the load is independent of the face sheet thickness and therefore the derivative vanishes. Solving equation (16) for ∂u ∂t f , j and inserting it into equation (17) yields As the face sheet thickness is assumed to be the same for all face sheet elements, and the derivative of the stiffness matrix are the entries related to each of the corresponding face sheet elements, equation (18) may be simplified to, where u j and K j are the element displacement vector and stiffness matrix of face sheet element j.
The sensitivity with regards to the face sheet thickness may then be acquired by summing over all face sheet elements, Finally, in order to evaluate the sensitivity, the derivatives of the element stiffness matrices are needed. Since the analysis was done using a commercial finite element solver, the stiffness matrix was not explicitly available as a function of the element thickness. To circumvent this, a surrogate model of the stiffness matrix was constructed. The element of a face sheet matrix was generated for a number of element thicknesses, and an expression on the form was fitted for each matrix entry, where b is the element entry to be evaluated and the a i are the coefficients to be computed for the surrogate model. This expression could then be differentiated, yielding,

Problem solution
In order to solve equation (5), an open-source solver called COIN-OR CBC 23 was used though the Python interface module cylp. This solver uses a branch-and-cut method, which first relaxes the integrality constraints in the problem and then solves this problem using the simplex method, giving a lower bound on the objective function. New constraints are then iteratively added to the relaxed problem to force the integer variables to take integer values at optimality. For details on the method, see the work of Conforti et al. 24 An important point is that if the difference between the best found integer solution and the relaxed problem is smaller than the mass of a single solid core element, the solving process is terminated and the current integer solution is accepted. The topology optimization is terminated once a solution is found that fulfills the constraints in equation (3), where the mass of the structure has converged, which is determined by the condition P k n¼kÀl m n À P kÀlÀ1 n¼kÀ2l m n P k n¼kÀl m n ≤ τ, where k is the latest iteration, l is the number of iterations for which the convergence is evaluated, m n is the structural mass at the n: th iteration and τ is the tolerance.

Model problem and numerical results
In this section the proposed optimization method is validated and employed to investigate the influence of the compliance constraints on the core topology and the resulting mass. First a model problem is introduced.

Modelling aspects
The proposed methodology was demonstrated on the clamped beam shown in Figure 1, the green part is the core and the red parts are the face sheets. The beam is 8 m long, has a core thickness of 20 cm and is 5 cm wide. A 2 MPa pressure was acting on a 40 cm wide section in the middle of the beam. Due to symmetry only the right half of the beam was modelled. The core of the beam was discretized into 800x40 4-node plane stress elements (Abaqus CPS4R), while each face sheet was discretized into 800 shear deformable beam elements (Abaqus B21) with a rectangular cross section having the same width as the core and a height according to the face sheet thickness.
The material properties of the face sheets and core are given in Table 1. The face sheet material properties are based on glass fiber reinforced plastic, while the core material properties are based on structural polymer foam.
For all optimizations, a filter radius of 4 cm was used, the convergence period was 15 iterations and the convergence tolerance was 0.001. Unless otherwise stated, the relaxation number ε j was 0.005, the maximum amount of elements allowed to be removed or added in each iteration was 160 (ϵ r = ϵ a = 0.005) and the optimization started from a fully populated core domain and a 7 mm face sheet thickness. For this starting configuration the compliance calculated from equation (2) is C full = 17,794 Nm.

Optimization of core with fixed face thickness
To verify the validity of the proposed method, a concurrent optimization run, assuming a relaxed compliance constraint of 28,000 Nm, was compared to a number of pure core topology optimizations. In these cases, the face sheet thickness was fixed at values between 4 and 9 mm. The relaxed constraint level was chosen to allow for the topology of the core to be modified. These core topology optimizations were performed using the TOBS method. The 28,000 Nm compliance constraint was chosen to be about 50% larger than the compliance of the initial structure, in order to give the algorithm some room to modify the structure. It should be noted that the compliance of the initial structure with a fully populated core domain was just above 28,000 Nm for the 4 mm face sheet thickness. This means that 4 mm is the minimum possible face sheet thickness, as the fully populated design domain is the core topology with the smallest possible compliance.
The results from this optimization may be seen below. In Table 2, the numerical results are given and the corresponding topologies of the core can be seen in Figure 2 for the different solutions. The mass variation for these fixed face thickness cases is shown in Figure 3.

Optimization of core topology and face thickness
With the proposed solution of equation (5) it is possible to compute solutions that correspond to an optimal balance between the face sheets and the core. Here two main cases will be presented, one keeping the same constraint level on the compliance as in the previous section and one, which consists of six different subcases in which the compliance constraint is successively relaxed.
Fixed compliance constraint. For the compliance constraint of 28,000 Nm, the resulting core topology is shown in Figure 4, and the face sheet thickness found as 5.17 mm, resulting in a total mass of 11.86 kg and a solid fraction of 0.55. These should be compared to the results in Figure 3, where the current result is indicated with a blue cross, and in Table 2.
Varying compliance constraint. To investigate how the optimal face sheet thickness and core topology is influenced by the constraint level, equation (5) was solved for a set of compliance constraints chosen as: 17,000, 20,000, 23,000, 25,000, 28,000 and 33,000 Nm. For each of these equation (5) was solved with the resulting core topologies shown in Figure 5 and numerical results in Table 3. Observe that the results from the 28,000 Nm constraint are repeated from the previous section for clarity.

Strain energy distribution
According to sandwich theory, 1 the deflection of a sandwich beam can be partitioned into shear and bending contributions. From equation (2) it then follows that the same partitioning can be made to the compliance. This allows for the distribution of strain energy between the core and face sheets for different compliance constraints to be studied. Assuming that the weak core and thin face sheet approximations hold also in the topologically optimized cores, the core carries all of the transverse force in the beam, while the face sheets carry all of the bending moment. Thus, the compliance contribution of shear and bending may be computed by summing up the strain energy in all of the core and all of the face sheet elements, respectively.
The results are shown in Figure 6 and 7, where the strain energies of the core and the face sheets, normalized with the total strain energy, are plotted against the compliance constraint. Interestingly, there seems to be an asymptotic behaviour as the compliance constraint is relaxed.

Mass distribution
A similar investigation may be made for the partitioning of mass between the face sheets and the core, see Figure 8 and 9. Also here there seems to be an asymptotic behaviour as the compliance constraint is relaxed.

Discussion
The numerical experiments presented in the previous section suggest that the concurrent optimization proposed in the current work is able to find a solution with low total mass. Obviously this could also have been achieved by running core topology optimization with a fixed face sheet thickness for a refined step. The final face sheet      thickness from the concurrent optimization is close to the 5 mm case, which had the lowest mass of all the solutions with fixed face sheet thicknesses.
As can be noted in Figure 5, the core structure seems to get thinner as either the face sheet thickness is increased or the compliance constraint is relaxed. In parallel the core structural members also tends to be inclined to the face sheets with an angle close to 45°. Furthermore the results in Figure 6 and 7 suggest that the relative contribution of shear and bending to the compliance seem to converge asymptotically as the compliance constraint is relaxed.
As shown in Appendix, assuming that deformation of the beam may be calculated as the sum of the bending and the shear deformations and that the beam is loaded by a concentrated point load in the middle, the compliance C ana may be written as, where u b and u s are the bending and shear deformations of the middle part of the beam, P is the magnitude of the point load applied in the middle of the beam and A and B are two constants. The current work is aimed at mass minimization, and the following compliance constrained mass minimization problem may now be formulated, where from equations (40) and (41), assuming that the core structural members have approximately a 45°angle to the face sheets and that the members only deform axially, the constants A and B are calculated as, and where L is the length of the beam. The optimal solution to this problem may be derived as, where m * f and m * c are the optimal mass of the face sheets and core. From this the ratio of the two masses is given as, which is independent of the magnitude of the compliance constraint itself. This would explain why the ratio of the face sheet and core mass converges asymptotically for the concurrent topology optimization as the compliance constraint is relaxed and the core topology converges to a core member configuration with 45°inclination with an optimal balance between the face sheets and the core.
Inserting the values of our problem from Figure 1 and Table 1 into equation (28) from, with p = 40,000 N, L = 8 m, ρ f = 1800 kg/m 3 , E f = 20 GPa, t c = 0.2 m, ρ c = 100 kg/m 3 , E c = 0.13 GPa and a compliance constraint C ¼ 66000 Nm (twice as large as in the topology optimization as the analytical model takes into account the entire beam and not just the right half), the resulting masses of the analytical solution can be seen in Table 4. The results from the concurrent topology optimization with 33,000 Nm constraint have been added for reference.
These results underline the proposed topology and face sheet thickness optimization method, indicating that the method is able to adequately balance the core member thickness against the face sheet thickness once it arrives at a core topology dominated by thin structural core members. The results are remarkably close to each other, recalling that the analytic solution is based on a point load and not a pressure acting in the middle of the beam as well as is assuming a highly simplified core.

Conclusions
The mixed-linear extension of the TOBS topology optimization method proposed in the current work finds a minimum mass of a clamped beam example, which is confirmed by a comparison to the results from doing topology optimization of the core with a priori fixed face sheet thickness. Varying the compliance constraint shows that the core topology seems to converge towards an idealized corrugated geometry as the constraint is successively relaxed. The strain energy in the face sheets, bending, and the core, shear, respectively, approaches an asymptotic partitioning. A similar observation may be made for the partitioning of the mass in between the face sheets and the core.