Enforcing a Force–Displacement Curve of a Nonlinear Structure Using Topology Optimization with Slope Constraints

: The objective of this study is to develop an optimization methodology to ﬁnd 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 di ﬀ erence 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 ﬁlter is employed to suppress the gray region observed between the solid and void regions, which generally occurs when using a density-based ﬁlter. 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 ﬁxed mesh in the topology optimization process. To validate the outlined methodology, several case studies with di ﬀ erent types of nonlinearity and structural features of the obtained layouts are investigated.


Introduction
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 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].

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: : (1) c is the objective function, and ρ is a vector of design variables, N t is the number of defined target points on the F-D curve, f i is the amplitude of the applied ith external force, u is the displacement vector in global coordinate, u * i is the prescribed value of the displacement with respect to the ith external force, and w c i is the weighting function. l 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. θ c i 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.
θ min i and θ max i are the upper and lower bounds of the slope constraint,ρ is a vector of physical density, and V 0 is the volume when the entire design domain is occupied by solid elements. V max is the upper bound of the volume fraction constraint, and r 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 S: The fourth-order tensor D is the elastic tensor, and E is the Green-Lagrangian strain tensor, which is represented as E = F T F − I /2. I is the identify matrix. F is the deformation gradient tensor defined as F = I + ∂u/∂X. X 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: f ext = f · p is the external force vector, with p 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 f is the forcing amplitude. f int is the internal nodal force vector, which is expressed by the summation of the elements of the internal nodal force vector (i.e., f int = e f int e ). The element nodal force vector is expressed as: where φ e is the stored elastic energy density in the element, and u e 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 K T 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].

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 ρ i is the filtered density of element i, ρ j is the design variable of element j, x j is the location of element j, υ j is the volume of element j, and N e,i is the number of the set of the neighborhood element i within a certain filter radius r. w x j is a weighting function specified by w x j = r − x j − x e . 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. ρ e 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 E 1 is the Young's modulus of the solid phase (the black regions), and E 0 is the Young's modulus of the void phase (the white regions), which generally have a relation of E 0 = 10 −9 ·E 1 . 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 φ L is the stored elastic energy density under small deformation. α e 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., φ e = φ(u e )E 1 ), whereas void elements are treated as linear elements (i.e., φ e = φ L (u e )E 0 ).
The interpolation factor α e is expressed similarly to the threshold projection filter in Equation (9) using a smoothed Heaviside function: where ρ 0 is a threshold between 0 and 1. If β 1 is large enough and the physical design density ρ e is close to 1, then the interpolation factor α e will be close to 1. Therefore, the nonlinear analysis is imposed on the corresponding element e. However, if ρ e is close to 0, then the interpolation factor α e 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:

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 P 1 , P 2 , and P 3 . First, the definable slope at the points is the tangential slope, which are denoted by K P1 T , K P2 T , and K P3 T . 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 θ P1 k , θ P2 k , and θ P3 k . 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 θ P k can be estimated at the same time as when the target points are defined. k k k 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   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 θ P1 c , θ P2 c , and θ P3 c . The slope θ P k is associated with the secant stiffness matrix, and the slope θ P c 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 K P T and θ P k , 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 θ P k and θ P c . 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 K is the secant stiffness matrix, which is composed of the linear stiffness matrix K L and the nonlinear stiffness matrix K NL . If a structure considered in the optimization problem is displaced by a concentrated force, (i.e., a single input case, in Equation (14), f ext i = f and f ext n = 0 ; n = 1 ∼ N , n i), then the slope θ P c ij 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 θ P c is associated with the inverse of the matrix K, and its sensitivity can be derived by the direct differentiation method (DDM): Equation (16) represents the sensitivity of the inverse of the stiffness matrix K 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).

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.

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. Let us first consider the case of a targeted softening nonlinearity. First, we measured 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 w c i 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 β max = 64 and an initial value of β = 1. The threshold η is fixed as η = 0.5. The threshold parameters of the energy interpolation scheme in Equation (12) are fixed as β 1 = 500 and ρ 0 = 0.01. The filter radius r is fixed as r filter = 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 (-Ο-) 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. . Four force-displacement curves measured from the initial layout (-Ο-), 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 ☆.
The curve measured from the initial design layout is depicted by a dash and circle line (-Ο-) 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. 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 (-Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 o Figure 3. Four force-displacement curves measured from the initial layout (-Ο-), 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 ☆.
The curve measured from the initial design layout is depicted by a dash and circle line (-Ο- Figure 3. The line seems straight, but it is slightly curved upward, meaning that the initial des behaves with hardening nonlinearity. As shown in Figure 3, on the one hand, the measured curv ☆-) from the optimization process, in which the slope constraint was incorporated, and prescribed target curve (-∇-) are in good agreement with each other. On the other hand, the measu curve (-□-) from the optimization process, in which the slope constraint was not incorporat obviously differs from that of the target curve. In addition, it seems straight as for the initial lay case but is also slightly curved upward, meaning that the inherent hardening nonlinearity that w in the initial layout remained during the optimization process even though the targeted nonlinea was softening. It can thus be seen that the optimization process has tried to reduce only the object value while maintaining the inherent nonlinearity in the initial layout. Therefore, it can be conclud that the slope constraint is necessary if different nonlinearity behaviors from the inherent nonlinea 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 a without slope constraints in the optimization process. Figure 4 shows the history plot for normalized objective function with respect to the objective function value on the initial design and volume fraction of each iteration in the optimization process; Figure 4a is the case where slope constraint is incorporated in the optimization process and Figure 4b is the case where the sl constraint is not incorporated. The total optimization process is repeated until 185 iterations and iterations, respectively.
-) 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 (c 0 ) 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. Appl. Sci. 2020, 10, x FOR PEER REVIEW  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  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 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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 22 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 N TP 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.
f N TP = 0 : 150 N TP : 150 , N TP = 3 : 3 : 15 (17) For N TP = 3, for instance, the forces f N TP 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 N TP = 6, the forces f N TP considered in the optimization problem were 0, 25, 50, 75, 100, 125, and 150 N. Thereafter, the displacements were measured by applying defined forces f N TP to the forcing position of the resultant layout in Figure 6a. The measured displacements were used as u * i 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. Appl. Sci. 2020, 10, x FOR PEER REVIEW  (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  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.
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 e  is less than 0.5 is set to 0 and whose design variable e  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 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 e ρ is less than 0.5 is set to 0 and whose design variable e ρ 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. 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 e ρ is less than 0.5 is set to 0 and whose design variable e ρ 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.
-), and measured curve from the optimal layout in Figure 6a (-).
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 ρ e is less than 0.5 is set to 0 and whose design variable ρ e 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. 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.

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 r filter = 2 × the element size for softening case and r filter = 6 × the element size for the hardening case.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 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. 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 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.

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.