All articles published by MDPI are made immediately available worldwide under an open access license. No special

permission is required to reuse all or part of the article published by MDPI, including figures and tables. For

articles published under an open access Creative Common CC BY license, any part of the article may be reused without

permission provided that the original article is clearly cited. For more information, please refer to

https://www.mdpi.com/openaccess.

Feature papers represent the most advanced research with significant potential for high impact in the field. A Feature

Paper should be a substantial original Article that involves several techniques or approaches, provides an outlook for

future research directions and describes possible research applications.

Feature papers are submitted upon individual invitation or recommendation by the scientific editors and must receive

positive feedback from the reviewers.

Editor’s Choice articles are based on recommendations by the scientific editors of MDPI journals from around the world.

Editors select a small number of articles recently published in the journal that they believe will be particularly

interesting to readers, or important in the respective research area. The aim is to provide a snapshot of some of the

most exciting work published in the various research areas of the journal.

The objective of this study is to develop an optimization methodology to find a layout that traces a prescribed force–displacement curve through a topology optimization approach. To this end, we propose an objective function to minimize the difference between a prescribed force–displacement curve and the curve calculated at each iteration of the optimization process. Slope constraints are introduced to solve issues encountered when using a small number of target points. In addition, a projection filter is employed to suppress the gray region observed between the solid and void regions, which generally occurs when using a density-based filter. A recently proposed energy interpolation scheme is implemented to stabilize the instability in the nonlinear analysis, which generally results from excessive distortion in the void region when the structure is modeled on a fixed mesh in the topology optimization process. To validate the outlined methodology, several case studies with different types of nonlinearity and structural features of the obtained layouts are investigated.

The generalized Hook’s elasticity law is not valid if a structure undergoes large deformation, for which the force–displacement (F-D) curve exhibits nonlinear behavior due to geometric nonlinearity [1,2,3]. Geometric nonlinearities can be the softening type, where variations of displacement increase for constant increments of the external force, or hardening types, where the variations decrease. For instance, a thin plate is generally governed by hardening nonlinearity [3]. Bensøe and Kikuchi [4] developed the topology optimization method, which has been applied to a wide range of areas over the last two decades, from mechanics to fluid dynamics, acoustics, and electrical dynamics. Most of these studies have been conducted using linear finite element analysis. However, Buhl et al. [5] dealt with an optimization problem concerning geometric nonlinearity due to large deformation, and many studies have since been introduced for topology optimization problems of mechanical structures associated with geometrically nonlinear modeling [6,7,8].

Seikomo and Noguchi [9] conducted one of the first studies of an optimal design problem that traces a prescribed force–displacement curve. They introduced one-step and multi-step approaches separately and they compared the results. Bruns et al. [10] and Bruns and Sigmund [11] conducted studies to trace a prescribed curve involving snap-through behavior. They used prescribed values on the curve as a constraint in the optimal design problem. In addition, to take into account snap-through and snap-back behaviors on the force–displacement curve, they built their method on a modified arc-length procedure instead of the Newton–Raphson procedure, which is normally used in nonlinear analysis.

Yoon et al. [12] developed a method to find an optimal layout tracing a prescribed curve defined on the reaction force and displacement plane through topology optimization. Recently, Mathias et al. [13] performed a study on stiffness optimization considering the secant matrix and tangent matrix, which can be defined from the F-D curve. The use of a small number of target points can be expected to decrease the computational cost. However, if the nonlinearity in the prescribed F-D curve differs from the inherent nonlinearity in the initial design, for example, the hardening nonlinearity for a thin plate, it can be impossible to find a layout that thoroughly follows the prescribed curve for the softening nonlinearity in the optimization problem with a small number of target points. For such cases, the slope constraint not only improves the convergence in the optimization design process but also decreases the computation cost in the Newton–Raphson stage. The details are discussed in the case study section.

Moreover, in topology optimization studies dealing with large deformation, low stiffness elements can undergo extensive distortion during the optimal design process, which leads to numerical instabilities in the nonlinear analysis when using the Newton–Raphson method [14,15,16]. To circumvent the numerical instability induced by low stiffness elements, Buhl et al. [5] and Pedersen et al. [17] excluded internal nodal forces that were surrounded by low stiffness elements during the optimization process. Yoon and Kim [18] parameterized the connectivity by considering every element to be connected by imaginary springs. Recently, Wang et al. [19] proposed an energy interpolation scheme using a threshold projection function.

In the present work, we propose a methodology to compute the layout of a structure that traces a predefined F-D curve in which softening, hardening, or both nonlinearities are involved through density-based topology optimization. We employ the energy interpolation scheme to mitigate the numerical instability. In addition, filtering is used to ensure a global convergent design in the topology optimization, through sensitivity filtering [20], density filtering [15], and the binary bat algorithmic filter [21] approaches. However, the filters inherently produce gray transition regions between solid and void regions. Recently, several projection schemes [22,23,24] have been proposed to alleviate the gray transition problem. In this study, we employ a threshold projection scheme expressed by a smoothed Heaviside function [25]. Finally, in order to deal with small numbers of target points along the F-D curve, we take into account the slope of the curve, as well as the values of the prescribed F-D curve, in order to adapt to the nonlinearities on the curve. Therefore, the slope is involved in the optimization problem as a constraint.

