1. Introduction
Increasing system complexity and functional integration lead to rising product weights. As a result, the energy efficiency of the products suffers. Lightweight design plays a decisive role in solving this design challenge. To face this challenge, various lightweight design strategies are available, which can also be combined. Particularly in material lightweight design, fibre reinforced polymers are often used. Their application is motivated by their exceptional mass-specific stiffness and strength. However, due to their anisotropic behaviour, existing classical design processes must be adapted to this class of materials, since the good properties are only found when the fibres are aligned in load direction. It is essential to consider this in the design process in order to unlock the full potential of these materials. A load-adapted design can be realized by continuous fibres, which are arranged along the load paths. But continuous fibres limit the design freedom and can cause high costs (Böhlke et al. Reference Böhlke, Hrymak, Kärger, Pallicity, Weidenmann, Wood, Böhlke, Henning, Hrymak, Kärger, Weidenmann and Wood2019). For these reasons, discontinuous long fibre reinforced polymers (with a fibre length of several millimetres) are used as an alternative. Due to their flowability, even complex geometries can be manufactured, which significantly increases the design freedom. In addition, they are notably cheaper than the continuous fibre reinforced alternatives and are suitable for industrial applications due to the possibility of high production volumes. However, the component geometry influences the material flow and thus the orientation of the fibres in the manufacturing process. And, as the resulting fibre orientations are among others responsible for the linear-elastic material properties, there is a direct dependency of design and manufacturing process, which has to be taken into account for the optimal load-adapted component design (Kärger Reference Kärger, Böhlke, Henning, Hrymak, Kärger, Weidenmann and Wood2019). This leads to a complex design process, since the product developer not only has a high degree of geometric design freedom, but also needs to consider the material properties resulting from the manufacturing process (Kärger Reference Kärger, Böhlke, Henning, Hrymak, Kärger, Weidenmann and Wood2019).
For the manufacturing of discontinuous long fibre reinforced polymer structures, primary forming processes such as compression molding are used, which requires cost-intensive molds. Optimizations of the component design can later on only be realized by experimental trial-and-error and expensive mold modifications (Rios, Davis, & Gramann Reference Rios, Davis and Gramann2001; Li et al. Reference Li, Chen, Xu, Dahl, Zeng, Mirdamadi and Su2017). In order to reduce the number of mold modifications and necessary prototypes and to be able to consider the influences of the manufacturing process already in early phases of product development, simulations of the compression molding process support the product developer. As a result, characteristics that are important for the structural performance such as stiffness or warpage and shrinkage can already be taken into account in the design process. The consideration of the anisotropic material properties resulting from the manufacturing process is therefore decisive for the validity of the simulation models. The use of isotropic material properties can lead to significant result differences (Davis, Gramann, & Rios Reference Davis, Gramann and Rios2002). Models for the prediction of fibre orientations have been the subject of extensive research for a long time (Folgar & Tucker III Reference Folgar and Tucker1984), although the use of fibre reinforced polymers in load-bearing components such as in the automotive industry is relatively new (e.g., BMW i3).
To further increase the stiffness of long fibre reinforced polymer components, constructional stiffening elements such as beads can be used (European Alliance for SMC/BMC 2016). For beaded structures, the stiffening effect is not only achieved by the used fibre reinforced material, but is significantly determined by the increase in the second moment of area due to the selected bead position and bead cross section (Revfi et al. Reference Revfi, Mikus, Behdinan and Albers2020). However, the bead position and the bead cross section affect the material flow and consequently the resulting component stiffness and strength due to local differences of the fibre orientations (Sanwald et al. Reference Sanwald, Henning, Bertram, Weidenmann, Brylka, Müller and Böhlke2013).
In the 1970s, design catalogs were developed to support product developers in designing beads in steel components (Oehler & Weber Reference Oehler and Weber1972), but these design proposals were only of limited applicability for complex load cases. In the 1990s, the design catalogs were gradually replaced by numerical optimization methods, since the simulations provided significantly stiffer results (Schwarz Reference Schwarz2002). Commercial optimization software tools such as Tosca Structure offer sensitivity-based and controller-based approaches for bead optimization (Dassault Systemes 2018). The controller-based approach dates back to Albers et al. (Reference Albers, Weiler, Emmrich and Lauber2005). By transferring the anisotropic material properties from the manufacturing simulation to a shell structure model, this model can be used for bead optimization. To generate the beads, the nodes of the shell model are moved perpendicular to the shell surface. The results of such bead optimizations for a model based on local anisotropic material properties for a one-sided clamped plate under bending load are exemplarily shown in Figure 1.
Figure 1. (a) Controller-based and (b) sensitivity-based bead optimization result for a one-sided clamped plate under bending load.
As can be seen from the sharp edges and strongly varying bead cross sections in Figure 1, none of the two approaches provides directly manufacturable bead structures without additional interpretation. Moreover, geometric restrictions such as minimum radius or flank angle cannot be defined. As a result, the design of an appropriate mold is considerably challenging, since no load-adapted bead cross section is obtained that contains specific values for all bead parameters (see Figure 2). In contrast to steel components, where the maximum bead cross section is limited by the residual formability of the sheet metal and for which simulative, manufacturing-integrated optimization approaches have already been developed (Majic et al. Reference Majic, Albers, Kalmbach and Clausen2013), the product developer has almost complete design freedom when designing the mold for beaded, long fibre reinforced polymer components. But, this design freedom is not only an advantage compared to classical manufacturing processes. It comes at the price of a huge variety in the design parameters which, in their interrelated complexity, are difficult to manage by the product developer.
Figure 2. Bead cross section and bead parameters according to Emmrich (Reference Emmrich2004).
For this reason, Revfi et al. (Reference Revfi, Mikus, Behdinan and Albers2020) developed a methodology for a manufacturing-based, load-adapted bead optimization of long fibre reinforced polymer components. A genetic algorithm is used to determine the stiffest bead cross section among the large number of parameter combinations (see Figure 2), taking into account a strength constraint. This methodology has been developed for sheet molding compound (SMC) where, due to its composition and process control during compression molding, negligible residual stresses at the component level and thus no warpage and shrinkage can be assumed (Atkins, Gandy, & Gentry Reference Atkins, Gandy and Gentry1982).
This is not the case with long fibre reinforced thermoplastics (LFT), which are processed by compression molding, too. Here, macro-mechanical residual stresses and warpage occur, which should be considered in the component design (Balaji Thattaiparthasarathy et al. Reference Balaji Thattaiparthasarathy, Pillay, Ning and Vaidya2008). The residual stresses present in the LFT component have two main causes. On the one hand, residual stresses are induced during mold filling due to the uneven solidification process from the surface to the core. Adhesion prevents the solid layer at the cold mold surface from moving and the newly solidified layers experience forces by the flowing melt. The melt pressure prevents shrinkage and longitudinal forces stretch the newly solidified layers (Fan et al. Reference Fan, Yu, Zuo and Speight2017). Thus, residual compressive stresses are generated on the component surface, while the slowly solidifying core causes tensile residual stresses (Sunderland, Yu, & Månson Reference Sunderland, Yu and Månson2001). These pressure-dependent stresses remain in the frozen component but can be partially reduced if the processing temperature is still high (Parlevliet, Bersee, & Beukers Reference Parlevliet, Bersee and Beukers2006). On the other hand, temperature-induced residual stresses arise during the cooling process and their extent is determined by the thermal expansion of the used LFT and the temperature difference. The residual stresses generated by the manufacturing process lead to warpage after ejection and therefore influence the component stiffness and strength, respectively. If residual stresses in the component are oriented in the same direction as externally induced stresses, they represent an additional load and have a negative impact on the component strength. If residual stresses and externally induced stresses are oriented in opposite directions, residual stresses can prevent component failure by reducing the locally occurring stresses.
For this reason, the present work analyses the influence of residual stresses on the stiffness-optimized design of bead cross sections in beaded LFT components. This is achieved by automatically coupling the commercial simulation software Moldflow 2019 and Abaqus 2019. Consequently, it is possible to provide the linear-elastic, local anisotropic material properties resulting from the fibre orientations as well as the residual stresses from Moldflow for the calculation of the resulting stiffness and stresses in a structural simulation in Abaqus. With the help of a genetic algorithm that varies all adjustable bead parameters (see Figure 2), the volume-specific stiffest bead cross section depending on the manufacturing process is determined under consideration of a strength constraint. This result is compared with the optimization result without consideration of the manufacturing-induced residual stresses and a result with an isotropic material modelling in order to derive recommendations for the product developer.
2. Methodology
The product development process for components made of LFT materials is characterized by a strong interdependency between structural design, material design and manufacturing process design. In theory, all parameters of the three domains have to be optimized simultaneously during product development. However, this is not yet possible today. Therefore, suitable models need to be developed. These models have to be able to represent the relevant parameters for the bead optimization to be investigated, while taking into account the residual stresses resulting from the manufacturing process. Due to the described interdependencies between design, material and manufacturing, it is indispensable to identify the relevant parameters of the respective domain and to simulate them for a reliable prediction of the stiffness. For this purpose, a simulation method has been developed that consists of a coupling of individual simulation models. The coupling diagram is shown in Figure 3. The advantage of coupling different models in one method is to overcome the isolated assessment of individual parameters and thus to be able to answer complex, cross-domain questions (Albers et al. Reference Albers, Spadinger, Serf, Reichert, Heldmaier, Schulz, Bursac, Schumacher, Vietor, Fiebig, Bletzinger and Maute2017). The reason for choosing commercial (simulation) software is their good availability and high robustness.
Figure 3. Coupling diagram for the optimization of manufacturing-based beads in fibre reinforced polymers considering residual stresses from the manufacturing process.
Starting with the first step of the workflow shown in Figure 3, the component to be beaded as well as the initial charge are generated in a CAD software, for example, Creo Parametric. The initial charge is the cut to shape and stacked semi-finished product. It has to be ensured that the volume of the initial charge is slightly higher than the component volume in order to avoid a short shot which means an underfilling of the mold. In the developed workflow, the height of the initial charge is adjusted automatically so that the volume of the initial charge is 105% of the final component. Afterwards, the manufacturing of the component by compression molding is simulated as a three-dimensional (3D) model in Moldflow with a defined number of element layers over the thickness. Within Moldflow, the Navier–Stokes equations for non-Newtonian viscosity are solved to determine the flow field. The 3D modelling allows the consideration of varying flow patterns over the thickness such as fountain flow, which is typical for thermoplastic polymers. The heated initial charge solidifies when in contact with the cold mold and thus, the flow velocity is zero at the surface and no slip occurs. At the inner region the velocity remains high and the material passes the solidified outer region.
For the calculation of the fibre orientations, there are different models which are developed for mold filling simulations. In the present work, the Anisotropic Rotary Diffusion – Reduced Strain Closure model model (ARD-RSC) is chosen. This orientation model, developed by Phelps & Tucker (Reference Phelps and Tucker2009), is particularly suitable for long fibres. However, as all one-phase approaches, this model can only calculate probabilities for the fibre orientation tensor and does not calculate the motion of each fibre itself. The calculation of the anisotropic material parameters resulting from the fibre orientation is performed using the micromechanics model of Mori & Tanaka (Reference Mori and Tanaka1973). In the work of Tucker III and Liang (Tucker & Liang Reference Tucker and Liang1999) different micromechanics models were compared and the Mori–Tanaka model was identified as the most suitable for fibre reinforced polymers. The calculation of local anisotropic material properties is performed in a two-step procedure (Eom & Veinbergs Reference Eom and Veinbergs2017). First, the fibres are assumed to be unidirectionally oriented, resulting in a transversely isotropic material. In the second step, the material properties are adapted to the simulated fibre orientation and local orthotropic material parameters are calculated. Five independent material parameters are required for this procedure: 
$ {E}_1 $, 
$ {E}_2 $, 
$ {\nu}_{12} $, 
$ {\nu}_{23} $ and 
$ {G}_{12} $.
When the simulation of the compression molding is completed, the fibre orientation tensors and the mechanical material properties are exported to provide them for the structural model in Abaqus (as shown in Figure 3). The fibre orientations are used to calculate the eigenvectors, which define the local coordinate systems for each element. Thus, each element has an individual orientation and individual mechanical properties. The MpCCI MapLib of Fraunhofer SCAI is used for the transfer. The Abaqus model consists of the same number of solid layers as chosen for the mold filling simulation. After mapping on the Abaqus solid model, the parameters are further transferred to a shell model (see Figure 4).
Figure 4. Two-stage mapping exemplary for three layers in thickness direction (Revfi, Spadinger, & Albers Reference Revfi, Spadinger and Albers2019).
This procedure is necessary because the tool used for the bead generation requires the description of the structural problem as a shell model. The external loads are applied on this shell model according to the load case to be investigated in order to derive the bead positions. To ensure that the calculated stress distribution corresponds to the load case, no residual stresses are used for the structural simulation of the component in this step.
Then, the stress results are transferred to the beading tool. The tool calculates the main bending stresses and the main bending stress directions for each element. The elements with the highest main bending stress values are potential trajectory starting points. Trajectories are generated along the highest bending stresses and their directions as long as they lie in a predefined design area (see Figure 5a). The trajectories are sorted in a list according to the bending stress value at the trajectory’s starting point. Each trajectory is a possible course of a bead. The beading tool allows a definition of course restrictions such as minimum bending radius and minimum length. Furthermore, a minimum distance between beads can be set. Beginning at the top of the trajectory list, the trajectories that fit the position constraints are added to the model and represent the centreline of the beads. Afterwards, the partition of the bead is created along the centrelines on the component (see Figure 5b). For a more accurate description of the bead geometry, the local element size has to be fine enough within the bead partition. The global element edge length can be chosen coarser. The only bead parameter that affects the partition and hence the mesh of the flat component is the bead width. A variation of the other bead parameters does not change the total number of elements. After meshing is complete, the nodes within the partition are moved perpendicular to the shell surface according to piecewise defined functions to create the beads (see Figure 5c). Based on this newly created beaded shell model, a solid mesh consisting of a defined number of element layers is created by a both-sided offset using the shell model as the midsurface.
Figure 5. (a) Design area (b) bead partitions and trajectories according to Revfi et al. (Reference Revfi, Spadinger and Albers2019) (c) schematic bead generation via perpendicular nodal displacement.
During the optimization process, several hundred beaded components are generated in this way. To generate different bead cross sections, the bead parameters (see Figure 2) are varied within a specified range. The upper limit of the flank angle and the lower limit of the radii depend on the mesh size. The perpendicular nodal displacement of the beading tool (see Figure 5c) can lead to a low number of elements in the flank, if the flank is designed too steep. This results in long stretched elements of poor quality. To avoid this, the maximum flank angle has to be chosen carefully. Every change in the geometry of the bead cross section leads to a change of the component volume. Since the initial charge of the original component cannot be used for the beaded components because of the volume difference, a corresponding initial charge is created for each bead geometry. To keep the volume of the initial charge at 105% of the final component, the initial charge’s height is automatically adjusted. All other parameters concerning the initial charge remain the same. For each generated bead geometry, the mold filling process and the following structural analysis are performed as previously described for the original component without beads, but this time residual stresses are considered.
In the present work, after mold filling, the component is cooled down to 25 °C (room temperature) and the in-cavity residual stresses are exported for this temperature. A thermo-viscoelastic model (see Eq. (1)), that describes the material behaviour under thermal and mechanical loads, allows the prediction of the in-mold stress distribution before ejection of the cooled component (Fan et al. Reference Fan, Yu, Zuo and Speight2017).
$$ {\sigma}_{ij}(t)={\int}_0^t{C}_{ij rs}\left({\xi}_{(t)}-{\xi}_{\left(\tau \right)}\right)\left(\frac{{\partial \varepsilon}_{rs}}{\partial \tau }-{\alpha}_{rs}\frac{\partial T}{\partial \tau}\right) d\tau $$ 
$ {C}_{ijrs} $ is the viscoelastic relaxation modulus tensor. 
$ {\sigma}_{ij}(t) $ and 
$ {\varepsilon}_{ij}(t) $ are the stress and strain tensor, respectively, where t is the time. T is the temperature and 
$ {\alpha}_{ij} $ is the tensor of coefficients of thermal expansion. 
$ {\xi}_{(t)} $ is a pseudo-time scale defined as 
$ \xi (t)={\int}_0^t\frac{1}{a_T} d\tau $ with 
$ {a}_T $ as the time temperature shift factor.
As thermo-viscoelastic material characterization leads to complexities, Moldflow uses the following simplified thermo-viscous-elastic stress-strain relationship for orthotropic material with the mechanical properties calculated from the mold filling simulation (Fan et al. Reference Fan, Yu, Zuo and Speight2017).
$$ \left[\begin{array}{c}\begin{array}{c}{\sigma}_{xx}\\ {}{\sigma}_{yy}\end{array}\\ {}{\sigma}_{xx}\\ {}\begin{array}{c}{\tau}_{xy}\\ {}{\tau}_{yz}\\ {}{\tau}_{zx}\end{array}\end{array}\right]=\left[\begin{array}{cccccc}\frac{1-{\nu}_{yz}{\nu}_{zy}}{E_y{E}_z\Delta }& \frac{\nu_{yx}+{\nu}_{zx}{\nu}_{yz}}{E_y{E}_z\Delta }& \frac{\nu_{zx}+{\nu}_{yx}{\nu}_{zy}}{E_y{E}_z\Delta }& 0& 0& 0\\ {}\frac{\nu_{xy}+{\nu}_{xz}{\nu}_{zx}}{E_z{E}_x\Delta }& \frac{1-{\nu}_{zx}{\nu}_{xz}}{E_z{E}_x\Delta }& \frac{\nu_{zy}+{\nu}_{zx}{\nu}_{xy}}{E_z{E}_x\Delta }& 0& 0& 0\\ {}\frac{\nu_{xy}+{\nu}_{xy}{\nu}_{yz}}{E_x{E}_y\Delta }& \frac{\nu_{zy}+{\nu}_{xz}{\nu}_{yz}}{E_x{E}_y\Delta }& \frac{1-{\nu}_{xy}{\nu}_{yz}}{E_x{E}_y\Delta }& 0& 0& 0\\ {}0& 0& 0& {G}_{xy}& 0& 0\\ {}0& 0& 0& 0& {G}_{yz}& 0\\ {}0& 0& 0& 0& 0& {G}_{zx}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\varepsilon}_{xx}\\ {}{\varepsilon}_{yy}\end{array}\\ {}{\varepsilon}_{xx}\\ {}\begin{array}{c}{\gamma}_{xy}\\ {}{\gamma}_{yz}\\ {}{\gamma}_{zx}\end{array}\end{array}\right] $$
$$ \Delta =\frac{1-{\nu}_{xy}{\nu}_{yx}-{\nu}_{yz}{\nu}_{zy}-{\nu}_{zx}{\nu}_{xz}-2{\nu}_{xy}{\nu}_{yz}{\nu}_{zx}}{E_x{E}_y{E}_z} $$The residual stresses include pressure-dependent stresses after freezing and temperature-dependent stresses resulting from thermal shrinkage during the cooling process. This is described in Eq. (4) (Autodesk Help 2020).
$$ \left\{{\sigma}_g\right\}=-\left[{D}_g\right]\left\{{\varepsilon}_{g0}\right\}+\left\{{\sigma}_{g0}\right\} $$ The index g refers to the global coordinate system. [D] is the stress–strain relationship matrix, [
$ {\varepsilon}_{g0} $] is the initial strain from zero pressure state or transitional temperature to room temperature and 
$ {\sigma}_{g0} $ is the initial stress.
For the structural simulations in the present work, the in-cavity residual stresses are introduced into the Abaqus model by an initial condition which defines the local stress tensor for each element. In the first step of the structural simulation, warpage is calculated and in the second step, the external load is applied.
For the optimizations carried out in this work, a genetic algorithm is used that consecutively creates generations of individuals to determine the stiffness-optimized solution. Genetic algorithms are suitable for parallelization and allow a wide investigation of the solution space reducing the risk of converging at a local optimum which is not the global optimum. The design variables are the bead parameters: height, width, flank angle, head radius and base radius. The bead parameters of each individual are automatically checked before each optimization loop if they meet the requirements for a sufficient number of elements in every segment of the bead cross section. The first generation consists of eight individuals. Every following generation will be generated based on the previous generation. The current optimal solution is kept and recombined with the other individuals that pass the later presented strength constraint. The probability of a recombination of each parameter is 60%. All individuals that do not pass the strength constraint are replaced by a mutation, here that means a slight variation of two bead parameters, of the current optimum. If all individuals pass the constraint, the two worst solutions are replaced by mutations of the current optimum in order to explore the solution space efficiently while maintaining a population size of eight individuals. Due to the limitation of parallel executable simulations, the number of eight individuals per generation is relatively low, resulting in more optimization loops. For the comparison of bead parameter sets and the according bead cross sections, a fitness value is used as objective function which needs to be minimized to determine the parameter set resulting in the highest specific stiffness. To assess the stiffness of the component, the work of the external forces (ALLWK in Abaqus) is used as energy quantity. This has the advantage compared to the Strain Energy (ALLSE) value, which is often used in topology optimization, that the residual stresses generated by the manufacturing process are not considered in the stiffness evaluation. In addition, the component volume is integrated into the fitness function. This is due to lightweight design reasons: the optimization method results in a component with a lower mass if there are individuals showing the same stiffness. The stiffness of a component depends on the two factors, geometry and material. The geometric stiffness results from the second moment of area, whereas the material stiffness is defined by the Young’s modulus. Accordingly, for direction-independent, isotropic materials, the stiffness can only be affected by changing the second moment of area. For these kinds of materials, the optimization of the objective function results in the bead geometry with minimum radii, maximum height and width and largest possible flank angle (Revfi et al. Reference Revfi, Mikus, Behdinan and Albers2020). However, for anisotropic materials, it can happen that due to a different flow behaviour, the material stiffness partly compensates the geometric stiffness (Revfi et al. Reference Revfi, Mikus, Behdinan and Albers2020). As a result, it is possible that two different bead cross sections result in the same stiffness. Therefore, the fitness value is calculated as given in Eq. (5). The objective is to optimize stiffness but at the same time to include the volume so that the optimization results in a lightweight solution if two individuals show a similar stiffness. The definition of the fitness function can be easily adjusted to other optimization problems such as minimum weight with a required stiffness.
$$ f=\frac{\mathrm{Volume}}{\mathrm{Work}} $$Additionally, a strength constraint is implemented. Since the presented optimization methodology serves the purpose of an initial design proposal and more accurate strength investigations are carried out at a later phase of the product development process, a simple, conservative maximum stress criterion is used. The general definition for a plane stress state is given in Eq. (6).
$$ 1\ge \max \left\{\frac{\sigma_1}{-{R}_{\parallel}^c},\frac{\sigma_1}{R_{\parallel}^t},\frac{\sigma_2}{-{R}_{\perp}^c},\frac{\sigma_2}{R_{\perp}^t},\left|\frac{\tau_{12}}{R_{\mathrm{shear}}}\right|\right\} $$ By replacing the longitudinal and transverse strengths by the smaller one, the formulation can be further simplified and only compressive strength, tensile strength and shear strength are necessary material parameters. In addition, the principal stresses are used for the strength criterion. With 
$ {\sigma}_I<{\sigma}_{II}<{\sigma}_{II I} $, the strength criterion for each element is given as in Eq. (7). Again, the definition of the constraint can be easily adjusted to other criterions.
$$ 1\ge \max \left\{\frac{\sigma_I}{-{R}^c},\frac{\sigma_{III}}{R^t},\left|\frac{\sigma_{III}-{\sigma}_I}{2\times {R}_{\mathrm{shear}}}\right|\right\} $$Based on the Eqs. (5) and (7) the optimization problem for this work can be mathematically formulated as shown in Eq. (8).
$$ {\displaystyle \begin{array}{c}\min f\left(\overrightarrow{x}\right)\\ {}\mathrm{with}\;\overrightarrow{x}={\left[\mathrm{bead}\hskip0.17em \mathrm{height},\mathrm{bead}\hskip0.17em \mathrm{width},\mathrm{flank}\hskip0.17em \mathrm{angle},\mathrm{head}\hskip0.17em \mathrm{radius},\mathrm{base}\hskip0.17em \mathrm{radius}\right]}^T\\ {}\mathrm{so}\ \mathrm{that}\ \mathrm{for}\ \mathrm{each}\ \mathrm{element}\\ {}{\sigma}_I\left(\overrightarrow{x}\right)\hskip0.20em \ge -{R}^c\\ {}{\sigma}_{III}\left(\overrightarrow{x}\right)\hskip0.20em \le \hskip0.30em {R}^t\\ {}{\sigma}_{III}\left(\overrightarrow{x}\right)-\sigma \left(\overrightarrow{x}\right)\hskip0.20em \le\ \left|2\times {R}_{\mathrm{shear}}\right|\end{array}} $$3. Model design
In this work, the analysis problem is a square plate under single bending. The problem is intentionally kept simple because the paper is focussed on the method development. Therefore, a simple example helps to evaluate the plausibility of the results. The square plate has an edge length of 200 mm and a height of 2 mm. A square initial charge is generated with an edge length of 130 mm. The meshed plate and the associated initial charge, which is positioned centrally, are shown in Figure 6.
Figure 6. Meshed plate (yellow) and initial charge (orange) (Revfi et al. Reference Revfi, Mikus, Behdinan and Albers2020).
The used thermoplastic material is a long glass fibre reinforced polypropyleneFootnote 1 from the Moldflow-integrated material database. The fibre length is 10 mm and the fibre content is 30 wt.%. The material data are applied as entered in the database. The model is meshed with tetrahedral elements with an edge length of 2 mm, which follows Moldflow’s recommendation of one to two times the thickness of the component (Autodesk Support 2020). The used meshing technique is Advancing Front and six element layers over the thickness are set. Important values and parameters needed for setting up the mold filling simulation are listed in Table 1. Moreover, the isotropic material properties calculated based on the values in Table 1 are provided.
Table 1. Additional modelling settings
To determine the trajectories for the beads, the two-stage mapping procedure is used (see Figure 4). For this, the Abaqus model in this contribution consists of six solid layers. After mapping on the Abaqus solid model the parameters are further transferred to a shell model consisting of seven layers.
The applied load case is shown in Figure 7. One side is clamped and the opposite edge is loaded displacement-controlled by moving its nodes 20 mm downwards. The structural simulation of the plate to derive the trajectories is performed using square elements with linear basis functions (S4) and an edge length of 2 mm (see Figure 4). This model does not consider residual stresses to optimize the trajectory position based on the defined load case. Having calculated the trajectories, the plate is partitioned (see Figure 5b). The global element edge length is set to 1.6 mm. For a more accurate description of the bead geometry, a local element size of 0.4 mm is applied within the bead partition. The solid mesh model is generated by a both-sided offset using the beaded shell model as the midsurface. The solid models always consist of six elements over the thickness. The total number of elements (C3D8R) is about 1.3 million with about 1.5 million nodes. For the used bead generation process, beaded geometries with the same bead width have the exact same number of elements.
Figure 7. Model with boundary conditions: one-sided clamped plate (Revfi et al. Reference Revfi, Mikus, Behdinan and Albers2020).
For the warpage simulation in the first step, which is only executed for beaded plates, the displacement is only restrained in out of plane direction (Z-direction) on the clamped edge and the loaded edge, respectively. Deformation in the XY-plane is possible. Additionally, a single node at the middle of the clamped edge is fully fixed in every direction. For these boundary conditions, an exemplary deformation state before applying the external load is shown in Figure 8. The purpose of this boundary condition definition is to model a realistic mounting of the plate in an assembly that requires a flat edge for the connection with adjacent components.
Figure 8. Residual stress-induced deformation state before the application of the external load, displacement in Z-direction (scale factor: 10).
In the second step, the external load is simulated. Here, the movement of the nodes at the clamped edge is completely restrained and the nodes of the opposite edge are moved 20 mm in negative Z-direction (see Figure 7). During the optimization process, sets of bead parameters within the chosen parameter range are generated (see Table 2). For the chosen local element size of 0.4 mm, which is a good tradeoff between accuracy and simulation time for this setup, these parameter limits ensure a sufficiently precise modelling of the bead geometry and a good mesh quality.
Table 2. Bead parameter range
To define the strength values of the investigated LFT material for the strength criterion, the compressive strength is assumed to be as high as the tensile strength. The strength values are given in Table 3. As elements being directly located at the clamped side are sensitive to numerical errors and thus overpredict occurring stresses, these elements are not considered for the strength constraint.
Table 3. Strength values
4. Results and discussion
In the following, the influence of the manufacturing-induced residual stresses on the stiffness-optimized bead cross section is investigated, taking into account the described strength constraint. For this purpose, three different optimization runs were carried out. In the first optimization run, the stiffness-optimized bead cross section is determined using isotropic material properties. These isotropic material properties were calculated on the basis of the selected LFT material and are shown in Table 1. The second optimization run considers the linear-elastic, local anisotropic material properties resulting from the manufacturing process but no residual stresses. The third optimization run is the most comprehensive optimization, which considers, in addition to the second one, the residual stresses resulting from the manufacturing process.
4.1. Isotropic optimization
Table 4 shows selected individuals (optimum and individuals in its direct vicinity) of the isotropic optimization. The stiffness-optimized bead cross section has a width of 25 mm, a height of 16 mm, a flank angle of 83°, a head radius of 2 mm and a base radius of 4 mm (highlighted in bold in Table 4). The maximum compressive stress occurring in this individual is −113 MPa. For the selected individuals in Table 4, the optimizer always chooses the largest bead width, since this bead parameter is not relevant to strength in this example, but it increases stiffness. An increase of the flank angle from 83° to 84° causes a higher second moment of area and thus a stiffer geometry, which is indicated by the lower fitness value, but at the same time leads to an exceedance of the maximum permissible strength value (114 > 113 MPa). A smaller base radius (3 mm), which also leads to a larger second moment of area, also causes invalid compressive stress values (123 > 113 MPa). The results with smaller flank angles or smaller bead height show that the optimization method works as intended. All these cross sections lead to a lower second moment of area. Since the stiffening effect of beads in structures with isotropic material properties is material-independent and is determined solely by the geometrical second moment of area, the optimization method is verified by these results.
Table 4. Selected individuals of the isotropic optimization – optimum highlighted in bold
4.2. Anisotropic optimization without residual stresses
In this section, the manufacturing-based bead cross sections are optimized. To generate these results, the method shown in Figure 3 is used. However, the residual stresses resulting from the manufacturing process are not considered in this optimization. Table 5 shows selected (optimum and individuals in its direct vicinity) individuals of this optimization run. The optimum is a bead width of 25 mm, a bead height of 13 mm, an 81° flank angle, a head radius of 2 mm and a base radius of 4 mm (highlighted in bold in Table 5). All selected individuals in Table 5 again have the maximum bead width, since this bead parameter is not relevant for strength in this example. However, it is striking that all individuals only have a bead height of 13 mm and are therefore 3 mm lower than the results of the isotropic optimization. This clearly shows the influence of the manufacturing process on the optimization result. By taking into account the local anisotropic material properties, stress states arise which prevent a stiffer geometry.
Table 5. Selected individuals of the anisotropic optimization without residual stresses – optimum highlighted in bold, further investigated individuals highlighted in italics
In the following, further individuals close to the optimum are examined in order to check the plausibility of the optimization result. A reduction of the base radius of the optimum cross section to 3 mm leads to a stiffer cross section, but exceeds the permissible stress (−119 MPa). The increase of the head radius to 3 mm produces a smaller second moment of area, hence this individual has, in combination with the resulting fibre orientation, a worse fitness value. The combination of a larger flank angle of 83° and at the same time an increase of the head radius (3 mm) leads to a lower specific stiffness as well as the choice of a smaller flank angle (78° and 79°) but the same radii as with the optimum cross section.
The two combinations that are marked in italics in Table 5 are particularly interesting. These two individuals have the same fitness value up to the second decimal place, although they have a difference of 4° in the flank angle. The geometry with a flank angle of 83° requires an external work of 3763.88 mJ for the specified displacement and is therefore stiffer than the geometry with a flank angle of 79° (ALLWK = 3723.27 mJ). This is to be expected from their second moment of area. However, the larger volume required for the cross section with a flank angle of 83° causes the fitness values to equalize, as the objective function relates the stiffness to the volume. To investigate this result further, Figure 9 shows the resulting probabilities for fibre orientations in Y-direction (A22 component of the fibre orientation tensor) of the three geometries with a width of 25 mm, 13 mm height, 2 mm head radius, 4 mm base radius and 81° (optimum), 79° and 83° flank angle. The results look very similar even if they show varying fibre orientations especially in the edge-near areas on the left and the right side.
Figure 9. Fibre orientation tensor: A22 component indicating vertical fibre orientation (Y-axis). (a) 81° (optimum), (b) 79° and (c) 83°.
From this, the angle φ between the fibre preferential direction and the isotropically determined maximum (magnitude) principal stress direction σP can be calculated (see Figure 10a). The isotropic model is used as reference because it provides information on how to optimally strengthen the component depending on the applied load case. Therefore, the angle φ serves as an indicator for the extent to which the fibres contribute to the stiffening of the component, since fibres that are oriented in the principal stress direction contribute more to stiffness.
Figure 10. (a) Fibre angle dependence of the Young’s modulus according to Roos, Maile, & Seidenfuß (Reference Roos, Maile and Seidenfuß2017). (b) Number of elements with fibres in the angle difference to the first principal stress for three geometries with different flank angles.
Figure 10b shows the number of elements in the component determined in this way in the specified angle ranges for the three geometries mentioned. The models all have the same number of elements of about 1.3 million as described in Section 2. Therefore, the number of elements shown in the specified angle range in Figure 10b is directly comparable. It can be seen that the stiffest geometry with the flank angle of 81° has more elements in all stiffness-relevant angle ranges between 0° and 45° than the bead cross section with 83° flank angle. This explains why the optimization gets a stiffer individual at 81° flank angle, although the flank angle of 83° is larger, which initially means a larger second moment of area for otherwise identical bead parameters. In this case, the material stiffness arising from the resulting fibre orientations obviously exceed the additional geometric stiffness generated by a steeper flank. At the same time, the 81° flank geometry has a smaller component volume compared to the 83° flank geometry. The combination of the two factors causes the 81° flank to have a better fitness value despite lower geometric stiffness. With a very similar argumentation it can be explained why the 79° flank is as good in fitness as the 83° flank, but worse than the 81° flank angle geometry. The 81° flank has consistently more elements in the 0°–15° range (see Figure 10b), which is extremely relevant for stiffness (see Figure 10a). Only in the areas between 15° and 35°, the 79° bead cross section shows better fibre orientations. In addition, the 81° bead has a larger second moment of area so that the stiffness is maximized in relation to the volume. The fact that the 79° flank has the same fitness value as the 83° cross section up to the second decimal place can be explained analogously. The bead with 83° flank angle is stiffer (ALLWK = 3763.88 mJ) than the 79° flank geometry (ALLWK = 3723.27 mJ). This is mainly due to the larger second moment of area. However, the fibres of the 79° bead are significantly better oriented for the present load case in all angular ranges than those of the 83° flank. Thus, the material stiffness compensates partly for the lower geometric stiffness. Related to the volume of the geometries, in which the 83° geometry has a significantly larger one, this results in an identical fitness value. In contrast to the isotropic material modelling, where the increase in geometric stiffness is never compensated through a higher volume, it can be clearly stated that the integration of the manufacturing process has a significant impact on the optimized bead cross section.
Of course, this statement must be evaluated against the background that no component tests are involved, but this is not possible due to cost reasons as each individual parameter combination would require its own mold. Therefore, the plausibility check has been carried out simulatively based on the verified optimization method (see Section 4.1). To further validate these results, the Moldflow simulations should be verified with comprehensive flow studies.
4.3. Anisotropic optimization with residual stresses resulting from the manufacturing process
The influence of the fibre orientations resulting from the manufacturing process on the stiffness-optimized bead cross section has already been shown in Section 4.2. To simulate the models with residual stresses resulting from the manufacturing process in an optimization loop, the workflow as shown in Figure 3 was followed. The valid boundary conditions and time steps (warpage and external load) are described in Section 3 and have to be kept in mind when interpreting the results as they have an impact on warpage and stresses, especially at the clamped side.
Table 6 shows the optimization result (highlighted in bold) and other selected individuals of the optimization run (individuals in direct vicinity of the optimum). It is directly evident that the residual stresses resulting from the process cause the stiffness-optimized bead cross sections to reach a bead height of only 12 mm while maintaining the strength constraint. Thus, the influence of the residual stresses is obvious. Without taking them into account, the product developer defines the bead geometry as 25 mm width, 13 mm height, 81° flank angle, head radius 2 mm and base radius 4 mm (see Table 5) or even worse when using isotropic material properties (see Table 4) and would consequently deviate significantly from the optimized cross section considering the residual stresses. As a result, the design would not meet the requirements and a trial-and-error procedure in production would ensue. Not only the height, but also the 5° smaller flank angle and 1 mm smaller base radius mean that extensive modifications of the mold would become necessary.
Table 6. Selected individuals of the anisotropic optimization with residual stresses – optimum highlighted in bold
The cross sections close to the optimum are analysed in the following. When comparing the external work (ALLWK) at the end of the load step, it is noticeable that the better fitness value of the geometry with 76° flank angle (ALLWK = 3351 mJ) compared to that with 78° flank angle (ALLWK = 3352 mJ) is due to the different volumes of the two components, since their ALLWK values are practically identical. Increasing the flank angle from 76° to 77° causes the component to fail due to too high compressive stresses (−116 MPa), which is initially surprising, since all other individuals in the immediate vicinity (75° and 78°) do not fail. For this reason, Figure 11 shows the residual stress states at the base radius of the centre bead after step one before the external load is applied and, marked in red, the element with the highest stresses in the loaded state for the respective geometry. In the case of the 77° flank angle, failure occurs at this element.
Figure 11. Residual stress states per flank angle (a) 76° (optimum), (b) 77°, (c) 78° and (d) 75°, and in red the element which has the highest stress in the loaded state and leads to failure at 77°.
It is noticeable that in the failed geometry with 77° flank angle, in contrast to the other geometries, the red highlighted element is located in an area with local residual stress concentration, which is subsequently increased by the stresses caused by the external load. The reason for this is the combination of the resulting fibre orientations and the bead cross section geometry, which obviously leads to an unfavourable warpage and contributes to an increase in local stress above the failure-critical compressive stress value of 113 MPa (see Table 3). The warpage in load direction (displacement in Z-direction) resulting from the geometry, fibre orientations and residual stresses is exemplarily shown in Figure 12 for the geometry with 76° flank angle. Moreover, the values of maximum displacement for 75°, 76°, 77° and 78° flank angles are presented. It is striking that the optimum shows the smallest residual stress-related deformation. However, this observation cannot be generalized, since the failed geometry with 77° flank angle causes the second smallest deformation in Z-direction. Consequently, in the present case it is not the maximum displacement in Z-direction that is decisive for the failure of the geometry, but the local superposition of the residual stresses with the stresses caused by the external load.
Figure 12. Warpage resulting from residual stresses as displacement in Z-direction (scale factor: 10) exemplarily shown for the 76° flank angle geometry (optimum).
Summing up this discussion, an initial design draft can only be created in consideration of the manufacturing process and provides much more information than if the manufacturing influences are not taken into account (see Section 4.1).
5. Conclusion
Lightweight solutions can only reach their full potential if the lightweight aspect is consistently integrated in the development process from the beginning on. Therefore, when using LFTs, it is essential to take into account the fibre orientations and residual stresses resulting from the manufacturing process right from the concept phase.
The simulation studies in this paper show the complexity of the design process for beads in long fibre reinforced thermoplastic components. The design guidelines for beads in steel components, which are a common instrument to transfer knowledge, are not transferable to components with discontinuous fibre reinforcement. This is already shown by the comparison of the isotropic optimization results with the manufacturing-based optimization results without residual stresses. In particular, the design parameter of the bead height, which is highly relevant for component stiffness, is significantly incorrectly predicted under the assumption of isotropic material properties, that is without considering the fibre orientations that occur in the manufacturing process.
If the residual stresses induced by the manufacturing process of LFT in compression molding are additionally taken into account, an even more complex design process results. Existing design rules for beads can no longer be used for this. The dependency of the bead geometry, the fibre orientations and the residual stresses resulting from the manufacturing process cannot be handled by the product developer without support. As a result, expensive modifications of the mold can become necessary if it is discovered during prototype construction that the components do not meet the requirements. Such a mentioned support of the product developer during the initial design process is the optimization method presented in this contribution. The iterative coupling of manufacturing and structural simulation in an automated optimization loop allows load-adapted designs to be generated in early phases of the product development process.
Even if the results within one optimization loop are very close and might not be completely resolvable by simulations, the comparison of isotropic, anisotropic and anisotropic with residual stresses modelling shows clear tendencies. The authors summarize that the results show the most stiffness-significant difference in the bead height. When being compared with optimization results based on isotropic or anisotropic material modelling without the consideration of residual stresses, it can be concluded that the bead height and the flank angle should be chosen smaller for LFT components. The strength constraint limits the maximum bead height. It is shown that the flank angle is highly sensitive to fibre orientations and residual stress distribution and has to be optimized individually. This bead parameter depends decisively on the defined boundary conditions, the initial charge geometry and its position in the mold since they have a direct impact on the resulting fibre orientations. However, these insights are only valid for the chosen boundary conditions (load case, initial charge geometry and position, fibre orientations, etc.) and the defined objective function. Nevertheless, due to the flexibility of the developed method, the objective function can easily be adapted to other problem formulations like optimizing the component weight under a stiffness constraint. Doing this might include the component thickness also to be a design variable.
Future research will focus on the transfer of the method on problems of industrial complexity. Moreover, the method will be extended on thermo-mechanical component load and its effect on the stiffness-optimized bead cross section. Besides, further geometric stiffening possibilities through additional ribs in the beads are investigated. This can be achieved, for example, through coupled topology optimization, which uses the volume defined by the bead cross section as design space. In addition, the LFT material system will be further characterized and the findings will be integrated into the mold filling simulations. By conducting flow studies to validate the Moldflow simulations and thus the input for the verified optimization method, the significance and transferability of the results can be further increased.
Acknowledgements
The research documented in this manuscript has been funded by the German Research Foundation (DFG) within the International Research Training Group ‘Integrated engineering of continuous-discontinuous long fibre reinforced polymer structures’ (GRK 2078). The support by the German Research Foundation (DFG) is gratefully acknowledged.
Open access
