Topology Optimization and Efficiency Evaluation of Short-Fiber Reinforced Composite Structures considering Anisotropy

1. Introduction

The aerospace and automotive industries have inherent requirements: high stiffness, high strength, and minimum mass are the primary concerns when designing structures [1,2,3,4]. For instance, high stiffness and minimal mass are required to limit the natural frequency of spacecraft to avoid resonance between the launch vehicle and itself [4,5,6], and in automobiles to improve performance while satisfying safety requirements [7,8,9]. Minimal mass is required for keeping the spacecraft structure lightweight to get to the desired orbital altitude by increasing the useful load fraction [10], lowering aircraft engineering costs by reducing airframe weight [11,12], and reducing aircraft and automobile operational costs by increasing fuel efficiency [11,12,13]. High strength is linked to the minimal mass requirement since materials with higher specific strength allow to keep mass to a minimum. The same applies to the high stiffness requirement by using materials with high specific stiffness.

The strength of the structure depends on the structure’s shape and material. Regarding materials, aluminum and titanium [14] are commonly used in the aerospace and automobile industries due to their high specific strength and stiffness. Moreover, new materials have been researched for automobile [15] and aerospace applications [16,17,18,19], including ceramics [20], composites [21], and nanocomposites [22]. Furthermore, nanocomposites [23,24], shape memory polymers [25,26], and short fiber-reinforced composites (SFRC) have experienced an increase in aerospace applications in recent years [20,27,28,29,30,31,32,33]. The principal metrics for comparing materials, rather than Young’s modulus or the ultimate tensile strength, are the specific stiffness and strength as both mechanical characteristics are normalized by material density. Current strategies for designing efficient structures can be categorized into guidelines and recommendations [34,35], methods and methodologies [36,37,38], and even parametric [39,40,41] and topology optimizations [40,42,43,44,45,46,47]. Topology optimization (TO) has proven to be an excellent mathematical method for designing efficient structures and optimizing material arrangement within a given design space for a set of loads, boundary conditions, and constraints with the goal of minimizing a specific objective function. This method, also known as the variable density model, was originally proposed in the work [48], where it was applied to the optimal design of two-dimensional structures to overcome the computational difficulties associated with stress analysis of variable thickness plates. Three-dimensional (3D) TO under strength and stiffness constraints, using material density as a design variable, was later presented in the studies [49,50,51]. Another approach to representing variable density using materials at the microstructure level and its application to TO problems can be found in the study [52], further explained in the work [53]. Designed structures are usually evaluated by considering maximum displacement, maximum strength, failure criteria values [54,55,56], or minimum mass of the structure, which are not reliable indicators since these values are local, except for mass, which is evaluated as an integral characteristic and does not provide comprehensive information about the structure other than its value. However, the load-carrying factor LCF G, proposed by Andrey Komarov in 1965 [48], allows us to relate both structural efficiency and weight requirements, facilitating objective comparisons between structures with relatively low levels of computational and modeling effort [57]. The fundamental physical meaning of the load-carrying factor is the integral characteristic of the structure, reflecting the internal forces in its elements (the mode of action of an external force on a structure up to its supports) and the extent of these internal forces (the extent of the external forces transmission paths). A disadvantage of the LCF is its dimensional dependence; if the geometrical dimensions of the structure and external loads are changed, the magnitude of the LCF also changes, making it difficult to compare different structural arrangements. Therefore, the usage of the LCF coefficient, as proposed by Valery Komarov in the work [58], is more appropriate. In the accompanying article [59], a proposal is made to employ the method of mixtures for assessing the weight efficiency of composite structures, necessitating a detailed examination of composite materials at the microstructure level. From a practical standpoint, the utilization of strength criteria at the representative volume mesolevel appears to be more advantageous when evaluating complex-shaped structures for weight efficiency. In the present work, the strength criterion for homogenized composite material will be considered.

In summary, current modern materials such as SFRC, TO, metrics for assessing the structure shape as failure criteria, and LCF have contributed to increase stiffness and strength while reducing the mass of structures. However, the interaction between them is often not carefully considered. The fiber orientation tensor is obtained by means of solving the Folgar-Tucker’s continuity equation, which is developed based on the fiber orientation kinetic theory of fiber suspensions [60,61,62,63]. SFRC is modelled as transversely isotropic material, even though its stiffness depends not only on the material’s mechanical characteristics but also on the fiber orientation tensor [64]; TO is typically performed in an isotropic medium [65,66], and LCF accounts for the equivalent stress of structures made of isotropic materials. In other words, anisotropy is either partially or completely ignored at various stages of the design process.