The objective of this study is to find a structural layout that follows a target nonlinearity through density-based topology optimization. In research carried out previously, a relatively large number of target points were used [12]. In this study, however, we propose a method that guarantees convergence in the optimization design process with fewer target points. To do so, the slope defined at the corresponding target points on the force–displacement curve is set as the constraint in the optimization problem. This slope is related to the secant matrix. Assuming that the external force applied to the structure is concentrated in one place and considering the force–displacement curve for displacement at any degree of freedom, i.e., for single-input multiple-output (SIMO), the sensitivity of the slope to the design variables can be derived by direct differentiation of the inverse of the secant matrix with respect to the design variables.

This study is organized as follows: In Section 2 the underlying theories are explained; in Section 2.1, we define the optimization problem and present the governing equations for a geometrically nonlinear finite element and sensitivity analysis for the objective function; in Section 2.2, we describe the density filter and parameterization for the energy interpolation; in Section 2.3, we describe concepts of the slope used in the optimization problem as a constraint and the sensitivity analysis; in Section 3, we demonstrate the performance of the proposed method using a thin plate; and in Section 4, we present the conclusions. This study dramatically improves the computational speed and the convergence of the optimization analysis by adding slope constraints as compared with similar researches conducted previously. In addition, the case studies in Section 3 were performed through in-house code using MATLAB. The topology optimization and optimal design in this study were basically constructed based on the 88 lines code [26].

2. Theory

2.1. Optimization Formulae

The objective function in the optimization problem is constructed using the relationship between the prescribed values on the F-D plane and the estimated values from the layouts in the optimization process. The displacements at the defined forces have to be calculated for each iteration of the optimization process. Structures show nonlinear behavior when they undergo large deformation and usually have inherent monotonous softening or hardening nonlinearities with respect to external forces. For instance, a thin rectangular plate exhibits monotonous hardening nonlinearity. In this work, we focus on monotonic softening and hardening nonlinearities, which allows us to use the Newton–Raphson method to calculate the displacement field. The solution does not converge when using a small number of target points on the F-D plane and a rectangular thin plate as the initial design model. This occurs because of the difference between the inherent nonlinearity in the initial design model and the desired nonlinearity. Further descriptions can be found in Section 2.3.

To help with the nonlinearity adaptation during the optimization process, the definable slopes at the target points are used as constraints in the optimization problem. Since we are dealing with a density-based topology optimization problem, an individual design variable can vary continuously in the interval [0, 1], where 0 and 1 imply the absence and presence of material, respectively. On the basis of these assumptions, the optimization problem can be defined as follows:

is the objective function, and is a vector of design variables, is the number of defined target points on the F-D curve, is the amplitude of the applied ith external force, is the displacement vector in global coordinate, is the prescribed value of the displacement with respect to the ith external force, and is the weighting function. is the position vector, in which a unit value is assigned at the corresponding degree of freedom where displacements are identified, whereas the other components are set to zero. The superscript T indicates the transpose of a vector. is the slope that can be defined on the F-D plane for the ith external force. Details about the slope and its derivative are given in Section 2.3.

and are the upper and lower bounds of the slope constraint, is a vector of physical density, and is the volume when the entire design domain is occupied by solid elements. is the upper bound of the volume fraction constraint, and is the residual force vector of the structure equilibrium. The residual force vector is iteratively estimated in the nonlinear analysis.

Large deformation with a small strain state is assumed in this study. This means that there is a linear relationship between the Green–Lagrangian strain tensor and the second Piola–Kirschhoff stress tensor :

The fourth-order tensor is the elastic tensor, and is the Green–Lagrangian strain tensor, which is represented as . is the identify matrix. is the deformation gradient tensor defined as . denotes the material coordinate in the undeformed geometry. The total Lagrangian finite element formulation is used.

The residual force vector in the equilibrium state is expressed as follows:

is the external force vector, with the loading vector in which a unit value is assigned at the corresponding degree of freedom where the loading is applied, whereas the other components are set to zero, and is the forcing amplitude. is the internal nodal force vector, which is expressed by the summation of the elements of the internal nodal force vector (i.e., ). The element nodal force vector is expressed as:

where is the stored elastic energy density in the element, and is an element nodal displacement vector.

The integration in Equation (4) is performed over the undeformed volume. Equation (3) is solved iteratively using the Newton–Raphson method. In this process, the tangent stiffness matrix has to be computed:

The details on how to determine the matrix are described in a previous study [27]. In this study, we use the adjoint variable method to calculate the sensitivity of the objective function:

where is the adjoint variable vector:

In the sensitivity-based topology optimization process, the design domain is iteratively updated by the method of moving asymptotes (MMA) [28]. In the present work we follow the standard optimization procedure with MMA [17,25,29].

