Narindra Ranaivomiarana
email: narindra.ranaivomiarana@onera.fr

Franc ¸ois-Xavier Irisarri
email: francois-xavier.irisarri@onera.fr

Dimitri Bettebghor
email: dimitri.bettebghor@onera.fr

Boris Desmorat
email: boris.desmorat@upmc.fr

Concurrent optimization of material spatial distribution and material anisotropy repartition for two dimensional structures

Keywords: Topology optimization, SIMP, Distributed orthotropy, Polar method, Material design

Concurrent optimization of material spatial distribution and material anisotropy repartition for two dimensional structures

Narindra Ranaivomiarana, François-Xavier Irisarri, Dimitri Bettebghor, Boris

Introduction

Weight and cost reductions of structures are permanent challenges for aeronautic industries. As a response to this concern, structural optimization is on the rise [START_REF] Haftka | Elements of Structural Optimization[END_REF], especially topology optimization. It consists of determining the best shape of a given structural component or the best layout of structures. An objective function (e.g. total mass, global stiffness) is considered and constraints are imposed. The problem of obtaining a targeted distribution of stiffness (or density) is today a crucial one in structural engineering, and in some cases the research on metamaterials addresses precisely this problem by means of a suitable arrangement of the microstructure (see e.g. [START_REF] Berrehili | The homogenized behavior of unidirectional fiber-reinforced composite materials in the case of debonded fibers[END_REF][START_REF] Boutin | Linear pantographic sheets: Asymptotic micro-macro models identification[END_REF][START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF][START_REF] Giorgio | Buckling modes in pantographic lattices[END_REF][START_REF] Melnik | Dynamic fiber reorientation in a fiber-reinforced hyperelastic material[END_REF][START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF]; a reference on fibrous structures is [START_REF] Dell'isola | Synthesis of fibrous complex structures: designing microstructure to deliver targeted macroscale response[END_REF]). For structures of isotropic material, methodologies in the literature are abundant: Sigmund et al. [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF] and Deaton et al. [START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF] survey developments of topology optimization. Interesting recent results are in [START_REF] Jia | Optimal structural topology of materials with micro-scale tension-compression asymmetry simulated using granular micromechanics[END_REF][START_REF] Lekszycki | A comparison between different approaches for modelling media with viscoelastic properties via optimization analyses[END_REF], and in particular in [START_REF] Andreaus | Optimal bone density distributions: numerical analysis of the osteocyte spatial influence in bone remodeling[END_REF][START_REF] Lekszycki | Functional adaptation of bone as an optimal control problem[END_REF] for applications to biomechanics. In order to reliably perform a topology optimization procedure, it is of crucial importance to have an insight into the effect of constitutive parameters onto the overall behaviour of the material. This can be achieved by means of suitable identification techniques, either experimental or numerical [START_REF] Han | Identification of elastic constants of alloys with sheet and fibre textures based on resonance measurements and finite element analysis[END_REF][START_REF] Placidi | Identification of Two-Dimensional Pantographic Structures with a Linear D4 Orthotropic Second Gradient Elastic Model Accounting for External Bulk Double Forces[END_REF]. Two main methods are used to parameterize the topology of the structure: boundary variation method and density based method. For boundary variation method, two main techniques are used in the literature: the level set method [START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF][START_REF] Allaire | Conception optimale des structures[END_REF][START_REF] Burger | A survey on level set methods for inverse problems and optimal design[END_REF] that uses a scalar function to represent the boundaries, and the phase-field method [START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF][START_REF] Bourdin | The phase-field method in optimal design[END_REF] that uses a phase field function over a domain composed of two phases. Density based methods use a density variable to define in each element of the structure if there is solid material [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] or void (0). The discrete problem is turned into a continuous one, conducting to intermediate densities difficult to interpret. Many penalization methods exist to steer the solution towards a discrete solid/void solution, among which: SIMP (Solid Isotropic Material Penalization) [START_REF] Bendsoe | Optimal shape design as a material distribution problem[END_REF], RAMP (Rational Approximation of Material Properties) [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF] and SINH [START_REF] Bruns | A reevaluation of the SIMP method with filtering and an alternative formulation for solidvoid topology optimization[END_REF]. Another technique of penalizing was proposed by Fuchs et al. using reciprocal variables (SRV) [START_REF] Fuchs | The SRV constraint for 0/1 topological design[END_REF]. Compared with density based methods, boundary variation methods give clear contour of the shape and give easy definition of geometrical parameters. Boundary variation methods do not require penalization method, but they need computation of topological gradients and construction of interpolation schemes that are complex [START_REF] Norato | A topological derivative method for topology optimization[END_REF]. In the present work, a density based method is used and the SIMP method is employed as penalization method. The SIMP method can be applied to a broad range of applications and physical models. A topology optimization technique was developed to solve thermal system in [START_REF] Gersborg-Hansen | Topology optimization of heat conduction problems using the finite volume method[END_REF]. Desmorat [START_REF] Desmorat | Structural rigidity optimization with an initial design dependent stress field. Application to thermo-elastic stress loads[END_REF] used the SIMP method for optimization with thermo-elastic stress loads. Elastic structures with frictionless unilateral contact have also been optimized in [START_REF] Desmorat | Structural rigidity optimization with frictionless unilateral contact[END_REF]. As this method reaches maturity, it is implemented in various commercial tools in the field of mechanics and multiphysics (e.g. fluid structure interaction). All the same, topology optimization methods are performed on a conceptual design level. The optimization methodology gives a first approach of the final manufacturable designed structure and a further step of reinterpreting the solution into a CAO model is needed. Afterwards, a parametric optimization is necessary for a more detailed design. An alternative solution for minimizing mass is the use of lighter materials. Composite structures have become popular in aerospace industries as they are lighter than metallic ones. Actually, composite materials have higher stiffness-to-weight and strength-to-weight ratios than metallic materials. Therefore, great research effort was made to optimize their use as their anisotropic behavior is different from that of metallic materials with isotropic behavior (a useful reference on the anisotropic behavior in the case of microstructured/micromorphic continua is [START_REF] Eremeyev | Material symmetry group and constitutive equations of micropolar anisotropic elastic solids[END_REF]). In particular, Ghiasi et al. reviewed the optimum stacking sequence of laminates by designing the anisotropy distribution with constant stiffness [START_REF] Ghiasi | Optimum stacking sequence design of composite materials part I: Constant stiffness design[END_REF] and variable stiffness [START_REF] Ghiasi | Optimum stacking sequence design of composite materials part II: Variable stiffness design[END_REF]. The main characteristic of the composite optimization is the design of anisotropy distribution. Composite optimization use two main methods to design the material anisotropy: direct approach and multilevel method. The direct approach consists of performing parametric optimization on stacking sequence design variables. Stacking sequence and orientation of the fibers are optimized in [START_REF] Irisarri | Multiobjective stacking sequence optimization for laminated composite structures[END_REF] using heuristic methods. Optimal thickness of laminated composites are sought in [START_REF] Sorensen | Topology and thickness optimization of laminated composites including manufacturing constraints[END_REF] using gradient-based methods. Nevertheless the direct approach is difficult to implement and requires high number of variables (varying with the number of plies). Furthermore, its main drawback is the lack of regularity in the mechanic response with respect to the local orientations in the material. In the multilevel method, first, the material anisotropy is optimized by considering an homogenized stiffness of the material, and second, a physical construction of the considered laminate is sought with the optimal design variables obtained in the first step (laminate retrieval) [START_REF] Ijsselmuiden | Optimal design of variable stiffness composite structures using lamination parameters[END_REF][START_REF] Bordogna | Aeroelastic optimization of variable stiffness composite wing with blending constraints[END_REF][START_REF] Meddaikar | Laminate optimization of blended composite structures using a modified Shepard's method and stacking sequence tables[END_REF]. The interest of this method is the possibility of taking advantage on the regularity of structural responses by using the homogenized stiffness. To parameterize the homogenized stiffness laminated composite plates, lamination parameters are commonly used [START_REF] Miki | Optimum design of laminated composite plates using lamination parameters[END_REF]. Jibawy et al. [START_REF] Jibawy | Hierarchical structural optimization of laminated plates using polar representation[END_REF] use polar invariants by change of frame to define the anisotropy of homogenized laminated plates. Examples of works using polar method to parameterize laminates can be found also in [START_REF] Julien | Conception Optimale de l'Anisotropie dans les Structures Stratifiées á Rigidité Variable par la Méthode[END_REF][START_REF] Vincenti | Optimal Orthotropy for Minimum Elastic Energy by the Polar Method[END_REF][START_REF] Catapano | Stiffness and strength optimization of the anisotropy distribution for laminated structures[END_REF]. Lamination parameters define only specific composite laminates. Rather, the polar method is used in this work to determine the material anisotropy as the polar invariants can parameterize a general anisotropic homogenized material. Composite optimization is generally done with a predefined shape of the structure, most of the time derived from the preexisting metallic part. Nevertheless, shapes of metallic structures are not appropriate for anisotropic composite structures. Indeed, the material anisotropy influences the optimal shape. Peeters et al. [START_REF] Peeters | Combining topology and lamination parameter optimisation[END_REF] obtained two different shapes in the case of a quasi-isotropic (QI) structure and in the case of a composite laminate structure obtained by a combination of topology optimization and composite optimization with lamination parameters. In the present study, results from initial work show systematically the influence of the orthotropic material orientations on the optimal shapes. Topology optimizations of an unidirectional composite made of long and straigth carbon fibers in an epoxy matrix were carried out for material orientations varying from 0 • to 180 • in steps of 1 • . A cantilever beam was optimized by minimizing its compliance and by imposing a volume constraint (50% of the total volume amount). The left side of the structure was clamped and a vertical load was applied in the middle of the right side. The optimizations give different optimal shapes with respect to the material orientation. Five different optimal shapes are shown in Figure 1 with orthotropic materials oriented in five different directions (0 • , 18 • , 45 • , 90 • ,-45 • ). Compared with the 0 • -case, the optimal shape in the 90 • -case has less material near the point load as the material is oriented as well as the load direction. Optimal shape in the 45 • -case is symmetric to the optimal shape in the -45 • -case: the orientations are symmetric beside x-axis. The compliance evolution in a polar representation is also shown in the figure with respect to the material orientation: the material oriented at 0 • has the lowest value and is therefore the stiffest. These results show the necessity of seeking the optimal shape and the optimal distribution of anisotropy of the considered structure all at once. In the literature, only few works seek simultaneously the optimal shape and the optimal anisotropy distribution. Irisarri et al. [START_REF] Irisarri | Optimisation of ply drop order in variable stiffness laminates[END_REF] used a shell structure with a predefined and fixed geometry to optimize a variable stiffness in the structure. The material anisotropy distribution was optimized on a wind turbine blade in [START_REF] Lund | Buckling topology optimization of laminated multi-material composite shell structures[END_REF]. Jibawy et al. [START_REF] Jibawy | Hierarchical structural optimization of laminated plates using polar representation[END_REF] sought the optimal distribution of anisotropy of a hole plate with fixed geometry. However, these cited methods seek the ply shape of a thin structure with a predefined and known shape of the shell (this last one does not change). Rather, Peeters et al. [START_REF] Peeters | Combining topology and lamination parameter optimisation[END_REF] simultaneously optimized the topology and the fiber angle distribution of two dimensional composite structures. In this last method, the optimal shape is sought at the same time as the optimal anisotropy distribution. The novelty of this work is the concurrent optimization of the topology and the material anisotropy of a two dimensional structure with general homogenized orthotropic material. The SIMP method is combined with the polar method to find the best spatial material distribution and material anisotropy repartition. A general orthotropic material is considered in this work. Polar method allows for representing any two dimensional fourth order elasticity tensor with invariants by change of frame [START_REF] Verchery | Les invariants des tenseurs dordre 4 du type de l'élasticité[END_REF]. In this way, thermodynamic bounds are established on the polar invariants to work with general orthotropic materials. Nonetheless, the method can be applied on laminates as their elastic bounds are known [START_REF] Vannucci | A Note on the Elastic and Geometric Bounds for Composite Laminates[END_REF]. The remainder of the paper is organized as follows. The problem of optimizing concurrently the material spatial distribution and the material anisotropy repartition is formulated in Section 2. The alternate directions algorithm used to solve iteratively and successfully local minimizations and global minimizations is presented in Section 3. An optimization of a cantilever beam is performed in Section 4 giving its optimal shape and optimal material anisotropy distribution. In Section 5, a comparison between sequential and concurrent optimization is made. Sequential optimization is made up of a topology optimization with a fixed isotropic material followed by an optimization of material anisotropy distribution with the obtained shape.

2 Problem formulation: simultaneous optimization of the material density and anisotropy

Parameterization of the distributed material density and anisotropy

Two main methods are used to parameterize the topology of the structure: boundary variation method (level set method [48,?,?], phase field method [START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF][START_REF] Bourdin | The phase-field method in optimal design[END_REF]) and density based method. Comparing to density based methods, boundary variation methods give clear contour of the shape and the definition of geometrical parameters is easy. However, the computation of topological gradients and the construction of interpolation schemes are complex. Thus, a density based method is used in this work. The shape of the structure is parameterized by a density field variable ρ(x). This density variable defines at each point x of the structure whether there is material (ρ(x) = 1) or void (ρ(x) = 0). ρ(x) takes any value in [ρ min , 1] while, in order to avoid singularity of elasticity tensors, the lowest admissible value ρ min is generally set to 10 -3 . Allowing ρ(x) to be valued in interval [ρ min , 1] involves intermediate densities appearing in the optimum topologies. Grey areas are difficult to interpret because they correspond to a mixture of void and material. To suppress grey areas, the so-called SIMP method [START_REF] Bendsoe | Optimal shape design as a material distribution problem[END_REF] is used. This method penalizes ρ(x) to be either 0 or 1 by using an exponant p. Optimized stiffness tensor C(x) and compliance tensor S(x) are expressed in function of the considered material stiffness tensor C 0 (x) and compliance tensor S 0 (x) (Eq. ( 1))

C(x) = ρ(x) p C 0 ⇔ S(x) = 1 ρ(x) p S 0 (1)
The anisotropic material can be parameterized by several methods. Allaire et al. [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] use homogenization method to parameterize sequential laminates made of a mixture of isotropic material and void. Peeters et al. use the lamination parameters to parameterize a laminated composite. Laminated structures are parameterized by lamination parameters in [START_REF] Peeters | Combining topology and lamination parameter optimisation[END_REF] and by polar invariants by change of frame in [START_REF] Jibawy | Hierarchical structural optimization of laminated plates using polar representation[END_REF][START_REF] Catapano | Stiffness and strength optimization of the anisotropy distribution for laminated structures[END_REF]. In the present work, the material is choosen to be general orthotropic (thermodynamically admissible). The stiffness properties of the considered anisotropic material is then represented by the elasticity tensor with homogenized coefficients. Spatial variations of the material anisotropy are allowed. A parameterization that allows to express the elasticity tensor in a general frame in a simple way is necessary.

Change of frame is cumbersome using the Cartesian representation. The polar method permits to write the elasticity tensor with its intrinsic properties using tensor invariants. By doing so, changing frame becomes simple as only a rotation of an angle with respect to the frame is needed. Thus, the polar invariants are chosen to express the stiffness tensor for an orthotropic material under assumption of plane stress. Equations [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF] shows the relation between the polar components (T 0 , T 1 , R 0 , R 1 , Φ 0 and Φ 1 ) and the Cartesian ones of the reduced stiffness tensor Q [START_REF] Julien | Conception Optimale de l'Anisotropie dans les Structures Stratifiées á Rigidité Variable par la Méthode[END_REF][START_REF] Vincenti | Optimal Orthotropy for Minimum Elastic Energy by the Polar Method[END_REF].

Q 1111 = T 0 + 2T 1 + R 0 cos4Φ 0 + 4R 1 cos2Φ 1 Q 1122 = -T 0 + 2T 1 -R 0 cos4Φ 0 Q 1112 = R 0 sin4Φ 0 + 2R 1 sin2Φ 1 Q 2222 = T 0 + 2T 1 + R 0 cos4Φ 0 -4R 1 cos2Φ 1 Q 2212 = -R 0 sin4Φ 0 + 2R 1 sin2Φ 1 Q 1212 = T 0 -R 0 cosΦ 0 (2)
Each Cartesian component of the reduced stiffness tensor is expressed with isotropic terms (T 0 , T 1 ) that do not depend on the orientation of the material, and anisotropic terms R 0 e 4iΦ 0 , R 1 e 2iΦ 1 that depend on the orientations Φ 0 , Φ 1 of the material. The polar invariants are the moduli T 0 , T 1 , R 0 , R 1 and the angle Φ 0 -Φ 1 . The isotropic parts do not influence the anisotropy of the material, thus (T 0 , T 1 ) are supposed to remain constant (in composite laminated plates made of identical unidirectional layers (with the same material and same thickness in each layer), the homogenized isotropic part (T 0 , T 1 ) of the laminate is equal to the isotropic part (T EL 0 , T EL 1 ) of the elementary layer (EL) [START_REF] Jibawy | Hierarchical structural optimization of laminated plates using polar representation[END_REF]). The material optimization is performed with respect to the anisotropic parts (R 0 , R 1 , Φ 1 ).

Optimization constraint: maximum volume and thermodynamical admissibility of the material

The optimization constraints are written in terms of the total volume amount of the structure and of the anisotropic part of the polar invariants by expressing their bounds. During the optimization, a target volume V 0 is defined for the structure. The volume V is equal to the material density ρ(x) integrated in the domain Ω . At each step of the optimization, the volume must satisfy the optimization constraint Eq. ( 3):

V = Ω ρ(x)dx = V 0 . (3)
The material to be designed is imposed to be orthotropic. For an orthotropic material:

Φ 0 -Φ 1 = K Π 4 with K = 0, 1. (4)
The orthotropic material used in this paper is taken to be as general as possible: the optimized orthotropic material is thermodynamically admissible, which means that the stiffness tensor is positive definite. The analytical bounds of the polar invariants are given in Eq. (3) [START_REF] Vannucci | Plane Anisotropy by the Polar Method[END_REF]:

           T 0 > 0 T 1 > 0 T 0 > R 0 T 0 T 1 > R 2 1 T 1 (T 2 0 -R 2 0 ) > 2R 2 1 (T 0 -R 0 cos4(Φ 0 -Φ 1 )) (5)

Double minimization of the complementary energy

In topology optimization, criteria such as buckling, frequency or compliance may be considered, for instance [START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF]. In this paper, we aim at maximizing the global structural stiffness measured by the compliance which is the external work. The criterion is written as:

Criterion = Ω f.udV + Γ 1 F.udS. (6)
The boundary of domain Ω is split into two parts: Γ 0 where a zero displacement is imposed, and Γ 1 where a surface load F is applied. f is the volume load and u, the displacement vector. The more the structure is rigid, the lower is the external work. Thus, maximizing the global structural stiffness means minimizing the compliance. Moreover, the compliance is equal to the double of the complementary energy. The optimization is made with respect to the density and the anisotropic part of the stiffness tensor polar invariants:

min {ρ,R 0 ,R 1 ,Φ 1 } Ω f.udV + Γ 1 F.udS = min {ρ,R 0 ,R 1 ,Φ 1 } Ω σ : C -1 : σ dV. (7)
The complementary energy can be written as the minimization of a positive quantity with respect to the statically admissible stress field τ which, for imposed displacements equal to zero, reads

Ω σ : C -1 : σ dV = min τSA Ω τ : C -1 : τdV. (8)
The stress field τ satisfies the elasticity problem (P), with assumption of small strains and small displacements:

       divτ + f = 0 in Ω τ.n = F on Γ 1 τ = C : ε(u) in Ω u = 0 on Γ 0 (P)
where ε(u) = 1 2 (∇ u + ∇u T ) is the strain tensor. By replacing the expression of the complementary in Eq [START_REF] Bendsoe | Optimal shape design as a material distribution problem[END_REF], the optimization problem is written as a double minimization with respect to the design variables {ρ, R 0 , R 1 , Φ 1 } and to the stress field τ. The density variable is subject to a maximal volume constraint and the polar invariants of the stiffness tensor are constrained by thermodynamic bounds:

min {ρ,R 0 ,R 1 ,Φ 1 } min τSA Ω τ : C -1 : τdV with            Ω ρ(x)dx = V 0 T 0 > R 0 T 0 T 1 > R 2 1 T 1 (T 2 0 -R 2 0 ) > 2R 2 1 (T 0 -R 0 cos 4(Φ 0 -Φ 1 )) Φ 0 -Φ 1 = K Π 4 (K = 0, 1) (9)
3 Complementary energy minimization using the alternate direction algorithm

Local minimizations of the complementary energy

Since the design variables {ρ, R 0 , R 1 , Φ 1 } are subjected only to algebraic constraints, the minimization with respect to them can be put inside the integral:

min τSA Ω min {ρ,R 0 ,R 1 ,Φ 1 } τ : C -1 : τdV (10)
The minimization of the complementary energy with respect to the design variables is solved locally in each point of the domain, for a fixed stress state. Since the density variable ρ(x) and the anisotropy variables {R 0 , R 1 , Φ 1 } are independent, the minimization is split in two steps. First the complementary energy is minimized with respect to the anisotropy variables. Second, the minimization with respect to the density variable is performed.

The complementary energy can be written as a simple function of the polar invariants of the stiffness tensor and the stress tensor. Calculating its derivative is then straightforward. Hence, the minimization of the complementary energy with respect to the anisotropy variables is done analytically. The optimal values of {R 0 , R 1 , Φ 1 } depending on the stress field are given in [START_REF] Julien | Conception Optimale de l'Anisotropie dans les Structures Stratifiées á Rigidité Variable par la Méthode[END_REF]. The volume constraint is taken into account in the minimization step with respect to the density variable ρ(x) through the introduction of a Lagrangian multiplier k:

min ρ 1 ρ(x) p τ : C -1 (R opt 0 , R opt 1 , Φ opt 1 ) : τ + kρ(x) (11)
The minimum of the local energy is attained by deriving the equation Eq [START_REF] Bourdin | The phase-field method in optimal design[END_REF] with respect to ρ.

ρ(x) = p τ : C -1 (R opt 0 , R opt 1 , Φ opt 1 ) : τ k 1 p+1 (12)
The value of k is calculated so that the volume constraint is satisfied.

Optimization Algorithm

Convex approximation methods such as MMA and GCMMA, (Globaly convergent) Method of Moving Asymptot [START_REF] Svanberg | The method of moving asymptotesa new method for structural optimization[END_REF][START_REF] Zillober | A globally convergent version of the method of moving asymptotes[END_REF] and descent algorithm method such as SQP, Sequential Quadratic Programming [START_REF] Arora | Structural optimization by mathematical programming methods[END_REF][START_REF] Schittkowski | Software for Mathematical Programming[END_REF], IPOPT [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] need the evaluation of the objective function as well as its gradient. Optimality criteria method [START_REF] Prager | Introduction to Structural Optimization[END_REF][START_REF] Kirsch | Optimum structural design : concepts, methods, and applications[END_REF][START_REF] Bendsoe | Optimization of Structural Topology, Shape, and Material[END_REF][START_REF] Rozvany | Layout optimization of structures[END_REF] computes the optimal values of design variables by expressing optimality conditions. Therefore, optimality criteria method is less expensive than the methods above in terms of numerical cost. This is reason why a method similar to optimality criteria is used in this work. The algorithm used to solve the numerical problem is an extension of the alternate directions algorithm [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF]. The advantage of the alternate direction algorithm is its simplicity and low numerical cost as the method iterates between local minimizations and finite element calculations of stresses. Moreover, the alternate direction algorithm is convergent, and it can take into account a large number of variables. The double minimization is solved with a fixed point method by considering the optimality conditions. At each iteration of the optimization, the minimization with respect to the design variables {ρ, R 0 , R 1 , Φ 1 } is first performed, then the minimization with respect to the stress field τ is operated. The minimizations are treated alternatively and separately. Thanks to the polar method, the local complementary energy is written with simple expressions. Hence, the local minimizations are solved analytically. The minimization with respect to the stress field τ corresponds to a finite element analysis (FEA). The cost of the algorithm is directly related to the finite element calculation cost. In addition to advantages cited earlier, the algorithm can manage an optimization only with respect to the density, only with respect to the anisotropy (polar parameters) or both all at once: the procedure changes in the local minimizations solved with explicit formulas. Figure 2 shows the flowchart of the algorithm:

1. The model and the design variables {ρ, R 0 , R 1 , Φ 1 } are initialized and a FEA is performed to compute the corresponding stress field. 2. Local minimizations with respect to the design variables are solved analytically, where the stress field is remained fixed. This step takes into account the volume constraint. Optimal polar parameters and density are obtained. 3. Global minimization is carried out, where the design variables are remained fixed. This step is a FEA to update the stress field corresponding to the obtained optimal design variables computed in the previous step. 4. A stop criterion is defined: if the objective function becomes stationary and the change in design variables is smaller than a preset threshold, convergence is achieved. If the stop criterion is satisfied, optimal parameters are found, otherwise, new local minimizations are performed (step 2). Numerical results are presented in this section to prove the efficiency of the method. The optimization is made for a two dimensional orthotropic linear elastic material. A Cantilever beam is optimized. The design domain is a rectangle clamped at the left side and loaded at the middle of the right side (Figure 3). The domain size is 40 mm x 20 mm discretized with a rectangular 80 x 40 mesh. The volume constraint is fixed at 50% of the total volume. The initial density is set to 1 in every element of the mesh. The initial material is an isotropic material where the values of T 0 and T 1 correspond to the isotropic part of a monolayered composite made of long and straigth carbon fibers in an epoxy matrix: T 0 = 26.88 GPa, T 1 = 24.74 GPa.

Strategy of penalization and convergence

The compliance and the volume are displayed in function of the iterations in Figure 4. The strategy of penalyzing the density is made in two steps where the exponent p in ρ(x) p is increased gradually. The exponent p is generally taken to be equal to 3 to get a 0/1 layout. A refinement is made afterwards where p is increased to 5 to suppress definitely intermediate density. Excepted for the first iteration, the volume does not change through the iterations because its total amount is imposed to be 50 % of the total feasible volume. The compliance decreases in each step (p = 3, p = 5). Global termination of each step is set to 0.1% on the compliance and to 0.01% on the density. The compliance increases when the value of p is increased as the structure becomes less stiff. In comparison with the initial isotropic cantilever with 100% volume amount, the optimal orthotropic one with 50% volume amount is 40% stiffer as the anisotropic part of this last one is

Optimal material spatial distribution and optimal orthotropy distribution

Figure 5 shows the optimal shape of the structure where black represents the presence of material and white its absence. The material is pictured when the density value is above 0.8. To avoid numerical problems such as checkerboard and mesh dependency, a filter is used: the density of an element depends on the density of its neighbors so that there is no brutal discontinuity on the density in the structure. The neighbour elements that influence the considered element are defined by a radius filter which can be interpreted as a bar width. In the optimal shape, the orthotropy is distributed: the material orthotropy changes continuously inside the structure. The optimal orthotropy direction Φ opt 1 is presented in Figure 6. It is aligned with the direction of the principal stress (highest absolute value). The direction changes continuously throughout the structure as the stress field is continuous, except on the areas that are solicited in shearing, where a bar intersect another one. In these areas, the optimized material is square symmetric (i.e. R 1 = 0). The apparent stiffness having the same value in Φ 1 modulo π 4 , it is continuous in space in the optimal design. The distribution of the moduli R 0 and R 1 are illustrated in Figure 7. The R 0 values are quasi-constant whereas the R 1 values vary from 0 to 25 GPa. In the areas where R 1 are minimum, the shearing is maximum. Actually, the value of R 1 depends on the deviatoric and spherics stress ratio: if the ratio is low, shearing is predominant and R 1 should be low; and vice versa. The optimal materials in these areas where R 1 = 0 are square symmetric materials, stiffened in two orthogonal directions. When R 1 is maximum, the optimal material is stiffened in one direction because it is solicited mostly in traction or compression. Sequential optimization consists of optimizing first the shape of a structure with a fixed isotropic material, and optimizing afterwards its material anisotropy with the obtained shape. Concurrent optimization involves finding the optimal shape with the optimal material anisotropy all at once. In this section, the sequential optimization of material density and anisotropy is compared with the concurrent one.

Cantilever beam optimization

The comparison between the two methods is first performed for the cantilever beam. The optimization initialization, for both methods, is the same as in the above section: the initial volume amount is 100% and the initial material is made of an isotropic material where the values of T 0 and T 1 correspond to the isotropic part of a monolayered composite made of long and straigth carbon fibers in an epoxy matrix. Figure 8 shows that the optimal shapes are different. Since the optimal shape is sought with a fixed isotropic material in the sequential optimization, whereas the optimal topology is sought simultaneously with the optimal anisotropy distribution in the concurrent optimization, it is coherent that the obtained shapes from the two methods are dissimilar. This result confirms that the optimal shape depends on the material anisotropy and that it is necessary to optimize them simultaneously. Table 1 gives the compliance values for both methods. There is two steps in the sequential optimization: first, the shape is sought with fixed isotropic material, second, the material anisotropy distribution is optimized. Between the first and the second step in the sequential optimization, the compliance drops by 60 %. The reason is that the material anisotropy is optimized in the second step. Hence, comparing to the initial isotropic structure, the anisotropic one is more stiff. Comparison between the sequential and the concurrent optimization gives a gap of 0.3% in terms of compliance, where the solution from the simultaneous optimization is slightly better. The reason of this small gap may be the elementary nature of the presented model, the simple kind of applied sollicitations and the study of the problem in two dimensions. It would be interesting to treat a more intricate problem with complex load and boundary conditions in three dimensions space. Fig. 8: Comparison between the optimal shapes of sequential and concurrent optimization for the cantilever beam

Bridge optimization

A more complex structure is now optimized: a bridge. The design domain is a rectangle clamped in some elements at the bottom (Figure 9), and a load pressure is applied on a deck. The domain size is 304 mm x 75 mm discretized with a rectangular 76 x 40 mesh. The elements pictured in black are not designed. The volume constraint is fixed at 20% of the total volume. The initial density is set to 1 in every element of the mesh. As for the cantilever beam, the initial material is an isotropic material where the values of T 0 and T 1 correspond to the isotropic part of a monolayered composite made of long and straigth carbon fibers in an epoxy matrix: T 0 = 26.88 GPa, T 1 = 24.74 GPa.

Figure 10 compares the optimal shapes of the bridge from sequential and concurrent optimization. As for the cantilever beam, the shapes resulting from the two methods are different. Furthermore, the concurrent optimized bridge has reduced size against the sequential one. This smaller bulk is an interesting result but explanation has not been found yet. The compliance of the bridge resulting from the concurrent optimization is 8% lower than that from the sequential optimization (Table 2). Comparing with the cantilever case, the loads and boundary conditions are more complex. Thus, the gap between the compliances is higher. The advantage of the concurrent optimization compared with the sequential one is that there is no a priori on the initialization: there is no need to guess the initial shape that should be used to optimize the anisotropic material. Indeed, the optimization process takes simultaneously into account the shape and the anisotropic material distribution. Moreover, the relevance of the concurrent optimization is proved when optimizing complex models (complex loads and boundary conditions). Actually, its resulting solution is stiffer than the resulting solution of the sequential optimization. And since the material distribution influences the optimal shapes, optimizing them all at once is necessary. The proposed methodology gives concurrently optimal material spatial distribution and optimal material anisotropy repartition by minimizing the compliance (maximizing the structure global stiffness). The optimization strategy lies on a similar method as optimality criteria where iterations are made successively and separately between local minimizations and finite element calculations. Parameterizing the shape of the structure with a density variable, and the material anisotropy with polar invariants by change of frame allows for solving the local minimizations analytically as the optimality conditions are explicit. This leads to a cost of the numerical algorithm relying only on the finite element calculations. Gradient methods are nevertheless favorable for more general objective functions than compliance. However, the simplicity and the convergence of the alternate directions algorithm makes its advantage. Furthermore, the algorithm permits the optimization only with respect to the density, to the material anisotropy or both all at once: the procedure changes in the local minimizations solved with explicit formulas. The method is straightforward, efficient and gives coherent results. The optimal distribution of the material anisotropy is continous throughout the structure as it depends on the stress field. The optimal material is stiffened in two orthogonal directions when solicited in shearing, and in one direction when solicited in traction or compression. The orthotropy direction is aligned with the direction of the principal stress (highest absolute value) and the apparent stiffness is countinous in space in the optimal design. Comparing to the sequential optimization (optimal shape sought with fixed isotropic material then optimal material anisotropy optimized with the obtained shape), the concurrent optimization is interesting as there is no a priori on the initialization. Indeed, the material anisotropy influences the optimal shape. Thus, the shape of isotropic material structure, used for optimizing the material anisotropy may be suboptimal in the sequential optimization.

Optimizing the shape and the material anisotropy does not exhibit this problem. Furthermore, when complex loadings are applied on the considered structure, the proposed solution from concurrent optimization is stiffer. This methodology is a very promising when considering real composite material distribution, as the only change to be performed will be on the admissible set of polar parameters that should take into account the feasibility of the considered composite material. Future work will aim at optimizing a three dimensional structure, by using the same method.

Fig. 1 :

 1 Fig. 1: Optimal compliance with respect to material orientation in a polar representation (middle, bottom) and five optimal shapes corresponding to 0 • (left, top), 18 • (middle top), 45 • (right, top), 90 • (left, bottom) and -45 • (right, bottom) material uniform orientations of the cantilever beam problem made of long and straigth carbon fibers in an epoxy matrix (left side of the structure is clamped, a vertical load is applied in the middle of the right side).

Fig. 2 :

 2 Fig. 2: Alternate directions algorithm used to solve the double minimization of the complementary energy

Fig. 3 :Fig. 4 :

 34 Fig. 3: Boundary conditions for the Cantilever beam optimization problem

Fig. 5 :Fig. 6 : 4 ( 1 Fig. 7 :

 56417 Fig. 5: Optimal topology of the cantilever beam with 50 % volume amount

  Optimal shape of the cantilever beam for concurrent optimization

Fig. 9 :Fig. 10 :

 910 Fig. 9: Load and boundary conditions for the bridge optimization problem

Table 1 :

 1 Compliances resulting from sequential (2 steps) and concurrent (1 step) optimization for the cantilever beam optimization

	Optimization Steps	Compliance [mJ/mm]
	Sequential	Shape with fixed isotropic material Material anisotropy distribution with fixed shape 13.75 35.89
	Concurrent	Shape and material anisotropy distribution	13.69

Table 2 :

 2 Compliances resulting from sequential (2 steps) and concurrent (1 step) optimization for the bridge optimization

	Optimization Steps	Compliance [mJ/mm]
	Sequential	Shape with fixed isotropic material Material anisotropy distribution with fixed shape 1.047 2.07
	Concurrent	Shape and material anisotropy distribution	0.96
	6 Conclusion

Acknowledgements The authors acknowledge STELIA Aerospace for funding this project.