Strength-constrained simultaneous optimization of topology and fiber orientation of fiber-reinforced composite structures for additive manufacturing

Abstract

Additive manufacturing (AM) enables flexible fabrication of lightweight fiber-reinforced composite (FRC) structures with topological optimized geometries and fiber orientations. However, premature failure may occur if the manufactured FRC structure is not properly designed against stress distribution. The present study develops a strength-constrained optimization algorithm by considering the 3D printed filament embedded with fiber to be orthotropic material. The Tsai–Wu criterion is incorporated in the topological and fiber orientation optimization process. The element density and the fiber angle are formulated as design variables to simultaneously optimize the distribution of matrix and fiber. Meanwhile, both design variables are filtered to prevent checkerboard patterns of matrix elements and to ensure the continuity of fiber angle. A P-norm approach combined with an adaptive normalization scheme is employed to achieve the global strength constraint and reduce the difference between the P-norm value and the actual maximum Tsai–Wu value. Sensitivity analysis of objective function and strength-constraint conditions with respect to design variables is carried out and it serves as a basis for the method of moving asymptotes. Topological optimization of three typical numerical examples under different loading conditions are conducted and the effectiveness of the proposed method is validated. 73.7% decrease in maximum Tsai–Wu value of the optimized L-shaped beam as well as a slight increase (∼0.2%) in compliance are obtained in comparison with the compliance-based design result.

Keywords

additive manufacturing fiber-reinforced composite topology optimization orientation design strength constraint

Highlights

• A strength-constrained simultaneous optimization algorithm for FRC structure is developed by incorporating the Tsai–Wu criterion.

• Multi-load MBB beam under different loading conditions is optimized to validate the proposed algorithm.

• Premature failure can be well eliminated without significantly increasing the compliance value.

Introduction

Fiber-reinforced composite (FRC) materials have been widely applied to manufacture lightweight structures due to their high strength, enough stiffness, favorable fracture resistance, etc. (Bunsell and Renard, 2005; Saroia et al., 2020). Optimal layout of fiber orientations and the reasonably minimized compliance of lightweight FRC structures under volume constraint conditions contribute most to increasing the overall mechanical properties of FRC structures. Topology optimization is an efficient and commonly used method for lightweight design of FRC structures to obtain an optimal distribution of matrix and orientations of fibers in a predefined domain (Rong et al., 2007; Sigmund and Maute, 2013; Zhu et al., 2020; Wein et al., 2020; Huang, 2021). While owing to the limitations of traditional manufacturing methods, the fibers of FRC structures are mostly restricted to predesigned fixed directions, for example, 0^°, ±45^°, 90^°, etc. (Setoodeh et al., 2005; Stegmann and Lund, 2005; Duan et al., 2015; Salas et al., 2018; Da Silva et al., 2020). Therefore, the positive role of fiber orientations in reinforcing the FRC structures cannot be fully performed.

Additive manufacturing (AM), as a layer-by-layer manufacturing process, has been widely used in the fields of aerospace, medical, construction, machinery, etc. due to the advantages of automation, flexibility, and low cost (Culmone et al., 2019; Wang et al., 2020, 2021a, 2021b). The advanced AM technologies enable the flexible fabrication of lightweight FRC structures with topological optimized geometries and fiber orientations. Particularly, the optimal layout of fiber orientations can effectively improve the mechanical properties of FRC structures. Tekinalp et al. (2014) found that the tensile strength and modulus of 3D printed short carbon fiber reinforced composites increased ∼115% and ∼700%, respectively, compared with traditional compression molded composites due to the fiber alignment. Anwer and Naguib (2018) have validated the effect of fiber alignment for hybrid carbon fiber reinforced 3D printing composites and found that the stiffness of 0° printed composites was increased by 34% as compared to the 45° printed composites. Hambach et al. (2016) stated that stress direction oriented carbon fibers can significantly improve the flexural properties of 3D printed cement paste and an ultimate flexural strength of about 120 MPa can be achieved. AM offers the opportunity to change fiber orientations during the printing process and presents a new technology for manufacturing the FRC structures with topological optimized configurations.

Due to the flexibility of this technology, the continuous fiber angle optimization (CFAO) approaches of lightweight FRC structures have been further developed to fully utilize the properties of fiber orientation and achieve better structural performances. Sudden orientation changes between adjacent element fiber angles would lead to stress concentration and reduce the structural reliability and safety, meanwhile making it difficult to manufacture optimized FRC components, the filter for fiber angle is therefore required. The CFAO approaches can be generally divided into two categories: multi-step methods and simultaneous optimization methods. For the multi-step methods, the stable treatment of multiple design variables was ensured by taking the model as isotropic materials in the first step iteration, and then the orientation of fibers was subsequently considered, while these methods ignore the effect of fiber orientation on the evolution of material distribution in the optimization process (Bahamonde et al., 2018; Lee et al., 2018). Compared with the multi-step methods, the simultaneous optimization methods can achieve the fiber orientation and topology design in each individual iteration. The simultaneous optimization to determine the optimal fiber orientation of orthotropic materials can be achieved by two different methods: analytical methods and mathematical programming methods.

For analytical methods, the orientation of fibers shall be prescribed in each iteration and realized by the stress-based method, strain-based method, principal strain method, principal stress method, and energy-based method (Pedersen, 1989; Díaz and Bendsøe, 1992; Cheng et al., 1994; Luo and Gea, 1998). These methods effectively seek optimal orientations of shear weak materials and strong materials (Völkl et al., 2018; Safonov, 2019; Yan et al., 2020; Caivano et al., 2020; Papapetrou et al., 2020). However, in the process of iteration, the preset fiber angle distribution is still obtained by finite element analysis based on the optimal structural form. This two-step method cannot achieve the synchronous optimization design of material and fiber layout. Meanwhile, these methods cannot adapt to the application under multiple loading conditions, due to the fiber angle corresponding to each iteration cannot be prescribed in advance.

Mathematical programming methods are preferred for the structural optimization in complex loading conditions or other physical situations by taking the element density and fiber orientation as design variables in each iteration to conduct the sensitivity analysis (Nomura et al., 2015). The optimization can be solved by the adoption of gradient-based optimization algorithms such as the interior-point optimization algorithm, the method of moving asymptotes (MMA), and the globally convergent version of MMA (GCMMA) (Svanberg, 1987, 2002). Most researches focus on the simultaneous optimization of 2D components for minimizing the compliance of structure (Hoglund and Smith, 2016; Tong et al., 2017; Boddeti et al., 2018; Kim et al., 2020, 2021; Ding and Xu, 2021; Li et al., 2021; Chu et al., 2021). There are also some topology optimization methods for 3D FRC structures (Jiang et al., 2019; Schmidt et al., 2020; Jantos et al., 2020; Fedulov et al., 2021; Jung et al., 2022).