2.2. Filtering and Parameterizations

Filtering enables the removal of checkerboard patterns [30,31] in optimal design results and the computation of layouts independent from mesh refinements [31]. Mesh-independent density filtering is also used [15,32]. The main concept of this filtering is defining the physical element density to be a weighted average of the neighboring design variables:

where is the filtered density of element i, is the design variable of element j, is the location of element j, is the volume of element j, and is the number of the set of the neighborhood element i within a certain filter radius r. is a weighting function specified by .

In order to mitigate the gray region between the solid and void regions, a threshold projection filter is applied using a smoothed Heaviside function [25]:

where is the threshold in the projection filter, and is a variable to control the sharpness of the projection filter; as →∞, Equation (9) becomes a discrete Heaviside projection. is a physical design density filtered by the threshold projection filter, which is used to represent the degree of occupancy of the material in element e.

The solid isotropic material with penalization (SIMP) method [33] is used to interpolate the Young’s modulus of each element, yielding the following expression for the interpolation:

where is the Young’s modulus of the solid phase (the black regions), and is the Young’s modulus of the void phase (the white regions), which generally have a relation of . The superscript p is a penalization factor (generally, p = 3). Wang et al. [19] proposed an energy interpolation scheme to mitigate the numerical instability due to the excessive distortion in low stiffness elements, which is inevitable in the topology optimization process. The energy interpolation form of element e is expressed as:

is the stored elastic energy density function, and is the stored elastic energy density under small deformation. is the interpolation factor for element e, which is 1 for solid elements and is 0 for void elements. Therefore, solid elements are treated as nonlinear elements in the numerical analysis (i.e., ), whereas void elements are treated as linear elements (i.e., ).

The interpolation factor is expressed similarly to the threshold projection filter in Equation (9) using a smoothed Heaviside function:

where is a threshold between 0 and 1. If is large enough and the physical design density is close to 1, then the interpolation factor will be close to 1. Therefore, the nonlinear analysis is imposed on the corresponding element e. However, if is close to 0, then the interpolation factor will be close to 0. Therefore, the linear analysis is imposed on the corresponding element e.

Using the density filter, threshold projection filter, and the energy interpolation scheme, the sensitivity form of the objective function in Equation (6) is updated as follows:

2.3. Slope Constraint

In this methodology, the slope of the F-D curve at the target points can be defined as constraints, irrespectively of the inherent nonlinearity in the initial design of the structure.

As an illustration, Figure 1 shows three target points defined on the F-D plane and the different slopes that can be associated with them. The target points are denoted by P1, P2, and P3. First, the definable slope at the points is the tangential slope, which are denoted by , , and . Another slope can be expressed by defining the straight line that connects the origin and the target points on the F-D plane and considering its angle with respect to the displacement axis. The associated slopes are denoted by , , and . In general, the tangent slopes introduced in the numerical analysis of nonlinear structures cannot be estimated until the values of the tangent slope at the points are calculated by numerical analysis, even if an optimal structure satisfying the target points is obtained. Therefore, it is not easy to include the tangent slope in the optimization design problem. However, the slopes represented by can be estimated at the same time as when the target points are defined.

If the three points are defined as the target points, as shown in Figure 1, and the optimization process is performed without the slope constraint, the topology optimization process maintains the inherent nonlinearity of the initial structure. For instance, a general rectangular thin plate maintains its hardening nonlinearity, and the resultant curve passes through the target points while retaining the hardening nonlinearity to minimize the objective function. The resultant curve is depicted by a dotted line (····). If one uses a large number of target points, this issue can be resolved without any constraints. However, when using only a small number of target points without the slope constraint, the optimization design process always falls into the local minimum, and the F-D curve behaves as shown by the dotted line in Figure 1.

Another slope can be defined at the target points by tracing the line that connects the origin and the target points and considering its angle with the force axis. These are denoted by , , and . The slope is associated with the secant stiffness matrix, and the slope is associated with the inverse of the secant stiffness matrix. In a previous study [13], the tangential stiffness and secant stiffness matrices were used as objective functions, which were associated with slopes concerning and , respectively, to minimize compliance for nonlinear structures. Even if they eventually maximize the slopes in the F-D plane, the sensitivities to the compliances according to the two stiffness matrices described in detail in the reference are not directly related to slopes and .

Let us now discuss the sensitivity analysis related to the slope constraints. In a similar way to a previous study [34], the principle of virtual work is applied to derive the finite element equation in a matrix-vector form for the in-plane motion of the thin rectangular plate:

where is the secant stiffness matrix, which is composed of the linear stiffness matrix and the nonlinear stiffness matrix . If a structure considered in the optimization problem is displaced by a concentrated force, (i.e., a single input case, in Equation (14), and ), then the slope at the target point P, which is defined for the external force at the ith degree of freedom and the displacement at jth degree of freedom, is represented by the inverse of the secant matrix.

The slope is associated with the inverse of the matrix , and its sensitivity can be derived by the direct differentiation method (DDM):

