Technique of considering the material anisotropy in topology optimization of short fibers composite structures

The use of composite materials in the aerospace industry has been increasing in the past 20 years. The application of short fiber reinforced composite as materials for aerospace components can solve the delamination problem by substituting the laminates, reduce the weight and the cost of the component. The design of short fiber reinforced composite structures requires obtaining the optimal structure geometry which is a function of the injection molding process. The novelty of this work is the development of a topology optimization algorithm that considers the material anisotropy by simulating the injection molding process at intermediate steps of topology optimization on a reduced mesh of the design area. It is shown that the consideration of the material anisotropy within topology optimization leads to a susceptible change in the topology of the structure and made it possible to increase the stiffness of the aerospace bracket molded from short fiber reinforced polymers.


Introduction
The purpose of this work is to develop an algorithm for topology optimization of structures made of short fiber reinforced polymers (SFRP) considering the anisotropy, which is determined by the hydrodynamic effect of the injection molding process.
Some applications examples of composite materials are the following presented: in the MC-21 aircraft, 31% of the total material distribution in the airframe structure is made up of composite materials [1]. The use of polyether ether ketone (PEEK) reinforced with short carbon fibers led to a reduction of the weight and the cost of manufacturing the Airbus A350-900 door hinge assembly by 40% [2]. Besides Joo SJ. et al. designed and manufactured a front bumper assemble made of a laminate of continuous glass fiber-reinforced thermoplastic tape and discontinuous chopped long glass fibers in a thermoplastic matrix [3]. SFRP can be used to create spatially-loaded structures of complex shape including tips of lattice sturctures [4]. Other applications of SFRP are presented in [5].
The relevance of the study is confirmed by contemporary works related to solving the problems of optimizing structures made of SFRP [6][7][8][9][10]. Park S W, Choi J H and Lee B C [6] performed a topology optimization of a vehicle front body structure made of SFRP in an isotropic statement, after which a three-dimensional geometric model is built and, additionally, a parametric optimization of wall thicknesses is performed based on an anisotropic solution. Even though Weiler H and Hensel T [7] simulated the injection and considered the anisotropy of the SFRP, the shape optimization of the structure is performed using a homogenization approach, based on morphing (displacement of the nodes of the existing mesh), which makes it possible to select the rounding radii, however, does not allow make significant changes in the structure. The solution of the problem of topology optimization for SFRP is known, presented by Ospald F and Herzog R [8], nevertheless in their work the calculation of the non-stationary problem of injection composites, for the sake of simplification, is replaced by the solution of the eikonal equation. In [9,10], a technique for topology optimization of structures made of short-reinforced materials is constructed, taking into account the anisotropy of the material, but a complete solution of the problem of thermoplastic injection molding is not performed. The prediction of the distribution of fiber orientation based on the geometric characteristics of the structure performed in [9,10] does not allow taking into account the welding lines and the associated changes in the stiffness of the structural elements.
Topological optimization allows to determine the optimal material layout of structures and greatly reduce their weight [11]. The works [12][13][14][15] provide an actual panorama of the effort in developing methods for topology optimization of structures made of composites. In the work [12] presented by Safonov A. topology optimization for continuous fiber-reinforced polymers is approached, simultaneously search for optimal topology density distribution and a fiber orientation tensor in spatially composite structures. Xu Y et al. [13] presented a topology optimization of ply orientation for multiple fiber-reinforced polymers. In the work [14] Lee J et al. performed a topology optimization for discrete orientation of fibers, even though the proposed procedure lacks a solution of topology for anisotropic materials. Hoang V N et al. [15] focused on the development of topology optimization of cellular composite structures. The efficiency of topological optimization can be estimated using the dimensionless criterion presented in the work of Komarov V A [16]. Estimation of weight efficiency of short-fiber reinforced structures presented in [17].
Stiffness of short fiber composite material can be predicted by the Mori-Tanaka homogenization scheme [18,19] with Advani-Tucker orientation averaged stiffness tensor procedure [20].
The novelty of this scientific work is the development of a topology optimization algorithm that considers the material anisotropy through the Mori-Tanaka homogenization scheme for each design domain finite element with taking into account actual fiber orientation field by simulating the injection molding process on a reduced mesh of the design area in intermediate steps of topology optimization. The goal is to obtain a topology in an anisotropy statement, which manifests an for molded short-fibers reinforced structures stiffness that exceeds the stiffness obtained using an isotropic topological optimization algorithm.
The algorithm allows solving a multidisciplinary problem at steps of the topology optimization cycle: simulate the injection of a reduced mesh area, determine the actual properties of the material at each point of the structure, and calculate the stress-strain state of a body of variable density. The solution to the topology optimization problem is performed in the Ansys Workbench software in combination with Autodesk Moldflow software for the simulation of the injection molding process. The design of an SFRP bracket under in-plane loading conditions was investigated. The goal was to minimize the compliance of the bracket of defined mass within the manufacturing constraints -"Pull out direction" in both directions and "Minimum member size" of the structural elements. It is shown that the consideration of the material anisotropy leads to a susceptible change in the topology of the structure and increasing its stiffness.

