Difference-based Equivalent Static Load Method with adaptive time selection and local stiffness adaption

Abstract

Structural optimization of crash-related problems usually involves nonlinearities in geometry, material, and contact. The Equivalent Static Load (ESL) method provides a method to solve such problems. It has previously been extended to employ an individual Finite Element model describing the deformed geometry at each considered time step under the name Difference-based Equivalent Static Load (DiESL) method. This paper demonstrates how an appropriate selection of the time steps in each cycle can further improve the convergence behavior of the DiESL method. It is shown that the adaptive selection of time steps leads to better objective values and more reliable convergence to the presumed global optimum. Furthermore, the DiESL extension enables the adaption of path-dependent structural properties of the original nonlinear problem like material stiffness in each linear auxiliary load case. In this paper, an adaption of the Young’s modulus on element level in the linear auxiliary problem corresponding to the local plasticization in the nonlinear dynamic problem is successfully implemented. Here, the test examples indicate that an observable improvement can only be obtained if neither the elements in the elastic nor in the plastic range are dominating the structure’s behavior.

1 Introduction

Linear static response structural optimization is highly efficient and is therefore embedded in a considerable number of applications commonly used during the design process in industry. Many commercial codes such as MSC NASTRAN, Altair OptiStruct, or VRAND GENESIS are available in this area, enabling sizing, shape, and topology optimization in an acceptable amount of time. The time saving can especially be attributed to the availability of (semi-)analytical sensitivity analysis enabling the application of very efficient gradient-based optimizers. In contrast to the well-established optimization based on linear analysis, the real challenge in optimization is the optimization of highly nonlinear dynamic systems. A prime example is the optimization of crash-related problems in automotive industry, which is also the main objective of this paper. Furthermore, we restrict ourselves to the use of commercial crash solvers for the nonlinear dynamic analysis to ensure that the proposed DiESL method can be applied to real automotive crash problems. Worldwide, only Finite Element Analysis (FEA) explicit solvers are used in automotive crash analysis. Here, the most dominating part of the nonlinearities does not result from material or geometric nonlinearities but from contact forces. They may not only change in magnitude but especially in location from time step to time step. As a consequence it is very difficult to apply sensitivity analysis which may be the reason that no commercial code in this area is able to compute sensitivities enabling gradient-based optimization for crash. Instead, metamodel-based methods are often used. Examples for such metamodels are polynomials, radial basis functions (Powell 1992; Schaback and Wendland 2001; Schaback 2002), Neural Networks (Hornik et al. 1990, 1992; Waszczyszyn 1999), and Kriging models (Matheron 1963; Cressie 1988, 1989, 1990). The drawback of the metamodel-based optimization is that the number of design variables is limited and applications with extremely high number of design variables such as topology optimization are impossible due to the large number of required designs needed to be evaluated in order to fit the metamodels.

A promising approach to circumvent the missing sensitivity problem of commercial crash solvers is to define linear auxiliary load cases enabling linear static response optimization. This means that the nonlinear dynamic optimization problem is solved by solving a sequence of linear static response optimization sub-problems. The ESL method provides a procedure to compute such auxiliary load cases. These auxiliary load cases are created by applying ESLs to linear statics.

The ESL method has been successfully applied to various kinds of optimization like sizing, shape, free sizing, and topology optimization (Choi and Park 2002; Park et al. 2005; Lee et al. 2013b; Shin et al. 2007; Lee et al. 2007, 2013a; Jeong et al. 2008; Jang et al. 2009; Hong et al. 2010; Kim and Park 2010; Park 2011; Lee and Park 2015; Karev et al. 2018, 2019; Choi et al. 2018). Nevertheless, the ESL method has some limitations and disadvantages. The main issues result from the fact that the ESLs are always calculated based on the undeformed initial geometry. To circumvent this drawback, a difference-based approach for calculating the ESLs—the DiESL method—has been introduced recently (Triller et al. 2021). It has been shown that the DiESL method enables a significant increase in approximation quality for displacements and strains from nonlinear dynamic problems while simultaneously providing faster convergence. The DiESL method splits the nonlinear displacement path into increments, for each of which one linear auxiliary load case is created. This enables the adaption of path-dependent structural properties, e.g., material stiffness, in the auxiliary load cases. The significantly lower tangent modulus of elements exceeding yield stress during the nonlinear dynamic analysis can be adopted in the corresponding linear auxiliary load cases by adapting the Young’s modulus on element level. Such an adaptation is tested in the following using a bilinear material model. In addition to the adaptation of the material properties, the question arises whether the approximation quality of the DiESL method can be improved by an appropriate selection of the time increments. Therefore, an adaptive selection of ESL-times is tested in this paper. The intention is to place ESL-times at points in time at which nonlinear structural changes are dominant.

This paper is structured in the following way: In Sect. 2, the DIESL approach is explained. Furthermore, the algorithm for adaptively selecting the ESL-times and the local adaption of the Young’s moduli is presented. In Sect. 3 both approaches are tested numerically using two different sizing optimization examples. Finally, a conclusion and an outlook is worked out.

2 The DiESL (Difference-based Equivalent Static Loads) method

The basic procedure of the DiESL method is very similar to the ESL method, introduced by Choi and Park (2002). The main idea is to create linear auxiliary load cases for nonlinear dynamic response optimization problems, which are used in optimization based on linear statics. This approach is advantageous for many reasons: Firstly, the missing sensitivity problem is solved. Secondly, well-developed commercial software systems can be used for analysis and optimization and no development of an own sensitivity analysis and optimization algorithm is necessary.

The original nonlinear dynamic response optimization problem can be stated as

$$\min f\left( {{\mathbf{x}},{\mathbf{u}}\left( {{\mathbf{x}},t} \right)} \right)$$

(1a)

$${\text{s.t.}}\quad g_{j} \left( {{\mathbf{x}},{\mathbf{u}}\left( {{\mathbf{x}},t} \right)} \right) \le 0;\quad j = 1, \ldots , m$$

(1b)

$${\mathbf{x}}^{{\text{L}}} \le {\mathbf{x}} \le {\mathbf{x}}^{{\text{U}}} ;\quad {\mathbf{x}} \in {\mathbb{R}}^{n} .$$

(1c)

Here, $f\left( {\mathbf{x}} \right)$ is the objective function, $m$ the number of constraints $g_{j}$, and ${\mathbf{x}}^{{\text{L}}}$ and ${\mathbf{x}}^{{\text{U}}}$ the lower and upper bounds of the design variables ${\mathbf{x}}$, respectively. The displacement vector ${\mathbf{u}}^{T} \left( t \right) = { }\left( {{\mathbf{u}}_{1}^{T} \left( t \right),{\mathbf{u}}_{2}^{T} \left( t \right),{ } \ldots ,{ }{\mathbf{u}}_{{{\text{n}}_{{\text{N}}} }}^{T} \left( t \right)} \right)$ is the solution of

$${\mathbf{M}}_{{{\text{NL}}}} \left( {\mathbf{x}} \right){\mathbf{\ddot{u}}}\left( t \right) + {\mathbf{C}}_{{{\text{NL}}}} \left( {\mathbf{x}} \right){\dot{\mathbf{u}}}\left( t \right) + {\mathbf{K}}_{{{\text{NL}}}} \left( {{\mathbf{x}},{ }{\mathbf{u}}\left( t \right)} \right){\mathbf{u}}\left( t \right) = {\mathbf{f}}\left( t \right)$$

(1d)

