Topology Optimization of Elastoplastic Behavior Conditions by Selectively Suppressing Plastic Work

Abstract

This work conducted topology optimization with an implicit analysis of elastoplastic constitutive equation in order to design supporting structures for unexpected heavy loading conditions. In this topology optimization model, plastic work was extracted from strain energy and selectively employed in the objective function according to deformation mode. While strain energy was minimized in elastic deformation areas, in elastoplastic deformation areas, the plastic work was minimized for the purpose of suppressing plastic deformation. This method can focus on suppressing plastic strain in the plastic deformation zone with maintaining elastic stiffness in the elastic deformation zone. These formulations were implemented into MATLAB and applied to three optimization problems. The elastoplastic optimization results were compared to pure elastic design results. The comparison showed that structures designed with accounting for plastic deformation had a reinforced area where plastic deformation occurs. Finally, a finite element analysis was conducted to compare the mechanical performances of structures with respect to the design method.

1. Introduction

Topology optimization is a scheme used to obtain a design with an optimized shape for a specific purpose, and is widely used in engineering problems [1,2,3,4,5,6,7,8,9,10,11]. The principle of topology optimization is to find an optimal distribution of material density within a design domain under constraint conditions. In structure design problems, strain energy is the first priority to be the objective function since minimizing strain energy can give a structure with a maximum stiffness. Under elastic deformation, calculating the strain energy is simple since stress and strain have a simple linear relation. In contrast, when a part of a material undergoes elastoplastic deformation, its material properties change according to the level of plastic strain, resulting in a nonlinear relation between stress and strain. This nonlinear relation leads to difficulties when attempting to determine stress and the material property matrix. Because of this, the majority of topology optimization applications are based on the assumption of elastic deformation [12,13,14,15,16,17,18]. Only a few researchers have tried to consider plastic deformation in topology optimization. Zhang et al. [19] employed enhanced assumed strain (EAS) elements to solve the locking issue, and Wallin et al. [20] made a multiplicative split of the deformation gradient. Schwarz et al. [21] conducted topology optimization for elastoplastic material by maximizing ductility. Blachowski et al. [22] recently proposed a stress intensity-driven elastoplastic topology optimization method to find a minimum weight structure that is able to carry a given load. Herfelt et al., [23] conducted a strength-based optimization resulting in a convex optimization which ensures that any local optimum is the global optimum. Alberdi and Khandelwal [24], Huang et al. [25], and Li et al. [26,27] have tried to design energy absorbing structures in order to absorb plastic energy. These studies are based on maximizing plastic work in order to dissipate plastic energy by using the designed energy absorbing structures.

This work studies a topology optimization method for designing supporting structures with considering elastoplastic deformation by suppressing plastic deformation. Since supporting structures need to withstand external loads for a long time, the structures are typically designed to be employed within a range of elastic deformation. However, when an unexpected load is applied to a part of the structure, that part can go through elastoplastic deformation. When the level of plastic deformation increases, stress is concentrated on the plastic deformation zone. This stress concentration accelerates the collapse of the structure by the plastic hinge phenomenon. For this reason, it is safer to design supporting structures with considering elastoplastic deformation in order to suppress the concentration of plastic deformation. For suppressing the stress concentration on elastoplastic deformation zones, this work separated plastic work from strain energy and selectively applied the plastic work in the objective function. The plastic strain and work are calculated based on the isotropic hardening model. In elastic deformation areas, strain energy was minimized as in general cases. In elastoplastic deformation areas, on the other hand, only the plastic work was employed as the objective function and minimized for the purpose of suppressing plastic deformation. This scheme is able to focus on suppressing plastic deformation in the plastic deformation zone while still maintaining elastic stiffness in the elastic deformation zone. For the optimization process in this work, the global equilibrium and constitutive condition were solved by implicit analysis based on the nonlinear weak form of the finite element method (FEM), which is widely used in structural analysis [28,29]. For topology optimization, the solid isotropic material with penalization (SIMP) [30,31,32] and Optimality Criteria (OC) methods were employed based on the reference of Sigmund’s scheme [18]. All of the formulations were implemented into MATLAB and applied to three optimization problems. Two different optimization methods, optimization under elastic deformation assumption and optimization accounting for elastoplastic deformation, led to different design shapes. A structural analysis with the FEM was also conducted to compare the mechanical performance of the designed shapes with respect to a design method.

