The Influence of Uncertain Loading on Topology-Optimized Designs

The design of structures using topology optimization can improve the structural performance and save material, in turn reducing costs. Using a framework of large-scale, three-dimensional topology optimization implemented by the authors in an open-source multiphysical software, we investigate the inﬂuence of uncertain loading on the optimized design. Direct diﬀerentiation is used to reveal the relationship between displacements and applied force, giving an eﬃcient and eﬀective tool to postprocess optimized topologies. The developed methodology for the assessment of the sensitivity with respect to applied forces is explored using two three-dimensional examples: the classic MBB cantilever and a cableway pylon. The advantages and limitations of this method are discussed.


Introduction
In the development of new structures and mechanical systems, the proper load collective is critical for correct design and dimensioning. Unforeseen usage, misuse and repurposing afflict the developed load collective with uncertainty and risk.
As topology optimization designs structures by finding the optimally allocated material for a specific load case or load collective, the uncertainty of loads is especially of interest. In this work, we pose the question, how "optimal" are designs with respect to the intrinsic uncertainty of load cases? Using the direct differentiation of the optimized compliance with respect to the applied loads, a framework is shown to quickly (and without further system evaluation) assess the influence of uncertain loading on the optimized design.
Topology optimization is a numerical method, which in concert with finite-element analysis, searches for the optimal allocation of material within a design domain (packaging) for a load case or a set of load cases (load collective). Topology optimization has seen recent attention as an essential tool for the Virtuous Circle of Lightweight Engineering Design [1] in which lighter systems result in lowered loading as well as reduced material use and lower costs and therefore resulting in even lighter designs. e advent of additive manufacturing has further strengthened the attention to this design tool. Additive technologies allow for a practical manufacturing method for the resulting optimized topologies, which in past required elaborate and skilled interpretation to derive manufacturable designs [2].
Uncertainty is a well-established challenge in topology optimization. To handle geometrical uncertainty associated with nano-manufacturing, a robust topology optimization method is developed in [3,4]. In [5], the latter methodology is extended to find not only geometrically robust designs, but also to improve convergence of the optimization algorithm when considering stress constraints. A robust method considering uncertain forces in a level-set framework is introduced in [6]. To consider uncertain loading direction, [7] introduces an interval-based robust topology optimization method.
Uncertainty in mechanical engineering design can be categorized into three categories: uncertain geometry, uncertain material and uncertain loading. Further, the modeling itself adds uncertainty via abstraction and simplification. e aim of this work is to assess the sensitivity of an optimum design with respect to uncertain loads and thereby giving feedback to the design engineer how susceptible the optimum design is to changes. e optimization will be carried out using Kratos Multiphysics and its Topology Optimization Application development by [8] as an extension to [9,10]. e developed methodology is based on the simple direct differentiation applied to the converged solutions and delivers meaningful results regarding the robustness with respect to the load case to the design engineer.

Topology Optimization
Topology optimization solves the design problem of the optimal material use by parameterizing the density, which is then scaled between zero and unity, i.e. no material and full material, where x i is the i th design variable. e design domain is discretized using finite elements and the density of each element is variable and therefore x i is the density of the i th element. In this work, topology optimization is implemented using the modified solid isotropic material with penalization (SIMP), first introduced by [11] and modified by [12], where E 0 is the Young's modulus for the base material. For this modification a minimal Young's modulus E min is introduced to avoid singularities in global stiffness matrix. Additionally, the penalty term p provides a more discrete solution of the problem. Within this work, a penalty term of p � 3 is used [13].

Optimization Formulation.
In topology optimization, the optimal material distribution is sought for a defined design domain and load case. For achieving efficiently stiff structures, the compliance C is minimized for a specified volume fraction v frac , resulting in the optimization formulation, h k x � 0, k � 1, . . . , n.
(3) g j (x) represents the inequality constraints of the optimization formulation, whereas h k (x) defines the equality constraints.
Compliance is defined here as the strain energy for a defined load. erefore, the minimization of compliance formulation can be interpreted as a maximization of the stiffness for a given load. e compliance of the given structure with the design variables x and the force F can be calculated as follows: Using the dependency of the displacement u and the global stiffness matrix k, we calculate the compliance on an element basis, shown here for the i th element, e element compliance is therefore calculated with the element displacement u i and the element stiffness matrix k i . e volume fraction v frac is the ratio of the material used to that of a full design domain and is defined as an inequality constraint with respect to a prescribed volume fraction v max , e modified SIMP approach is used for the element compliance in equation (6) and summing these up over the elements of the design domain based on the elemental stiffness matrix of the base material of the i th element k 0,i , the modified optimization problem reads as min