This work aims to develop a methodology for obtaining topology-optimal structures made of SFRC. The objectives of this work are as follows: to make the stiffness matrix dependent on the fiber orientation tensor, to obtain the fiber orientation tensor by performing injection molding simulations, to formulate a metric for assessing SFRC structures that simultaneously evaluates the shape and material of the structure, and to experimentally verify the developed methodology. It is assumed that by considering the anisotropy during TO, the resulting topology will perform better than a topology obtained in an isotropic medium. Moreover, the anisotropy is directly related to the arrangement of fibers within the structure, specifically to the fiber orientation tensor, which is obtained by solving Folgar-Tucker’s continuity equation [67] through modeling the injection molding process within the designated design region in Autodesk Moldflow. The consideration of anisotropy induced by the fiber arrangement is achieved by implementing the Advani-Tucker orientation averaging procedure on the material stiffness matrix [68] and attaching the resultant tensor to each mesh element. This procedure was realized with AnisoTopo [69]. The attachment of the stiffness matrix was accomplished by interpolating the meshes of the injection and structural models with Digimat MAP. Finally, TO was performed without altering the stiffness arrangement within the design region. The resultant topology can be assessed using a metric such as a modified LCF, which is described in terms of failure criteria for composite materials, such as the Tsai-Hill failure criterion.

The numerical and experimental results of this study confirm that considering the anisotropy during TO increases the stiffness of the resultant topologies. This is achieved by obtaining the fiber orientation tensor by solving the Navier-Stokes equations for non-Newtonian fluids and making the stiffness matrix dependent on the obtained fiber orientation tensor.

2. Materials and Methods

2.1. Materials and material models

Three different materials and their non-linear models were considered (see Figure 1 and Table 1). These materials include two anisotropic materials: a 50% glass fiber reinforced polyamide 6, denoted as PA6 50GF [70], along with its associated model as presented in [71], and a 30% carbon fiber reinforced polyamide 6, referred to as PA6 30CF [72] with its corresponding model [73]. Additionally, an isotropic material, the D16T aluminum alloy [74], which shares similarities with the 2024 -T4 ASTM aluminum alloy. It’s worth noting that no specific preparation steps or treatments were applied to these materials before testing. Mechanical characteristics of these materials were determined through tensile tests conducted on samples cut from injected-molded plates with dimensions 200 x 150 x 4 mm, following ISO 527-2:2012(en) standard [75]. While the primary mechanical characteristics relevant to this study are summarized in Table 1, more comprehensive details on the material characteristics can be found in their respective documentation.

The material models for PA6 50GF and PA6 30GF were calibrated with experimental data as described in [73]. In summary, this calibration method involves modeling the material microstructure as a two-phase material and taking into account the fiber orientation in the calculation of mechanical properties. The material’s mechanical properties were obtained by homogenizing the fibers and matrix using a second-order Mori–Tanaka homogenization scheme in Digimat MF. As the mechanical properties correspond to a unidirectional fiber-reinforced composite, it is modeled as a transversely isotropic material. Accounting for the fiber orientation was achieved by applying Tucker’s averaging procedure. The matrix stress–strain state is described using the J2 plasticity model based on von Mises’ equivalent stress. When the equivalent stress σ_eq exceeds the yield stress σ_y, the reaction becomes nonlinear, and plastic deformation occurs. The plastic strength is then defined as

σ_plastic = σ_y + R(ε_p),

(1)

where σ_y is the yield stress; R(ε_p)=kε_p + R_∞[1 - e^-mεp] represents the isotropic strain exponential and linear hardening law, and ε_p is the accumulated plastic strain. Here, k is the linear hardening modulus in MPa; m is the hardening exponent; and R_∞ is the hardening modulus, MPa. The material model parameters were adjusted by minimizing the difference between the tensile strain-stress curves of the composite material and the experimental results.

2.2. Methods

2.2.1 Topology optimization

Topology optimization (TO) represents an effective approach for optimizing material distribution within a structure to enhance the transfer of internal forces from areas of load-carrying areas to support regions. Since the strain energy quantifies the energy stored in a body due to deformation, minimizing this response leads to a reduction in the body’s compliance [51]. Strain energy is mathematically defined as W = 0.5 Fu, where F = Ku is the force expressed as the product of the stiffness K and deformation u. The stiffness is influenced by anisotropy resulting from the arrangement of fibers within the body, specifically the fiber orientation tensor.

In this study, the objective function of TO was to minimize the structure compliance by reducing the total strain energy of an anisotropic composite material. This was achieved by adjusting the topology density ρ within the design region Ω while adhering to the given design region volume to retain g₁(x) and constraints such as the “Member size” and the “Pull-out direction”. The formulation was defined as follows:

minimize f(x) = 0.5u^TK(ρ(x), A(x))u,

by varying x ∈ Ω,

subject to h(x) = K(ρ(x), A(x))u - F = 0,

g₁(x) = ∫ ρ(x)dΩ - V_ret ≤ 0,

g₂(x) = ρ(x) > 0,

g₃(x) = ρ(x) ≤ 1,

g₄(x) = ∫|∇ρ(x)|dΩ - δ ≤ 0,

g₅(x) = ρ_i - ρ_k ≤ 0         ∀ x_i = x_k, y_i = y_k, |z_i| ≥ |z_k|,

Here, K represents the global stiffness matrix, u denotes the nodal displacement vector, F is the nodal force vector, and x is the vector containing design domain elements with coordinates x, y, z. A represents the fiber orientation tensor, and δ is related to the minimum structural member size. TO, in this work, was performed using the “Sequential Convex Programming” solver within the Ansys Mechanical Workbench software.