and contains the nonlinear displacements of all $n_{{\text{N}}}$ nodes at the time $t$. The general procedure of the DiESL method to create the linear auxiliary load cases is the following (Fig. 1): First, a nonlinear dynamic analysis is performed. For specified time steps $t^{i}$, $i = 1, \ldots , n_{{\text{T}}}$ the displacement fields ${\mathbf{u}}\left( {t^{i} } \right)$ are derived from nonlinear dynamic analysis afterward. The basic difference between ESL and DiESL is depicted in Fig. 2: the standard ESL method uses the undeformed geometry to compute the nodal displacements ${\mathbf{u}}^{i}$ for any given time $t^{i}$ (Fig. 2, middle). Consequently it falls short of following the nonlinear displacement path (Fig. 2, left). In contrast, DiESL is able to follow the nonlinear displacement path by splitting it into linear increments. Each increment ${\Delta }{\mathbf{u}}\left( {t^{i} } \right)$ is computed based on the linear submodel at time $t^{i}$ (Fig. 2, right), which we call ${\text{LSM}}^{i}$ in the following. The i-th ${\text{LSM}}^{i}$ is defined by the coordinates of all nodes at time step $t^{i}$. The vector

$${\mathbf{r}}^{T} \left( {{\text{t}}^{i} } \right) = { }\left( {{\mathbf{r}}_{1}^{T} \left( {t^{i} } \right),{ }{\mathbf{r}}_{2}^{T} \left( {t^{i} } \right),{ } \ldots ,{\mathbf{r}}_{{n_{{\text{N}}} }}^{T} \left( {t^{i} } \right)} \right)$$

(2)

contains the coordinates ${\mathbf{r}}_{j} \left( {t^{i} } \right)$ of all $n_{{\text{N}}}$ nodes. Accordingly, ${\mathbf{r}}\left( {t^{0} } \right)$ describes the coordinates of all nodes of the undeformed model.

Fig. 1

Optimization process of the DiESL method for a nonlinear dynamic problem

Fig. 2

Displacement path of an arbitrary node during the deformation of a structure (left) and the corresponding displacement ${\mathbf{u}}\left( {t^{i} } \right)$ (middle) and ${\Delta }{\mathbf{u}}\left( {t^{i} } \right)$ (right) used for the computation of the ESLs and DiESLs at time steps $t^{i}$, respectively.

The coordinates of a linear submodel ${\text{LSM}}^{i}$ describing the deformed geometry at time $t^{i}$ can therefore be calculated by

$${\mathbf{r}}\left( {t^{i} } \right) = {\mathbf{r}}\left( {t^{0} } \right) + {\mathbf{u}}\left( {t^{i} } \right).$$

(3)

The incremental nonlinear displacements leading from ${\mathbf{r}}\left( {t^{i} } \right)$ to ${\mathbf{r}}\left( {t^{i + 1} } \right)$ are calculated by

$${\Delta }{\mathbf{u}}\left( {t^{i} } \right) = {\mathbf{u}}\left( {t^{i + 1} } \right) - {\mathbf{u}}\left( {t^{i} } \right).$$

(4)

Note that all ${\text{LSM}}^{i}$ have the same mesh topology and that they differ only in the coordinates ${\mathbf{r}}\left( {t^{i} } \right)$. In each ${\text{LSM}}^{i}$ we then calculate the loads ${\Delta }{\mathbf{f}}_{{{\text{DiESL}}}}^{i}$ yielding the incremental displacement ${\Delta }{\mathbf{u}}\left( {t^{i} } \right)$ in linear statics leading from the structure ${\mathbf{r}}\left( {t^{i} } \right)$ to the subsequent deformed structure ${\mathbf{r}}\left( {t^{i + 1} } \right)$.

To determine the incremental equivalent static loads ${\Delta }{\mathbf{f}}_{{{\text{DiESL}}}}^{i}$ in the linear submodel ${\text{LSM}}^{i}$, the corresponding stiffness matrix ${\mathbf{K}}^{i} = {\mathbf{K}}\left( {{\mathbf{x}},{\mathbf{r}}\left( {t^{i} } \right)} \right)$, which depends on the design variables ${\mathbf{x}}$ and the nodal coordinates ${\mathbf{r}}\left( {t^{i} } \right)$ of ${\text{LSM}}^{i}$, has to be multiplied by the vector of the incremental nonlinear displacements ${\Delta }{\mathbf{u}}\left( {t^{i} } \right)$

$${\mathbf{K}}^{i} {\Delta }{\mathbf{u}}\left( {t^{i} } \right) = { }\Delta {\mathbf{f}}_{{{\text{DiESL}}}}^{i} { };{ }i = 0,{ } \ldots ,{ }n_{{\text{T}}} - 1.$$

(5)

Using the incremental equivalent static loads $\Delta {\mathbf{f}}_{{{\text{DiESL}}}}^{i}$, gradient-based linear static response optimization can now be performed. In contrast to the ESL method, which requires only a single FE model representing the undeformed structure of the initial model with the coordinates ${\mathbf{r}}\left( {t^{0} } \right)$, $n_{{\text{T}}}$ FE models have to be considered in one optimization run in the DiESL method. Consequently, Multi-Model Optimization (MMO) has to be applied where more than one FE model is taken into account simultaneously in one linear static response optimization run. Here, the FE equation of linear statics

$${\mathbf{K}}^{i} {\Delta }{\mathbf{u}}^{i} = { }\Delta {\mathbf{f}}_{{{\text{DiESL}}}}^{i}$$

(6)

is solved for each ${\text{LSM}}^{i}$, which yields the incremental linear displacements ${{\Delta }}{\mathbf{u}}^{{iT}} = \left( {{{\Delta }}{\mathbf{u}}_{1}^{{iT}} ,{{\Delta }}{\mathbf{u}}_{2}^{{{\text{i}}T}} , \ldots ,{{\Delta }}{\mathbf{u}}_{{n_{{\text{N}}} }}^{{iT}} } \right)$. The displacement fields ${\Delta }{\mathbf{u}}\left( {t^{i} } \right)$ and ${\Delta }{\mathbf{u}}^{i}$ of the nonlinear and the linear analysis are identical only at the beginning of the linear static response optimization (inner loop). At the end of the inner loop the linear and the nonlinear dynamic responses no longer match because the linear auxiliary load cases are only an approximation of the nonlinear behavior of the system. Therefore, nonlinear dynamic analysis has to be applied again with the updated design variables. If the difference is too high the process is iterated (outer loop) until the difference is small enough and additional termination criteria (Shin et al. 2007; Jeong et al. 2010; Kim and Park 2010; Park 2011) are fulfilled. It can be expected that the difference between the linear and nonlinear dynamic responses increases with the length of the inner loop optimization path. If this difference is too high, the update in the outer loop may result in huge changes in the search directions and convergence is slowed down or even not guaranteed anymore. For this reason commercial ESL codes, such as VRAND ESLDYNA, limit the number of iterations and offer the application of move limits to restrict the length of the inner loop optimization path in each cycle. In the following the outer loop iterations are called cycles to distinguish them from the iterations of the inner loop in the linear static response sub-problem.

Summarizing, the DiESL approach solves the following optimization problem:

$$\min f\left( {{\mathbf{x}},{\Delta }{\mathbf{u}}^{0} \left( {\mathbf{x}} \right), \ldots , {\Delta }{\mathbf{u}}^{{n_{{\text{T}}} - 1}} \left( {\mathbf{x}} \right)} \right)$$

(7a)

$${\text{s.t.}}\,\quad g_{j} \left( {{\mathbf{x}},{\Delta }{\mathbf{u}}^{0} \left( {\mathbf{x}} \right), \ldots , {\Delta }{\mathbf{u}}^{{n_{{\text{T}}} - 1}} \left( {\mathbf{x}} \right)} \right) \le 0;\quad j = 1, \ldots , m$$

(7b)

$${\mathbf{x}}^{{\text{L}}} \le {\mathbf{x}} \le {\mathbf{x}}^{{\text{U}}} ;\quad {\mathbf{x}} \in {\mathbb{R}}^{n} ,$$

(7c)