Sensitivity Analysis.
e sensitivities, i.e. gradients, of the design variables with respect to the objective and constraint functions are computed and provided to the optimization algorithm. With this information, the design variable values are updated for the next iteration. As topology optimization typically carries more design variables than objective and constraint functions, the adjoint method is used in order to efficiently calculate the sensitivities. 2 Mathematical Problems in Engineering e compliance sensitivity based on equation (6) is calculated for all elements. Applying the adjoint method and the simplifications as well as the modified SIMP approach, the objective function gradient of every element can be calculated as e constraint sensitivity too can be calculated for every element. Due to the fact that every element has a constant volume, the constraint function sensitivity can be written as: 2.3. Sensitivity Filtering. In order to avoid common problems with topology optimization, such as checkerboard patterns, different filtering methods can be used. In the case of this work, the sensitivity filtering is implemented in the optimization process in order to avoid such difficulties [14].
Checkerboard patterns refer to areas in the design domain that have full and zero density in neighboring elements, giving rise to a structure that looks like a checker-or chessboard (checkerboard patterns can be seen in Figure 1, where a too low filter radius was used for the filtering process). Checker boarding results in an unmanufacturable design as well as artificially (numerically) high stiffness [13]. e filtering of the sensitivities considering the surrounding elements alleviates this issue. Based on the value of these sensitivities and a weighting factor H i,j , the gradient of the given element is altered. e weighting factor H i,j evaluates the distance from the i th element to the j th neighboring element. e weighting of the surrounding elements of a given element can be seen in Figure 1, where the distance determines the influence of the neighboring element properties. Depending on the filter radius, a certain number of neighboring elements are considered. Mathematically, this filtering of the sensitivities is given by e weighting factor is defined by [12].
where r min represents the filter radius. With the difference of the radius value and the distance between the given element i and the element from which the sensitivity is considered (described by the index j), the factor can be calculated.

Optimization Algorithm.
e optimization problem is solved in this work using the Optimality Criteria (OC) algorithm, which is commonly used in topology optimization [15]. e OC algorithm searches iteratively for the optimal design using the sensitivities of objective and constraint functions. A Lagrange multiplier is introduced to handle the constraint function, here volume fraction. is is done by using a bisectional algorithm. e update of the design variable values is given by [16].
For this iterative process, the tuning parameters m and η are needed, where m represents the move limit of the design variable and set here to 0.2. e tuning parameter η is introduced for numerical damping purposes and is set here to 0.5 [13].
To evaluate the design variable, B i is defined from the optimality condition as: with zf/zx i being the objective function sensitivity defined as the compliance sensitivity zC/zx i and zg/zx i the constraint sensitivity, which has the value of unity for every element. e former is introduced in section 2.2 of which the filtered sensitivities are used (section 2.3). e volume constraint sensitivity is constant throughout the optimization and is simply the identity matrix for uniform meshes, as used here [17].
To determine if convergence has been reached, the relative change in the objective function is observed, For the following observations, the convergence criteria was set to ϵ stop � 10 − 6 . Further, there is a heuristic stopping criteria of maximum iterations, which was set to 200 iterations in this investigation.

Compliance Sensitivity and Robustness with respect to
Force. In this work, the variation of the force and the resulting impact on the compliance is investigated. is is carried out in a post-optimization step to assess the consequence in the applied load. By changing the force vector on the structure and optimizing the design domain, the uncertainty in load cases and the resulting optimized topologies are investigated.
Directly differentiating (5) with respect to the i th force and solving for displacement sensitivity with respect to the force vector reveals Mathematical Problems in Engineering e compliance sensitivity with respect to force is found by applying the chain rule to (4), e sensitivity of the stiffness matrix with respect to force is zero in linear-elastic analysis. e sensitivity of the force vector with respect to a single force (or multiple forces) is a vector with ones at the point of interest and zero else, denoted here as e i . e compliance sensitivity with respect to force simplifies to zC zF i � e i u + u e i , Summarizing, the compliance sensitivity with respect to force is twice the resulting displacement of the degree of freedom of the force (or forces) in question.
Although design robustness can have a wide range of understandings in engineering design, in this work it is defined in line with the following definitions: Mathematically, this is therefore defined in this work as the inverse of the sensitivity with respect to the uncertainty or noise factors, i.e. here compliance with respect to force, It should be noted that the absolute value of the compliance sensitivity is used as robustness measure of how sensitive a structure is and not how it changes. Specifically, the higher the sensitivity (gradient) of the compliance with respect to the force, the higher the resulting variation and the lower the robustness [21]. For the extreme case of infinite robustness, the sensitivity is zero while for zero robustness, the sensitivity is infinity.