Methodology
The goal of topology optimization is to find the structure layout in the design domain that satisfy the minimum compliance (maximum stiffness) for a given loading condition and various constraints such as a given amount of material, manufacturing requirement, etc. The mechanical properties of a product made of SFRP depend on the injection molding process. The objective of topology optimization performed in this work is the structure compliance minimization by minimizing the total strain energy calculated for an anisotropic composite material with respect to the topology density inside the design 3 region subject to given amount of material and constraints such as "Pull-out direction" (for paper case in both directions relative to the plane XY) and "Member size": where W -total strain energy;  -topology density;  -design region; K -stiffness matrix; u -deformation vector; F -forces vector; x-design domain elements vector with coordinates x, y, z; aij -fiber orientation tensor. Short fiber composite material stiffness can be predicted by the Mori-Tanaka model [18]. The five effective elastic constants of this model depend in fiber and matrix elasticity parameters, fiber size and volume fraction, and are calculated using the equations from [19]. Unidirectional short-fibers composite material can be considered as transversely isotropic material. To tacking into account the fiber orientation we use the Advani-Tucker orientation averaged stiffness tensor procedure with hybrid approximation of fourth-order orienatation tensor [20]. The approach presented in the works [18][19][20] implemented in AnisoTopo [21] C-code for prediction stiffnes matrix on each mesh element in APDL format.
The mechanical characteristics of a short-reinforced composite material determination at each step of topological optimization requires to know the fiber orientation tensor at each design domain element. To calculate the orientation tensor on internal iterations of topological optimization, we use Autodesk Moldflow which has been run in batch mode and managed via Synergy API. To determine the boundaries of the structure at intermediate iterations of topological optimization we use the reduced mesh, containing the design domain elements with a topological density value higher than the threshold value. Reduced mesh can be exported at intermediate iterations of topological optimization and usable for Autodesk Moldflow filling calculation. For the fiber orientation tensor extrapolation at elements of design domain outside the reduced mesh we use Digimat MAP.
Consider the steps of the topological optimization algorithm taking into account the material anisotropy in more detail (figure 1). Step 1: Algorithm initialization consists in determining the design region and boundary conditions. Design domain mesh is creating in ANSYS Workbench Mehsing and exporting in ANS format (designRegionMesh.ans) using an APDL script. The definition of unchanging topology areas is carried out by specifying files design.txt and frozen.txt with desing and frozen elements numbers and allows us to ensure the preservation of elements in the boundary conditions (fixings, loads and melt entry points) places. The designRegionMesh.ans, design.txt, frozen.txt files is placed in the AnisoTopo directory and will not change during optimization.
Step 2: Injection molding filling process is calculated at designRegionMesh.ans mesh in Autodesk Moldflow.
Step 3: Based on the calculated fiber orientation tensor, the values of the stiffness matrices in each element of the designRegionMesh.ans mesh determined in APDL format using AnisoTopo code (apdl_pre.txt, figure 2). Step 4: The technique can be implemented in two ways -by molding calculating at each intermediate iteration or by sequentially performing a series of topological optimization calculations with refined material properties. The first case allows us to quickly take into account the change in fiber orientation when changing the topology, the second way is more robust. To take into account the orientation values at each iteration, it is necessary to import stiffness matrices into the calculation tree of ANSYS Workbench Mechanical (figure 3a), while for sequential performing a series of topology optimization calculations with refined material properties, it is necessary to import stiffness matrices into the geometry tree (figure 3b).
(a) (b) Figure 3. Import of stiffness matrices for calculating SFRP anisotropic properties: (a) at each topology iteration and (b) sequential topology optimizations.
Step 5: Solving the couple task of topology optimization and composite molding requires to the topological algorithm interruption and export the reduced by topological density threshold desing space mesh at intermediate iterations of topological optimization.
In the case of successive approximations of topology optimization iteration, it is sufficient after solving the next task of topology optimization iteration, export the model in STL format. In the case of simulating the injection molding process at each step of topology optimization, the cycle is interrupted through command insertion, which starts the next iteration in the presence of a flag file ( figure 4).
Density distribution at intermediate iterations can be obtained from TOPO files (file.topo) which has Hierarchical Data Format (HDF) and exported in TXT format with the help of software for browsing and editing HDF files. To obtain a reduced mesh, an APDL script is used (figure 5), which removes elements from the mesh with a density value less than the threshold. The resulting reduced mesh is written in ANS format (reducedMesh.ans) that is available for import into Autodesk Moldflow. In the case of a not fully connected mesh, the separated from molding gate elements do not participate in the filling calculation. The fiber orientation on such elements is obtained by extrapolation.   Step 6: The reduced mesh from Ansys (reducedMesh.ans) is imported into Autodesk Moldflow. If necessary, it is made mesh refinment and removing disconnected regions. Injection molding filling process is calculated ( figure 6) and exported in PAT and XML file formats (meshMoldflow.pat and fiberMoldflow.xml), containing a refined reduced mesh and values of the fiber orientation tensor. Step 7: For the topological optimization of structures from a short-reinforced composite, taking into account its anisotropy, it is necessary to know about the fiber orientation at each element of the design domain. Digimat-MAP allows us to transfer the values of the fiber orientation tensor calculated on the Extrapolation of calculated fiber orientation tensor values is done with a user-defined tolerance ( figure 8). Increasing the tolerance increases the time it takes to extrapolate the results. The necessary requirement for the extrapolation tolerance is that the boundaries defined by this tolerance capture the elements of the design domain, which are included in the subsequent iterations of topological optimization. This approach allows, on the one hand, to take into account at molding calculation the sizes and cross-sections of structural elements at the current iteration of topological optimization, and on the other hand, to correctly and smoothly predict the mechanical characteristics of the material in adjacent elements of topological optimization. On elements outside the tolerance, the tensor of orientation of the reinforcing fibers is set corresponding to the isotropic material (1/3, 1/3, 1/3). The low value of the topological density in elements outside the tolerance eliminates the lack of information about the fiber orientation.
The result of the fiber orientation transfer is the file fiberOrientAnalisys.xml containing the fiber orientation tensor values on all elements of designRegionMesh.ans mesh at the current iteration of topological optimization, which allows us to calculate the stiffness of all design domain elements using the Mori-Tanaka and Advani-Tucker models. Step 8: After the refining the values of the fiber orientation tensor the algorithm evaluates the convergence criterion and, if necessary, continues with step 3. Solution convergence is confirmed when the residuals between the topology density value in the current iteration and the previous iteration are below a given limit (in paper case is 0.05%). If the convergence condition is fulfilled then the algorithm continues with step 9.
Step 9: Finally, when the solution converges the model is exported in STL format.