2.2.1.1 Unidisciplinary topology optimization considering constant anisotropy

Anisotropy was modeled by adjusting the elastic constants of the stiffness matrix. Since composites reinforced with 30–50 wt.% fiber along the flow direction exhibit double the stiffness compared to that across the flow direction [76], the anisotropy was incorporated by modifying stiffness matrix components. The mechanical characteristics of the anisotropic material were aligned with the global coordinate axis and remained constant throughout the TO, allowing for the optimization of structure arrangement based on a predefined stiffness distribution.

2.2.1.2 Multidisciplinary topology optimization considering variable anisotropy

As explained in the previous section, modeling the anisotropy during structural analysis involves incorporating the relevant material characteristics of an anisotropic material. However, this approach is suitable for cutting structures from injection-molded plates (two-dimensionally) but not for molded structures (three-dimensionally). In molded structures, the fiber orientation within the channels, which are structural elements from the perspective of solid mechanics, cannot be defined until TO is performed.

The technique outlined in [76] was applied to account for fiber arrangement within the channels and dynamically adjust the mechanical characteristics of each structural element by recalculating the element’s fiber orientation tensor. Algorithm 1 presents the pseudocode for multidisciplinary TO considering variable anisotropy, with detailed information on files and program execution available in [77].

The process begins with the introduction of initial values, including boundary conditions for structural analysis BC_struct and injection molding simulation BC_injMold, material properties of the matrix and fiber MP_matrix+fiber, injection molding material MP_injMold, TO parameters OP_topoOpt, topology density threshold th, design region volume percentage to retain V_ret, interpolation tolerance between the injection molding and structural analysis meshes δ_map, and an optional geo.stp_file for geometry. MP_matrix+fiber and OP_topoOpt are saved as materialProperties.txt and topoParameters.txt, respectively, for correct utilization in AnisoTopo [69]. The structural domain region mesh, elements participating in TO, and elements related to boundary conditions are obtained from Ansys Workbench and saved as designRegionMesh.ans, design.txt, and frozen.txt, respectively. TO with variable fiber orientation proceeds after specifying the convergence criterion for TO algorithm, which is the relative difference ε between the previous and current objective function value, not exceeding the objective relative difference ε_obj. Key moments in TO include:

First, the structural mesh is exported to Autodesk Moldflow. If it’s the first iteration, reducedMesh.ans will be the same as designRegionMesh.ans as it hasn’t been previously reduced. Otherwise, the structural mesh is reduced using the “delete_elements” algorithm and saved as reducedMesh.ans (details in [77]).

Second, the reducedMesh.ans, BC_injMold, and MP_injMold are introduced to Moldflow for calculating the fiber orientation tensor A, which is exported along with the injection molding mesh mesh_injMold and saved in the file meshMoldFlow.pat and fiberOrientMoldFlow.xml, respectively. Fiber orientation mapping from injection molding to structural analysis mesh is performed in DigimatMAP, and the mapped fiber orientation A’ is stored in the file fiberOrientAnsys.xml.

Third, in AnisoTopo [69], A’ is introduced to calculate the anisotropic stiffness matrix (K) using the Advani-Tucker orientation averaging technique applied to the material stiffness matrix [4]. AnisoTopo’s code employs the Mori-Tanaka homogenization method to compute mechanical properties.

The composite strain and stress are contingent on the strain and stress of both the matrix and fiber in proportion to their volume fractions:

$$\epsilon =\left(1-{\varphi }_f\right){\epsilon }_m+{\varphi }_f{\epsilon }_f,$$

$$\sigma =\left(1-{\varphi }_f\right){\sigma }_m+{\varphi }_f{\sigma }_f,$$

where ${\varphi }_f$ is the fiber volume fraction, and the subscripts m denote the strain and stress of matrix values, while the subscripts f correspond with the strain and stress of the fiber.

The unidirectional short-fiber reinforced composite material is modeled as transversely isotropic. The elastic moduli, as introduced by Tandon and Weng [78], were used to calculate the elastic coefficients:

$$\frac{E_{11}}{E_m}=\frac{1}{1+\frac{{\varphi }_f\left(A_1+2{\upsilon }_m{A}_2\right)}{A_6}}$$

$$\frac{E_{22}}{E_m}=\frac{1}{1+\frac{{\varphi }_f\left[-2{\upsilon }_m{A}_3+\left(1-{\upsilon }_m\right)A_4+\left(1+{\upsilon }_m\right)A_5{A}_6\right]}{2A_6}},$$

where $E_m$, ${\upsilon }_m$ are the Young’s modulus and Poisson ratio of the matrix, respectively. The parameters $A_i$ are the functions of Eshelby’s tensor and can be found in [79]. In this particular study, we empoly Eshelby’s tensor for an elliptical inclusion, which depends on the fiber’s aspect ratio. The Tucker’s averaging procedure is used to account for the fiber orientation tensor, described as follows:

$$C_{ijkl}=B_1{a}_{ijkl}+B_2\left(a_{ij}{\delta }_{kl}+{\delta }_{ij}{a}_{kl}\right)+B_3\left(a_{ik}{\delta }_{jl}+a_{il}{\delta }_{jk}+a_{jl}{\delta }_{ik}+a_{jk}{\delta }_{il}\right)$$

$$+B_4\left({\delta }_{ij}{\delta }_{kl}\right)+B_5\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jl}\right),$$

where a_ijkl is the fourth-order fiber orientation tensor, δ_ij is the second-order unit tensor, and the coefficients B are related to the components of the stiffness matrix of the transversely isotropic unidirectional composite [80]. This resulting tensor is linked to each mesh element and is exported as apdl_pre.txt to Ansys Mechanical Workbench.

Fourth, the structural analysis is carried out. In all iterations but the first, the structural analysis is performed on the current topology determined by the topology density ρ (details in [76]). Subsequently, ε is calculated, and the topology is extracted based on V_def and exported as density.topo. If the current ε is less than ε_obj for the number of times specified by the k-iterations criterion K_ε [81] (in this work, three times in a row), TO is stopped, and the last topology is exported as topoOptStruct.stl. Otherwise, the counter g is incremented by 1, and the loop repeats until convergence is achieved.

Algorithm 1. Multidisciplinary topology optimization.
Input: BCstruct, BCinjMold, MPinjMold, MPmatrix+fiber, OPtopoOpt, th, δmap, Vdef, εobj, geo.stp (Optional)
Output: topoOptStruct.stl
write materialProperties.txt ← MPfiber+matrix
write topoParameters.txt ← OPtopoOpt
meshstruct, designelements, frozenelements = AnsysWorkbench_Mesh(BCstruct, geo.stp);
write designRegionMesh.ans ← meshstruct, design.txt ← designelements, frozen.txt ← frozenelements;
g = 1;
counter_epsilon = 0;
while (counter_epsilon < Kε) do
	if g == 1 then
		reducedMesh.ans = designRegionMesh.ans;
	else
		domain_mesh_reduced = delete_elements(designRegionMesh.ans, th, (density.txt)g-1);
		write reducedMesh.ans ← domain_mesh_reduced
	end if
	meshinjMold, A = AutodeskMoldFlow(reducedMesh.ans, BCinjMold, MPinjMold);
	write meshMoldFlow.pat ← meshinjMold, fiberOrientMoldFlow.xml ← A;
	A’ = DigimatMAP(fiberOrientMoldFlow.xml, meshMoldFlow.pat, designRegionMesh.ans, δ);
	write fiberOrientAnsys.xml ← A’;
	KEL = AnisoTopo(materialProperties.txt, fiberOrientAnsys.xml, topoPararmeters.txt, designRegionMesh.ans, design.txt, frozen.txt);
	write apdl_pre.txt ← KEL;
	if g == 1 then
		Wg = AnsysWorkbench_StructuralAnalysis(designRegionMesh.ans, apdl_pre.txt, BCstruct);
	else
		Wg = AnsysWorkbench_StructuralAnalysis(designRegionMesh.ans, apdl_pre.txt, BCstruct, ρ);
		εg = |(Wg - Wg-1)/Wg-1|
		if εg <= εobj then
			counter_epsilon ++
		else
			counter_epsilon = 0
		end if
	end if
	ρ = AnsysWorkbench_TopologyOptimization_Iteration(designRegionMesh.ans, design.txt, frozen.txt, Vdef, topoPara, apdl_pre.txt, Wg)
	write density.topo ← ρ;
	Convert density.topo to (density.txt)g with HDFView();
	g++
end while
topoOptStruct = delete_elements(designRegionMesh.ans, th, (density.txt)g-1);
write topoOptStruct.stl ← topoOptStruct

2.2.2 Metric for evaluating structure design quality of composite materials

The load-carrying factor (LCF), denoted as G, along with its coefficient C_K, is employed to assess the quality of a structural arrangement [57,58]. Their typical formulations are defined as follows:

$$\mathrm{G}={\int }_{\mathrm{V}}{\mathsf{\sigma }}_{\mathrm{e}\mathrm{q}}\mathrm{d}\mathrm{V},$$

(2)

$${\mathrm{C}}_{\mathrm{K}}=\frac{\mathrm{G}}{\mathrm{F}\mathrm{l}}$$

(3)

Here, σ_eq represents the equivalent stress, V is the volume of the structure, F is the characteristic load in N, and l is the characteristic linear dimension in m (l represents the distance between areas where loads are applied to the locations of the supports). To properly evaluate the quality of the structural arrangement in composite materials, the LCF has been redefined based on stress criteria.

For structures made of isotropic materials, the LCF is expressed as

$$\mathrm{G}={\mathsf{\sigma }}^{\mathrm{U}\mathrm{T}\mathrm{S}}{\int }_{\mathrm{V}}{\mathrm{F}}_{\mathrm{e}\mathrm{q}}\mathrm{d}\mathrm{V},$$

(4)