Equation (16) represents the sensitivity of the inverse of the stiffness matrix with respect to the physical design density vector . The superscript -1 indicates the inverse of a matrix. Assuming that a single external force is applied and that a single response is measured, the element at the corresponding degree of freedom for the external force and the response is a source of sensitivity of the slope to the physical design variable. The sensitivity of the slope is represented by only the single element of the inverse secant stiffness matrix to the corresponding degree of freedom. This is the main reason to employ the inverse secant stiffness matrix instead of the secant matrix. Sensitivity to the design density vector can be induced from the chain rule, as shown in Equation (13).

3. Case Studies

3.1. Clamped-Clamped Case

For validation purposes, let us consider the clamped-clamped rectangular structure presented in Figure 2. The tested structure is 0.6 m long (L = 0.6 m), 0.05 m wide (B = 0.05 m), and 1e-3 m thick (t = 1e-3). Young’s modulus and Poisson’s ratio are 20 GPa and 0.33, respectively. An external force f is exerted on two-thirds of the total length L from the left end in the upward direction. The structure is discretized using a 20 by 240 mesh of four-node elements.

Four amplitudes, namely 0 N, 50 N, 100 N, and 150 N, were defined as the external forces of interest on the F-D curve. Target points were defined by setting displacements at these forces on the external force location. For the parameters in the optimization problem, the volume is set as 70% as compared with the design where all the elements are solid. In the initial design model, the design variable was set as 0.7 over the entire design domain. The slopes are set to the minimum and maximum values with a marginal value of ±10%. The weighting function is set to the corresponding external force levels (i.e., 50 N, 100 N, and 150 N). The value for the projection parameter in Equation (9) is doubled at every 50 iterations of topology optimization or at convergence, with a maximal value of and an initial value of . The threshold is fixed as . The threshold parameters of the energy interpolation scheme in Equation (12) are fixed as = 500 and = 0.01. The filter radius r is fixed as rfilter = 2 × the element size.

Let us first consider the case of a targeted softening nonlinearity. First, we measured displacements at the external force levels on the initial design layout. On the basis of the measured F-D curve from the initial design layout, the prescribed target points for the softening nonlinearity on the F-D plane are arbitrary determined around the F-D curve for the initial design layout, as shown in Figure 3. In order to study the influence of the slope constraint proposed in this work, we apply the optimization processes with and without the slope constraints.

The curve measured from the initial design layout is depicted by a dash and circle line (–O–) in Figure 3. The line seems straight, but it is slightly curved upward, meaning that the initial design behaves with hardening nonlinearity. As shown in Figure 3, on the one hand, the measured curve (–☆–) from the optimization process, in which the slope constraint was incorporated, and the prescribed target curve (–∇–) are in good agreement with each other. On the other hand, the measured curve (–□–) from the optimization process, in which the slope constraint was not incorporated, obviously differs from that of the target curve. In addition, it seems straight as for the initial layout case but is also slightly curved upward, meaning that the inherent hardening nonlinearity that was in the initial layout remained during the optimization process even though the targeted nonlinearity was softening. It can thus be seen that the optimization process has tried to reduce only the objective value while maintaining the inherent nonlinearity in the initial layout. Therefore, it can be concluded that the slope constraint is necessary if different nonlinearity behaviors from the inherent nonlinearity in the initial layout is desired in the prescribed target curve.

In addition, we found it interesting to analyze the process of obtaining results based on with and without slope constraints in the optimization process. Figure 4 shows the history plot for the normalized objective function with respect to the objective function value on the initial design (c0) and volume fraction of each iteration in the optimization process; Figure 4a is the case where the slope constraint is incorporated in the optimization process and Figure 4b is the case where the slope constraint is not incorporated. The total optimization process is repeated until 185 iterations and 268 iterations, respectively.

On the one hand, it is observed that the slope constraint helped to approach the optimal result layout faster in the optimization process. The objective function curve is mostly smooth except for a specific iteration range (120 ~ 170 iteration). On the other hand, when there is no slope constraint, it can be seen that the objective function curve fluctuates in most of the iteration range. In general, the Newton–Raphson method is based on the previous iteration results, therefore, for a smooth curve case, the desired convergence result in the Newton–Raphson process is close to the values of the previous iteration. If the fluctuation is severe in the objective and volume fraction curves, however, a large number of analyses should be performed internally in the Newton–Raphson process to find the convergent point. According to Figure 4, it was observed that, for the case without slope constraint, the Newton–Raphson process takes 1.5~2 times longer for a single force level as compared with the case with slope constraint. In conclusion, although an additional constraint condition was added to the optimization problem, the total analysis time required for the optimization design process did not differ between the two cases, and in some cases, it took a shorter time in the case of constraint. The same tendency was also observed in the hardening nonlinearity case.

As illustrated in Figure 5, the outlined optimization procedure could also be validated for a hardening nonlinearity. As for the softening case, target points are arbitrary decided based on the measured F-D curve from the initial design layout. In order to verify the nonlinear trend of the obtained designs for forcing values other than the targeted one, the F-D curves in the range of 0 to 150 N were measured at intervals of 2 N for both hardening and softening layouts.

