Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2024.48.0022 Acta Polytechnica CTU Proceedings 48:22–26, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague A GEOMETRY PROJECTION METHOD FOR DESIGNING AND OPTIMIZING ADDITIVELY MANUFACTURED VARIABLE-STIFFNESS COMPOSITE LAMINATES Yogesh Gandhia,∗, Julián Noratob, Ana Pavlovica, Giangiacomo Minaka a Alma Mater Studiorum – Università di Bologna, Department of Industrial Engineering, Via Fontanelle 40, Forlì-47121, Italy b University of Connecticut, Department of Mechanical Engineering, 191 Auditorium Road, U-3139 Storrs, CT 06269, United States ∗ corresponding author: yogesh.gandhi@unibo.it Abstract. A method for designing laminates is presented using geometry projection to optimize the layout of additively manufactured variable-stiffness composite laminates. By considering fiber- reinforced bars as geometric primitives, the geometry projection methodology is extended to include optimizing regions with intersecting load paths. This is achieved by utilizing a dual representation of bars, which considers the geometric parameters and the element-wise density field representation. The dual representation enables the combining and overlapping of bars, resulting in a localized orthotropic material response at overlapping regions that mitigates the transverse compliant response of fiber- reinforced components. The proposed method’s effectiveness is demonstrated through minimizing the compliance of the Messerschmitt-Bölkow-Blohm beam problem, a well-known benchmark problem in topology optimization. Keywords: Topology optimization, geometry projection, continuous fiber-reinforced polymers, variable- stiffness laminates. 1. Introduction Our study focuses on developing and optimiz- ing variable-stiffness composite laminates (VSCLS) that can be produced using additive manufactur- ing. These VSCLs consist of layers with differ- ent orientations of fibers, which are customized to achieve specific mechanical properties, such as maximizing the stiffness-to-weight ratio. Design- ing VSCLs involves considering various factors, in- cluding fiber orientation, stacking order, and layer thickness. Unlike constant stiffness composites [1], optimizing VSCLs, as discussed elsewhere [2], is challenging due to the absence of analytical for- mulations. Therefore, designing VSCLs necessi- tate the use of discretization methods, such as the finite element method (FEM), to represent vari- ations in material constituents. This leads to increased design variables, making the optimiza- tion process computationally intensive. As a re- sult, computational design tools like topology op- timization (TO) have become crucial for achiev- ing the desired mechanical performance. These tools can also accommodate modifications result- ing from advanced manufacturing processes [3], such as continuous fiber-fused filament fabrication (CF4). Topology optimization (TO) is an iterative pro- cess for finding the optimal material distribution in a design domain, minimizing an objective function under constraints. Various methods include homog- enization, solid isotropic material with penalization (SIMP) [4], level-set [5], evolutionary structural opti- mization, phase field, and feature mapping [6]. De- tails on these approaches have been reviewed else- where [7, 8]. Various techniques have also been ad- vanced to determine the optimal distribution of fiber orientations within the design domain, as reviewed elsewhere [9]. In this paper, a geometry projection (GP) based [10] procedure is formulated to enable the design of VS- CLs. While using fiber-reinforced bars (FRBs) as features, we modify the GP methodology in multiple ways. Initially, we leverage the dual nature of the ge- ometry projection method to define overlapping FRBs in the design domain. These overlapping FRBs are then retained and modeled using composite laminate theory to calculate homogenized stiffness matrices, resulting in a local orthotropic material response. Ad- ditionally, forming overlapping FRBs in the design domain reduces the high strain energy density at in- tersecting load paths, which is further amiable for optimization. The component-wise GP formulation enables the seamless printing of VSCLs with the clear manifestation of overlapping components. Nonethe- less, the proposed method to design VSCLs is limited to a single layer, which we call GP-AM. 22 https://doi.org/10.14311/APP.2024.48.0022 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 48/2024 GP-AM 2. Geometry projection formulation for variable-stiffness composite laminates The GP procedure begins by mapping the design variables of the bar to a density field ρb (x; zb), where x represents any point in the design region. The dual geometric parameters/density representation can consider individual bars as high-level geometric objects or field variables in component densities. The design in this work involves the combination of multiple FRBs. Each bar b ∈ B, where B denotes the set of all bar indices, is represented as a rectangle with semicircular caps whose medial axis is a line segment, which occupies a region Ωb ∈ R2. The medial axis is characterized by its two endpoints (x1b,x2b), and the offset distance corresponds to the radius rb of the bar. A membership variable αb ∈ [0, 1] is assigned to each bar and penalized similarly to density-based methods, allowing the optimizer to remove or include it in the design. The design variable vector zb for bar b is therefore expressed as: zb := (x1b,x2b, rb, αb) . (1) The projected density at a location x is determined by the intersection of a ball with a radius r and cen- tered at x with Ωb, i.e.: ρb(x; zb) := |Br x ∩ Ωb (zb)| |Br x| . (2) In 2D, under the assumption that r is significantly smaller than the dimensions of the bar, the intersection of Br x and ∂Ωb can be approximated as a line segment. As a result, the proportion of the projected density can be determined by calculating the proportion of the circular segment with a height of h = r − ϕb, where ϕb represents the signed distance from x to ∂Ωb. The projected density for bar b is a uniquely determined function of ϕb, effectively functioning as a regularized Heaviside function: ρb(x; zb) := H̃ ( ϕb(x; zb) r ) . (3) The expression for H̃ is given. H̃(x) =  0, if x ≤ −1 1 + 1 π ( x √ 1 − x2 − arccosx ) , if |x| < 1 1, if x ≥ 1 dH̃ dx = { 2 √ 1 − x2/π, if |x| < 1 0, otherwise . A penalized density is determined for each bar which is utilized to compute its elastic stiffness tensor. This is similar to the SIMP technique employed in density- based topology optimization [11]. The penalized den- sity is expressed as: ρ̆eff b (x; zb) := (αbρb(x; zb))q , (4) where we recall q is the penalization factor. 2.1. Combining components Eq. (2) pertains to the projection of a single bar. This study aims to focus on continuous fiber reinforcement, particularly in areas where multiple reinforcements overlap. This requires the material interpolation to accurately represent the stiffness that arises from over- lapping bars. A straightforward approach is to com- pute the element’s elasticity tensor as the sum of all elasticity tensors: Ce = Cv + nb∑ b=1 ρ̆eff be (Cb − Cv) . (5) The elasticity tensor interpolation Ce given by (5) can be used to interpolate between the solid material (bar material Cb) and void material Cv. 2.2. Elasticity tensor for overlapping components The bar’s material coordinate system (MCS), repre- sented by {ê1b, ê2b, ê3b}, is typically not the same as the laminate coordinate system (LCS), {e1, e2, e3}. We assume that the fiber reinforcement aligns with ê1b, and we define ê2b so that it is orthogonal to ê1b. Furthermore, we define e3 = ê1b × ê2b such that it corresponds to the global out-of-plane axis. For plane stress of the laminate, the components of the elasticity tensor Cb of bar b in LCS are given by Cp b = T ⊤ 1 Ĉ p bT 1 and Cs b = T ⊤ 2 Ĉ s bT 2, (6) with T 1 = [ c2 s2 cs s2 c2 −cs −2cs 2cs c2 − s2 ] , T 2 = [ c −s s c ] , (7) where, c and s denote the cosine and sine of the angle θ respectively, where θ represents the angle between the axes ê1 and x. The 3×3 matrix Ĉ p b establishes the relationship between the in-plane strains {εx, εy, γxy} and the stresses {σx, σy, τxy} in LCS. Similarly, the 2 × 2 matrix Ĉ s b governs the relationship between the out-of-plane shear strains {γxz, γyz} and the stresses {τxz, τyz} in LCS. 2.3. Elemental laminate elasticity matrix For a laminate that is symmetrical with respect to the middle plane (z = 0) and consists of 2nb layers, the membrane-bending coupling matrix Dmd is equal to zero. We utilize the first shear deformation theory (FSDT) to analyze the mechanical behavior of VSCLs. In the FSDT, we relax Kirchhoff’s theory by consid- ering that the transverse normal no longer remains 23 Y. Gandhi, J. Norato, A. Pavlovic, G. Minak Acta Polytechnica CTU Proceedings Material E1 E2 v12 G12 G13 G23 Carbon epoxy AS4/3501-6 113.6 GPa 9.650 GPa 0.334 6.0 GPa 6.0 GPa 3.1 GPa Table 1. Material properties used for all the examples. perpendicular to the mid-plane (z = 0) after deforma- tion. This implies that the elastic displacement field within the laminate is expanded by as: u(x, y, z) = u(x, y) + zψx(x, y) v(x, y, z) = v0(x, y) + zψy(x, y) w(x, y, z) = w0(x, y), (8) where u0(x, y) and v0(x, y) represent in-plane displace- ments, while w0(x, y) represents the out-of-plane dis- placement. The subscript (·)0 indicates the displace- ment field of the reference plane. Using four-node, bilinear, quadrilateral, and plane- stress elements, we construct element stiffness matri- ces based on the abovementioned assumptions. The stiffness matrix for the eth element can be represented as: Dme = nb∑ i=1 [hi+1 − hi] Cp ei; Dse = nb∑ i=1 [hi+1 − hi]κCs ei; Dde = nb∑ i=1 1 3 [ h3 i+1 − h3 i ] Cp ei. (9) Each bar is defined by the top and bottom planes hi+1 and hi. The shear correction factor is κ = 5/6. Additionally, we assume that the layer thickness H is uniform, which means that (hi+1 − hi) = H, and we also assume that the effect of stacking on the bending stiffness is disregarded, with hb = 0, leading to (h3 b+1 − h3 b) = H3 for all bars. 3. Optimization problem and sensitivity analysis We address the minimization of compliance for a spec- ified volume fraction limit. We formulate the opti- mization problem as follows: min{zb} f := log(c+ 1) subject to: v ≤ v̄ KU = f zi ≤ zi ≤ z̄i, i = 1, 2, . . . , nz, (10) where compliance, denoted by c = U⊤f , and U and f represent the global displacement and force vectors, respectively. The value v̄ sets an upper limit on the volume fraction, while K represents the global stiff- ness matrix. The interval [zi, z̄i] defines the lower and upper bounds on the ith design variable. Figure 1. Design region, boundary conditions, and initial design for beam in 3−point in-plane bending. The volume fraction, v, can be defined as v(e) := 1∑ e ∣∣Ω(e) ∣∣ ∣∣∣Ω(e) ∣∣∣ ∑ b ρeff be . (11) Finally, for the design-independent loading, the sensi- tivity of Eq. (10) and Eq. (11) is as follows: ∂zic = − ∑ e u⊤ ( ∂ziK (e) ) u ∂ziv = ∑ e ∣∣Ω(e) ∣∣ ∂ziv (e)∑ e ∣∣Ω(e) ∣∣ . (12) The stiffness matrices design sensitivities Da, where {a = m, d, s} in Eq. (9) can be ex- pressed in matrix form as follows: ∂ziDe =  ∂ziDme 03 03×2 03 ∂ziDde 03×2 02×3 02×3 ∂ziDse  . (13) For the matrix 0n, it’s a zero matrix of size n× n, and for 0m×n, it’s a zero matrix of size m× n. The chain rule is used to compute the design sensitivity of Eq. (13). This requires finding the derivative of the transformation in Eq. (6) and the penalized effec- tive density. The computation of these derivatives is straightforward but has been left out for brevity. More detailed information can be found elsewhere [10]. 4. Examples We consider bars made of CFRP; fiber orientation is continuously aligned to the bar’s axis. Table 1 lists unidirectional carbon-epoxy AS4/3501-6 material properties used for the bars. The following settings are considered until men- tioned otherwise. The method-of-moving-asymptotes (MMA) [12] optimizer is used for the optimization rou- tine, with the default parameters described in here [13]. The void material assigned with a Young’s modulus of Evoid = 0.001 GPa and a Poisson’s ratio of vvoid = 0.3. In the initial design (see Figure 1), bars have a radius 24 vol. 48/2024 GP-AM Design Combined density 0 1 2 3 Table 2. Optimal designs and combined density for the MBB beam. The first column corresponds to the compliance value obtained by design, i.e., c = 0.0905 kN · mm. In the second column, transparency is used in the color of the bars for the GP-AM designs to help distinguish overlapping components. The color of the bars for the optimal designs in the last column corresponds to their penalized membership variable value αq b, and the bars are plotted with transparency to facilitate the visualization of all bars. 0 100 200 300 10−1 100 iteration f f0 0 100 200 3000.2 0.3 0.4 0.5 iteration v v0 Figure 2. Objective (left) and volume constraint (right) history for the MBB beam design of Tab. 2 corresponding to the GP-AM design. equal to the average of their upper and lower bounds during initialization. The sizing variable, α, is set to 0.5, and the move limit, m, is fixed and set to 0.02 throughout the optimization process. The procedure of optimization entails three criteria for stopping. The first criterion is achieved when the 2-norm of the alteration in the vector of design vari- ables is smaller than 0.002. The second criterion is attained when the KKT (Karush-Kuhn-Tucker) op- timality condition norm is below 0.002. The third criterion is fulfilled when the modification in the ob- jective function is less than 10−9. The optimization process is stopped if any of these criteria are met. 4.1. Beam in 3-point bending The example considers a 150 mm × 50 mm simply sup- ported beam in 3-point bending, as illustrated in Figure 1. The maximum volume fraction is 0.5. Be- cause the problem is symmetrical, only the right side of the structure is modeled. Initially, the entire plate is designed with 27 bar, and the design variables are constrained within the following bounds.{ (0, 0) 2 0 } ≤ { (xb1, xb2) rb αb } ≤ { (150, 50) 5 1 } (14) In this situation, the highest bending stress occurs in the upper left corner where the load is applied, and the normal stresses cause tension on the bottom edge and compression on the top edge, with significant shear stresses near the neutral axis. Additionally, the optimal design for a single load application aligns bars’ fiber reinforcement axis with the principal stress at each location [14]. In Table 2, it can be observed that there are hori- zontal bars at the top and bottom edges to represent the normal stress distribution, while inclined bars are used for the shear stress distribution. The overlap- ping bars contribute to a larger second-moment area of inertia, resulting in higher bending stiffness, which optimally supports the applied load. Furthermore, the method proposed in this study yields an optimal solu- tion comparable to the findings of other density-based methods for the minimum compliance problem. Figure 2 illustrates that the optimization process fol- lows a usual convergence pattern for compliance; there is a significant decrease in the initial iterations, accom- panied by minor adjustments to the design variables in the subsequent optimization iterations. Additionally, it displays a standard smooth convergence behavior. 5. Conclusions The numerical experiment results demonstrate the effectiveness of the GP-AM in generating optimal so- lutions for the minimum compliance problem. From a manufacturing standpoint, the designs obtained with the proposed methods are amenable to manufacturing techniques for VSCL. Although not presented here for brevity, the proposed method produces good designs for several design regions, for example, problems: rect- angular plate in pure torsion; cantilever beam under in- and out-of-plane bending; and square membrane in out-of-plane bending. 25 Y. Gandhi, J. Norato, A. Pavlovic, G. Minak Acta Polytechnica CTU Proceedings 6. Acknowledgement Financed by the European Union-NextGenerationEU (Grant No 2022NW83RE_002). The opinions ex- pressed are those of the authors only and should not be considered representative of the European Union or the European Commission’s official position. Neither the European Union nor the European Commission can be held responsible for them. References [1] H. Ghiasi, D. Pasini, L. Lessard. Optimum stacking sequence design of composite materials Part I: Constant stiffness design. Composite Structures 90(1):1–11, 2009. https: //doi.org/10.1016/j.compstruct.2009.01.006 [2] H. Ghiasi, K. Fayazbakhsh, D. Pasini, L. Lessard. Optimum stacking sequence design of composite materials Part II: Variable stiffness design. Composite Structures 93(1):1–13, 2010. https: //doi.org/10.1016/j.compstruct.2010.06.001 [3] J. Plocher, A. Panesar. Review on design and structural optimisation in additive manufacturing: Towards next-generation lightweight structures. Materials & Design 183:108164, 2019. https://doi.org/10.1016/j.matdes.2019.108164 [4] G. I. N. Rozvany. A critical review of established methods of structural topology optimization. Structural and Multidisciplinary Optimization 37:217–237, 2009. https://doi.org/10.1007/s00158-007-0217-0 [5] N. P. Van Dijk, K. Maute, M. Langelaar, F. Van Keulen. Level-set methods for structural topology optimization: A review. Structural and Multidisciplinary Optimization 48:437–472, 2013. https://doi.org/10.1007/s00158-013-0912-y [6] F. Wein, P. D. Dunning, J. A. Norato. A review on feature-mapping methods for structural optimization. Structural and Multidisciplinary Optimization 62:1597–1638, 2020. https://doi.org/10.1007/s00158-020-02649-6 [7] O. Sigmund, K. Maute. Topology optimization approaches. Structural and Multidisciplinary Optimization 48(6):1031–1055, 2013. https://doi.org/10.1007/s00158-013-0978-6 [8] J. D. Deaton, R. V. Grandhi. A survey of structural and multidisciplinary continuum topology optimization: Post 2000. Structural and Multidisciplinary Optimization 49:1–38, 2014. https://doi.org/10.1007/s00158-013-0956-z [9] Y. Gandhi, G. Minak. A review on topology optimization strategies for additively manufactured continuous fiber-reinforced composite structures. Applied Sciences 12(21):11211, 2022. https://doi.org/10.3390/app122111211 [10] H. Smith, J. A. Norato. Topology optimization with discrete geometric components made of composite materials. Computer Methods in Applied Mechanics and Engineering 376:113582, 2021. https://doi.org/10.1016/j.cma.2020.113582 [11] M. P. Bendsøe, O. Sigmund. Material interpolation schemes in topology optimization. Archive of Applied Mechanics 69:635–654, 1999. https://doi.org/10.1007/s004190050248 [12] K. Svanberg. The method of moving asymptotes – a new method for structural optimization. International Journal for Numerical Methods in Engineering 24(2):359– 373, 1987. https://doi.org/10.1002/nme.1620240207 [13] K. Svanberg. MMA and GCMMA, versions September 2007. Optimization and systems theory 104, 2007. [14] P. Pedersen. Some general optimal design results using anisotropic, power law nonlinear elasticity. Structural optimization 15:73–80, 1998. https://doi.org/10.1007/BF01278492 26 https://doi.org/10.1016/j.compstruct.2009.01.006 https://doi.org/10.1016/j.compstruct.2009.01.006 https://doi.org/10.1016/j.compstruct.2010.06.001 https://doi.org/10.1016/j.compstruct.2010.06.001 https://doi.org/10.1016/j.matdes.2019.108164 https://doi.org/10.1007/s00158-007-0217-0 https://doi.org/10.1007/s00158-013-0912-y https://doi.org/10.1007/s00158-020-02649-6 https://doi.org/10.1007/s00158-013-0978-6 https://doi.org/10.1007/s00158-013-0956-z https://doi.org/10.3390/app122111211 https://doi.org/10.1016/j.cma.2020.113582 https://doi.org/10.1007/s004190050248 https://doi.org/10.1002/nme.1620240207 https://doi.org/10.1007/BF01278492 Acta Polytechnica CTU Proceedings 48:22–26, 2024 1 Introduction 2 Geometry projection formulation for variable-stiffness composite laminates 2.1 Combining components 2.2 Elasticity tensor for overlapping components 2.3 Elemental laminate elasticity matrix 3 Optimization problem and sensitivity analysis 4 Examples 4.1 Beam in 3-point bending 5 Conclusions 6 Acknowledgement References