where ${\Delta }{\mathbf{u}}^{i}$ is the solution of

$${\mathbf{K}}\left( {{\mathbf{x}},{\mathbf{r}}\left( {t^{i} } \right)} \right){\Delta }{\mathbf{u}}^{i} = { \Delta }{\mathbf{f}}_{{{\text{DiESL}}}}^{i} { };{ }i = 0,{ } \ldots ,{ }n_{{\text{T}}} - 1$$

(7d)

in ${\text{LSM}}^{i}$, the difference-based equivalent static loads ${\Delta }{\mathbf{f}}_{{{\text{DiESL}}}}^{i}$ are computed as

$${\Delta }{\mathbf{f}}_{{{\text{DiESL}}}}^{i} = {\mathbf{K}}\left( {{\mathbf{x}},{\mathbf{r}}\left( {t^{i} } \right)} \right){\Delta }{\mathbf{u}}\left( {t^{i} } \right){ };{ }i = 0,{ } \ldots ,{ }n_{{\text{T}}} - 1$$

(7e)

with ${\Delta }{\mathbf{u}}\left( {t^{i} } \right) = {\mathbf{u}}\left( {t^{i + 1} } \right) - {\mathbf{u}}\left( {t^{i} } \right)$ in which ${\mathbf{u}}\left( {t^{i} } \right)$ is the solution of

$${\mathbf{M}}_{{{\text{NL}}}} \left( {\mathbf{x}} \right){\mathbf{\ddot{u}}}\left( t \right) + {\mathbf{C}}_{{{\text{NL}}}} \left( {\mathbf{x}} \right){\mathbf{\dot{u}}}\left( t \right) + {\mathbf{K}}_{{{\text{NL}}}} \left( {{\mathbf{x}}, {\mathbf{u}}\left( t \right)} \right){\mathbf{u}}\left( t \right) = {\mathbf{f}}\left( t \right)$$

(7f)

for the selected time steps $t^{i}$. The objective function and the constraints used in the linear static response optimization and in the nonlinear dynamic response optimization are identical. The only difference is that the nonlinear dynamic responses are approximated by linear static responses. For dynamic responses such as section forces, velocities, and accelerations, which are not defined in linear statics, it is necessary to find proper approximations. For velocity and acceleration, for instance, simple finite forward (Jeong et al. 2010; Triller et al. 2021) or central (Karev et al. 2018, 2019) differences between adjacent load cases (corresponding to adjacent time steps) can be used.

2.1 Computation of displacements

During optimization of all LSMs using MMO, each LSM analysis yields incremental displacements. The total linear displacement of a node at time $t^{i}$ is used as an approximation of the respective nonlinear displacement ${\mathbf{u}}^{i} \approx {\mathbf{u}}\left( {t^{i} } \right)$. It can be computed recursively as

$${\mathbf{u}}^{i} = {\mathbf{u}}^{i - 1} + {\Delta }{\mathbf{u}}^{i - 1} ,$$

(8)

where ${\Delta }{\mathbf{u}}^{i - 1}$ is the solution of (7) or (8b) for ${\text{LSM}}^{i - 1}$. Accordingly, the cumulated displacements can be calculated as

$${\mathbf{u}}^{i} = { }{\mathbf{u}}^{0} + \mathop \sum \limits_{j = 0}^{i - 1} {\Delta }{\mathbf{u}}^{j} .$$

(9)

In general ${\mathbf{u}}^{0} = 0$ applies. The accumulation is processed by the MMO, which handles all LSMs and accumulates their solutions ${\Delta }{\mathbf{u}}^{i}$.

2.2 Computation of adaptive ESL-times

Instead of using equidistant time steps (Triller et al. 2021), we now propose how to select the value of all $n_{{\text{T}}}$ ESL-times ${\mathbf{t}}$ adaptively in each cycle. The ESL-times should be placed at points in time at which nonlinearities are dominant. Therefore, a time-dependent function indicating the nonlinearities is needed. In this paper, we focus on crash problems, where a rigid impactor collides with a structure. In this case the contact force ${\text{f}}$ may be used to detect nonlinearities. For example in Fig. 3 the contact force curve and three different states of deformation for a crash box are illustrated (refer to Sect. 3.2 for a description of the model). The oscillating course of the contact force reflects the creation of plastic hinges and subsequent contacting.

Fig. 3

Deformation of a crash box model at three different times and corresponding contact force curve between impactor and crash box (times of deformations indicated with red circles)

Thus, it can be expected that especially at the extrema of the contact force curve nonlinearities are dominating. The idea is now to fit the contact force curve $f\left( t \right)$ by a polygonal line $l\left( t \right)$ with ESL-times ${\mathbf{t}} = \left\{ {t_{0} , t_{1} , t_{2} , \ldots , t_{{n_{{\text{T}}} }} } \right\}$ as breakpoints, directly after the nonlinear dynamic analysis in each cycle. To do so, the results of the nonlinear dynamic analysis (i.e., contact force and displacement field) must be stored for all $n_{{\uptau }} + 1$ possible ESL-times ${{\varvec{\uptau}}} = \left\{ {\tau_{0} ,\tau_{0} + \Delta \tau ,{ }\tau_{0} + 2\Delta \tau , \ldots ,\tau_{{n_{{\uptau }} }} - \Delta \tau ,\tau_{{n_{{\uptau }} }} } \right\}$. Due to storage requirements $n_{{\uptau }}$ should be limited. In this publication, ${\Delta }\tau = 1 {\text{ms}}$ is selected leading to $n_{\tau } = 100$ if the overall simulation time is 100. Based on the stored data, we approximate the function $f\left( t \right)$ as a piecewise linear function:

$$f\left( t \right) = f_{i - 1} + \frac{{f_{i} - f_{i - 1} }}{{\tau_{i} - \tau_{i - 1} }}\left( {t - \tau_{i - 1} } \right);{ }\tau_{i - 1} \le t < \tau_{i} ;i = 1, \ldots ,n_{{\uptau }}.$$

(10)

The number of ESL-times $n_{{\text{T}}}$ should also be limited, because each ESL-time corresponds to an LSM and each LSM increases the computational effort in the design domain (inner loop). In case of limited resources this can become prohibitive for very large models. However, commercial optimizers like OptiStruct employ (MPI) in their MMO implementation which permits the analysis of each LSM on an individual host. Consequently, the set of ESL-times ${\mathbf{t}}$ is a sub-set of all possible ESL-times ${{\varvec{\uptau}}}$. To determine the optimal selection of ${\mathbf{t}}$, the piecewise linear function $l\left( t \right)$, defined by the ESL-times ${\mathbf{t}}$

$$l\left( t \right) = f\left( {t_{j - 1} } \right) + \frac{{f\left( {t_{j} } \right) - f\left( {t_{j - 1} } \right)}}{{t_{j} - t_{j - 1} }}\left( {t - t_{j - 1} } \right);{ }t_{j - 1} < t \le t_{j} ; j = 1, \ldots ,n_{{\text{T}}}$$

(11)

is fitted to the piecewise linear function $f\left( t \right)$ defined by ${{\varvec{\uptau}}}$ (Fig. 4). Therefore, we minimize the sum of squared residuals between $l\left( t \right)$ and $f\left( t \right)$:

$${\text{min }}SSR\left( {\mathbf{x}} \right); SSR\left( {\mathbf{x}} \right) = \frac{1}{{ n_{{\uptau }} }}\mathop \sum \limits_{i = 0}^{{n_{{\uptau }} }} \left( {f\left( {{\uptau }_{i} } \right) - l\left( {{\uptau }_{i} } \right)} \right)^{2} ; {\mathbf{x}} \in {\mathbb{R}}^{{n_{{\text{D}}} }}$$

(12a)