Both curves increase monotonically, which means that they do not exhibit a snap-through or snap-back phenomenon in the region of the external force level of interest. In the softening case Figure 5a, the F-D curve shows a softening trend and the optimal curve passes through the target points, which confirms the effectiveness of the proposed method. In hardening cases, as shown in Figure 5b, it does not pass through the target points exactly. This means that the objective function is located at a local minimum, and the proposed method does not always guarantee to give perfectly optimized results for all types of defined target points. A similar case was introduced in a previous study [12].

Figure 6 shows the optimal layouts without deformation at 0N and with deformations at the target points in Figure 5. The upward arrows ↑ indicate the forcing location. Gray regions are hardly seen between the solid and void regions, which shows that the projection filter worked well. In the layout for the softening case in Figure 6a, a specific region leading to the overall nonlinear features of the structure is highlighted by a dashed box, and a schematic representation is provided. The actual layout in the box can be divided into two parts that are denoted by A and B. In the simplified layout, part A is represented by an arched beam, and part B is the foundation to support the part A through clamped boundary conditions. This type of structural configuration like an arched beam has a geometrically similar configuration to the softening spring [35,36]. As stated in the reference, a softening spring is induced by combination of sources of positive and negative stiffness provided by a bistable mechanism. On the one hand, a bistable mechanism is characterized by a snap-through phenomenon. Therefore, it can be interpreted that the desired softening nonlinearity in the optimization process is induced by the snap-through phenomenon occurring in the dotted region. Here, the region in the dash-dot box acts as a prismatic joint in the reference. In addition, it is well known that an arched beam is related to the softening nonlinearity by the buckling phenomenon [37,38]. On the other hand, in the hardening nonlinear case, it is not easy to find distinctive features in the configuration of the resultant layout, and it is considered that the entire structure is associated with the hardening nonlinearity. As shown in Figure 6, the structures have small displacements but a large strain. In this case, the void regions do not have a large deformation, and eventually the use of the energy interpolation method has little effect on the optimization process and the results. The case of large deformation is covered in the next section. Next, in order to investigate the effectiveness of the slope constraint proposed in this study, we proceeded with the optimal design process while increasing the number of target points without the slope constraint. An investigation of the effect of the slope constraint was done in the softening case. First, the resultant layout on the softening case, Figure 6a, was set as the reference layout. Then, the number of target points NTP was increased from three to 15 in steps of three; therefore, five additional optimization design processes were performed. The magnitude of the force considered in the performed individual optimization process is as follows.

For NTP = 3, for instance, the forces considered in the optimization problem were 0 N, 50 N, 100 N, and 150 N. These forces are equal to the force magnitudes in the optimization processes performed in the preceding case studies for softening and hardening. For NTP = 6, the forces considered in the optimization problem were 0, 25, 50, 75, 100, 125, and 150 N. Thereafter, the displacements were measured by applying defined forces to the forcing position of the resultant layout in Figure 6a. The measured displacements were used as in Equation (1), and the optimization design processes were performed without the slope constraints. The optimization design parameters were the same as those used previously.

Figure 7 shows the relationship between the force–displacement curves obtained in the optimization design processes while increasing the number of target points without the slope constraints, and the force–displacement curve (—) measured from the optimal layout (i.e., the same as in Figure 5a (—)). We can verify the performance of the proposed slope constraint by confirming how the optimal F-D curves approach the reference curve (—) when the number of target points is increased. As seen in the figure, the trend of the F-D curve is close to the reference curve as the number of target points increases. This trend is also evident at around 150 N external force. Furthermore, the softening tendency of the force–displacement curve in the cases where the number of target points increases is shown to be strong. In the case of 15 target points, the resultant layout with the desired nonlinear force–displacement characteristic defined as the reference curve was obtained. The CPU of the computer used in this study was Intel (R) Core i9-9900k and the memory (RAM) was 64 GB. Figure 5a shows the procedure for the optimization design process performed considering three target points and slope constraints. Figure 6a shows the resultant layout. The analysis time required was 5 h and 45 min. The time duration for optimization design process to obtain the resultant F-D curve shown in Figure 7 was 15 days and 22 h, which was obtained by considering only the objective function using 15 target points without slope constraints. In conclusion, the proposed slope constraint not only shortens the optimization process time by using only a small number of target points in the optimization problem but also improves the convergence in the optimization design process.

In summary, in the absence of a slope constraint, a large number of target points are required to obtain a layout that follows the prescribed target points. Of course, however, the computational analysis takes a long time in this case. In order to overcome this, this study proposed a slope constraint, which makes it possible to obtain the layout of a structure with a nonlinear F-D curve of the desired with just three target points.