Here, σ^UTS is the ultimate tensile stress, F_eq = σ_V / σ^UTS is the maximum stress criterion, defined as the ratio of von Mises stress σ_V to the material’s ultimate tensile stress σ^UTS. The LCF coefficient remains the same as in Equation 3.

In the case of structures made of anisotropic materials, the LCF is expressed as

$${\mathrm{G}}_{\mathrm{T}\mathrm{H}}={\mathsf{\sigma }}_0^{\mathrm{U}\mathrm{T}\mathrm{S}}{\int }_{\mathrm{V}}{\mathrm{F}}_{\mathrm{T}\mathrm{H}}\mathrm{d}\mathrm{V},$$

(5)

Here, σ₀^UTS is the ultimate tensile strength either along the longitudinal direction (along the fiber), and F_TH is the average Tsai-Hill criterion, determined using Advani-Tucker’s averaging procedure [68]. It is defined as follows:

$${\mathrm{F}}_{\mathrm{T}\mathrm{H}}=E_1{a}_{ijkl}+E_2\left(a_{ij}{\delta }_{kl}+a_{kl}{\delta }_{ij}\right)+E_3\left(a_{ik}{\delta }_{jl}+a_{il}{\delta }_{jk}+a_{jl}{\delta }_{ik}+a_{jk}{\delta }_{il}\right)+E_4\left({\delta }_{ij}{\delta }_{kl}\right)+E_5\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right),$$

(6)

where E₁ = F₁₁₁₁^ud - 2 F₁₁₂₂^ud + F₂₂₃₃^ud - 4 F₁₂₁₂^ud + F₂₃₂₃^ud; E₂ = F₁₁₂₂^ud - F₂₂₃₃^ud; E₃ = F₁₂₁₂^ud - F₂₃₂₃^ud; E₄ = F₂₂₃₃^ud; E₅ = F₂₃₂₃^ud. The values of the Tsai-Hill criteria tensor are determined using the following expression:

$$F_{TH}^{ud}=\frac{{\mathsf{\sigma }}_{11}^2}{{\mathrm{X}}^2}-\frac{{\mathsf{\sigma }}_{11}\left({\mathsf{\sigma }}_{22}+{\mathsf{\sigma }}_{33}\right)}{{\mathrm{X}}^2}+\frac{{\mathsf{\sigma }}_{22}^2+{\mathsf{\sigma }}_{33}^2}{{\mathrm{Y}}^2}+\left(\frac{1}{{\mathrm{X}}^2}-\frac{1}{{\mathrm{Y}}^2}\right){\mathsf{\sigma }}_{22}{\mathsf{\sigma }}_{33}+\frac{{\mathsf{\sigma }}_{12}^2+{\mathsf{\sigma }}_{13}^2}{{\mathrm{S}}^2}+\left(\frac{4}{{\mathrm{Y}}^2}-\frac{1}{{\mathrm{X}}^2}\right){\mathsf{\sigma }}_{23}^2,$$

(7)

Here, σ_ij represents the components of the stress tensor (component 11 corresponds to the fiber’s longitudinal axis, etc.), X is the longitudinal strength limit, Y is the transverse strength limit, S is the transverse shear strength. The LCF coefficient, accounting for the anisotropy, can be written as

$${\mathrm{C}}_{\mathrm{K}}^{\mathrm{T}\mathrm{H}}=\frac{{\mathrm{G}}_{\mathrm{T}\mathrm{H}}}{\mathrm{F}\mathrm{l}}$$

(8)

2.2.3 Brackets manufacturing and load testing technique

The two- and three-dimensional brackets were manufactured by injection molding using a Negri Bossi VE 210-1700 injection molding machine. The filling parameters were as follows: melt temperature 230 °C, mold temperature 80 °C. The pellets were dried before the injection at a temperature of 90 C for 4 h in a plastic pellet dryer. The two-dimensional brackets were cut from the injected molded plates, while the three-dimensional brackets were molded. Mold for three-dimensional brackets injection molding was manufactured from St-3 steel plates underwent CNC machining using a 4-flute 4 mm AlTiN D4x50x4Dx4F coated carbide endmill.

A static test was conducted to verify the developed methodology on both brackets. Before conducting the static test, the brackets were weighed using an electronic balance with a resolution of 0.01 g. Both brackets are loaded along the Y-axis until they fail. Restricting plates were used to ensure that the brackets failed due to fracture and not by buckling.

The mechanical testing equipment for both brackets consisted of an MTS 322 testing machine with mechanical grips, an MTS 793 controller, and an MTS 661.20F-03 force sensor, as well as special equipment shown in Figure 2, which included auxiliary rods and plates. Before performing the tensile test, the bracket support 1 and loading rod 2 were positioned. A mounting tool was used to check the collinearity between the bracket’s support and loading rod axes, and a construction level was used to ensure the vertical straightness of the loading scheme. Subsequently, the bracket was positioned between the bracket support plates and secured with M6 bolts. Displacement-restricting plates 3 were attached to the bracket support 1 using M6 bolts. Linking rod 4 is attached to the loading rod 1 with M6 bolts. The alignment of the bracket and linking-rod axes was achieved by adjusting the height of the loading rod. Finally, the bracket and linking rod were fixed in place using M6 bolts. After positioning all the components and ensuring that the bracket was not pre-stressed, the gauge sensor was set to zero.