Even though various compliance-based simultaneous optimization methods have been applied to the design of FRC structures to acquire higher stiffness, while most of them failed to incorporate the strength constrained criterion, and may result in a total loss of performance if premature failures occur. To address this problem, the von Mises yield criterion has been generally used for isotropic materials in stress-constrained topology optimization methods to predict material failure (Moon and Yoon, 2013; Kiyono et al., 2016; Yang et al., 2018; Zhao et al., 2020; Xia et al., 2020). Considering 3D printing with continuous fiber reinforcement in the filament causes orthotropy, the Tsai–Wu strength criterion is more suitable for evaluating the failure of FRC additively manufactured. The Tsai–Wu strength criterion has strong applicability and high prediction accuracy than the maximum stress and maximum strain criteria, and has been widely applied to the design and evaluation of composite structures (Nikbakt et al., 2018; Mirzendehdel et al., 2018; Xu et al., 2021). The application of this criterion can be classified into two categories: as an objective function and as a constraint function. As a constraint function, the number of strength constraints increases with the element number in the finite element model, resulting in a large computational burden in the optimization process. A common strategy to this dilemma is to approximate the maximum value via global measures, such as the P-norm and Kreisselmeier–Steinhauser (KS) functions. And a P-norm function with an adaptive normalization scheme was proposed to reduce the difference between the maximum value and P-norm value (Le et al., 2010). Despite the advantages of the Tsai–Wu failure criterion, it is not used for the simultaneous optimization of FRC structures to reduce the risk of structural failure and the material failure occurred in the compliance-based design of FRC structures (Papapetrou et al., 2020). Besides, the existing compliance-based simultaneous optimization methods focus on the single constraint and single loading condition, which cannot achieve universal applicability.

In this paper, a strength-constrained criterion is introduced to simultaneously obtain the optimal topology and fiber orientation of FRC structures to address the premature failure problem, where the element density and fiber angle of design variables are filtered. Moreover, apart from the single loading condition, the proposed topology optimization method is also applied in the multiple loading condition which is more consistent with the engineering practice. The Tsai–Wu criterion is used to predict the failures of FRC materials and the sensitivity computations of strength constraint functions with respect to design variables are provided. By transforming the multi-constraint problem into a mathematical model and calculating the sensitivity value, the optimal design can be achieved effectively. MMA for solving the multi-variable and multi-constraint problem is used to obtain optimal topology and fiber orientation of FRC structures considering strength constraints. Three numerical examples under a single loading condition or multiple loading conditions are optimized to validate the effectiveness of the proposed method by comparing them to compliance-based optimized designs. The presented method and results are expected to provide guidance to practical additively manufacturing of FRC structures.

Method

Topology optimization formulation

In the present topology optimization process, element density and fiber angle are two design variables. Meanwhile, the strength is considered as a constraint condition. The strength-constrained simultaneous method obtaining the optimal topology and fiber orientation of 3D printed FRC structures is stated as follows

\begin{array}{l} find ρ = {ρ_{1}, ρ_{2}, \dots \dots, ρ_{n}}^{T}, θ = {θ_{1}, θ_{2}, \dots \dots, θ_{n}}^{T} \\ \min C (ρ, θ) = U {(ρ, θ)}^{T} K (ρ, θ) U (ρ, θ) = \sum_{e =1}^{n} ρ_{e}^{q} {u_{e}}^{Τ} k_{e} (θ_{e}) u_{e} \\ s .t . K (ρ, θ) U (ρ, θ) = F \\ V= \sum_{e =1}^{n} ρ_{e} v_{e} \leq f V_{0} \\ g_{e} (ρ, θ) \leq 1 \\ ρ_{\min} \leq ρ_{e} \leq 1 \\ θ_{\min} \leq θ_{e} \leq θ_{\max}, (e = 1,2, \dots, n) \end{array}

(1)

where ρ and θ are the element density and fiber angle vectors, respectively; C( ρ , θ ) is the structural compliance and it is equal to the sum of the element strain energy; The equilibrium equations are solved by the finite element method as K ( ρ , θ ) U ( ρ , θ ) = F , where K , U , and F are the global stiffness matrix, displacement vector, and applied load vector, respectively; q is the penalization power; u _e is the displacement vector associated with the e^th element; k _e(θ_e) is the e^th element stiffness matrix of a solid element; V and V₀ are the material and design domain volumes, respectively; v_e is the volume of the e^th element and f is the prescribed volume fraction; g_e( ρ , θ ) is the Tsai–Wu value of the e^th element; ρ_min is non-zero to avoid singularity of the stiffness matrix and is set as 10⁻⁸; θ_min and θ_max are the lower and upper bound of the design variable θ_e, respectively; n is the number of elements used to discretize the design domain.

The flowchart of the proposed algorithm is represented in Figure 1. And computational steps are summarized as follows:

Figure 1.

Flowchart of the strength-constrained simultaneous optimization method.

Step 1

Initialize material parameters, the design domain, load, support, finite element analysis mesh, and optimization parameters.

Step 2

Filter the design variables, perform a static finite element analysis, and obtain the node displacement of each discrete element to calculate the sensitivity values.

Step 3

Check the convergence: if yes, output the final result; if no, update design variables and repeat the steps described in Step 2 until the constraints are satisfied and the allowed convergence error is reached. The following equation (2) is used to check the convergence

| \frac{C^{I} - C^{I - 1}}{C^{I - 1}} | \leq ε

(2)

where C^I and C^I−1 are objective values corresponding to the I^th and (I-1)^th iterations, respectively, and ɛ is the allowed convergence error which is set as a fixed small value.

Filtering the design variables

To obtain the feasible solution and avoid the formation of checkerboard patterns, a common approach is to apply a filter to the sensitivity variables and design variables (Andreassen et al., 2011). Herein, a filter to design variables is introduced. The modified design variables in element e can be expressed as

{\bar{x}}_{e} = \frac{1}{\sum_{i \in N_{e}} H_{e i}} \sum_{i \in N_{e}} H_{e i} x_{i}

(3)

where

{\bar{x}}_{e}

is the filtered design variable;

x_{i}

is the original design variable. The neighborhood of element e, here termed N_e, is generally specified by the elements that have centers within a particular filter radius (r_min) of the center of element e; H_ei is a weight factor defined as H_ei= max (0, r_min-△(e, i)), where △(e, i) is the center-to-center distance of elements i and e. Because the e^th finite element has two unique design variables, that is, density (ρ_e) and fiber angle (θ_e), and the filtering technique is applied to them.

Stiffness and stress penalization

To create black-and-white structures, a penalization factor (q) of stiffness is introduced to reduce the number of intermediate density elements. There is a singularity problem in the optimization process that the low-density element still maintains a high stress value, which shall be completely removed to ensure the accuracy of the optimized results. Stress penalization method (Le et al., 2010) is therefore introduced in the optimization process to avoid the singularity problem.

The stiffness penalization function ( ${\bar{ρ}}_{e}^{q}$ ) is inserted when the global stiffness matrix is assembled from the solid material element stiffness matrices ( $k_{e} ({\bar{θ}}_{e})$ ). For the plane examples in the present study, quadrilateral plane stress elements and a fixed uniform mesh are used to discretize the design domain for structural analysis (Jiang et al., 2019). Under this circumstance, the stiffness matrix can be calculated as

\begin{matrix} K (\bar{ρ}, \bar{θ}) = \sum_{e =1}^{n} {\bar{ρ}}_{e}^{q} k_{e} ({\bar{θ}}_{e}) \\ k_{e} ({\bar{θ}}_{e}) = \int_{Ω e} B^{T} D ({\bar{θ}}_{e}) B d Ω \end{matrix}