Lastly, we tried to verify if the desired nonlinear behavior would indeed appear in the fabricated structure based on the resultant layout. To do this, every element whose design variable is less than 0.5 is set to 0 and whose design variable is larger than 0.5 is set to 1. The reference value 0.5 is decided based on the threshold value which was used in the projection filter. Next, the void elements with zero design variables are removed and the structure layout is composed solely of the solid elements with the design variable of 1. The F-D curve for modified layout is finally computed using the commercial software COMSOL, at the same external force levels as those used in the optimization problem.

Figure 8 shows that the F-D curves, which are measured from the modified layout in which only the solid elements are composed, exhibit the same nonlinearities as the curves obtained in the optimization process, and as the curves obtained using COMSOL. Interestingly, one would expect the layout (—□—) in which there are void elements to be stiffer than the layout (—○—) in which the void elements are removed, but the opposite observation is made. This is because the void elements also affect the stiffness of the entire structure. However, the results show that the effect of stiffness reinforced in the solid region is greater than the stiffness reduced by eliminating the void elements. This is natural because in the optimization design, the void region represents a low sensitivity region and the solid region represents a high sensitivity region. Therefore, it is obvious that the stiffening effect reinforced in the high sensitivity region is larger than the stiffness effect removed in the low sensitivity region.

3.2. Clamped-Free Case

In this section, we consider the clamped-free rectangular structure presented in Figure 9. The tested structure is the same as the structure tested in the previous section, and external force f is exerted on the upper side of the free end in the upward direction. The structure is also discretized using a 20 by 240 mesh of four-node elements.

Similar to the previous section, the design variable was set as 0.7 over the entire design domain. Let us first examine the F-D curve in the initial design model. The curve measured from the initial design layout is depicted by a dotted line (∙∙∙∙) in Figure 10. The line seems straight, but it is slightly curved upward, meaning that the initial design behaves with hardening nonlinearity. Three cases were tested with different types of nonlinearity, which are prescribed by defining target points at 0, 5, 100, and 150 N. We first consider softening nonlinearity. As for the previous case, the displacements in each external force level, i.e., the target points, are arbitrary decided based on the measured F-D curve from the initial design layout. The circles (○) in the figure represent the target points in the optimization design process. All the parameters in the optimization problem are set to the same values used in the previous section, except the filter radius r which is fixed as rfilter = 2 × the element size for softening case and rfilter = 6 × the element size for the hardening case.

Figure 10 also shows the F-D curves obtained with the optimal layouts, for forcing amplitudes other than the targeted ones. All optimal curves increase monotonically, which means that they do not exhibit a snap-through and snap-back phenomenon in the region of the external force level of interest.

In both cases, as shown in Figure 10, the optimal curves pass through the target points, which confirms the effectiveness of the proposed method.

Figure 11 shows the optimal layouts with deformations at the target points in Figure 10. The upward arrows ↑ indicate the forcing location. Gray regions are hardly seen between the solid and void regions, which shows that the projection filter worked well. In the results of the softening cases in Figure 11a, the outer frame is composed of a relatively thick layer, but the inside is composed of thin layers. Magnified parts of these thin layers in the softening results are shown in the dashed boxes. When a large external force is imposed (i.e., 150 N), buckling occurs. Hence, it can be concluded that softening nonlinearity is induced in the thin layers by buckling, which show macroscopic bending behavior when a large force is applied.

In contrast, there is only a relatively thick frame without internal or external distinction in the case of hardening. In our experience, the optimization process does not converge when using a large filter radius for the softening cases and when using a small filter radius for the hardening case. This is probably due to the different features of the required structural configuration depending on the nonlinear properties of the structure, and the filter size is considered to be sensitive to reflect the properties.

The difference between the case studies in the previous section and the cases covered in this section is that this section deals with large-strain and large-deformation cases. As stated in the previous section, for a small deformation but large strain case in the clamped-clamped boundary condition, almost the same results were obtained regardless of whether the energy interpolation was used or not. However, when we do not use energy interpolation for the large-deformation case discussed in this section, the Newton–Raphson process performed to measure the displacement during the optimization process often fails to converge by exceeding the maximum number of iteration set. Finally, the F-D curve result at the final iteration of the optimization process does not converge to the target points.

In conclusion, energy interpolation is a necessary tool to mitigate the effect caused by large deformation in the void region. Moreover, as shown in Figure 11, the effect of the energy interpolation technique can be confirmed from the fact that every void region is smooth and without excessive distortion. This feature can be identified at 150 N in Figure 11b, especially where the largest displacement region has been zoomed in with a dotted box.

4. Conclusions

A topology optimization method has been proposed to find the layout of a structure that traces the prescribed nonlinearity on the F-D curve. The desired nonlinearity was realized using a small number of target points on the F-D plane as compared with previous studies. To this end, the slope at the target points was included as a constraint. This serves as a mediator that helps to break the inherent nonlinearity in the initial design model and to adapt to the nonlinearity described by the target points during the optimization process. The slope constraints were also shown to improve the convergence and decrease the computation cost in the optimization process. The necessity of the slope constraint was verified by comparing the optimization design results obtained from increasing the number of target points excluding the slope constraint with those from a smaller number of target points including the slope constraint.