3. Results

TO was performed on two case studies to obtain the structural layout in two- and three-dimensional spaces. In the two-dimensional case, the topologies obtained with and without considering constant anisotropy were compared. In the three-dimensional case, the solutions obtained using fixed and variable fiber orientation were compared. The materials used were PA6 50GF, PA6 30CF, and D16T. The aluminum structures, along with the LCF coefficient, serve as a control for comparing the resultant topologies.

3.1. Topology optimal two-dimensional structures

3.1.1. Topology optimization and structural arrangement quality assessment

In the two-dimensional case, the design space is a rectangle with dimensions 75 x 50 x 4 mm, as shown in Figure 3, where frozen elements near the loads and supports are marked in grey. The structured mesh consists of 47862 elements, with an edge size of 0.25 mm. This was assessed by element quality mesh metric, with a minimal value of 0.4009, a maximum value of 0.9949, an average value of 0.8834, and a standard deviation of 6.4556 x 10^-2.

The force F was applied to the elements located within the ring at the center-right side of the region, while the displacement on the elements within the rings on the left side was set to 0 mm. The design region volume percentage to retain m_ret was set to 12%. The objective relative difference ε_obj value was established as 0.1%. The material properties were defined as follows: for isotropic material - Young’s modulus of 8 GPa and a Poisson’s ratio of 0.25; for orthotropic material - the following elastic constants E_X = 13 GPa, E_Y = 7 GPa, E_Z = 6.5 GPa, υ_XY = 0.272, υ_YZ = 0.365, υ_XZ = 0.254, G_XY = 1.979 GPa, G_YZ = 1.639 GPa, G_XZ = 1.763 GPa (the X-axis correspond to the bracket symmetry axis, the Y-axis corresponds to force direction, the Z axis is determined by the right-hand rule).

The convergence results are plotted in Figure 4. A converged solution was achieved after 44 and 34 iterations for TCA (Topology Constant Anisotropy) and TCI (Topology Constant Isotropy) cases, respectively. The resultant topologies are shown in Figure 5; the topology considering isotropy and constant anisotropy will be further referred to as TCI and TCA, respectively (where C stands for constant, even though isotropy is constant by definition).

The von Mises-based and Tsai-Hill failure criterion fields of the resultant two-dimensional topologies, loaded at a force corresponding to the relation ^F/_m =70 N/gr, are presented in Figure 6. It can be observed that the von Mises failure criterion underestimates the strength of structural members loaded transverse to the fiber direction. Meanwhile, the Tsai-Hill failure criterion allows for a more accurate estimation of the strength of structures made of short-reinforced composite materials.

The failure criteria and volume of each element in equations (3) and (7) were multiplied and summed to calculate the LCF and LCF coefficients. The LCF and LCF coefficient of each topology and failure criterion combination are presented in Table 2. The relative percentage difference was calculated between the LCF coefficients of TCA and TCI.

3.1.2. Experimental verification

The TCA and TCI brackets (30 samples in total) were cut using a milling machine from plates molded from a short-reinforced composite materials PA6 50GF, PA6 30CF and, for comparison with isotropic material, from an aluminum D16T plate. Mechanical tests of the brackets were carried out on MTS 322 machine (Figure 7). The results are presented in Figure 8 in terms of specific force and normalized deformation by the characteristic dimension (see Section 2.4). The normalized specific stiffness was calculated as the ratio of F/m to /l in the elastic zone in all cases, where F is loading force, m is bracket mass, is the displacement of the bracket lug along the line of force, l is the distance from the embedment to the axis of application of the force.

Non-linear structural analysis was performed on the obtained topologies to compare the numerical and experimental results (Figure 9). The experimental lines are summarized in the form of average values over the samples and the scatter field by the value of the standard deviation. The normalized specific stiffness of the numerical results was calculated from 25 to 50% of the maximal specific force for each TCA and TCI and material combination. For experimental results, it was calculated from 2.07 to 4.12% for PA6 50GF, from 4 to 5.74% for PA6 30CF, and from 1.5 to 3% for D16T. The integral characteristics of each structure as well as their normalized specific stiffness are presented in Table 3.

The statistics of the experimentally obtained specific stiffness of TCA and TCI are shown in Table 4. The high coefficient of variation of D16T topologies results from the machining precision during manufacturing of these topologies.

3.2. Topology-optimal three-dimensional structures

3.2.1. Topology optimization and topology assessment

In Figure 10, the design region and boundary conditions for the TO are presented. The design region has dimensions of 105 x 60 x 10 mm. The optimization was carried out on a mesh of 163325 tetragonal elements, each with a size of 1.5 mm. The force F was applied to the center-right side element of the region, while the elements located at the left-side corners had a fixed displacement of 0 mm. The target volume fraction to be retained was set to 20%, and the minimum size of the structural elements was 6 mm. Anisotropic material, PA6 50 GF and PA6 30CF, as well as isotropic material D16T, were considered for the analysis.