2. Formulations

2.1. Implicit Analysis of Global Equilibrium

The weak form of the structural analysis in the topology optimization is given as below:

$$\begin{array}{l}\mathit{G}\left(\mathit{u}\right)={\mathit{F}}_{int}\left(\mathit{u}\right)-{\mathit{F}}_{ext}=\mathbf{0},\\ \mathrm{where} {\mathit{F}}_{int}\left(\mathit{u}\right)={\int }^{ }{\mathit{B}}^T\mathbf{\sigma }d\nu \mathrm{and} {\mathit{F}}_{ext}={\int }^{ }{\mathit{N}}^T\overline{\mathit{t}}ds+{\int }^{ }{\mathit{N}}^T\mathit{b}d\nu .\end{array}$$

(1)

G(u) is the global equilibrium function and should be zero to satisfy the equilibrium. F_int(u) is the internal force vector and F_ext is the external force vector. u presents the displacement tensor, B is the strain-displacement matrix, so-called the B matrix in the FEM, and N is the shape function matrix. σ is the stress tensor, $\overline{\mathit{t}}$ is the traction tensor and b is the body force vector. Since the body force can be negligible in many engineering problems, this work does not account for the body force. In an elastic problem, the mechanical property can be assumed to be a constant and F_int(u) can be simplified as below:

$${\mathit{F}}_{int}\left(\mathit{u}\right)={\int }^{ }{\mathit{B}}^{\mathit{T}}\mathbf{\sigma }d\nu =\left[{\int }^{ }{\mathit{B}}^{\mathit{T}}\mathit{D}\mathit{B}d\nu \right]\mathit{u}={\mathit{K}}_e\mathit{u}.$$

(2)

D means the elastic modulus matrix and K_e is the stiffness matrix of elastic deformation. In this case, G(u) can be easily solved as a linear equation. However, in elastoplastic problems, G(u) is a nonlinear equation because the material property changes with the level of plastic deformation. In this case, a linearization method is required to solve G(u) and this work employs the Newton–Rapson scheme as below:

$$\mathit{G}{\left(\mathit{u}\right)}_{i+1}=\frac{\partial \mathit{G}}{\partial \mathit{u}}\left({\mathit{u}}_i\right)d{\mathit{u}}_{i+1}+\mathit{G}\left({\mathit{u}}_i\right)=\mathbf{0}.$$

(3)

i denotes iteration number in the Newton–Rapson scheme. Substituting Equation (1) into Equation (3) leads to Equation (4) for the implicit analysis of global equilibrium.

$$\begin{array}{l}\mathit{G}{\left(\mathit{u}\right)}_{i+1}={\mathit{K}}_{ep}\left({\mathit{u}}_i\right)d{\mathit{u}}_{i+1}-{\mathit{F}}_{ext}+{\mathit{F}}_{int}\left({\mathit{u}}_i\right)=\mathbf{0}, \\ \mathrm{where} {\mathit{K}}_{ep}\left({\mathit{u}}_i\right)={\int }^{ }{\mathit{B}}^T{\mathit{C}}^t\mathit{B}d\nu +{\int }^{ }{\mathit{H}}^T\mathit{\sigma }\mathit{H}d\nu .\end{array}$$

(4)

K_ep(u_i) is the stiffness matrix of elastoplastic deformation, ${\mathit{C}}^t$ represents the consistent tangent modulus which changes according to plastic strain. H is the partial derivative matrix of the shape function. In order to solve Equation (4), stress and the consistent tangent modulus at i-th iteration should be determined.