An energy interpolation scheme was used to suppress the instability in the nonlinear analysis caused by the void regions, which inevitably appears in the topology optimization process. In addition, a projection filter with a density filter was used to mitigate the gray region between the solid and void regions. We performed test studies with several types of nonlinearity to verify the method. Moreover, we examined the configuration features of a plane structure that exhibits softening and hardening nonlinearities through the results obtained from the case studies.

Author Contributions

Conceptualization, J.L. and G.K.; methodology, J.L.; writing—original draft preparation, J.L.; writing—review and editing, T.D. and G.X.; supervision, G.K.; All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This work was supported by the Dong-A University research fund.

Conflicts of Interest

The authors declare no conflict of interest.

References

Bathe, K.J.; Bolourchi, S. A geometric and material nonlinear plate and shell element. Comput. Struct.1980, 11, 23–48. [Google Scholar] [CrossRef]

Reddy, J.N. An Introduction to Nonlinear Finite Element Analysis; Oxford University Press: New York, NY, USA, 2004. [Google Scholar]

Reddy, J.N. Theory and Analysis of Elastic Plates and Shells, 2nd ed.; CRC Press: New York, NY, USA, 2007. [Google Scholar]

Bendsøe, M.P.; Kikuchi, N. Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Eng.1988, 71, 197–224. [Google Scholar] [CrossRef]

Buhl, T.; Pedersen, C.; Sigmund, O. Stiffness design of geometrically nonlinear structures using topology optimization. Struct. Multidiscip. Optim.2000, 19, 93–104. [Google Scholar] [CrossRef]

Wang, F.; Sigmund, O.; Jensen, J. Design of materials with prescribed nonlinear properties. J. Mech. Phys. Solids2014, 67, 156–174. [Google Scholar] [CrossRef]

Jansen, M.; Lombaert, G.; Schevenels, M. Robust topology optimization of structures with imperfect geometry based on geometric nonlinear analysis. Comput. Methods Appl. Mech. Eng.2015, 285, 452–467. [Google Scholar] [CrossRef]

Luo, Y.; Wang, Y.M.; Kang, Z. Topology optimization of geometrically nonlinear structures based on an additive hyperelasticity technique. Comput. Methods Appl. Mech. Eng.2015, 286, 422–441. [Google Scholar] [CrossRef]

Sekimoto, T.; Noguchi, H. Homologous topology optimization in large displacement and buckling problems. JSME Int. J. Ser. A Solid Mech. Mater. Eng.2001, 44, 616–622. [Google Scholar] [CrossRef][Green Version]

Bruns, T.E.; Sigmund, O.; Tortorelli, D.A. Numerical methods for the topology optimization of structures that exhibit snap-through. Int. J. Numer. Methods Eng.2001, 55, 1215–1237. [Google Scholar] [CrossRef]

Bruns, T.E.; Sigmund, O. Toward the topology design of mechanisms that exhibit snap-through behavior. Comput. Methods Appl. Mech. Eng.2004, 193, 3973–4000. [Google Scholar] [CrossRef]

Yoon, G.; Noh, J.; Kim, Y. Topology optimization of geometrically nonlinear structures tracing given load-displacement curves. J. Mech. Mater. Struct.2011, 6, 605–625. [Google Scholar] [CrossRef]

Bruns, T.E.; Tortorelli, D.A. An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms. Int. J. Numer. Methods Eng.2003, 57, 1413–1430. [Google Scholar] [CrossRef]

Pedersen, C.B.W.; Buhl, T.; Sigmund, O. Topology synthesis of large-displacement compliant mechanisms. Int. J. Numer. Methods Eng.2001, 50, 2683–2705. [Google Scholar] [CrossRef]

Yoon, G.; Kim, Y. Element connectivity parameterization for topology optimization of geometrically nonlinear structures. Int. J. Solids Struct.2005, 42, 1983–2009. [Google Scholar] [CrossRef]

Díaz, A.; Sigmund, O. Checkerboard patterns in layout optimization. Struct. Optim.1995, 10, 40–45. [Google Scholar] [CrossRef]

Sigmund, O.; Petersson, J. Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct. Optim.1998, 16, 68–75. [Google Scholar] [CrossRef]

Bourdin, B. Filters in topology optimization. Int. J. Numer. Methods Eng.2001, 50, 2143–2158. [Google Scholar] [CrossRef]

Bendsøe, M.P.; Sigmund, O. Material interpolation schemes in topology optimization. Arch. Appl. Mech.1999, 69, 635–654. [Google Scholar] [CrossRef]

Ribeiro, P.; Petyt, M. Multi modal geometrical nonlinear free vibration of fully clamped composite laminated plates. J. Sound Vib.1999, 225, 127–152. [Google Scholar] [CrossRef]

Rivin, E.I. Stiffness and Damping in Mechanical Design; Marcel Dekker, Inc.: New York, NY, USA, 1999. [Google Scholar]