Results of Numerical Investigations to Uncertain Loading
In this study, uncertain loading is explored in topology optimization with three-dimensional domains on two examples. e first example is the classic MBB cantilever, e.g. [22], in which the out-of-plane dimension is much smaller than the in-plane dimensions, resulting in a relatively thinwalled cantilever. e second example, developed in [8] finds the optimized topology of a cableway pylon.

Moderately
in-Walled MBB Cantilever. e Messerschmitt-Bölkow-Blohm (MBB) cantilever is a classic benchmark example in topology optimization. It stems from a structural problem of the company Messerschmitt-Bölkow-Blohm (now integrated into Airbus SE) in the design of a floor structure of the Airbus passenger aircraft [23]. Considering the MBB cantilever, the change of the compliance in respect of the changing force is investigated. For this investigation, a three-dimensional MBB cantilever of the dimensions 60 × 20 × 4 mm (see Figure 2) is used. e given cantilever beam is considered having a generic isotropic material with a Young's modulus of 1 Pa. e finite-element model has a regular mesh using 1 mm element length, resulting in 4800 hexahedral elements and 6408 nodes.
Two base load cases will be investigated below: (1)

MBB Cantilever under Downward Loading Condition.
e first load case foresees a downward (negative y) load. In the following, this load is varied in x, y and z and the compliance sensitivities are calculated.
To investigate the optimized design, the load is varied. In different optimization processes the force in y-direction is increased or decreased, whereas in the xand z-direction a force is added in the positive or negative direction. To avoid to large step sizes and to evaluate also small variations, the change of the force or the added load was 0.1, 0.2 and 0.3 N in both positive and negative directions. e resulting optimized topologies can be seen in Figure 3. On the other hand, when investigating the sensitivities, the step sizes where chosen smaller. Due to the fact that a sensitivity is an information of the given condition, the best way to numerically calculate these values is to choose a very low step size. Additionally, the MBB beam is thin walled in the z-direction, which could lead to a sensitive behavior and therefore a smaller step size should improve the results. Figure 3 shows the optimized compliance with respect to the force in x, y and z, respectively. For both variation of force in x and y, the relationship is nearly linear. For the force in z, the compliance increases in both negative and positive directions as this is a increase in the load. In fact this compliance sensitivity shows a quadratic form with respect to force (see Figure 4). Further of interest is the compliance sensitivity in x: the compliance increases with a compression of the structure, but with a load in positive x-direction, the compliance decreases. is means that although an increase in the magnitude of the load, the stiffness of the structure can be increased. e analytical sensitivities of the compliance with respect to the forces used in the diagrams in Figure 3 can be calculated. For this, the displacement is taken at the node where the force is acting on the structure. erefore, the displacement of the degree of freedom of the applied force. e sensitivity is defined simply by twice the displacement (20).
For the calculation of the compliance sensitivity, the displacements are considered where the structure is loaded (here a point load), e compliance sensitivities calculated via (20) are given by It should be noted that force in y-direction is negative and therefore an increase of the force is negative, which changes the sign of the sensitivity from negative to positive.
e numerical values of the sensitivities given by the optimization tool and the different force variation can be seen in Figure 5 and they are listed in Table 1. e aim of this investigation was to see how the analytical compliance sensitivities with respect to the force vary, when changing the stepsize in which the force is altered. It can be seen, that the analytical and numerical results are very close together. Only in the z-direction the sensitivities vary notably.
From the compliance sensitivity, we can then calculate the robustness of compliance with respect to A change in the defined load case, From this result and comparing it with Figure 3 as well as Figure 5, we can see the weakness of this method. e compliance sensitivity with respect to the out-of-plane force F z is orders of magnitude smaller than the other forces. In fact, the analytical sensitivity is numerically zero, showing Box and Fung (1994): Robustness means insensitivity to variation [18]. Arvidsson and Gremyr (2008): Robust Design methodology means systematic efforts to achieve insensitivity to noise factors [19]. Taguchi (1986): Robustness is the state where the technology, product, or process performance is minimally sensitive to factors causing variability [20].   no change. As can be seen in the right graph of Figure 3, this value is correct, though only locally at 0 N. is can be therefore especially problematic in the case, when there is no load present in the direction of the sensitivity. To explore this further, we now look at a second load case for which there are non-zero applied forces for the respective compliance sensitivities.