$${\text{s.t.}}\quad \left| {t_{{i,{\mathbf{fix}}}} - x_{j} } \right| \ge \Delta \tau ;\quad j = 1, \ldots , n_{{\text{D}}} ,\;i = 0, \ldots , n_{{{\text{fix}}}} ,$$

(12b)

$$x_{j} - x_{j - 1} \ge \Delta \tau ;\quad j = 1, \ldots , n_{{\text{D}}}$$

(12c)

Fig. 4

Piecewise linear fit of contact force curve $f\left( t \right)$ by polygonal line with ESL-times ${\mathbf{t}}^{*}$ as breakpoints before (left) and after (right) optimization

where the number $n_{{\text{D}}}$ of discrete design variables ${\mathbf{x}}$ is not equal $n_{{\text{T}}}$ because we allow the definition of $n_{{{\text{fix}}}}$ prescribed ESL-times ${\mathbf{t}}_{{{\text{fix}}}} = \left\{ {t_{{0,{\text{fix}}}} = \tau_{0} , \ldots ,t_{{n_{{{\text{fix}}}} }} = t_{{n_{{\text{T}}} }} } \right\}$. The latter must contain the first and the last ESL-time but could also contain ESL-times in between, e.g., the maximum of $f\left( t \right)$. Consequently, the number of discrete design variables ${\mathbf{x}}$ is given as

$$n_{{\text{D}}} = n_{{\text{T}}} + 1 - n_{{{\text{fix}}}}.$$

(13)

Note that this optimization problem must not be mixed up with the inner loop in Fig. 1. It has to be applied between nonlinear analysis and the computation of ESLs (step 7 in overall flow scheme in chapter 2.4).

With this approach times corresponding to certain curve characteristics can be prescribed, e.g., the maximum contact force. Furthermore, the last ESL-time $t_{{n_{{\text{T}}} }}$ can be set such that only the relevant part of the deformation process is captured, which for example can be specified by the time of maximum intrusion. In this case, the sum of squared residuals ($SSR)$ in Eq. (12a) is modified as follows: The value $n_{{\uptau }}$ is replaced by $\widetilde{{n_{{\uptau }} }}$ such that $\tau_{{\widetilde{{n_{{\uptau }} }}}} = t_{{n_{{\text{T}}} }}$. To ensure that each possible ESL-time $\tau_{i}$ is only taken once, the distance between two adjacent ESL-times or design variables ${\mathbf{x}}$ must be at least the time increment $\Delta \tau$. This optimization problem is solved using sequential least square programming in Python. To reduce computational effort, the original discrete optimization problem is solved by the relaxed continuous auxiliary problem and the results are rounded afterward according to the stored data at the times ${{\varvec{\uptau}}}$:

$$x_{{{\text{r}},i}}^{*} = \frac{{x_{i}^{*} }}{\Delta \tau } + 0.5\Delta \tau ; i = 1, \ldots , n_{{\text{D}}}.$$

(14)

The fitted ESL-times ${\mathbf{t}}^{*}$ thus follow from the rounded results $x_{{{\text{r}},i}}^{*}$ and all prescribed times ${\mathbf{t}}_{{{\text{fix}}}}$. Furthermore, the initial values ${\mathbf{x}}^{\left( 0 \right)}$ are determined by a successive removal strategy: the number of design variables $n_{{\text{D}}}$ is initially set to the largest value possible such that each of the breakpoints $t_{j}$ of $l\left( t \right)$ is located on a time $\tau_{i}$. In this case, $SSR\left( {\mathbf{x}} \right) = 0$. Then one breakpoint is removed, namely the breakpoint with the smallest impact on SSR. This process of removing one breakpoint is successively repeated until the desired number of design variables $n_{{\text{D}}}$ is obtained. This strategy of successive breakpoint removal can be beneficial especially for contact force curves where the extremes are unevenly distributed. In this case, the likelihood of the optimizer terminating in a local optimum decreases. This advantage is illustrated in Fig. 5, where $SSR\left( {{\mathbf{x}}_{{\text{r}}}^{*} } \right)$ is given for different numbers of ESL-times $n_{{\text{T}}}$: initialing optimization problem (13a) with equidistant breakpoints fails to find the optimum solution for all shown values of $n_{{\text{T}}}$. It is important to note that the number of ESL-times $n_{{\text{T}}}$ remains constant throughout the entire optimization (outer loop) in this implementation.

Fig. 5

$SSR\left( {{\mathbf{x}}_{{\text{r}}}^{*} } \right)$ over $n_{{\text{T}}}$ for different choices of initial designs ${\mathbf{x}}^{\left( 0 \right)}$: equidistant versus successive breakpoint removal

2.3 Local adaption of Young’s modulus

The DiESL approach also enables the adaptation of path-dependent structural properties in the LSMs corresponding to the associated increment or times. In the course of the nonlinear dynamic analysis, plasticization occurs, for example, and the stiffness of some elements is drastically re. To adopt stiffness changes to the linear auxiliary problem, a bilinear material model is introduced to the design domain, defined by the Young’s modulus ${\text{E}}$ and the tangent modulus ${\text{E}}_{{\text{T}}}$. Depending on the existence of a plastic strain ${\upvarepsilon }_{{j,{\text{pl}}}} \left( {t^{i} } \right)$ within an element in the analysis domain, the element’s material property is adapted to a smaller value ${\text{E}}_{{\text{T}}}$ in the corresponding ${\text{LSM}}^{i}$. Additionally, the Poisson’s ratio is set to $0.49$ to approximate the incompressibility in the plasticized area. This is a very simplified and crude approach. It cannot describe elements that start in the elastic domain and enter the plastic domain within one time increment. Nevertheless, it should be sufficient to study the beneficial effects on the approximation quality of DiESL. Following this approach the stiffness matrix

$${\hat{\mathbf{K}}}^{i} = {\mathbf{K}}\left( {{\mathbf{x}},{\mathbf{r}}\left( {t^{i} } \right),{ }{{\varvec{\upvarepsilon}}}_{{{\text{pl}}}} \left( {t^{i} } \right)} \right)$$

(15)

now also depends on the vector

$${{\varvec{\upvarepsilon}}}_{{{\text{pl}}}}^{T} \left( {t^{i} } \right){ } = { }\left( {{\upvarepsilon }_{{1,{\text{pl}}}}^{T} \left( {t^{i} } \right),{\upvarepsilon }_{{2,{\text{pl}}}}^{T} \left( {t^{i} } \right),{ } \ldots ,{\upvarepsilon }_{{n_{{\text{E}}} ,{\text{pl}}}}^{T} \left( {t^{i} } \right)} \right)$$

(16)

containing the plastic strains of all $n_{{\text{E}}}$ elements in the analysis domain. In Fig. 6 this procedure is illustrated for a deformed state of the side impact model, which will be introduced in Sect. 3. The determination of the ESL-times as well as the local adaption of Young’s moduli has to be performed right before the computation of ${\Delta }{\mathbf{f}}_{{{\text{DiESL}}}}^{i}$. The computation of the DiESLs and the linear static response optimization then is performed using the adopted stiffness matrices ${\hat{\mathbf{K}}}^{i}$ instead of ${\mathbf{K}}^{i}$.

Fig. 6

Plastic strain of nonlinear dynamic analysis of the side impact (Sect. 3) for all $x_{i} = 0.8{\text{ mm}}$ at time $t = 82{\text{ ms}}$ and corresponding linear submodel

2.4 Implementation

The overall program flow of the DiESL method has been implemented in Python. It uses the commercial solver LSTC LS-DYNA (LSTC 2015) for nonlinear dynamic analysis and OptiStruct (HyperWorks 2019) for the computation of the ESLs and for the linear static response optimization. Two termination criteria have been used. The first criterion is that the design has to be feasible. Since the linear static sub-problem only provides estimations for the actual problem, the optimized responses of the linear static analysis often differ from those of the nonlinear dynamic analysis. Thus, the convergence check has to be performed after the nonlinear dynamic analysis and therefore a small constraint violation is tolerated. Hence, the implemented first convergence criterion is that the maximum normalized constraint violation $g_{{{\text{max}}}}$ has to be smaller than a specified limit $\delta > 0$:

$$g_{\max } \le \delta .$$

(17)

For both test problems $\delta = 0.01$ has been used, which is sufficiently small such that all converged designs still can be considered feasible. The second criterion is that the relative change of the objective function between two subsequent cycles $k - 1$ and $k$ has to be smaller than a given value $\varepsilon$:

$$\frac{{\left| {f\left( {{\mathbf{x}}^{\left( k \right)} } \right) - f\left( {{\mathbf{x}}^{{\left( {k - 1} \right)}} } \right)} \right|}}{{\left| {f\left( {{\mathbf{x}}^{\left( k \right)} } \right)} \right|}} \le \varepsilon.$$

(18)

Again for both test problems the value $\varepsilon = 0.01$ has been used. If this criterion is satisfied, the objective hardly changes and continuing the optimization is not worthwhile in most cases. Note that time-dependent responses (e.g., intrusions) are only checked at the selected ESL-times. Furthermore, the optimization is terminated with no convergence after 40 cycles, because the move limits $D$ become too small to allow a considerable change of the design thereafter.

To attribute the fact that the linear static sub-problem is only an approximation of the actual nonlinear dynamic response optimization problem, the length of the optimization path of each inner loop is limited. For this purpose, the same strategy is used as in the commercial code ESLDYNA (VRAND 2012) based on two parameter sets. The first set controls the move limits of the current cycle $k$:

$$\begin{gathered} \hat{x}_{j}^{{{\text{U}},\left( k \right)}} = {\text{min}}\left( {x_{j}^{{\text{U}}} ,{ }x_{j}^{{\left( {k - 1} \right)}} + D\left| {x_{j}^{{\left( {k - 1} \right)}} } \right|} \right) \hfill \\ \hat{x}_{j}^{{{\text{L}},\left( k \right)}} = {\text{max}}\left( {x_{j}^{{\text{L}}} ,{ }x_{j}^{{\left( {k - 1} \right)}} - D\left| {x_{j}^{{\left( {k - 1} \right)}} } \right|} \right)\quad D \in \left[ {0,{ }1} \right]. \hfill \\ \end{gathered}$$

(19)

Parameter $D$ controls the size of the current move limits in cycle $k$ by means of two parameters $D_{{{\text{ini}}}}$ and $\beta$. The parameter $D_{{{\text{ini}}}}$ corresponds to the initial value of $D$ in (20). The reduction factor $\beta_{{{\text{red}}}} \in { }\left( {0,{ }\left. 1 \right]} \right.$ is used to control how $D$ changes from cycle to cycle (outer loop)

$$\begin{array}{*{20}c} {D^{\left( 1 \right)} = D_{{{\text{ini}}}} { };} & {k = 1} \\ {D^{\left( k \right)} = D^{\left( 1 \right)} \beta^{k - 1} { };} & {k > 1} \\ \end{array}$$

(20)

In this publication, the parameter set $D_{{{\text{ini}}}} { } = { }0.2$ and $\beta = 0.9$ has been used in all test problems. Since the move limits only limit the length of the optimization path of each iteration a second restriction is needed. The second parameter set controls the number of iterations (inner loop) by the parameter $max_{{{\text{iter}}}}$ defining the maximum number of iterations per cycle (outer loop). For all examples shown below, the parameter $max_{{{\text{iter}}}} = 2$ has been used. This means that each cycle contains 2 iterations.

Due to the fact that DiESL does not start from a single undeformed initial model but from deformed structures related to the different times, it may happen that the linear static response optimization terminates with an error during the calculation of the ESLs due to excessively deformed elements and thus poor element quality. In order to realize a robust application of DiESL, an automated repair mechanism for the mesh of the affected ${\text{LSM}}$ s has been developed by deleting the failed elements. The deleted elements of the previous cycle are added back at the beginning of each cycle and it is checked again if there are failed elements.

The overall flow scheme explained before is implemented as follows:

Step 1:	Set initial design variables and parameters ($k = 0,{\mathbf{x}}^{{\left( {k = 0} \right)}}$, $\delta$,$\epsilon$, $D_{{{\text{ini}}}}$, $\beta ,$ $max_{{{\text{iter}}}}$)
Step 2:	Perform nonlinear dynamic analysis with ${\mathbf{x}}^{\left( k \right)}$
Step 3:	If adaptive time selection is used, determine ESL-times ${\mathbf{t}}^{*}$ by fitting appropriate response, e.g., contact force
Step 4:	If $k > 0$, check convergence criteria: If Eqs. (17) and (18) are satisfied then terminate the process
Step 5:	Calculate the incremental displacements ${\Delta }{\mathbf{u}}\left( {t^{i} } \right)$ and the node coordinates ${\mathbf{r}}\left( {t^{i} } \right)$ of all LSMs for all selected time steps $t^{i}$
Step 6:	Check the element quality of each LSM FE mesh. If check was not successful, delete failed elements in respective LSM and repeat Step 6 for the remaining mesh
Step 7	If local adaption of Young’s moduli is used, check maximum plastic strain of each element: If ${\upvarepsilon }_{{j,{\text{pl}}}} \left( {t^{i} } \right) > 0$ then adapt Young’s modulus and Poisson’s ratio of element $j$ in ${\text{LSM}}^{i}$ to ${\text{E}}_{{\text{T}}}$
Step 8:	Calculate the incremental equivalent static loads $\Delta {\mathbf{f}}_{{{\text{DiESL}}}}^{i}$
Step 9:	Update the move limits $D$ according to Eq. (20)
Step 10:	Solve the linear static response optimization problem with the incremental equivalent static loads $\Delta {\mathbf{f}}_{{{\text{DiESL}}}}^{i}$
Step 11:	Update the design variables in the nonlinear model, set $k = k + 1$ and go to step 2

3 Test problems

Two test problems for the examination of the influence of the adaptive selection of ESL-times and the local adaption of Young’s moduli on the convergence of the DiESL method are presented in this chapter. The first test problem is a simplified model of a crash side impact problem with 7 design variables. The second test problem is a crash box defined by 12 design variables colliding with a rigid impactor. During the deformation process of the crash box, a lot of buckling and self-contacting can be observed. Therefore, this highly nonlinear example is well suited for the examination of the influence of the local adaption of ESL-times and the local Young’s modulus adaption on the convergence behavior of the DiESL method. Both test problems represent sizing optimization problems in which the mass of the structure is to be minimized, while the intrusion of an impactor is constrained. A separate contact formulation for the contact between impactor and structure is used in LS-DYNA for both test problems. For the side impact problem, a contact between impactor and structure is defined in the OptiStruct model such that the displacement of the impactor can be constrained directly. For the crash box, no contact is defined in the OptiStruct model. Instead, the displacement of a structural node which is known to be in permanent contact with the impactor is constrained instead (Fig. 14). For both test problems, the adaptive time selection is compared to an equidistant distribution of ESL-times. Furthermore, we compare the results of both problems with and without local adaption of the Young’s moduli in the design domain. Finally, the combination of both extensions is tested.

3.1 Side impact