Boehler, Q.; Vedrines, M.; Abdelaziz, S.; Poignet, P.; Renaud, P. Parallel singularities for the design of softening springs using compliant mechanisms. In Proceedings of the ASME 2015 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Boston, MA, USA, 2–5 August 2015. [Google Scholar] [CrossRef]

Timoshenko, S. Theory of Elastic Stability; International Student Edition; McGraw-Hill: New York, NY, USA, 1936. [Google Scholar]

Figure 1.

Introduction of definable slopes in the force–displacement (F-D) curve and estimated F-D curves with and without a slope constraint in the optimization process.

Figure 1.

Introduction of definable slopes in the force–displacement (F-D) curve and estimated F-D curves with and without a slope constraint in the optimization process.

Figure 2.

Clamped-clamped test problem.

Figure 2.

Clamped-clamped test problem.

Figure 3.

Four force–displacement curves measured from the initial layout (–O–), prescribed in the optimization problem as target points (–∇–), measured from the optimal layout without slope constraints (–□–) and measured from the optimal layout with slope constraints (–☆–).

Figure 3.

Four force–displacement curves measured from the initial layout (–O–), prescribed in the optimization problem as target points (–∇–), measured from the optimal layout without slope constraints (–□–) and measured from the optimal layout with slope constraints (–☆–).

Figure 4.

History plots of objective function (–✽–) and volume fraction (–O–) for the optimization procedure. (a) With slope constraint; (b) Without slope constraint.

Figure 4.

History plots of objective function (–✽–) and volume fraction (–O–) for the optimization procedure. (a) With slope constraint; (b) Without slope constraint.

Figure 5.

Force–displacement curves measured from the initial design domain (⋯) and optimal layouts (—) with circles (○) representing the target points. (a) Softening; (b) Hardening.

Figure 5.

Force–displacement curves measured from the initial design domain (⋯) and optimal layouts (—) with circles (○) representing the target points. (a) Softening; (b) Hardening.

Figure 6.

Deformed layouts of optimal topology at target points. (a) Softening; (b) Hardening.

Figure 6.

Deformed layouts of optimal topology at target points. (a) Softening; (b) Hardening.

Figure 7.

Force–displacement curves measured from the optimization processes with increasing the number of target points without the slope constraints. For three target points case (–□–), for six target points case (–ж–), for nine target points case (–□–), for twelve target points case (–○–), for fifteen target points case (–✡–), and measured curve from the optimal layout in Figure 6a (—).

Figure 7.

Force–displacement curves measured from the optimization processes with increasing the number of target points without the slope constraints. For three target points case (–□–), for six target points case (–ж–), for nine target points case (–□–), for twelve target points case (–○–), for fifteen target points case (–✡–), and measured curve from the optimal layout in Figure 6a (—).

Figure 8.

Force–displacement curve comparisons between optimal layouts from the optimization process (—□—) and commercial software COMSOL (—○—). (a) Softening; (b) Hardening.

Figure 8.

Force–displacement curve comparisons between optimal layouts from the optimization process (—□—) and commercial software COMSOL (—○—). (a) Softening; (b) Hardening.

Figure 9.

Clamped-free test problem.

Figure 9.

Clamped-free test problem.

Figure 10.

Force–displacement curves measured from the initial design domain (∙∙∙∙) and optimal layouts (—) with circles (○) representing the target points. (a) Softening; (b) Hardening.

Figure 10.

Force–displacement curves measured from the initial design domain (∙∙∙∙) and optimal layouts (—) with circles (○) representing the target points. (a) Softening; (b) Hardening.

Figure 11.

Deformed layouts of optimal topology at target points. (a) Softening; (b) Hardening.

Figure 11.

Deformed layouts of optimal topology at target points. (a) Softening; (b) Hardening.

Lee, J.; Detroux, T.; Kerschen, G.

Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints. Appl. Sci.2020, 10, 2676.

https://doi.org/10.3390/app10082676

AMA Style

Lee J, Detroux T, Kerschen G.

Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints. Applied Sciences. 2020; 10(8):2676.

https://doi.org/10.3390/app10082676

Chicago/Turabian Style

Lee, Jongsuh, Thibaut Detroux, and Gaëtan Kerschen.

2020. “Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints” Applied Sciences 10, no. 8: 2676.

https://doi.org/10.3390/app10082676

Find Other Styles

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

No

No

Article Access Statistics

For more information on the journal statistics, click here.

Multiple requests from the same IP address are counted as one view.

Lee, J.; Detroux, T.; Kerschen, G.

Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints. Appl. Sci.2020, 10, 2676.

https://doi.org/10.3390/app10082676

AMA Style

Lee J, Detroux T, Kerschen G.

Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints. Applied Sciences. 2020; 10(8):2676.

https://doi.org/10.3390/app10082676

Chicago/Turabian Style

Lee, Jongsuh, Thibaut Detroux, and Gaëtan Kerschen.

2020. “Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints” Applied Sciences 10, no. 8: 2676.

https://doi.org/10.3390/app10082676

Find Other Styles

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely

those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or

the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas,

methods, instructions or products referred to in the content.