MBB Cantilever under Diagonal Loading Condition.
As the lack of out-of-plane loading shows nonlinear relationship when a force is added in this direction, a diagonal load is investigated. e optimized topology of this base load case can be seen in Figure 6. Again a load variation is applied and these results are seen in Figure 7. As this time a force in z-direction is present, the results of the variation in this direction are notably better as they are with the MBB beam in the previous example, were the single load was applied.
is time numerical and analytical results in all directions lie close together. Even though the relationship shows some nonlinearity when varying the loads up to ± 60% of the baseline, the trend is well represented.
We calculate the compliance sensitivity with respect to the force again based on the displacement.
is shows agreeance with the numerical results, cf. Figure 7. Based on the sensitivity, we find the robustness, We now see that the out-of-plane force is better represented by the analytical value, despite the relatively thinwalled nature of the design space in this direction. As expected, the compliance in out-of-plane load is orders of magnitude more sensitive with respect to in-plane force.

MBB Cantilever with Variation of the Sensitivity Filter
Radius. e filter is used on the sensitivities in order to get a better structure and to avoid i.e. checkerboard patterns (cf. Section 2.3). e size of the radius defines the area from which surrounding elements are considered (see Figure 1). e impact of the filter radius is investigated on the optimized compliance and its sensitivities with both the downward and diagonal loading conditions. is was done by optimizing the MBB beam with the different filter radii and comparing the results in terms of geometry, compliance value and compliance sensitivity.
In Figure 8, the impact of the filter radius on the optimized topology is seen for the downward loading condition. With a filter size of 1 mm, equivalent to the element size, a significant checkerboarding effect is seen. As expected, with increasing filter radius, the structure becomes more discrete. Figure 9(a) shows that this results in a significant increase in compliance. It also shows, that up to the radius size of r min � 1mm the compliance does not change. is is due to the element size, the elements used in the MBB example have the size of 1mm × 1mm × 1mm and so no neighboring element is detected by the filter. As it is expected (and can be seen in Figure 8) the compliance increases with an increasing filter radius. is can be explained by the resulting topology, as can be seen in Figure 8, the topology gets less branched and more compact, this though means that a certain degree in stiffness is lost. is can also be seen in Figure 9(b), where the compliance sensitivities w.r.t the force change over the variation of the filter radius are presented. When the filter radius is increased the topologies have growing sensitivity to force change. is leads to the assumption that a more manufacturable and discrete topology is paid for with a less stiff structure. Additionally, the structure possesses a lower robustness with respect to changes in the load.
When investigating the impact of the filter radius on the MBB beam with the diagonal load case from Figure 7, the results are very similar to the ones from the MBB beam with the singular load case. e results of the diagonal load case can be seen in Figure 10. Although the results are very   similar, the sensitivity with respect to the out-of-plane load F z is much higher, as expected. is can be explained with the thin-walled nature of the MBB beam in the z-direction. e increasing compliance with the increasing filter radius can be explained as follows: We are paying for a more continuous structure (i.e. lack of checkerboarding) with a higher compliance. With an ever larger filter radius, this phenomenon is further increased. e resulting topology is less branched and more compact. But the price is that the structure loses stiffness. e above effect can be seen analogously for the sensitivities. As seen in section 2.5, the compliance sensitivity of the structure is computed with the displacement in the degree of freedom of the load in question (or the sum therefore). With less structural stiffness, the point will be further deformed, leading to an increase in its sensitivity.