Figure 7 shows a rigid pole (${\text{mass}} = 44.67 {\text{kg}}$) with initial speed $v_{{\text{y}}} = - 8{\text{ ms}}$ colliding with a frame structure. Symmetry conditions with respect to the YZ plane are inserted to reduce the computational effort and to improve numerical stability. The translational degrees of freedom of the pole in $x$- and $z$-direction as well as all its rotations are locked. The structure is clamped along the distant edge in all 6 degrees of freedom using single-point constraints (SPC). The frame structure is made of steel (Young’s modulus $E = 210 {\text{GPa}}$, density $\rho = 7850 {\text{kg/m}}^{3}$, fully integrated shell elements, number of nodes = 6412) and bilinear plastic material behavior is applied (tangent modulus $E_{{\text{T}}}$ = 0.6 GPa, yield stress = 0.25 GPa). In this case, the material properties of the linear auxiliary models can be adapted correctly because the linear models use the same two stiffness values. The simulation end time is $100 {\text{ms}}.$ This captures the pole’s maximum intrusion and part of its rebound. Consequently, for the equidistant ESL-time distribution (ET) the last ESL-time step is selected as 100 ms. In contrast to this, the last ESL-time is always set to the time of maximum intrusion when using adaptive selection of ESL-times (AT). Furthermore, the time corresponding to the maximum crash force is prescribed in each DiESL cycle.

Fig. 7

FE model of the pole impact and labeling of design variables (left). Contact force curve $f\left( t \right)$ and corresponding optimized fit $l\left( t \right)$ for $n_{{\text{T}}} = 20$ and the initial design $x_{j} = 0.8 {\text{mm}}$ (right)

As shown in Fig. 7 the design of the FE model is specified by seven design variables each corresponding to a single sheet metal thickness (except $x_{1}$ representing all 12 facets of both crossbars and $x_{7}$ representing both end caps of rocker profile). The same mesh is used for both linear and nonlinear analysis to avoid results mapping. The objective is to minimize the mass of the frame structure while constraining the maximum intrusion of the pole in negative y-direction $d\left( {\varvec{x}} \right)$. Mathematically, the problem is given as follows:

$$\begin{gathered} \min mass\left( {\mathbf{x}} \right) \hfill \\ {\text{s.t.}}\quad d\left( {\mathbf{x}} \right) \le 200{\text{ mm}} \hfill \\ 0.5{\text{ mm}} \le x_{j} \le 3{\text{ mm }}; j = 1, \ldots ,7. \hfill \\ \end{gathered}$$

First the adaptive selection of ESL-times (AT) is compared to an equidistant distribution of ESL-times (ET). A multi-start optimization study was conducted in order to evaluate the results independently from the initial values: A DOE was created to generate 20 configurations of uniformly distributed design variables (space filling design) used as initial values. The DOE was created as a Strength Two Orthogonal Array. Furthermore, the number of ESL-times $n_{{\text{T}}}$ is varied between 5, 10, and 20. For evaluation, two criteria are used: the number of cycles* necessary for convergence and the resulting objective value ($mass^{*} )$. Both are related to the approximation quality of the linear auxiliary problem: both faster convergence and smaller $mass^{*}$ suggest a better quality of the linear approximation model in the design domain. Figure 8 displays the resulting $mass^{*}$ against the cycles necessary for convergence for each multi-start optimization run. The results are lumped in two clusters that can be distinguished by their mass. They represent two different local optima in the design space. Table 1 shows the averaged results over all configurations for both approaches. Furthermore, the averaged results for the two clusters are given. It can be seen that the averaged $mass^{*}$ is similar for ET and AT independent of the observed scenario (all runs or individual cluster). The average number of cycles* until convergence is slightly better for AT. This advantage increases if the number of ESL-times is reduced. This can be confirmed by comparing individual results of both approaches in Fig. 8.

Fig. 8

Resulting masses and corresponding cycles of all multi-start optimization runs, varying $n_{{\text{T}}}$, and the selection of ESL-times

Table 1 Multi-start results for equidistant distribution of ESL-times (ET) and adaptive time selection (AT)

The superior convergence behavior of the AT approach with respect to cycles* for smaller $n_{{\text{T}}}$ may be attributed to the fact that the last ESL-time is set as the time of maximum intrusion of the impactor. Hence, AT has a better utilization of the small number of ESL-times which concentrate on the relevant deformation process.

Now we want to examine the influence of the local adaption of Young’s moduli (LA) with and without adaptive time selection of ESL-times on the DiESL method’s approximation quality. Again the averaged results are presented in tabular format (Table 2) and the individual results are plotted in Fig. 9. It turned out that a considerable number of runs with local adaption of Young’s moduli (LA) converged to a new optimum with $mass^{*} \le 7.0 {\text{kg}}$. Based on the huge amount of design points evaluated in the course of this and other studies and from examination of contour plots we conclude that this is the global optimum. However, in average, the runs with LA need more cycles to converge than the corresponding runs without adaption (NLA). These convergence issues can be attributed to the oscillation of some design variables over cycles in the region of the global optimum.

Table 2 Multi-start results for equidistant distribution of ESL-times (ET) and adaptive time selection (AT) with (LA) and without (NLA) local adaption of Young’s moduli

Fig. 9

Resulting masses and corresponding cycles of all multi-start results using $n_{{\text{T}}} = 20$ for equidistant ESL-times (left) and with adaptive ESL-times selection (right)

In order to understand why LA converges to a better optimum, we want to compare the constraint’s contour lines for the nonlinear dynamic problem and the corresponding DiESL approximation for the selected design points. To facilitate visualization of the contour lines suitable regions have to be identified in which only 2 design variables are relevant. We start by comparing the optima resulting from all runs summarized in Table 2. Depending on the values of the design variables, ${\mathbf{x}}^{*}$ and the resulting $mass^{*}$ four different groups can be observed, which are illustrated on the left side of Fig. 10 as parallel coordinates. Here, the average value of the resulting designs ${\mathbf{x}}^{*}$ and the corresponding standard deviations are plotted for each optimum. It can be seen that all design variables except for $x_{1} , x_{4}$, and $x_{5}$ are at their lower bound. Furthermore, it can be concluded that the four groups can be classified by their value $x_{4}$.On the right side of Fig. 10$x_{1}^{*}$ and $x_{5}^{*}$ are therefore plotted over $x_{4}^{*}$. In the range $x_{4}^{*} \le 1.1{\text{ mm}}$, all values $x_{5}^{*}$ are at their upper bound. Therefore, we conclude that $x_{5}$ is irrelevant in this range. On the other hand for $x_{4}^{*} > 1.1 {\text{mm}}$, all values $x_{1}^{*}$ are at their lower bound. Thus, we end up with two ranges, each spanned by only two design variables.

Fig. 10

Averaged ${\mathbf{x}}^{*}$ with corresponding standard deviation for all runs using $n_{{\text{T}}} = 20$ split into 4 groups: a $mass^{*} \ge 7.2 {\text{kg}}$, b $7.0 {\text{kg }} \le mass^{*} < 7.2 {\text{kg}}$ and $x_{4} > 1.1 {\text{mm}}$ runs, c $7.0 {\text{kg }} \le mass^{*} < 7.2 {\text{kg}}$ and $x_{4} \le 1.1 {\text{mm}}$, and d $mass^{*} < 7.0 {\text{kg}}$ (left) and thicknesses $x_{1}$(cross), $x_{5}$ (circle) over $x_{4}$ (right)

For each range, a DiESL approximation is created in one representative design point, and the intrusion contour lines of these approximations are compared to those from the nonlinear dynamic problem. The contour lines of both design points are illustrated in Fig. 11 for NLA and LA each combined with ET and AT. In general, all DiESL approximations fit the overall trend of the original problem for both design points. The noisy contour line of the nonlinear dynamic problem is smoothed by the DiESL approximation, which is worthwhile when applying gradient-based optimization. The local adaption of Young’s moduli improves the approximation quality of DiESL in both design points. This improvement can be seen most clearly in Fig. 11 right, which represents the global optimum with $mass^{*} \le 7.0{\text{ kg}}$. The contour lines illustrate why this optimum is found only with LA: The angle between the objective’s and constraint’s contour line changes its sign. While the NLA DiESL approximations suggest that increasing $x_{4}$ will yield a smaller mass at constant intrusion, both LA approximations indicate the opposite. Furthermore, with LA both mass and intrusion contour lines are almost parallel near $x_{4} = 0.6 {\text{mm}}$, which may explain the observed oscillations of design variables in the region of the global optimum: Two nearby design points may lead to angles between mass and intrusion with opposite sign, this will cause the optimizer to invert its search direction.