Results and discussion
The work of the algorithm for short fiber reinforced structures topological optimization, taking into account the anisotropy of the material, will be considered using the example of aerospace bracket. The fixing and loading scheme of the bracket is shown in figure 9. The design area is 105x60x10 mm. The problem is considered on a mesh of 163 325 tetragonal elements 1.5 mm in size. Topology optimization subject to 25% volume of material retainment, the "Minimum member size" equal to 5.8 mm and "Pull out direction" in both directions relative to the plane XY.
Considering the solution of the problem in the isotropic and the anisotropic formulations -for materials with characteristics corresponding to those of polyamide-6, reinforced with 50% weight glass fibe (volume fraction is 35%). The isotropic material has an elastic modulus of 8 GPa and a Poisson's ratio of 0.25. The SFRP consists of a matrix with a density of 1.4 g/cm 3 and a short glass fiber with a density of 2.6 g/cm 3   The "Sequential Convex Programming" optimization method based on a gradient-based procedure was used; iterations were stopped when the convergence accuracy was reached equal to 0.05%. Convergence plot are shown at figure 10. The results of topology optimization with and without consideration of the anisotropy are shown in figure 11. It can be seen that consideration of the anisotropy led to a change in the topology -the appearance of three diagonal structural elements within the bracket instead of two elements of a larger diameter. Optimization result exported in STL format. Volume of bracket is 14483 mm3 for isotropy statement and 14269 mm3 for anisotropy statement. The field of equivalent stresses at 1500 N force load with taking into account fiber orientation and material anisotropy (figure 12) is mainly uniform, To estimate the effectiveness of anisotropy account in topology optimization, let us compare the deformation of the brackets at 1500 N load obtained for different bracket shapes in the reference case of their manufacture from a linear isotropic material with 8GPa Young's modulus and 0.25 Poisson's ratio ( figure 13) and in case of molding from the PA-6 with 50% weight of short glass fibers, the anisotropy of which is determined by molding ( figure 14). The reference case (figure 13) shows us that if the bracet is manufactured from an isotropic material, the isotropic topology shape is not inferior to the anisotropic one. In case of the molding and anisotropy of the composite is taking into account (figure 14) anisotropy topology shape is stiffer than isotropic one by 6.55%.

Conclusion
An algorithm for topology optimization of structures manufactured by the injection molding process taking into account the anisotropy of SFRP was developed. The results of topology optimization with and without consideration of the anisotropy demonstrated that the consideration of the anisotropy led to a change in the topology of the structure. Taking into account the material anisotropy in the topological optimization algorithm made it possible to increase the stiffness of the aerospace bracket by 6.55%. During the optimization slightly improvement on the topology after a few topology iterations allowed to simulate the injection molding one time for several topology iterations. The optimization process could be further improved by finding the best ratio between numbers of injection molding and topology interations. The current optimization was performed at a constant gate position. The optimal gate position to minimize the negative effect of the weld lines in the strength of the structure is a topic of future research.
Contribution E I Kurkin -general approach, algorithms development, ANSYS Workbench topology optimization interruption, example task statement; E A Kishov -C and APDL code implementation, example task statement; O E Lukyanov -3-d modeling; O U Espinosa Barcenas -example task calculations. All authors -text of the article.