Cableway Pylon.
e investigation of the force variation was extended to the large-scale example of a cableway pylon.
Topology optimization of such tower-like structures can be performed via ground structure [24] and the above described SIMP approach. e former method has been successfully applied to tower structures, e.g. [25][26][27], though dictates a design of lattice like members. e SIMP approach has seen use for tower-like structures [26,28] and will be used here.
An example of such a cableway pylon is shown in Figure 11(a). e according parts of the cableway pylon can be seen in Figure 11(b). e model was constructed such that the area on the ground is fixed and therefore no displacements in this area are allowed. e design domain allows for the gondola to pass the pylon without any problems.
e forces onto the pylon are spread over the area on top of the pylon (see red area in Figure 11(c)), where the design domain is connected with the yoke (see Figure 11(b)). For the investigation of the cableway pylon, steel was used and the Young's modulus was set to 206900 MPa.
Due to different load cases and varying situations (gondola ascent or descent, wind in the driving direction or perpendicular to the driving direction etc.), which can occur within the lifetime of the cableway pylon, the resulting force working onto the structure has components in every direction. A simplified load case was developed, evaluating the different load cases and choosing a worst case (in every direction). ese loads were set for the optimization process in order to cover all possible load combinations:    e force variation for the investigation is realized by varying the existing forces in the main directions. For this the forces are increased or decreased by 1% or 0.5%. e variation resulted in different compliance values of the structure. In Figure 12 the resulting changes of the forces can be seen. e numerical results show excellent agreement with the analytical sensitivities. is is the case as we are not assessing the sensitivities at unloaded degrees of freedom as with the MBB cantilever. Further, no dimension can be considered relatively thin-walled, in contrast to the MBB case.
e analytical sensitivities of the optimized cableway pylon can be calculated with the displacement at the nodes where the force is working on the structure. In the case of the pylon, the force is working on 25 nodes spread over the area on the top. To get the sensitivity, the displacements in the force direction have to be added up and then multiplied by two. is results in the following sensitivities: In Figure 12 the compliance change calculated with the sensitivities is compared with the analytical results. It  can be seen that with an increase in the force, the compliance increases too. is means higher forces lead to a structure with lower stiffness. e graphs in this figure show, that the analytical and numerical results are very similar. Additionally it also shows, that in the case of the cableway pylon, which does not possess a thin-walled nature, no notable difference of analytical or numerical values can be seen. e robustness of the compliance with respect to Loading is the inverse of the compliance sensitivity resulting in 11.9242 134.3129 rough this assessment we can easily see that the resulting design for the cableway pylon is most robust with respect to a downward load, i.e. an axial load, while it is more sensitive in bending, as one may expect. erefore, especially these loads must be carefully derived in the planning of such a structure.

Conclusion
e definition of loads, i.e. the loads that a structure will see in its life, is afflicted by uncertainty. As topology optimization designs a structure based on its loading, the proper definition of load collectives is therefore essential to gain useful results. In this work, we derive a methodology to calculate the sensitivity of the resulting compliance with respect to the change in the load case without any further evaluations.
is simple, yet useful method gives the designing engineer feedback to the optimized design and the sensitivity to changes in loading.
is method though shows limitation when the structure is relatively thin walled, i.e. one dimension is of a smaller order of magnitude and unloaded in this direction. As is shown with the MBB cantilever, when unloaded, the out-ofplane force shows zero sensitivity. In the numerical studies, a new out-of-plane load greatly changes the resulting optimized topology and the sensitivity is by no means zeros. With the cableway pylon, this aspect is not a problem, showing that the thin-walled case is a special case and in need of further investigation. e investigation of the filter radius shows that with an increasing filter radius a more manufacturable topology is achieved. e more concrete the structure gets the lower the stiffness will be. Higher compliance sensitivities with respect to the load are to be expected with a higher filter radius. is means that even though the structure has an easier topology we lose structural properties. erefore a compromise has to be made when choosing the filter radius.
is method is not only applicable to complianceminimizing topology optimization resulting in stiff structures, but to the design of compliant mechanisms, which are designed to have large deformations in certain degrees of freedom (e.g. [5,29,30]). With such designs, the design of the force and therefore the compliance sensitivity with respect to force is of great concern.
In future investigations, we plan on extending the optimization with stress constraints and to quantify the uncertainty of the force when these are enforced. In this case not only the sensitivities of the compliance objective change would be observed, also the sensitivities of the stresses constraints could be investigated-so called shadow prices-giving even more information to the engineer. e robustness of force magnitude and direction are also of concern of compliant mechanisms in which the compliance is maximized.
Summarizing, this work provides a method requiring no further computational effort to assess the compliance sensitivity of topology-optimized designs. Although, weaknesses exist in regards to compliance sensitivity with respect to loads that were not considered, its is meaningful to assess the robustness, to compare sensitivities of different loads and the sensitivity of a structure to loading in general.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
A preprint of this article has previously been published in https://engrxiv.org/preprint/view/2001.

Conflicts of Interest
e authors declare that they have no conflicts of interest.