(4)

where K is the global stiffness matrix;

k_{e} ({\bar{θ}}_{e})

and

D ({\bar{θ}}_{e})

are the stiffness matrix and elastic matrix of the e^th solid element, respectively, and both depend on the fiber angle θ_e of the e^th element; B is the strain–displacement matrix;

Ω_{e}

represents the region occupied by the e^th finite element; q is the penalization factor of stiffness; here, q is set as a fixed value of 3 in the entire optimization process.

The stresses at the center of each element are calculated as

σ_{e} ({\bar{ρ}}_{e}, {\bar{θ}}_{e}) = {\bar{ρ}}_{e}^{s} D ({\bar{θ}}_{e}) B u_{e}

(5)

where u _e is the displacement vector associated with the e^th element, and s is the penalization factor of the stress. Unlike the penalization factor q, the recommended value range of s is [0.5, 1) (Yang et al., 2018). For purpose of avoiding the singularity problem, case tests with different values have been verified and s is set as a fixed value of 0.8 in the entire optimization process. The effective stress vector of the e^th element under the local coordinate system is expressed as

\begin{matrix} σ_{e 1,2} ({\bar{ρ}}_{e}, {\bar{θ}}_{e}) = T^{- 1} σ_{e} ({\bar{ρ}}_{e}, {\bar{θ}}_{e}) = {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B u_{e} \\ σ_{e 1,2} ({\bar{ρ}}_{e}, {\bar{θ}}_{e}) = [\begin{matrix} σ_{1} \\ σ_{2} \\ σ_{6} \end{matrix}] \end{matrix}

(6)

where D ₀ is the original elastic matrix under the local coordinate system, and T is the transform matrix which is used to conduct the matrix coordinate transform (Chu et al., 2021); σ₁, σ₂, and σ₆ are the normal stresses and shear stress under the local coordinate system, respectively.

Global P-norm strength measure

The Tsai–Wu failure criterion is adopted in the 2D plane stress condition (Mirzendehdel et al., 2018). The Tsai–Wu value of the e^th element is expressed as

g_{e} = F_{1} σ_{1} + F_{2} σ_{2} + F_{11} σ_{1}^{2} + F_{22} σ_{2}^{2} + 2 F_{12} σ_{1} σ_{2} + F_{66} σ_{6}^{2}

(7)

where the coefficients F_i and F_ij can be obtained in terms of compressive strength values X_c and Y_c, tensile strength values X_t and Y_t, and shear strength S, which are expressed as

\begin{array}{l} F_{1} = \frac{1}{X_{t}} - \frac{1}{X_{c}}, F_{2} = \frac{1}{Y_{t}} - \frac{1}{Y_{c}} \\ F_{11} = \frac{1}{X_{t} X_{c}}, F_{22} = \frac{1}{Y_{t} Y_{c}} \\ F_{12} = - \frac{1}{2} \sqrt{F_{11} F_{22}}, F_{66} = \frac{1}{S^{2}} \end{array}

(8)

The number of strength constraints is related to the number of finite elements used to discretize the design domain, and therefore it is too large for computations and practical applications. Especially for the strength-constrained optimization under a multiple loading condition, the number of strength constraints is multiplied, which is more likely to lead to the problem of difficult convergence of iteration results. This computational burden is often reduced using a single global strength measure, which makes the adjoint sensitivity analysis computationally efficient. These n constraints can be expressed in terms of a single maximum strength constraint

g_{\max} (ρ, θ) = \max_{e = 1,2, ..., n} (g_{e} (ρ, θ)) \leq 1

(9)

However, the maximum function is not differentiable and must be smoothed. The P-norm aggregation function is introduced to approximate the maximum strength by clustering all constraints into a global constraint. It is expressed as

g_{pn} = {(\sum_{e =1}^{n} g_{e}^{P})}^{\frac{1}{P}}

(10)

where g_pn is the clustered value that approximates the maximum strength, g_e is the Tsai–Wu value at the centroid of the e^th element, and P is the P-norm factor. The clustered value g_pn is known to approach the maximum strength (g_max) when P→∞ and a higher value of P is preferable as it provides a more accurate approximation of g_max. However, excessively high values of P result in numerical problems. Consequently, the value of the P-norm factor must be appropriately selected. Because the Tsai-Wu values can be negative and odd values of P-norm factor produces wrong sensitivities, the value of P must be even. Based on previous research (Kiyono et al., 2016; Yang et al., 2018) and numerical tests with different values, P is set as a fixed value of 8 in the entire optimization process to ensure numerical stability. To reduce the difference between g_pn and g_max, the correction factor introduced by Le et al. (2010) is adopted and formulated as

g_{\max} \approx c g_{pn}

(11)

where c is the correction factor that is calculated at each optimization iteration and is expressed as

c^{I} = α^{I} \frac{g_{\max}^{I - 1}}{g_{p n}^{I - 1}} + (1 - α^{I}) c^{I - 1}

(12)

where c^I and c^I−1 are the correction factors of the I^th and (I-1)^th iterations, respectively; α^I controls the variations between c^I and c^I−1 and is set as a fixed value of 0.5.

Sensitivity analysis

Compared with the compliance-based simultaneous optimization methods, the introduction of strength constraints is the key point in this study. Therefore, the sensitivity analysis of the strength constraint with respect to design variables is more crucial for improving the algorithm. The sensitivity of a function φ with respect to the design variable (x_j) is obtained using the chain rule

\frac{\partial φ}{\partial x_{j}} = \sum_{e \in N_{j}} \frac{\partial φ}{\partial {\bar{x}}_{e}} \frac{\partial {\bar{x}}_{e}}{\partial x_{j}} = \sum_{e \in N_{j}} \frac{1}{\sum_{i \in N_{e}} H_{e i}} H_{j e} \frac{\partial φ}{\partial {\bar{x}}_{e}}

(13)

where φ represents either the objective function or the constraint function; the original design variable (x_j) represents the element density or the element fiber angle;

{\bar{x}}_{e}

is the filtered design variable; N_j is the neighborhood of element j, which is generally specified by the elements that have centers within a particular filter radius (r_min) of the center of element j; H_je is a weight factor defined as H_je = max (0, r_min-△(j,e)), where △(j,e) is the center-to-center distance of elements j and e.

According to the definition in equation (1), the sensitivity of the objective function with respect to the filtered design variable is

\begin{array}{l} \frac{\partial C}{\partial {\bar{ρ}}_{e}} = - q {\bar{ρ}}_{e}^{q - 1} {u_{e}}^{Τ} k_{e} ({\bar{θ}}_{e}) u_{e} \\ \frac{\partial C}{\partial {\bar{θ}}_{e}} = - {\bar{ρ}}_{e}^{q} {u_{e}}^{Τ} \frac{\partial k_{e} ({\bar{θ}}_{e})}{\partial {\bar{θ}}_{e}} u_{e} = - {\bar{ρ}}_{e}^{q} {u_{e}}^{Τ} \int_{Ω e} B^{T} \frac{\partial D ({\bar{θ}}_{e})}{\partial {\bar{θ}}_{e}} B d Ω u_{e} \end{array}

(14)

According to the definition in equation (1), the sensitivity of the volume constraint with respect to the filtered design variable is

\frac{\partial V}{\partial {\bar{ρ}}_{e}} = v_{e}, \frac{\partial V}{\partial {\bar{θ}}_{e}} = 0

(15)

In terms of the chain rule, the sensitivity of the strength constraint with respect to the filtered design variable is

\frac{\partial g_{\max}}{\partial {\bar{x}}_{i}} = c \sum_{e = 1}^{n} \frac{\partial g_{pn}}{\partial g_{e}} {(\frac{\partial g_{e}}{\partial σ_{e 1,2}})}^{T} \frac{\partial σ_{e 1,2}}{\partial {\bar{x}}_{i}}

(16)

where the derivative of the P-norm function with respect to the Tsai–Wu value can be obtained

\frac{\partial g_{pn}}{\partial g_{e}} = {(\sum_{e =1}^{n} g_{e}^{P})}^{\frac{1}{P} - 1} g_{e}^{P - 1}

(17)

By the definition in equation (7), the derivatives of the Tsai–Wu value with respect to its stress components are

\begin{array}{l} \frac{\partial g_{e}}{\partial σ_{1}} = F_{1} + 2 F_{11} σ_{1} + 2 F_{12} σ_{2} \\ \frac{\partial g_{e}}{\partial σ_{2}} = F_{2} + 2 F_{22} σ_{2} + 2 F_{12} σ_{1} \\ \frac{\partial g_{e}}{\partial σ_{6}} = 2 F_{66} σ_{6} \end{array}

(18)

According to the definition in equation (6), the derivative of the effective stress vector under the local coordinate system with respect to the filtered design variable can be expressed as

\begin{array}{l} \frac{\partial σ_{e 1,2}}{\partial {\bar{ρ}}_{i}} = \frac{\partial {\bar{ρ}}_{e}^{s}}{\partial {\bar{ρ}}_{i}} D_{0} T {({\bar{θ}}_{e})}^{T} B u_{e} + {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B \frac{\partial u_{e}}{\partial {\bar{ρ}}_{i}} \\ \frac{\partial σ_{e 1,2}}{\partial {\bar{θ}}_{i}} = {\bar{ρ}}_{e}^{s} D_{0} \frac{\partial T {({\bar{θ}}_{e})}^{T}}{\partial {\bar{θ}}_{i}} B u_{e} + {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B \frac{\partial u_{e}}{\partial {\bar{θ}}_{i}} \end{array}

(19)

where

\frac{\partial {\bar{ρ}}_{e}^{s}}{\partial {\bar{ρ}}_{i}} \neq 0

and

\frac{\partial T {({\bar{θ}}_{e})}^{T}}{\partial {\bar{θ}}_{i}} \neq 0

only for i=e,

\frac{\partial u_{e}}{\partial {\bar{ρ}}_{i}}

and

\frac{\partial u_{e}}{\partial {\bar{θ}}_{i}}

can be obtained through the finite element method as KU = F.

The matrix L _e gathers the nodal displacements of the e^th element from the global displacement vector satisfying u _e = L_eU. Assuming that the external load F is independent of the design variables, then

\begin{array}{l} \frac{\partial U}{\partial {\bar{ρ}}_{i}} = - K^{- 1} \frac{\partial K}{\partial {\bar{ρ}}_{i}} U \\ \frac{\partial u_{e}}{\partial {\bar{ρ}}_{i}} = - L_{e} K^{- 1} \frac{\partial K}{\partial {\bar{ρ}}_{i}} U \\ \frac{\partial U}{\partial {\bar{θ}}_{i}} = - K^{- 1} \frac{\partial K}{\partial {\bar{θ}}_{i}} U \\ \frac{\partial u_{e}}{\partial {\bar{θ}}_{i}} = - L_{e} K^{- 1} \frac{\partial K}{\partial {\bar{θ}}_{i}} U \end{array}

(20)

The sensitivity of the strength constraint with respect to the filtered design variable can be further expressed as

\begin{array}{l} \frac{\partial g_{\max}}{\partial {\bar{ρ}}_{i}} = c \sum_{e = 1}^{n} \frac{\partial g_{pn}}{\partial g_{e}} {(\frac{\partial g_{e}}{\partial σ_{e 1,2}})}^{T} (\frac{\partial {\bar{ρ}}_{e}^{s}}{\partial {\bar{ρ}}_{i}} D_{0} T {({\bar{θ}}_{e})}^{T} B u_{e} - {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B L_{e} K^{- 1} \frac{\partial K}{\partial {\bar{ρ}}_{i}} U) \\ = c \frac{\partial g_{pn}}{\partial g_{i}} {(\frac{\partial g_{i}}{\partial σ_{i 1,2}})}^{T} s {\bar{ρ}}_{i}^{s - 1} D_{0} T {({\bar{θ}}_{i})}^{T} B u_{i} - c \sum_{e = 1}^{n} \frac{\partial g_{pn}}{\partial g_{e}} {(\frac{\partial g_{e}}{\partial σ_{e 1,2}})}^{T} {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B L_{e} K^{- 1} \frac{\partial K}{\partial {\bar{ρ}}_{i}} U \\ \frac{\partial g_{\max}}{\partial {\bar{θ}}_{i}} = c \sum_{e = 1}^{n} \frac{\partial g_{pn}}{\partial g_{e}} {(\frac{\partial g_{e}}{\partial σ_{e 1,2}})}^{T} ({\bar{ρ}}_{e}^{s} D_{0} \frac{\partial T {({\bar{θ}}_{e})}^{T}}{\partial {\bar{θ}}_{i}} B u_{e} - {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B L_{e} K^{- 1} \frac{\partial K}{\partial {\bar{θ}}_{i}} U) \\ = c \frac{\partial g_{pn}}{\partial g_{i}} {(\frac{\partial g_{i}}{\partial σ_{i 1,2}})}^{T} {\bar{ρ}}_{i}^{s} D_{0} \frac{\partial T {({\bar{θ}}_{i})}^{T}}{\partial {\bar{θ}}_{i}} B u_{i} - c \sum_{e = 1}^{n} \frac{\partial g_{pn}}{\partial g_{e}} {(\frac{\partial g_{e}}{\partial σ_{e 1,2}})}^{T} {\bar{ρ}}_{e}^{s} D_{0} T {({\bar{θ}}_{e})}^{T} B L_{e} K^{- 1} \frac{\partial K}{\partial {\bar{θ}}_{i}} U \end{array}

(21)

In this paper, the adjoint method is employed for sensitivity analysis because only a global strength aggregation constraint exists. The adjoint vector ( λ ) can be defined according to the following equation

Kλ= \sum_{e = 1}^{n} \frac{\partial g_{pn}}{\partial g_{e}} {\bar{ρ}}_{e}^{s} L_{e}^{T} B^{T} T ({\bar{θ}}_{e}) D_{0}^{T} \frac{\partial g_{e}}{\partial σ_{e 1,2}}

(22)

where λ is obtained by solving the adjoint equation.

Substituting equation (22) into equation (21), the sensitivity of the strength constraint with respect to the filtered design variable is calculated as

\begin{array}{l} \frac{\partial g_{\max}}{\partial {\bar{ρ}}_{i}} = c \frac{\partial g_{pn}}{\partial g_{i}} {(\frac{\partial g_{i}}{\partial σ_{i 1,2}})}^{T} s {\bar{ρ}}_{i}^{s - 1} D_{0} T {({\bar{θ}}_{i})}^{T} B u_{i} - c λ^{T} \frac{\partial K}{\partial {\bar{ρ}}_{i}} U \\ \frac{\partial g_{\max}}{\partial {\bar{θ}}_{i}} = c \frac{\partial g_{pn}}{\partial g_{i}} {(\frac{\partial g_{i}}{\partial σ_{i 1,2}})}^{T} {\bar{ρ}}_{i}^{s} D_{0} \frac{\partial T {({\bar{θ}}_{i})}^{T}}{\partial {\bar{θ}}_{i}} B u_{i} - c λ^{T} \frac{\partial K}{\partial {\bar{θ}}_{i}} U \end{array}

(23)

Recalling the stiffness matrix in equation (4), the derivative of the global stiffness matrix with respect to the filtered design variable is

\begin{array}{l} \frac{\partial K}{\partial {\bar{ρ}}_{i}} = q {\bar{ρ}}_{i}^{q - 1} L_{i}^{T} k_{i} ({\bar{θ}}_{i}) L_{i} \\ \frac{\partial K}{\partial {\bar{θ}}_{i}} = {\bar{ρ}}_{i}^{q} L_{i}^{T} \int_{Ω i} B^{T} \frac{\partial D ({\bar{θ}}_{i})}{\partial {\bar{θ}}_{i}} B d Ω L_{i} \end{array}

(24)

The sensitivity of the strength constraint with respect to the filtered design variable can eventually be evaluated as

\begin{array}{l} \frac{\partial g_{\max}}{\partial {\bar{ρ}}_{i}} = c \frac{\partial g_{pn}}{\partial g_{i}} {(\frac{\partial g_{i}}{\partial σ_{i 1,2}})}^{T} s {\bar{ρ}}_{i}^{s - 1} D_{0} T {({\bar{θ}}_{i})}^{T} B u_{i} - c λ_{i}^{T} {q {\bar{ρ}}_{i}^{q - 1} k}_{i} ({\bar{θ}}_{i}) u_{i} \\ \frac{\partial g_{\max}}{\partial {\bar{θ}}_{i}} = c \frac{\partial g_{pn}}{\partial g_{i}} {(\frac{\partial g_{i}}{\partial σ_{i 1,2}})}^{T} {\bar{ρ}}_{i}^{s} D_{0} \frac{\partial T {({\bar{θ}}_{i})}^{T}}{\partial {\bar{θ}}_{i}} B u_{i} - c λ_{i}^{T} {\bar{ρ}}_{i}^{q} \int_{Ω i} B^{T} \frac{\partial D ({\bar{θ}}_{i})}{\partial {\bar{θ}}_{i}} B d Ω u_{i} \end{array}

(25)

where λ _i is the vector of the adjoint nodal values of the i^th element.

Case study and simultaneous optimization results

Three commonly adopted numerical examples are used to validate the proposed method. The optimized results of strength-constrained and compliance-based simultaneous optimization are analyzed with compliance values and maximum Tsai–Wu values, respectively, and the proposed method is extended to the multi-load scenario. In all numerical examples, the material properties and strengths are set as E₁ = 39 GPa, E₂ = 8.4 GPa, G₁₂ = 4.2 GPa, v₁₂ = 0.26, X_t = 1062 MPa, X_c = 610 MPa, Y_t = 31 MPa, Y_c = 118 MPa, S = 72 MPa, respectively. To avoid oscillation of the objective function value in the iterative process, the move limits of the element density and the element fiber angle are set to 0.1 and π/9, respectively. And all the initial values of the element density design variable are set to 1. The parameter values of θ_min and θ_max are referred to the work of Xia and Shi (2017). The filter radius is twice the element size. In the multi-variable and multi-constraint situations, the selection of MMA optimizer parameters is particularly important, which directly affects the convergence and stability of iteration results. The values of parameters such as asyinit, asyincr, raa0, asydecr, epsimin, and albefa are used to update design variables. Parameter values of the MMA optimizer are listed in Table 1.

Table 1.

Parameter values of the MMA optimizer.

Asyinit	Asyincr	Raa0	Asydecr	Epsimin	Albefa
0.5	1.5	0.01	0.8	10⁻⁸	0.4

L-shaped beam

The first case study is an L-shaped beam. Because the design domain contains a re-entrant corner with an initial geometric stress singularity at which the highest stress is located, the material failure cannot be avoided through traditional compliance-based optimization design. Thus, design results must be modified significantly. Relevant research has demonstrated that structures with rounded corners are preferred to prevent stress concentration at the corners for isotropic materials.

The dimensions and boundary/load conditions are shown in Figure 2, where L = 100 mm and the thickness = 1 mm. The top edge is fixed and a vertical load of F = 200 N is applied at the top point of the right vertical edge. The vertical load is distributed on five nodes to avoid stress concentration. The mesh consists of 6400 elements, and the volume fraction of this example is 0.4. The initial fiber orientations and the set value of the convergence condition are π/2 and 1×10⁻⁵, respectively.

Figure 2.

Dimensions and boundary/load conditions of L-shaped beam.

Figure 3 shows the convergence histories of the L-shaped beam of strength-constrained simultaneous optimization. The iterative curve of the maximum Tsai–Wu value tends to increase first and then decrease, and converges to the constant value of approximately 1. The initial compliance value is 223.87 N·mm, and it attains the maximum value of 2492.12 N·mm at the 10^th iteration. Thereafter, the iterative curve of the compliance exhibits a descending trend. Finally, results satisfy the constraints and convergence conditions, and the iteration terminates at the 616^th iteration. As Figure 3 shows, all element density design variables begin to transform from the initial solid material to the optimized result. The fiber angles are represented by different colors and all of the fiber angles gradually change from the initial π/2 to even distribution. Finally, a reasonable distribution of the material and fiber orientation is realized. The final compliance value is 239.33 N·mm. The compliance value increases by 6.9% from 223.87 N·mm to 239.33 N·mm when the volume fraction of the scenario decreases by 60%. From optimized results of the L-shaped beam, not only the structural strength is significantly improved, but also the loss of structural stiffness can be negligible, although the structural material is reduced. It is proved that the proposed method is conducive to the realization of the structural design on lightweight high-performance 3D printed FRC structures.

Figure 3.

Convergence histories of strength-constrained simultaneous optimization.

Figure 4 shows the optimized results of the L-shaped beam considering strength constraints. The fiber orientation is basically consistent with the structural form and maintains good continuity in Figure 4(a). And the corner of the L-shaped beam is rounded, which is similar to the stress-constrained topology optimization of isotropic materials. The optimized Tsai–Wu value distribution is shown in Figure 4(b) and the absolute value is selected for plotting. The figure shows that the maximum value appears at the L-shaped beam corner, which is 0.9992. The Tsai-Wu values of the optimized components’ connection positions are slightly larger than that of other positions. This is due to the greater variation of the fiber orientation in these locations. The Tsai–Wu values for the scenario are smaller than 1, which avoids causing structural damage.

Figure 4.

Optimized results of strength-constrained simultaneous optimization (a) density and fiber distribution; (b) Tsai–Wu value distribution.

To further prove the necessity of considering strength constraints, the optimized results of the L-shaped beam without strength constraints are shown in Figure 5. Comparing the result of compliance-based simultaneous optimization (Figure 5(a)) and the result of strength-constrained simultaneous optimization (Figure 4(a)), fiber orientations of the two results are basically distributed along with structural forms and have good continuity. However, the structural forms are considerably different, particularly in the corner positions. The corner of the L-shaped beam considering strength constraints is rounded, and the optimized result without strength constraints is still a re-entrant corner. Comparing Figure 4(b) and Figure 5(b), the maximum Tsai–Wu values of optimized results still appear at the L-shaped beam corner. But the maximum Tsai–Wu value without strength constraints is larger than 1 which is 3.7959, and there is a noticeable failure at the corner. The result considering strength constraints avoid material failure because of satisfying strength constraints. A comparison of two optimized results shows that the maximum Tsai–Wu value decreases from 3.7959 to 0.9992 and is reduced by 73.7% through the addition of strength constraints, which illustrates the necessity of considering strength constraints.

Figure 5.

Optimized results of compliance-based simultaneous optimization (a) density and fiber distribution; (b) Tsai–Wu value distribution.

Figure 6 shows the values of compliance and Tsai-Wu in different scenarios. A represents the initial design without simultaneous optimization, B represents the compliance-based simultaneous optimization, and C represents the strength-constrained simultaneous optimization. Through a comparison with the objective function value in different scenarios, the compliance value of C is similar to the initial design without simultaneous optimization, and the difference is only 15.46 N·mm. The maximum Tsai–Wu values of A and B are 4.7737 and 3.7959, respectively, which are significantly larger than 1, indicating that the material failure has occurred. However, the maximum Tsai–Wu value of C is less than 1 and the material failure is avoided, which further proves the necessity of considering strength constraints during the simultaneous optimization of the topology and fiber orientation.

Figure 6.

Compliance and maximum Tsai–Wu values in different scenarios.

Cantilever beam

This example is a cantilever beam with a length of L = 100 mm. The dimensions and boundary/load conditions of the cantilever beam problem are shown in Figure 7. The vertical load of F = 400 N is distributed on five nodes to avoid stress concentration. The mesh consists of 100 × 50 elements, and the volume fraction is 0.5. The initial fiber orientations and the set value of the convergence condition are 0 and 1 × 10⁻⁵, respectively.

Figure 7.

Dimensions and boundary/load conditions of a cantilever beam.

The convergence histories of the cantilever beam considering strength constraints are shown in Figure 8. The iterative curve of the maximum Tsai–Wu value tends to increase first and then decrease and converges to a constant value of approximately 1. The initial compliance value is 236.35 N·mm. At the 175^th iteration, the results satisfy the constraints and convergence conditions, and the iteration terminates. The final compliance value is 289.34 N·mm. When the volume fraction of the cantilever beam is decreased by 50%, corresponding compliance value is observed to increase up to 22.4% from 236.35 N·mm to 289.34 N·mm. As shown in Figure 8, the objective function curve is provided with good convergence and stability. In addition, attributed to the variation of Tsai-Wu value with the change of design variables, the maximum Tsai–Wu value curve fluctuates slightly, which can be neglected in the final design result.

Figure 8.

Convergence histories of strength-constrained simultaneous optimization.

Figure 9 shows the results of the simultaneous topology and fiber orientation optimization of the cantilever beam. By comparing Figures 9(a) and (c), the fiber orientation is basically consistent with the structural form. But the optimized result of strength-constrained simultaneous optimization is significantly asymmetric, which is considerably different compared with the compliance-based simultaneous optimization. Precisely because of the introduction of strength constraints, the force transmission path of the structure is changed and thus generates the satisfying design result. The compliance values of strength-constrained and compliance-based methods are 289.34 N·mm and 285.14 N·mm, respectively. The difference is only 4.2 N·mm. Figure 9(b) shows the Tsai–Wu value distribution of the final topology optimization results considering strength constraints. Compared with Figure 9(d), the maximum Tsai–Wu value is 0.9158 and reduced by 60.4%, which does not reach the failure value. Because there is good fiber continuity at the intersection nodes, the occurrence of failure is avoided, which confirms that the constraints are well satisfied. The comparative analysis shows that the addition of strength constraints could greatly change the optimal design result, rather than fine-tuning the structural form obtained by the compliance-based method.

Figure 9.

Optimized results of simultaneous optimization (a) density and fiber distribution; (b) Tsai–Wu value distribution of strength-constrained method; (c) density and fiber distribution; (d) Tsai–Wu value distribution of compliance-based method.

Multi-load scenario

Notably, in practical application, the loading condition of the engineering structure is usually complicated, and most of them are under the multi-load situation. The final example is a Messerschmitt–Bolkow–Blohm (MBB) beam design problem in which two situations are considered. The difference between a single loading condition and a multiple loading condition is discussed. The two loads are loaded at the bisection points, and each load is 200 N. The dimensions and boundary/load conditions are shown in Figure 10, where L = 120 mm and the thickness = 1 mm. The mesh consists of 120 × 60 elements, and the volume fraction of this example is 0.3. The initial fiber orientations and the set value of the convergence condition are 0 and 4 × 10⁻⁵, respectively.

Figure 10.

Dimensions and boundary/load conditions of an MBB beam.

For the single loading condition, both loads are applied simultaneously and the finite element analysis is performed only once. The optimal topology and fiber orientation in the strength-constrained case has been obtained and it is presented in Figure 11(a). The optimal result in the unconstrained case is presented in Figure11(c). The compliance values are 97.18 N·mm and 93.33 N·mm, respectively. By applying strength constraints, the width and tilt angle of internal components change significantly to meet the design criteria. As shown in Figure 11(d), the maximum Tsai–Wu value of compliance-based method exceeds the threshold of material failure. On the other hand, in Figure 11(b), the maximum Tsai–Wu value is given by 0.9948 and the strength is under control clearly.

Figure 11.

Optimized results of the single loading condition (a) density and fiber distribution; (b) Tsai–Wu value distribution of strength-constrained method; (c) density and fiber distribution; (d) Tsai–Wu value distribution of compliance-based method.

The multiple loading condition is used to calculate the compliance and Tsai-Wu values of optimized results of the single load case. The compliance value is 119.18 N·mm and Figure 12 shows the Tsai–Wu value distribution under two load conditions. The figure shows that the maximum Tsai–Wu value is larger than 1 and optimized results of the single-load case are not applicable to multiple-load cases. Therefore, the method under the single loading condition is far from meeting the requirements on structural design optimization. It is necessary to realize a more reasonable optimization designing scheme under the multiple loading situation.

Figure 12.

Tsai–Wu value distribution of optimized results under two load conditions.

For the multiple loading condition, the optimal topology is not the superposition of topological forms under the single loading condition. Each of multi-loads acts on the structure and the finite element analysis is performed separately. The objective compliance is calculated as the sum of compliances corresponding to each load and the number of strength constraints after P-norm aggregation is 2.

Figure 13 shows the convergence histories of the MBB beam considering strength constraints under multiple loading. Element density design variables begin to transform from the initial solid material, and fiber-angle design variables gradually change from the initial 0 to be consistent with the structural form. Finally, a reasonable distribution of the material and fiber orientation is realized. The iterative curve of the maximum Tsai–Wu value tends to decrease and converge to a constant value of approximately 1. The initial compliance value is 74.17 N·mm, and it attains the maximum value of 849.87 N·mm at the 11^th iteration. Thereafter, the iterative curve of the compliance exhibits a descending trend. When the iteration number is 100, the iterative change is slow indicating that the structural form becomes clear. The final compliance value is 68.80 N·mm. The compliance value decreases by 7.2% when the volume fraction decreases by 70%, which proves that optimizing the fiber orientation can effectively increase the structural stiffness. Therefore, the reduction of the quantity of materials does not necessarily have an adverse effect on the mechanical properties. By virtue of the strengthening effect of fiber orientation and the introduction of the reasonable constraint conditions, the optimal design of 3D printed FRC structures with the characteristics of high strength, high stiffness could be realized under the multi-load condition.

Figure 13.

Convergence histories of strength-constrained simultaneous optimization under multiple loading.

The optimized results of the MBB beam under multiple loading are presented. In Figures 14(a) and (b), the fiber orientation is basically consistent with the structure form and the optimized result is almost symmetric. Owing to the addition of strength constraints, the optimized layout of material and fiber must be changed. The structure forms have changed dramatically to avoid material failure. As shown in Figures 14(c) and (d), the maximum Tsai–Wu value is 0.9806, which does not reach the failure value. Compared with Figures 14(e) and (f), the design of strength-constrained simultaneous optimization can be guaranteed without material failure, indicating that the proposed method is applicable for multiple loading scenarios. Commonly, the introduction of constraints could lead to an increase in objective functions. The compliance shows an increase of 3.82%, in comparison with 66.27 N·mm obtained by the compliance-based simultaneous optimization. As shown in the result, the proposed method can reduce the increase of objective functions as well as the stiffness loss under the same volume constraint conditions.

Figure 14.

Optimized results of simultaneous optimization under multiple loading (a) density and fiber distribution of strength-constrained method; (b) density and fiber distribution of compliance-based method; (c) Tsai–Wu value distribution in case 1; (d) Tsai–Wu value distribution in case 2 of strength-constrained method; (e) Tsai–Wu value distribution in case 1; (f) Tsai–Wu value distribution in case 2 of compliance-based method.

Conclusions

In this paper, a strength-constrained simultaneous optimization method is presented to obtain the optimal topology and fiber orientation of 3D printed FRC structures. Moreover, apart from the single loading condition, the proposed topology optimization method is also extended to the application of multi-load situation. Three typical numerical examples are simulated to validate the feasibility and effectiveness of the proposed method. Conclusions can be drawn as follows:

(1) The proposed method is able to avoid premature failures of FRC structures with topological optimized configurations compared with the compliance-based optimization design method. The maximum Tsai–Wu value of the optimized L-shaped beam design is reduced by 73.7% compared with the compliance-based design. The strength of structure could be increased remarkably in the case of minimizing the stiffness loss.

(2) The P-norm approach combined with an adaptive normalization scheme was adopted in the strength-constrained simultaneous optimization and the computational burden problem with many strength constraints was well solved. Especially in multi-constraint optimization cases, this approach exhibits great advantages, which could realize the effective control of strength constraint and minimum the computational cost.

(3) Desired convergence of the optimization process for different cases are obtained. The fiber orientation maintains good continuity and is consistent with the topological configuration, which greatly facilitates the additive fabrication of FRC structures.

(4) Element density and fiber angle are selected as two variables in the iterative process, which achieves the mathematicalization of the mechanical problem, and avoids the shortcomings of human intervention on the analytical methods. Furthermore, by utilizing the interaction between structural topology and fiber angle, the free optimization design of 3D printed FRC structures could be realized.

The proposed strength-constrained method can be extended to structures with multi-objective optimization, and meanwhile can be extended to the structural optimization of 3D structures. AM constraints, such as nozzle size restrictions and overhang constraints shall be included in the calculation process, to realize the actual fabrication of the optimized FRC structures. Next-step work shall focus on the integrated design of structural optimization and printing path planning to facilitate the smooth additive fabrication of lightweight FRC structures.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: National Natural Science Foundation of China (grant numbers 51878241, 52178198, 51808183), and Natural Science Foundation of Hebei Province (grant numbers E2021202039 and E2019202484).

ORCID iD

Guowei Ma

References

Anwer

Naguib

(2018) Multi-functional flexible carbon fiber composites with controlled fiber alignment using additive manufacturing. Additive Manufacturing 22: 360–367

Andreassen

Clausen

Schevenels

, et al. (2011) Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization 43: 1–16

Bahamonde

Tatting

Harik

, et al. (2018) An analysis framework for topology optimization of 3D printed reinforced composites. In: 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Kissimmee, Florida, 8–12 January 2018

Boddeti

Ding

Kaijima

, et al. (2018) Simultaneous digital design and additive manufacture of structures and materials. Scientific Reports 8: 15560

Bunsell

Renard

(2005) Fundamentals of Fibre Reinforced Composite Materials. Boca Raton: CRC Press

Caivano

Tridello

Paolino

, et al. (2020) Topology and fibre orientation simultaneous optimisation: A design methodology for fibre-reinforced composite components. Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials: Design and Applications 234: 1267–1279

Cheng

Kikuchi

(1994) An improved approach for determining the optimal orientation of orthotropic material. Structural Optimization 8: 101–112

Chu

Xiao

Gao

, et al. (2021) Robust topology optimization for fiber-reinforced composite structures under loading uncertainty. Computer Methods in Applied Mechanics and Engineering 384: 113935

Culmone

Smit

Breedveld

(2019) Additive manufacturing of medical instruments: A state-of-the-art review. Additive Manufacturing 27: 461–473

10.

Da Silva

ALF

Salas

Silva

ECN

, et al. (2020) Topology optimization of fibers orientation in hyperelastic composite material. Composite Structures 231: 111488

11.

Díaz

Bendsøe

(1992) Shape optimization of structures for multiple loading conditions using a homogenization method. Structural Optimization 4: 17–22

12.

Ding

(2021) A novel discrete-continuous material orientation optimization model for stiffness-based concurrent design of fiber composite. Composite Structures 273: 114288

13.

Duan

Yan

Zhao

(2015) Integrated optimization of the material and structure of composites based on the Heaviside penalization of discrete material model. Structural and Multidisciplinary Optimization 51: 721–732

14.

Fedulov

Fedorenko

Khaziev

, et al. (2021) Optimization of parts manufactured using continuous fiber three-dimensional printing technology. Composites Part B: Engineering 227: 109406

15.

Hambach

Möller

Neumann

, et al. (2016) Portland cement paste with aligned carbon fibers exhibiting exceptionally high flexural strength (> 100MPa). Cement and Concrete Research 89: 80–86

16.

Hoglund

Smith

(2016) Continuous fiber angle topology optimization for polymer fused filament fabrication. In: Proceedings of the 26th Annual International Solid Freeform Fabrication Symposium – An Additive Manufacturing Conference, Austin, Texas, USA, 8–10 August 2016, pp. 1078–1090

17.

Huang

(2021) On smooth or 0/1 designs of the fixed-mesh element-based topology optimization. Advances in Engineering Software 151: 102942

18.

Jantos

Hackl

Junker

(2020) Topology optimization with anisotropic materials, including a filter to smooth fiber pathways. Structural and Multidisciplinary Optimization 61: 2135–2154

19.

Jiang

Hoglund

Smith

(2019) Continuous fiber angle topology optimization for polymer composite deposition additive manufacturing applications. Fibers 7: 14

20.

Jung

Lee

Nomura

, et al. (2022) Inverse design of three-dimensional fiber reinforced composites with spatially-varying fiber size and orientation using multiscale topology optimization. Composite Structures 279: 114768

21.

Kim

Lee

Yoo

(2021) Machine learning-combined topology optimization for functionary graded composite structure design. Computer Methods in Applied Mechanics and Engineering 387: 114158

22.

Kim

Lee

Nomura

, et al. (2020) Topology optimization of functionally graded anisotropic composite structures using homogenization design method. Computer Methods in Applied Mechanics and Engineering 369: 113220

23.

Kiyono

Vatanabe

Silva

ECN

, et al. (2016) A new multi-p-norm formulation approach for stress-based topology optimization design. Composite Structures 156: 10–19

24.

Norato

Bruns

, et al. (2010) Stress-based topology optimization for continua. Structural and Multidisciplinary Optimization 41: 605–620

25.

Lee

Kim

Nomura

, et al. (2018) Topology optimization for continuous and discrete orientation design of functionally graded fiber-reinforced composite structures. Composite Structures 201: 217–233

26.

Gao

, et al. (2021) Full-scale topology optimization for fiber-reinforced structures with continuous fiber paths. Computer Methods in Applied Mechanics and Engineering 377: 113668

27.

Luo

Gea

(1998) Optimal orientation of orthotropic materials using an energy based method. Structural Optimization 15: 230–236

28.

Mirzendehdel

Rankouhi

Suresh

(2018) Strength-based topology optimization for anisotropic parts. Additive Manufacturing 19: 104–113

29.

Moon

Yoon

(2013) A newly developed qp-relaxation method for element connectivity parameterization to achieve stress-based topology optimization for geometrically nonlinear structures. Computer Methods in Applied Mechanics and Engineering 265: 226–241

30.

Nikbakt

Kamarian

Shakeri

(2018) A review on optimization of composite structures Part I: laminated composites. Composite Structures 195: 158–185

31.

Nomura

Dede

Lee

, et al. (2015) General topology optimization method with continuous and discrete orientation design using isoparametric projection. International Journal for Numerical Methods in Engineering 101(8): 571–605

32.

Papapetrou

Patel

Tamijani

(2020) Stiffness-based optimization framework for the topology and fiber paths of continuous fiber composites. Composites Part B: Engineering 183: 107681

33.

Pedersen

(1989) On optimal orientation of orthotropic materials. Structural Optimization 1(2): 101–106

34.

Rong

Jiang

Xie

(2007) Evolutionary structural topology optimization for continuum structures with structural size and topology variables. Advances in Structural Engineering 10(6): 681–695

35.

Safonov

(2019) 3D topology optimization of continuous fiber-reinforced structures via natural evolution method. Composite Structures 215: 289–297

36.

Salas

Ramírez-Gil

Montealegre-Rubio

, et al. (2018) Optimized dynamic design of laminated piezocomposite multi-entry actuators considering fiber orientation. Computer Methods in Applied Mechanics and Engineering 335: 223–254

37.

Saroia

Wang

Wei

, et al. (2020) A review on 3D printed matrix polymer composites: its potential and future challenges. The International Journal of Advanced Manufacturing Technology 106: 1695–1721

38.

Schmidt

Couret

Gout

, et al. (2020) Structural topology optimization with smoothly varying fiber orientations. Structural and Multidisciplinary Optimization 62: 3026–3105

39.

Setoodeh

Abdalla

Gürdal

(2005) Combined topology and fiber path design of composite layers using cellular automata. Structural and Multidisciplinary Optimization 30: 413–421

40.

Sigmund

Maute

(2013) Topology optimization approaches. Structural and Multidisciplinary Optimization 48: 1031–1055

41.

Stegmann

Lund

(2005) Discrete material optimization of general composite shell structures. International Journal for Numerical Methods in Engineering 62: 2009–2027

42.

Svanberg

(1987) The method of moving asymptotes - A new method for structural optimization. International Journal for Numerical Methods in Engineering 24: 359–373

43.

Svanberg

(2002) A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM Journal on Optimization 12(2): 555–573

44.

Tekinalp

Kunc

Velez-Garcia

, et al. (2014) Highly oriented carbon fiber–polymer composites via additive manufacturing. Composites Science and Technology 105: 144–150

45.

Tong

Zhang

(2017) Optimal fiber orientation and topology design for compliant mechanisms with fiber-reinforced composites. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 231: 2302–2312

46.

Völkl

Klein

Franz

, et al. (2018) An efficient bionic topology optimization method for transversely isotropic materials. Composite Structures 204: 359–367

47.

Wang

Jiang

, et al. (2020) Mechanical behaviors of 3D printed lightweight concrete structure with hollow section. Archives of Civil and Mechanical Engineering 20: 16

48.

Wang

Liu

Yang

, et al. (2021a) Bonding performance of 3D printing concrete with self-locking interfaces exposed to compression–shear and compression-splitting stresses. Additive Manufacturing 42: 101992

49.

Wang

Liu

, et al. (2021b) Interlayer reinforcement of 3D printed concrete by the in-process deposition of U-nails. Cement and Concrete Research 148: 106535

50.

Wein

Dunning

Norato

(2020) A review on feature-mapping methods for structural optimization. Structural and Multidisciplinary Optimization 62: 1597–1638

51.

Xia

Shi

(2017) Optimization of composite structures with continuous spatial variation of fiber angle through Shepard interpolation. Composite Structures 182: 273–282

52.

Xia

Wang

Liu

(2020) Uncertainty-oriented topology optimization of interval parametric structures with local stress and displacement reliability constraints. Computer Methods in Applied Mechanics and Engineering 358: 112644

53.

Gao

, et al. (2021) On design of carbon fiber reinforced plastic (CFRP) laminated structure with different failure criteria. International Journal of Mechanical Sciences 196: 106251

54.

Yan

Hua

, et al. (2020) Concurrent topology optimization of structures and orientation of anisotropic materials. Engineering Optimization 52: 1598–1611

55.

Yang

Liu

Zhang

, et al. (2018) Stress-constrained topology optimization based on maximum stress measures. Computers & Structures 198: 23–39

56.

Zhao

Han

, et al. (2020) Continuum structural topological optimization with dynamic stress response constraints. Advances in Engineering Software 148: 102834

57.

Zhu

Zhang

, et al. (2020) Design of compliant mechanisms using continuum topology optimization: a review. Mechanism and Machine Theory 143: 103622