As in the previous case of study, the resultant topology, considering anisotropy and isotropy, will be further referred to as TVA (Topology Variable Anisotropy) and TVI (Topology Variable case with Isotropy material), respectively. The convergence progress is illustrated in Figure 11, where a converged solution was achieved after 71 and 66 iterations, respectively. The resultant topologies are displayed in Figure 12. A comparison of the results in Figure 12 reveals that considering the material’s anisotropy in topology optimization results in changes to the structural layout. In Figure 12b, the connections of the main structural members become more rounded, and some rods are added in the central area of the part.

3.1.2. Topology reconstruction

The resultant TO topologies from Section 3.2.1 are idealized results and are not final geometries ready for verification and manufacturing. They were manually reconstructed in Siemens NX. Figure 13 presents the baseline and reconstructed TVA and TVI topologies.

Verification was performed by comparing the fiber orientation tensor at three points within the topology and evaluating the LCF coefficient of the topologies in different reconstruction stages. Since injection molding simulation during TO was performed on a tetrahedral mesh, this mesh was evaluated as well. The internal-exported, baseline, and reconstructed topologies, where the fiber orientation tensor was extracted, and the components of their tensors are presented in Figure 14.

The internal-exported, baseline, and reconstructed topologies were subject to lineal structural analysis using linear PA6 50GF, PA6 30CF and D16T. Figure 15 illustrates the von Mises-based and Tsai-Hill failure criteria fields for the baseline, smoothed TVA and reconstructed TVA and TVI when the topologies were subjected to loads corresponding to ^F/_M = 70 N/gr. Table 5 displays the corresponding C_K values of both topologies and provides a comparison between them.

The failure criteria F and volume dV of each element, as described in equations (3) and (7), were multiplied and summed to calculate the LCF coefficient. The LCF coefficient of each topology and material combination is presented in Table 5. The characteristic linear dimension l for the baseline topologies is 105 mm, whereas for the reconstructed topologies, it is 102.55 mm. The ultimate tensile stress σ^UTS for PA6 30CF topologies is 169.35 MPa, while for D16T topologies, it is 476 MPa. The percentage change is calculated with respect to the correspondent baseline topology.

3.2.3. Experimental validation

TVA and TVI brackets are made by injection molding (Figure 16). To strengthen the lugs, aluminium embedded elements (bushings) were added. The simulation of the injection molding process was qualitatively validated by comparing the real and simulated flow at a specific filling time (Figure 17). The used at design stage flow model makes it possible to predict the material distribution during the molding process with good accuracy.

Fiber orientation was assessed by comparing the numerical fiber orientation against the real fiber orientation, which was observed under an electronic microscope, Tescan Vega 3T, at the fracture zone. Figure 18 shows the points at which the fiber orientation was evaluated.

Mechanical tests of the TVA and TVI brackets from PA6 50GF, PA6 30CF and D16T (51 samples in total) were carried out on MTS 322 machine (Figure 19). The experimentally obtained specific force-normalized deformation curves for TVA and TVI brackets are presented in Figure 20.

Non-linear structural analysis of TVA and TVI brackets was performed on the reconstructed topologies to compare the numerical and experimental results (Figure 21). The experimental lines are summarized in the form of average values over the samples and the scatter field by the value of the standard deviation.

The normalized specific stiffness was calculated in the elastic zone in all cases with respect to normalized deformation. For the numerical results, the normalized specific stiffness ranged from 5.66 to 13.17%, while for the experimental results, it ranged from 7.5 to 15.8% for PA6 50GF and from 8.5 to 17.3% for PA5 30CF and D16T. The integral characteristics of each structure, as well as their normalized specific stiffness, are presented in Table 6.

The statistics of the experimentally obtained specific stiffness for TVA and TVI are shown in Table 7. The high coefficient of variation of D16T topologies is attributed to the machining precision during the manufacturing of these topologies.

4. Discussion

Assessment of the topologies by the LCF coefficient (see Section 3.1.1 and Section 3.1.2) reveals that the classic coefficient C_K is more suitable for structures made of isotropic materials. The classic C_K does not convey any information regarding the anisotropy, making it incapable of accurately estimating the quality of structures made from anisotropic material. For example, the classic C_K indicates the advantage of TCI made from anisotropic materials over TCA, despite the latter having lower specific stiffness (see Table 2 and Table 3). The C_K^TH incorporates the effect of anisotropy on the stress-state of the structure by optimizing the placement of structural elements to achieve more efficient mold filling. The minimal C_K^TH values of the TCA and TCI correspond to the maximal specific stiffness of these topologies. However, a significant drawback arises from the fact that the integral over the volume of the failure criteria is the product of the average criteria values over the volume, making LCF insensitive to both under and over-stressed elements. We recommend evaluating the structure’s topology not only using any of the possible formulations of the LCF coefficient but also employing a metric based on the coefficient of variance of the failure criteria (where the average and deviation are calculated from the failure criteria value over the volume). Moreover, reformulating the objective function to minimize both the average and deviation of the total strain energy of the topology should lead to a more equally strong structure.