Fig. 11

Intrusion contour lines of objective and constraints for the nonlinear dynamic problem and the DiESL approximation for range 1: design point $x_{4} = 1.98 {\text{mm}}$; $x_{5} = 2.73 {\text{mm}}$; $d\left( {\mathbf{x}} \right) = 199.1 {\text{mm}}$ (left) and range 2: design point $x_{1} = 0.8 {\text{mm}}$; $x_{4} = 0.6 {\text{mm}}$; $d\left( {\mathbf{x}} \right) = 200.6 {\text{mm}}$ (right). Notes: (1) All remaining design variables are set to their respective bound ($0.5 {\text{mm}}$ or $3.0 {\text{mm}}$), (2) two mass contour lines have been added to each plot

It can therefore be concluded that LA seems to generally improve DiESL approximation quality. The delayed convergence caused by oscillations around the global optimum which were only observed when using LA cannot be attributed to LA as an inherent weakness, instead it can be explained with the nature of the examined example. The fact that only LA managed to find the global optimum is clearly a product of the improved approximation quality.

In order to better understand why the DiESL method reflects the original problem’s global trend, we plot the contour lines of the design point in range 2 (Fig. 11 right) for all selected ESL-times ${\mathbf{t}}$ in Figs. 12 and 13. Since the ESL-times differ, Fig. 12 contains the equidistant ESL-times and Fig. 13 the ones resulting from AT. For sake of clarity, the legend is omitted as it is identical to that in Fig. 11. As explained earlier, the total displacements in DiESL are the cumulated sum of displacement increments of all previous LSMs. Consequently, the sensitivities or the orientation of the contour lines are also sums of all previous LSMs. The orientation of the contour lines or direction of sensitivities can therefore also be understood as the weighted average of the LSMs sensitivities, where the weight factors are the sensitivity magnitudes. Examining Figs. 12 and 13 the previous findings can be reconfirmed: both extensions “adaptive time selection” (AT) as well as “local stiffness adaption” (LA) improve the approximation quality of the standard DiESL approach using equidistant times (ET) and no local stiffness adaption NLA. It can also be seen that poor DiESL approximations at early time steps get propagated to subsequent times and have a negative impact, although this may get mitigated.

Fig. 12

Contour lines of intrusion constraint for the nonlinear dynamic problem and the ET DiESL approximation for range 2: design point $x_{1} = 0.8 {\text{mm}}$; $x_{4} = 0.6 {\text{mm}}$ for all equidistant ESL-times and respective intrusions. Note: the legend is given in Fig. 11

Fig. 13

Contour lines of intrusion constraint for the nonlinear dynamic problem and the AT DiESL approximation for range 2: design point $x_{1} = 0.8 {\text{mm}}$; $x_{4} = 0.6 {\text{mm}}$ of all adaptively selected ESL-times and respective intrusions. Note the legend is given in Fig. 11

3.2 Crash box

Figure 14 shows a crash box being crushed by a rigid impactor (${\text{mass}} = 622 {\text{kg}}$) with initial speed $v_{{\text{y}}} = 4.167{\text{ ms}}$. This model has been published originally by (Ma et al. 2020) but several changes have been made. For example, symmetry conditions with respect to the YZ plane have been inserted in order to reduce the computational effort and to improve numerical stability. The impactor’s translational degrees of freedom in $x$- and $z$-direction as well as all its rotations are locked. The structure is clamped along the distant edge in all 6 degrees of freedom using single-point constraints (SPC). The crash box structure is made of steel (Young’s modulus $E = 210 {\text{GPa}}$, density $\rho = 7850 {\text{kg/m}}^{3}$, fully integrated shell elements, number of nodes = 5842). In contrast to the side impact model, a piecewise linear material behavior is applied here (Fig. 14, right). This means the material behavior cannot be adapted correctly by using a bilinear material model in the design domain. For sake of simplicity the stiffness of elements in the plastic range is globally approximated using the tangent modulus $E_{{\text{T}}}$ = 0.4 GPa, this corresponds to the slope of the dashed yield curve in Fig. 14, right. This is considered to be accurate enough, since as illustrated in Fig. 15 most parts of the structure are in the plastic range even at the very beginning of the deformation process.

Fig. 14

FE model of the crash box and labeling of design variables (left); applied material model and corresponding approximation for local adaption of Young’s moduli (right)

Fig. 15

Plastic strain in nonlinear dynamic analysis of the crash box model for an initial design of the multi-start study at an early time step $t = 8{\text{ ms}}$ (left) and corresponding linear submodel (right)

As shown in Fig. 14 the design of the crash box is specified by twelve design variables, each corresponding to a single sheet metal thickness. As before, the same mesh is used for both linear and nonlinear analysis to avoid results mapping. A contact is defined between the crash box structure and the impactor in the nonlinear model only. The linear model does not include a contact definition, nor does it include the impactor. The objective of the optimization problem again is to minimize the mass of the structure while constraining the maximum intrusion in y-direction $d\left( {\varvec{x}} \right)$. The impactor’s intrusion in the design domain is approximated by using the displacement of a structural node in the impact zone. Mathematically the problem is given as follows:

$$\begin{gathered} \min mass\left( {\mathbf{x}} \right) \hfill \\ {\text{s.t.}}\quad d\left( {\mathbf{x}} \right) \le 160{\text{ mm}} \hfill \\ 0.5{\text{ mm}} \le x_{j} \le 2.5{\text{ mm }}; j = 1, \ldots ,12. \hfill \\ \end{gathered}$$

Figure 16 shows the contact force curve for a very soft design. This soft design with consequently high intrusion has been chosen deliberately in order to determine the largest ESL-time to be considered for the equidistant ESL-time distribution (ET): the last equidistant ESL-time is selected as $t_{{n_{{\text{T}}} }} = 98 {\text{ms}}$, such that the maximum intrusion of the soft design is captured. Since the maximum intrusion exceeds the constraint of 160 mm in this case, feasible designs are expected to reach the maximum intrusion earlier. As before, with the adaptive selection of ESL-times (AT), the last ESL-time is always prescribed as the time at which maximum intrusion occurs for any given design. Additionally, the time corresponding to the maximum crash force is prescribed in each DiESL cycle as illustrated in Fig. 16, right.

Fig. 16

Contact force curve for $x_{j} = 1.2 {\text{mm}}$ and corresponding piecewise linear fit using $n_{{\text{T}}} = 14$ equidistant ESL-times (left) and $n_{{\text{T}}} = 14$ adaptive ESL-times (right)

For the crash box, the evaluation of both extensions AT and LA introduced in this paper model follows the same procedure as explained before for the side impact example. Now, the number of ESL-times varied between 7, 14, and 31. Furthermore, the number of multi-start optimization runs has been increased to 40 uniformly distributed initial designs, since the number of design variables is higher in this case.

First we want to check the influence of the adaptive selection of ESL-times on the approximation quality of the DiESL method. In Table 3, the averaged multi-start results are given for both settings ET and AT. We introduce one subgroup of all runs, namely, the runs converged to the best optimum found which is characterized by 2 features: firstly, all design variables cluster in a small region of the design space and secondly they have a mass smaller than 0.55 kg. Filtering optimal multi-start designs by the second criterion always yields designs located in the first cluster criterion, hence the second is used to define the subgroup. The cluster outlining this best optimum is illustrated in Fig. 17 (right) using $n_{{\text{T}}} = 14$. Here, the average value for each design variable of the resulting design ${\mathbf{x}}^{*}$ and the corresponding standard deviation is plotted for all 40 runs as well as for the best optimum. The standard deviation of all runs is very high which means that many local optima were found. This can be explained by the nature of this crash box example: it is a highly nonlinear example with a lot of local optima, and we did not employ a global optimizer. Therefore, it is rather less surprising that the DiESL method yields several different local minima. Nevertheless, it can be seen in Fig. 17 (left) and in Table 3 that with AT the method converges more often to the best optimum (e.g., 4 ET vs. 19 AT for $n_{{\text{T}}} = 14$). This trend holds if the number of ESL-times is varied (Table 3, Fig. 18 right).