2.2. Implicit Stress Integration

In order to determine $\mathbf{\sigma }$ in F_int(u_i) and K_ep(u_i) in Equation (4), the increment of stress ($d\mathbf{\sigma }$) is given as below:

$$d\mathbf{\sigma }=\mathit{D}d{\mathbf{\epsilon }}^e, \mathrm{where} d{\mathbf{\epsilon }}^e=d\mathbf{\epsilon }-d{\mathbf{\epsilon }}^p.$$

(5)

$d{\mathbf{\epsilon }}^e$ is the increment of elastic strain, $d\mathbf{\epsilon }$ is the increment of total strain, and $d{\mathbf{\epsilon }}^p$ denotes the increment in plastic strain. In the numerical analysis, Equation (5) is accumulated with the return mapping method as below:

$${\mathbf{\sigma }}_{n+1}={\mathbf{\sigma }}_n+d{\mathbf{\sigma }}_{n+1}={\mathbf{\sigma }}^T-\mathit{D}d{\mathbf{\epsilon }}_{n+1}^p, \mathrm{where} {\mathbf{\sigma }}^T={\mathbf{\sigma }}_n+\mathit{D}d{\mathbf{\epsilon }}_{n+1}.$$

(6)

${\mathbf{\sigma }}^T$ is the trial stress, and n is the step number in the return mapping method. In order to determine $d{\mathbf{\epsilon }}_{n+1}^p$, this work employs the associated rule, as below:

$$d{\mathbf{\epsilon }}_{n+1}^p={\left(\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}\right)}_{n+1}d{\overline{\epsilon }}_{n+1}^p,$$

(7)

where $\overline{\sigma }$ represents the effective stress, and $d{\overline{\epsilon }}^p$ denotes the increment of effective plastic strain. In order to determine $d{\overline{\epsilon }}_{n+1}^p$, the following three conditions (yield condition, stress update, and work-hardening) should be satisfied in elastoplastic deformation [33,34].

$${\mathrm{g}}_1=\overline{\sigma }\left({\mathbf{\sigma }}_{n+1}\right)-{\phi }_{n+1}\left({\overline{\epsilon }}_{n+1}^p\right)=0,$$

(8)

$${\mathbf{g}}_2={\mathit{D}}^{-1}\left({\mathbf{\sigma }}_{n+1}-{\mathbf{\sigma }}^T\right)+{\left(\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}\right)}_{n+1}d{\overline{\epsilon }}_{n+1}^p=0,$$

(9)

$${\mathrm{g}}_3=d{\phi }_{n+1}{H}^{-1}-d{\overline{\epsilon }}_{n+1}^p=0.$$

(10)

where $\phi$ is the flow stress and H is a work-hardening parameter in the stress–strain curve. This work employs the isotropic hardening law defining H as the slope in the stress–strain curve at the respective time. Solving the above equations leads to the equation below.

$$\begin{array}{l}d{\overline{\epsilon }}_{n+1}^p=\sum \frac{{\mathrm{g}}_1-\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}\left(\mathit{E}\right){\mathbf{g}}_2+{\mathrm{g}}_{3}H}{\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}\left(\mathit{E}\right)\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}+H},\\ \mathrm{where} \mathit{E}={\left[{\mathit{D}}^{-1}-\frac{{\partial }^2\overline{\sigma }}{\partial {\mathbf{\sigma }}^2}\left(d^2{\overline{\epsilon }}_{n+1}^p\right)\right]}^{-1}.\end{array}$$

(11)

With Equations (6) and (11), σ can be determined and substituted into F_int(u) and K_ep(u_i) terms in Equation (4).

2.3. Consistent Tangent Modulus

The ${\mathit{C}}^t$ term in K(u_i) should also be determined in the implicit analysis. The consistency condition brings the following equation:

$$\frac{\partial \overline{\sigma }}{\partial \mathit{\sigma }}d\mathbf{\sigma }-h^{\prime }d{\overline{\epsilon }}^p=0, \mathrm{where} h^{\prime }=\frac{d\phi }{d{\overline{\epsilon }}^p}.$$

(12)

With Equations (6), (7) and (12), it results in the following equation:

$$d{\overline{\epsilon }}^p=\frac{\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}E}{\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}E\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}+h^{\prime }}d\mathbf{\epsilon }.$$

(13)

With Equations (5) and (13), the relation of the consistent tangent modulus, ${\mathit{C}}^t$, can be given:

$$d\mathit{\sigma }={\mathit{C}}^{t}d\epsilon , \mathrm{where} {\mathit{C}}^t=\left[\mathit{E}-\frac{\mathit{E}\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}⨂\mathit{E}\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}}{\frac{\partial \overline{\sigma }}{\partial \mathit{\sigma }}\mathit{E}\frac{\partial \overline{\sigma }}{\partial \mathbf{\sigma }}+h^{\prime }}\right].$$

(14)

2.4. Hardening Model and Material Parameters

An appropriate hardening model should be employed to determine $\varphi$ in Equations (8) and (10). This work selected a modified Hockett–Sherby hardening model [35], which has been validated in some engineering problems [36,37,38]:

$$\phi \left({\overline{\epsilon }}^p\right)=A-Bexp\left(-C{\left({\overline{\epsilon }}^p\right)}^b\right)+D{\overline{\epsilon }}^p,$$

(15)

where A, B, C, b, and D are the constants of the modified Hockett–Sherby model. In this work, a mild steel was selected as the material. The values of the constants in the Hockett–Sherby hardening model for the material are shown in

Table 1. Material parameters of the modified Hockett–Sherby model.

Parameters of the Modified Hockett–Sherby Model
A (MPa)	B (MPa)	C	b	D (MPa)
568.74	368.19	1.81	0.60	78.1

and the stress-strain curve of the Hockett–Sherby hardening model for the material is shown in Figure 1. Young’s modulus was 200 GPa and Poisson’s ratio was 0.3 in this work.

2.5. Optimization

Minimizing strain energy is given as below:

$$\begin{array}{}\mathrm{(16)}& min :c\left(x\right)={\sum }_{e=1}^m{\left(x_e\right)}^{pp}{\iint }^{ }\left(\mathbf{\sigma }d\mathbf{\epsilon }\right)d{\nu }_e,& subject to :d=\frac{{\sum }_{e=1}^m{x}_e}{m}& {\mathrm{g}}_1=0& \mathit{G}{\left(\mathit{u}\right)}^k=\mathbf{0}& x_{min}\le x_e\le 1.\end{array}$$

In the optimization, the SIMP [30,31,32] was applied and a standard OC method [18,39] was used to solve Equation (16). c(x) is the objective function, $x_e$ represents the density of an element in the topology optimization; the subscript e means an element. $x_{min}$ is a minimum value of $x_e$ to prevent singularity problem ($x_{min}=0.0001$ in this work). m is the total number of elements and pp is the penalization power ($pp=3$ in this work). d is a volume fraction, the ratio of the final volume to the initial volume ($d=0.5$ in this work). Equation (16) is optimized based on the final strain of each time step and repeated until the target time step (details are explained in Figure 2). G(u)^k is the global equilibrium equation in the implicit analysis of Equation (4) and the superscript k is the time step number of the global loading time. ${\mathrm{g}}_1=0$ means satisfying the constitutive equation shown in Equation (8). In order to concentrate on suppressing plastic deformation, the strain energy can be separated into two parts, elastic energy ($W^e$) and plastically dissipated work ($W^p$) as below:

$$\mathbf{\sigma }d\mathbf{\epsilon }=\mathbf{\sigma }\left(d{\mathbf{\epsilon }}^e+d{\mathbf{\epsilon }}^p\right)=d{W}^e+d{W}^p.$$

(17)

Then, the objective function can be modified as below:

$$min :c\left(x\right)={\left(x_e\right)}^{pp}{\sum }_{e=1}^m{\iint }^{ }\left[\mathbf{\sigma }d{\mathbf{\epsilon }}^e\cdot w\left(x\right)+\mathbf{\sigma }d{\mathbf{\epsilon }}^p\right]d{\nu }_e,$$

(18)

$$\begin{array}{l}subject to :d=\frac{{\sum }_{e=1}^m{x}_e}{m}\\ {\mathrm{g}}_1=0\\ \mathit{G}{\left(\mathit{u}\right)}^{\mathit{K}}=\mathbf{0}\\ x_{min}\le x_e\le 1 ,\\ \mathrm{where} \{\begin{array}{c}w\left(x\right)=1if d{W}^p=0\\ w\left(x\right)=0otherwise\end{array}.\end{array}$$

In elastic deformation elements, minimizing strain energy is applied as in general topology optimization cases. On the other hand, in elastoplastic deformation elements, only plastic work is employed as the objective function in order to concentrate on minimizing plastic work. Then, the sensitivity of the objective function can be determined as below (refer to [18,21]):

$$\begin{array}{l}\frac{\partial c\left(x\right)}{\partial x_e}=\frac{\partial }{\partial x_e}\left({\iint }^{ }\left[{\left(x_e\right)}^{pp}\mathbf{\sigma }d{\mathbf{\epsilon }}^e\cdot w\left(x\right)+{\left(x_e\right)}^{pp}\mathbf{\sigma }d{\mathbf{\epsilon }}^p\right]d{\nu }_e\right),\\ \mathrm{where} \{\begin{array}{c}w\left(x\right)=1 if d{W}^p=0\\ w\left(x\right)=0otherwise\end{array}.\end{array}$$

(19)

For the mesh-independency filter [18], the sensitivity is modified as below:

$$\begin{array}{l}{\frac{\partial c\left(x\right)}{\partial x_e}|}_{modified}=\frac{1}{x_e{\sum }_{b=1}^m{T}_b}{\sum }_{b=1}^m{T}_b{x}_b\frac{\partial c\left(x\right)}{\partial x_b},\\ \mathrm{where} T_b=\left(r_{min}-l\left(e,b\right)\right)\cdot q\left(l\left(e,b\right)\right),\\ \{\begin{array}{c}q\left(l\left(e,b\right)\right)=1if l\left(e,b\right)\le r_{min}\\ q\left(l\left(e,b\right)\right)=0otherwise\end{array}.\end{array}$$

(20)

l(e,b) is the distance between two elements e and b. r_min is the radius of the filter and the radius is 1.5 in this work. The optimization problem is solved with the OC method to find a new density [39]. The formulations were implemented into MATLAB and Figure 2 presents a flow chart of the MATLAB code. At each time step, the FEM simulation is performed first to solve elastoplastic behavior based on Equations (1)–(15). From the FEM result, the optimization process is conducted. In the optimization process, the objective function is obtained and the mesh-independency filter is applied. Then the density of element $x_e$ is updated at each iteration, as shown in Figure 2. The steps of the optimization are based on Equations (18)–(20), and details are explained in [18,30,31,32]. When the optimization process converges, it moves to the next time step in order to start the FEM simulation for the increased loading condition. This process is repeated until the target time step.

3. Problems and Results

Figure 3 presents the definition of three design problems. Figure 3a is the well-known Messerschmitt–Bölkow–Blohm (MBB) beam problem, which is problem A in this work, and Figure 3b shows a tip loaded cantilever beam problem—problem B. Figure 3c presents a structure with a spring supporter—problem C. The load increases with the loading time step and reaches the boundary condition shown in Figure 3 at the last time step 1.0. The time increment is set Δt = 0.1 in all the examples. Figure 4, Figure 5 and Figure 6 show the results of the design problems. In Figure 4, Figure 5 and Figure 6, each subfigure (a) is the result based on the assumption of elastic deformation, and (b) represents the result of the presented method, the selectively applied plastic work in elastoplastic deformation. In problem A, the designed structure based on elastic deformation, design A–A in this work, has the well-known shape [18] as shown in Figure 4a. The design with elastoplastic deformation (design A–B) shows a slightly different design with a thicker top part, marked A in Figure 4b, since the part A has more chances to go through plastic deformation. In problem B, the result under elastic deformation (design B–A), shown in Figure 5a, represents the well-known shape [40]. The design with the presented method (design B–B) brings a strengthened design to the area, marked B, as shown in Figure 5b. In problem C, each design has a quite different shape. The design under elastoplastic deformation—design C–B in Figure 6b—strengthens the stick part which adjoins the spring (denoted C in Figure 6b) when it is compared with the design under elastic deformation, design C–A in Figure 6a. In addition, design C–B includes a horizontal bar to make D part stronger, as shown in Figure 6b.

4. Discussion

In order to analyze the mechanical performance of the designed shapes, an FEM model of each design was generated with the four-node linear elements, as shown in Figure 7. The boundary condition was the same for each problem shown in Figure 3. The FEM simulation was conducted with the MATLAB code as mentioned above. Figure 8, Figure 9 and Figure 10 show the plastic strain distribution and the displacement field in each model. Figure 8a, Figure 9a and Figure 10a show that the designed structures, under the elastic deformation assumption, lead to large plastic strain zones in the three problems. Since the plastic deformation zones are generated on narrow load paths, the plastic deformation zone could easily result in a plastic hinge, which is a cause of structure collapse. On the other hand, the designed structures in the elastoplastic deformation can drastically reduce the plastic strain, as shown in Figure 8c, Figure 9c and Figure 10c. In terms of the magnitude of displacement, designs A–A and A–B lead to a similar contour of displacement field in Figure 8b,d; however, the maximum magnitude of displacement of design A–B decreases to less than half compared to that of design A–A, as shown in Figure 11. In the same manner, the displacement magnitude of design B–B is about half of that of design B–A as shown in Figure 11, although the displacement contours are very similar in Figure 9b,d. In problem C, a different tendency appears in the displacement analysis. Design C–B not only has a lower value of displacement (about 97% reduction as shown in Figure 11), but also the contours of the displacement field are different when compared to design C–A, as shown in Figure 10b,d. In design C–A, the maximum displacement is generated on the stick part connected with the spring. However, design C–B has maximum displacement on the top and left part. This is because the stick part of design C–A is strengthened by design C–B.

In order to clearly present the elastic and elastoplastic zones, Figure 12 separates the two zones by two colors for the three examples. The green and red colors present elastic and elastoplastic zones, respectively. In all the examples, the designs based on the elastoplastic deformation, shown in Figure 12b,d,f, clearly reduce the elastoplastic zones compared to the designs based on the pure elastic deformation, as shown in Figure 12a,c,e.

For more discussion, the distribution of plastic strain is analyzed in detail in Figure 13. For the analysis, each section was set on a plastic deformation zone in each problem: section O-O’ for problem A in Figure 8a,c, section P-P’ for problem B in Figure 9a,c, and section Q-Q’ for problem C in Figure 10a,c. As shown in Figure 13, the new designs with elastoplastic deformation result in drastically reduced plastic strain in all of the problems. Through each section line, the level of plastic strain decreases on average more than 80%, and the curve of the plastic strain distribution becomes flatter in the structures designed with the selectively employed plastic work.

Figure 14 shows the dissipated plastic work in the structures caused by plastic deformation according to analysis time. The time is a normalized time, when the time of the final step in the implicit analysis is 1. In all of the problems, the designed structures with elastoplastic deformation can delay and reduce the dissipation of plastic energy. In design C–B, the plastic work in the structure is almost zero. These results mean that the structures designed by selectively applying plastic work have less chances of collapsing by a plastic hinge. In short, the designed structures which account for elastoplastic behavior can suppress plastic deformation and reduce the chances of being collapsed.

In Figure 15, with the MMB design based on the elastoplastic deformation (design A–B), it is checked whether the designed structure behaves well under other loading conditions, than the design condition, by the FEM simulation. In this study, a critical collapse loading condition (f_collapse = 360 N) is defined because this condition leads to that the entire central cross section of the beam undergoes plastic deformation, as shown in Figure 15a. Note that the load is distributed within 3 mm as in the original example in Figure 3. Figure 15a–d present that as the applied load decreased, the elastoplastic zone decreases. When the load is below 0.55f_collapse, almost all regions of the beam are in elastic condition, as shown in Figure 15d. The results show that the designed structure can be used in other loading conditions. Note that the condition of f = 0.83f_collapse (in Figure 15c) is the original condition, as shown in Figure 3.

Finally, the stability of the optimization process is discussed. Since the plastic work is separated from the strain energy in elastoplastic deformation zone, the objective function is not perfectly continuous. The plastic work calculation is affected by the time increment and loading boundary condition, so that the stability of the numerical procedure is also affected by them. Figure 16 compares the change in the normalized objective function for different time increments and loading conditions in the MMB example. In this study, the original time increment used in the MMB example (in Figure 3) is set as the reference time increment (Δt = Δt_ref.). For the load condition, the collapse loading condition in Figure 15 (f_collapse = 360 N) is also set as the reference loading; the load is distributed within 3 mm as in the original example in Figure 3. Three cases are compared in Figure 16. The first case is the same to the original example (Δt = Δt_ref. and f = 0.83f_collapse) of Figure 3. The second case has the condition of Δt = 1.5Δt_ref. and f = 0.83f_collapse to check the time increment effect. The third case is designed by the condition of Δt = Δt_ref. and f = 0.9f_collapse. The objective function is normalized by the initial value of it for each case, and the normalized time 1.0 is the time to reach the target loading condition. The maximum number of iteration at each time step is limited to 100 in this study. As shown in Figure 16, the different conditions do not lead to significant change in design results. Because the volume fraction constraint and the position of elastoplastic zone are the same. For the stability, at the beginning of the optimization for elastic deformation, the three curves are similar. The optimization of the larger time increment (blue curve) case shows a tendency to be minimized, but has unstable fluctuations during elastoplastic deformation. The orange curve for the larger load condition, close to the collapse load, shows more unstable fluctuations during minimizing the objective function. The black curve of the original condition is reliably optimized in the elastoplastic design process. Based on these results, it is recommended to use the appropriate time step and loading condition in optimization processes.

5. Conclusions

This work conducted topology optimization to design supporting structures using an implicit analysis of elastoplastic deformation. The details are explained below:

(1)

The optimization algorithm was combined with the nonlinear weak form of the finite element method.

(2)

In the objective function, plastic work was separated from stain energy and the separated plastic work was selectively applied according to a deformation mode. In elastic deformation areas, strain energy was minimized, as in general cases. In elastoplastic deformation areas, only the plastic work was minimized for the purpose of suppressing plastic deformation. This method is able to focus on suppressing plastic deformation in the elastoplastic deformation areas while still accounting for elastic stiffness in the elastic deformation areas.

(3)

The structures designed while considering elastoplastic deformation strengthened the areas where plastic deformation occurs. The reinforced structures drastically reduced plastic strain and the magnitude of displacement when compared to shapes designed with only elastic deformation.

(4)

The structures designed with elastoplastic deformation can delay the beginning of plastic deformation and reduce plastically dissipated energy in the structures. This reduction and delay of the plastic deforming can decrease chances of plastic collapse.

(5)

The time step and loading condition affects the stability of the optimization. It is recommended to use appropriate time step and loading condition in optimization processes.

Based on the above results, this study shows that the selectively applied plastic work in the objective function can result in the design of safer supporting structures which reduce plastic deformation under conditions of large deformation.