The results obtained in Section 3 confirm that considering anisotropy during TO increases the normalized specific stiffness of the resultant topologies. Specifically, the normalized specific stiffness increased numerically by 6.06–10.18% and 0.77–0.88% in the two- and three-dimensional cases of study, respectively (see Section 3.1.2 and Section 3.2.3). The impact of considering anisotropy during TO is more pronounced in the 2D case. It can be assumed that the contrast in the normalized specific stiffness between the two- and three-dimensional cases is due to the regular orientation of the flow—and the fibers within it—along the structural elements in the three-dimensional case, causing the mechanical characteristics in these elements to be closer to those of the 0° oriented material. The increase in the change in the normalized specific stiffness of the reconstructed topologies compared to the baseline topologies (see Section 3.2.3) can be attributed to inaccuracies introduced by the CAD engineer during reconstruction and the replacement of fixed supports with bolted joints. Moreover, an increment is observed when comparing the C_K^TH of the baseline and reconstructed topologies, suggesting that C_K^TH and the normalized specific stiffness are correlated and that C_K^TH serves as an indicator of the degree of pristine (flawless) structure compared to the baseline topology. The experimental results showed that the normalized specific stiffness increased by 4.84–5.63% and 5.66–7.30% in the two- and three-dimensional cases, respectively (see Section 3.1.2 and Section 3.2.3). This indicates that two important points in this work: making the stiffness matrix dependent on the fiber orientation tensor, and obtaining the fiber orientation tensor by solving the Navier-Stokes equations along with the Folgar-Tucker’s continuity equation allow us to achieve a stiffer structure (since the 2D case does not account for these points). The deviation of the numerical from the experimental results can be attributed to differences in the loading scheme, the material model, and the omission of the weld lines during the numerical calculation. During the experiment with TVA and TVI brackets, it was noted that the poor adhesion between the aluminum bushing and the SFRP structure generated a weak joint interfacial strength, which could be the reason for the lower normalized specific stiffness than that predicted by the numerical model. Therefore, the investigation of the adhesion of aluminum and titanium alloys with SFRP can be part of future studies.

Regarding the topology, it is evident that TO attempts to orient the structural elements along the fiber direction when accounting for anisotropy, in contrast to TO in an isotropic medium, where the elements are positioned further from the symmetry line to reduce bending moments. The same topology distribution can be observed in the results of the work [63], where fiber orientation was simultaneously and sequentially considered. Therefore, is confirmed that for better structural performance of the structural elements, they should be oriented in the direction of the fiber, in other words, in the direction of maximum stiffness. Moreover, this last statement also demonstrates the minimal difference between the topologies in the three-dimensional case.

We have several recommendations for improving the repeatability of this work. First, the maximum strain rate at which the topologies are loaded during the tensile test should be 0.01 s^-1. Second, the effect of boundary conditions should be minimized to achieve higher contrast between solutions, alternatively, different boundary conditions should be investigated. For instance, increasing the design region, placing the loads further from the support, or analyzing complex loading schemes such as the geometry presented in the work [82]. Third, changing the type of material model used during TO from linear to non-linear should be considered, as previous works [83] have demonstrated that this leads to stronger structures. Finally, solving for the fiber orientation tensor during TO should be performed on a smooth surface (see Figure 14).

5. Conclusion

In this work, a methodology for obtaining topology-optimal structures made of short fiber-reinforced composites (SFRC) while considering the material’s anisotropy was presented. The calculation of the composite material’s stiffness matrix was performed using the fiber orientation tensor, which was obtained by solving the plastic molding equations for non-Newtonian fluids. Both calculations and experimental results confirmed that accounting for the material’s anisotropy can significantly enhance the rigidity of structures made from short-reinforced composite materials.

In instances where predetermined anisotropy was taken into consideration, the orientation of the reinforcing fibers led to stiffness increases ranging from 5 to 10%. In the case of molding truss-type structures, the stiffness increase was somewhat lower, ranging from 0.8 to 7.3%. This discrepancy arises from the fact that, in truss-type structures, the majority of the fibers align with the load-bearing elements of the truss, while the truss contributes less to the overall rigidity of the structure.

Furthermore, the results obtained using this technique were compared with those of topology-optimal structures. Thus, the increase in stiffness is also expressed as a percentage. In specific applications, especially within the aerospace industry, such an increase in stiffness can prove to be significantly advantageous.

This work demonstrates that the modified C_K^TH, in terms of the Tsai-Hill failure criterion, can be effectively employed in the design of fiber-reinforced polymer-based composite structures. It’s worth noting that the C_K approach can be formulated using other failure criteria as well. For instance, the Tsai–Wu failure criterion may be more appropriate when distinguishing between tension and compression strengths is crucial.

Currently, various options exist for assessing weight efficiency using dimensionless criteria. Nevertheless, this field has yet to achieve a comprehensive form, and the approach proposed in this work contributes to the advancement of dimensionless methods for assessing the quality of structures made from short-reinforced composite materials.