Table 3 Averaged multi-start results for equidistant distribution of ESL-times (ET) and adaptive time selection (AT)

Fig. 17

Resulting masses and corresponding cycles of multi-start study for $n_{{\text{T}}} = 14$ (left) and averaged ${\mathbf{x}}^{*}$ with corresponding standard deviation for the two groups: a all runs and b runs converged to the best optimum for $n_{{\text{T}}} = 14$ (right)

Fig. 18

Average mass against number of ESL-times for all runs (left) and runs converged to the best optimum subgroup (middle); number of runs converged to the best optimum subgroup (right)

Additionally, the average $mass^{*}$ of both, all runs and subgroup runs, is smaller if adaptive time selection is applied when compared to equidistant times, cf. Figure 18 (left, middle). It can also be seen that the benefits of AT become more pronounced for small numbers of ESL-times $n_{{\text{T}}}$: Fig. 18 (left) shows that the averaged $mass^{*}$ of all runs increases and Fig. 18 (right) illustrates that the number of runs converging to the subgroup goes down as $n_{{\text{T}}}$ is reduced. This holds for both methods ET and AT, but ET suffers much more. In contrast, the average number of cycles necessary for convergence for all runs does not reveal a consistent trend, see Table 3.

Summarizing for the crash box example, the adaptive selection of ESL-times leads to several advantages: lower objective value, faster convergence to the best optimum subgroup and even better objective values for runs converged to the subgroup. These observations lead us to the conclusion that adaptive time selection indeed improves the approximation quality of DiESL in this example.

Finally, the local adaption of Young’s moduli and its combination with AT is investigated. The number of ESL-times is set to $n_{{\text{T}}} = 31$ for all runs. In this case we limit ourselves to the tabular comparison in Table 4. Evaluating all runs for ET, the average $mass^{*}$ increases if LA is applied. Also, the number of runs converged to the best optimum reduces. In this case, the local adaption of Young’s moduli seems to reduce the approximation quality of the DiESL method. This may be explained by the following consideration: with the local adaption of Young’s moduli the tangent stiffness at the start of the increment is employed for the whole increment represented by an LSM. With ET these increments may be relatively big. If the increments are chosen regardless on the nonlinear behavior of the structure, the tangent stiffness may be a worse approximation for the whole increment.

Table 4 Averaged multi-start results for local adaption of Young’s moduli in each LSM (LA), no local adaption (NLA) in combination with both equidistant (ET) and adaptive (AT) selection of ESL-times using $n_{{\text{T}}} = 31$

For AT the combination with LA has negligible impacts on the number of cycles and $mass^{*}$ and leads to no considerable advantage. There is one conceivable explanation for the missing benefits obtained with LA: In the crash box example there is a huge amount of elements which enter the plastic domain very early during the simulation. Consequently, nearly all LSMs predominantly consist of elements in the plastic range. In that case, elements in the elastic range have only a small impact on the overall behavior. Therefore, the improvement due to the LA implementation can be neglected. Note that in the extreme case where all elements are in the plastic domain, the LA approach has no impact on the optimization results at all: for a reduced global stiffness in the linear models, the ESL forces scale accordingly, but the linear model yields identical results from the linear optimizer’s perspective. This means the LA approach can only improve the approximation quality of the DiESL models if there is a reasonable mix of plastic and non-plastic elements.

4 Conclusion

In this paper, two extensions of the DiESL method are presented: First the adaptive selection of ESL-times (AT) and secondly a local adaption of Young’s moduli (LA) in the design domain. Both extensions are evaluated on the basis of two representative sizing optimization examples by minimizing mass with a displacement constraint on an impactor.

The evaluation is conducted by means of several DOE-based multi-start optimization studies with varying initial designs. Three criteria are used to evaluate the influence of the extensions on the DiESL method’s performance: the average number of cycles required for convergence, the number of multi-start runs converging to the global optimum, and the resulting average objective value.

The adaptive selection of ESL-times yields different results depending on the example considered. For the side impact example no significant influence on the chosen criteria can be observed, as long as the number of ESL-times is sufficiently high. As the numbers of ESL-times decrease the DiESL methods converge faster if the ESL-times are selected adaptively. Observing the crash box example, the adaptive selection of ESL-times leads to a considerable improvement of all criteria: It is shown that on average the extension converges both faster and more often to the global optimum, and the average objective is better for the overall multi-start runs. These advantages become more dominant if the number of ESL-times is reduced. It is plausible that the difference between both approaches is less dominant for bigger $n_{{\text{T}}}$ because in that case ${\Delta }t$ is small even for equidistant times (ET). It is worth noticing that AT enables the reduction of $n_{{\text{T}}}$ without drastically reducing the approximation quality compared to ET. As the number of $n_{{\text{T}}}$ reduces, less LSMs have to be assembled and the computational effort decreases, which is especially beneficial when handling large models or if computational resources are limited. However, in some cases the usage of ET may also be beneficial, for example if structural responses like velocities or accelerations are involved in the optimization problem, because they require the usage of finite differences. The major differences between the crash box and the side impact example are the number and intensity of structural nonlinearities forming in the course of the simulation. For example the crash box develops more strong contacts and plastic hinges than the side impact. This difference is reflected in the shape of the contact force curve and it may explain the different results of the two examples.

Due to the fact that the optimal number of ESL-times $n_{{\text{T}}}$ is unknown at the beginning of the optimization it may be beneficial if not only the locations of breakpoints are optimized but also $n_{{\text{T}}}$ itself instead of defining it. A crude approach is to compute the SSR values for a predefined range of $n_{{\text{T}}}$ values and to select the smallest $n_{{\text{T}}}$ still fulfilling the SSR accuracy requirements. An alternative may be to use a more sophisticated approximation of the response curve and an optimization algorithm wherein both the number and the position of ESL-times are optimized. Such methods are described for instance in Costa et al. (2018), Bertolino et al. (2021) and they could be the basis for improving the AT algorithm in future.

The local adaption of Young’s moduli also yields different results depending on the observed example: In the case of the side impact, LA yielded a new optimum. The advantages of LA are clearly illustrated by comparing the intrusion constraint’s contour lines of the original nonlinear dynamic problem with the DiESL approximation. However, for the crash box less runs yield the best optimum when combining LA and ET. This disadvantage could be compensated by the combination of LA with AT, but still it does not lead to an improvement. The reason may be related to the ratio of elements in the elastic and plastic domain: Most elements of the crash box plasticize already at very early stages of the deformation process, whereas in the case of the side impact only some parts of the structure do. We therefore conclude that the local adaption of Young’s moduli only leads to improved approximation quality of the DiESL method if neither the plasticized nor the elastic element portion dominates the ratio.

Summarizing, the applicability of both extensions depends on the problem at hand. The adaptive selection of ESL-times (AT) leads to a considerable improvement of the DiESL method’s approximation quality if the contact force clearly indicates the presence of intense nonlinearities. The local adaption of Young’s moduli (LA) only yields an observable improvement if neither the elements in the elastic nor in the plastic domain are dominating the structure’s behavior. In general, the DiESL method leads to sufficiently good approximations even without either extension. In order to better assess the influence of both extensions in practice, further examples should be optimized.