text:
\section*{Abstract}
\re{T}opology optimization is an important basis for the \re{design} of components\re{.
Here,} the optimal structure is found \re{within a design space subject to boundary conditions} as well as the material law.
\re{Additionally, the specific material law has a strong impact on the final design.
Even more: a, for instance, linear-elastically structure is not optimal if plastic deformation will be induced by the loads.}
\re{Hence,} a physically correct and resource-efficient \re{inclusion of plasticity modeling is needed}.
In this contribution, we present an extension of the thermodynamic topology optimization that accounts for \re{the} non-linear material behavior due to the evolution of plastic strains.
For this purpose, we develop a novel surrogate plasticity model that allows to compute the \re{correct} plastic strain tensor corresponding to the current structure design.
We show \re{the} agreement \re{of the model} with the classic plasticity model without dissipation and that the interaction of the topology optimization with plastic material behavior results in structural changes.

\section{Introduction}
Engineers are always looking for structures that meet the specific requirements in an optimal way.
One possibility for finding these structures is provided by optimization schemes which are classified as follows: i) improving the principal idea, ii) modifying the material, iii) thickness dimensioning, iv) optimization of shape and v) optimization of topology \cite{schumacher2013Optimierung,harzheim2008strukturoptimierung}.
Herein, the optimization scheme that demands the minimum amount of restrictions is given by topology optimization.
The consideration of the real materials properties offers additional potential for the optimal design of components.
Therefore, it is important to account for the physical material behavior even during the process of topology optimization.

There are various variants of topology optimization available as, e.\,g.{}, the optimization for temperature evolution, noise reduction, dynamic response, or structural stiffness.
All of these approaches have in common that the related physical balance laws, in most cases the balance of linear momentum, are solved along with a mathematical optimization problem which is given in terms of an objective function.
The most common objective is the minimization of compliance, i.\,e.{} the maximization of stiffness according to a target structure volume \cite{sigmund2013topology}.
Therefore, topology optimization determines the position and arrangement of material within a given design space and boundary conditions such that the objective function is minimized.
The topology of a structure can be parameterized via different approaches during the specific numerical investigation.
For the numerical solution of the physical balance law, e.\,g.{}, the balance of linear momentum, usually the finite element method (FEM) is employed.
Consequently, the finite elements introduce a discretization of the design space, and it is thus most common to assign a density value for each discrete subvolume, i.\,e.{} for each finite element.
For this assignment, a variety of different numerical schemes has been developed among which the probably most popular is given by ``Solid Material with Penalization'' (SIMP) proposed by Bendsøe and Sigmund in \cite{Bendsoe1989Optimal,bendsoe2003Topology}.
The fundamental idea of SIMP is the introduction of a non-linear interpolation function between void and full material such that a black and white design is aspired due to the inherent non-convex total potential.
Further popular developments are overviewed in \cite{sigmund2013topology,Deaton2014ASurvey}.

In a series of papers, we aimed at contributing to the problem of topology optimization: by using thermodynamic extremal principles, evolutionary access to the problem of topology optimization has been presented for which we referred our method to as thermodynamic topology optimization (TTO).
It routes back to \cite{junker2015variational} while further important developments have been presented for the numerical treatment in \cite{jantos2019accurate} and for hyperelastic material in \cite{Junker2021hyperelastic}.
This topology optimization makes use of an extended Hamilton principle which is well-known in the context of material modeling, cf. \cite{junker2021extended}.
Therefore, the extended Hamilton functional is formulated and its stationarity conditions serve as update procedure for the evolution of the topology.
In this manner, no classical optimization problem is formulated.
Since the free energy function is part of the extended Hamilton functional, the result is very similar to classical schemes for topology optimization \cite{jantos2019comparison}.
The manipulation of topology is managed by the density as state variable which is defined for each discretized area.
One advantage of this thermodynamic topology optimization is that no further optimization algorithm is needed.
In contrast, the relative density is described by a transient partial differential equation (PDE) in which the local free energy density serves as source term.
Consequently, the material optimization problem is converted to an evolutionary problem.
The field equation for the topology results from the stationary condition of the extended Hamilton functional.
Additionally, evaluation of the functional results in the field equations for displacement and internal (state) variable which accounts for the (local) microstructure of the material.
From this follows that the extended Hamilton functional according to topology optimization also features to take any physically dissipative material behavior into account.

In context of \re{accounting for a} real material behavior during \re{the} optimization, the plastic material behavior plays a major role which requires a thermodynamically and mathematically rigorous treatment due to its complexity.
Direct access of plastic material behavior within topology optimization might be given by using classical plasticity models with the characteristic stress/strain diagram resulting in a hysteresis curve in cyclic loading.
Several successful examples are provided in the literature: a first approach to account for a classical elasto-plastic material model within an adaptive material topology optimization was proposed by Maute \textit{et al.}{} \cite{Maute1998Adaptive}.
Approaches to determine plasticity by homogenization strategies are also possible, cf. \cite{Yuge1995Optimization}.
This is particularly interesting for plastic parts of composites \cite{Swan1998VoigtReuss}.
Furthermore, topology optimization based on finite strains can be considered for plasticity \cite{Wallin2016Topology}.
Nakshatrala and Tortorelli \cite{Nakshatrala2015Topology} optimize dynamically loaded structures while accounting for plastic material behavior.
A different option was proposed by the consideration of damage, cf. \cite{Li2017Topology,Alberdi2017Topology}.
For all such non-linear stress/strain relations, the optimization routine usually demands an additional algorithm for convergence.
Here, one prominent possibility is provided by the ``method of moving asymptotes'' (MMA).

Unfortunately, the strategy of considering physical material models usually renders such optimization rather time-consuming: due to the local path-dependence, the physical loading process needs to be discretized with several time steps each of \re{which} demands the solution of the physical finite element problem.
The nested finite element simulations for the physical process and the topology optimization problem \re{demand} a remarkably higher amount of computation time.
To compensate this drawback, several strategies  can be found which aim to directly include plasticity into the optimization process: one approach is to couple an elastic model with local stress constraint as mentioned e.\,g.{} by \cite{Duysinx1998Topology,Duysinx1999Topology,Bruggi2012Topology,Luo2012Topology}.
Another idea by Amir \cite{Amir2016StressConstrained} is to define a single global stress constraint within the formulation of optimization problem to bypass the local calculation for each material point. Bruns \textit{et al.}{} \cite{bruns2002numerical} propose to constrain the load-displacement by limit points in case of damage.
Another way is to account for the nonlinear material behavior on a second, microscopic scale by developing a new model reduction technique which is proposed by Fritzen \textit{et al.}{} \cite{Fritzen2016Topology} and extended by Xia \textit{et al.}{} \cite{Xia2017Evolutionary}.
A special characteristic of this approach is the use of an evolutionary optimization method on the macroscopic scale.
Furthermore, surrogate models have been developed which avoid the need of solving physical finite element problems.
Zhao~et~al.~\cite{Zhao2019Material,Zhao2020Topology}, for instance, developed a surrogate model as a fictitious non-linear elastic material model which asymptotically approximates a perfect elasto-plastic behavior.
They accounted for the von Mises criterion in \cite{Zhao2019Material} and also developed an approach valid for the Drucker-Prager criterion in \cite{Zhao2020Topology}.
Due to the absent path-dependence, the computation of the sensitivity is straight forward and only one finite element analysis needs to be computed for each iteration step.
Therefore, this approach has \re{a remarkable} disadvantage that the resulting stress/strain curve matches the similar classical elasto-plastic curve even at a material point level only at the limit points.
Furthermore, there is no possibility to compute the plastic strain tensor which serve as thermodynamic state variable.

Hence, a topology optimization method accounting for plastic material behavior in a resource-efficient manner is missing.
In this contribution, we aim at expanding the thermodynamic topology optimization such that it can be applied to elasto-plastic materials with a novel 3D surrogate material model.
Therefore, the surrogate model is based on a classical plastic material model whereby the mentioned disadvantages of non-linearity and path-dependence should be prevented by modifications resulting from the \re{nature} of optimization: we determine topology optimization results for the (maximal) external loading; unload\re{ing} and cyclic load curves are not \re{investigated here}.
During an evolutionary optimization process the topology \re{and} the resulting structural stiffness evolve, respectively.
Different\re{ly} stiff structures loaded with the same external loading result in different displacement fields.
Therefore, the displacement field \re{also} evolves during the optimization process.
For instance, high strains might be present in the beginning of the optimization process with associated high plastic strains.
However, the evolution of local stiffness result\re{s} in reduced strains and consequently \re{reduced} plastic strains.
Since this \re{``}unloading\re{''} process does not correspond to the physical evolution of plastic strains but occurs due to the evolution of the topology optimization we denote for this as ``virtual unloading\re{''}.
\re{For a} classical \re{plasticity} model\re{, the virtual unloading} evoke\re{s} dissipation which results in the typical hysteresis curve.
However, we are interested in a material model that computes the plastic strains for each displacement state as it results from physical loading.
Therefore, the surrogate model \re{must reflect} the physically correct \re{stress/strain} curve \re{during loading, i.\,e.{} without any effects due to virtual unloading} -- the unloading process must not correspond to the physical evolution of plastic strains.
In the case of virtual unloading the loading branch in the stress/strain curve needs to be followed back\re{wards}.
To this end, we propose a hysteresis-free behavior for the surrogate model by suppressing the dissipative character.
Finally, this results in the following benefits \re{of} the surrogate model: i) the surrogate model match\re{es} the results of plastic material models at the material point level, ii) the physical state variable can be measured in terms of the plastic strain tensor and iii) the total number of finite element simulations can be reduced while \re{the material behavior is still} physically accurate.
\re{Consequently, our} surrogate model in terms of the thermodynamic topology optimization allows finding optimal structures if plastic material behavior occurs.

Th\re{e} paper is structured as follows: first, we recall the basics of the thermodynamic topology optimization by use of Hamilton's principle and complement our previous approaches by the inclusion of plasticity.
To this end, we develop a surrogate material model for our topology optimization approach that accounts for plasticity without consideration of dissipation-related hysteresis effects.
Afterwards we present a suitable strategy for numerical implementation.
\re{Finally}, the functionality of the proposed approach is tested and analyzed by means of computing topology optimizations for several boundary value problems.

\section{Surrogate model for plasticity developed in thermodynamic topology optimization}

The thermodynamic topology optimization is based on Hamilton's principle which is usually a variational strategy for material modeling \cite{junker2021extended}.
Assuming stationary of an extended Hamilton functional follows the thermodynamic and physical laws and yields field equations for all state variables i.\,e.{} displacements, temperature and internal variables.
Expanding the functional for topology optimization provides the benefit that the optimization problem can be tackled by a system of partial differential equations.
Consequently, the mathematical optimization problem is converted into an evolutionary problem.
At the same time, the stationarity of the extended Hamilton functional comprises the evolution of microstructural material behavior which affects the evolution of topology.
Furthermore, constraints on the topology design as well as on the material behavior can be considered easily by taking use of Lagrange or Karush Kuhn Tucker parameters.
It is worth mentioning that no classical optimization problem is solved in thermodynamic topology optimization.
In contrast, the stationarity condition of the Hamilton functional with respect to the density variable serves as update scheme for the topology.

We use the following notation for tensor operations:
the single contraction is noted as ``$\cdot$'' reading $\boldface{a} \cdot \boldface{b} = c \Leftrightarrow a_i b_i = c$ when applied to two vectors $\boldface{a}$ and $\boldface{b}$, while it results in $\boldface{A} \cdot \boldface{b} = \boldface{c} \Leftrightarrow A_{ij} b_j = c_i$ when applied to a vector and a second-order tensor $\boldface{A}$.
Moreover, the double contraction is denoted as ``$:$''.
It results in $\boldface{A} : \boldface{B} = c \Leftrightarrow A_{ij} B_{ij} = c$ when applied to two second-order tensors while it results in $\mathbb{A} : \boldface{B} = \boldface{C} \Leftrightarrow \re{\mathbb{A}}_{ijkl} B_{kl} = C_{ij}$ when applied to a fourth-order tensor $\mathbb{A}$ and a second-order tensor.
Finally, the tensor product, i.\,e.{} the dyadic product, is noted as ``$\otimes$'' and reads $\boldface{a} \otimes \boldface{b} = \boldface{C} \Leftrightarrow a_i b_j = C_{ij}$ when applied to two vectors and $\boldface{A} \otimes \boldface{B} = \mathbb{C} \Leftrightarrow A_{ij} B_{kl} = \re{\mathbb{C}}_{ijkl}$ when applied to two second-order tensors.

In this contribution, the approach of topology optimization does not account for dynamic effects and therefore, we consider quasi-static loading.
\re{Here}, the extended Hamilton functional \cite{junker2021extended} for a quasi-static and isothermal case reads
\begin{equation}
  \bar\mathcal{H} = \bar\mathcal{H} \left[\boldface{u},\boldface{y} \right] := \mathcal{G}\left[\boldface{u}, \boldface{y} \right] + \mathcal{D}\left[\boldsymbol{\alpha} \right]
\end{equation}
and sums the Gibbs energy $\mathcal{G}$ and the dissipation-related work $\mathcal{D}$.
This functional depends on the displacements $\boldface{u}$ and the state variable $\boldface{y} = \{\boldsymbol{\alpha}, \chi\}$.
The state variable is decomposed into the vectorial quantity $\boldsymbol{\alpha}$ collecting all internal variables which describe the physical material behavior in terms of the microstructural state.
In our case of an elasto-plastic material, we thus chose $\boldsymbol{\alpha}=\boldsymbol{\varepsilon}^{\mathrm{p}}$ where $\boldsymbol{\varepsilon}^{\mathrm{p}}$ denotes the plastic part of the strain and $\boldsymbol{\varepsilon}^\mathrm{e}$ the elastic part, i.\,e.{} $\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^\mathrm{e}+\boldsymbol{\varepsilon}^{\mathrm{p}}$.
The quantity $\chi$ denotes the density variable for defining the topology.
Here, the density variable $\chi \in [\chi_{\mathrm{min}},1]$ with $\chi_{\mathrm{min}} > 0$ represents void ``white'' material for $\chi = \chi_{\mathrm{min}}$, the full ``black'' material for $\chi = 1$, and a mixed ``gray'' phase for $\chi \in ]\chi_{\mathrm{min}},1[$.
The relative density is then \re{modeled} via the SIMP approach \cite{bendsoe2003Topology} by the interpolation function
\begin{equation}
\label{eq:chiInterpolation}
\rho(\chi)=\chi^3 \, ,
\end{equation}
for instance.
Other approaches are also possible, see~\cite{Junker2021hyperelastic} where a sigmoid function has been used.

According to Hamilton's principle the stationary condition of the functional is provided as
\begin{equation}
 \bar\mathcal{H} = \bar\mathcal{H}[\boldface{u}, \boldsymbol{\varepsilon}^{\mathrm{p}}, \chi] := \mathcal{G}\left[\boldface{u}, \boldsymbol{\varepsilon}^{\mathrm{p}}, \chi \right] + \mathcal{D}\left[\boldsymbol{\varepsilon}^{\mathrm{p}} \right] \rightarrow \underset{\boldface{u},\boldsymbol{\varepsilon}^{\mathrm{p}}, \chi}{\text{stat}} \ .
\end{equation}
Therein, $\mathcal{G}$ is defined as difference between the energy stored in the body with volume $\Omega$ and the work due to external forces.
It hence reads
\begin{equation}
  \label{eq:GibbsEnergy}
  \mathcal{G} \left[\boldface{u}, \boldsymbol{\varepsilon}^{\mathrm{p}}, \chi \right] := \int_{\Omega} \Psi(\boldsymbol{\varepsilon}^\mathrm{e}, \chi) \ \mathrm{d} V - \int_{\Omega} \boldface{b}^{\star} \cdot \boldface{u} \ \mathrm{d} V - \int_{\Gamma_{\sigma}} \boldface{t}^{\star} \cdot \boldface{u} \ \mathrm{d} A
\end{equation}
with the Helmholtz free energy $\Psi$, the body forces $\boldface{b}^{\star}$ and the traction vector $\boldface{t}^{\star}$.
The boundary conditions are defined as Dirichlet conditions for $\boldface{u}^{\star}$ on $\Gamma_{u}$ and as Neumann conditions for $\boldface{t}^{\star}$ on $\Gamma_{\sigma}$.
Hence, the complete boundary $\partial\Omega$ of the body is given by $\partial\Omega = \Gamma_{u} \cup \Gamma_{\sigma}$ and $\Gamma_{u} \cap \Gamma_{\sigma} = \varnothing$.
Furthermore, the dissipation-related work is defined by
\begin{equation}
  \mathcal{D} := \int_{\Omega} \bfp^\diss : \boldsymbol{\varepsilon}^{\mathrm{p}} \ \mathrm{d} V
\end{equation}
with the non-conservative force $\bfp^\diss$ which can be derived from the dissipation function $\Delta^{\mathrm{diss}}$ by
\begin{equation}
\label{eq:pdiss}
  \bfp^\diss := \pf{\Delta^{\mathrm{diss}}}{\dot{\boldsymbol{\varepsilon}}^\mathrm{p}} \ .
\end{equation}
More details on the thermodynamic basis are provided in~\cite{junker2021extended}.
According to \cite{Junker2021hyperelastic}, the physically motivated Hamilton functional $\bar{\mathcal{H}}$ can be extended for thermodynamic topology optimization by adding
\begin{equation}
  \label{eq:HamiltonFunctinalTTO}
  \mathcal{H}\left[\boldface{u}, \boldsymbol{\varepsilon}^{\mathrm{p}}, \chi \right] := \bar\mathcal{H} \left[\boldface{u}, \boldsymbol{\varepsilon}^{\mathrm{p}}, \chi \right] - \mathcal{R} \left[\chi \right] + \mathcal{C} \left[\boldsymbol{\varepsilon}^{\mathrm{p}}, \chi \right] \ .
\end{equation}
where additional constraints are included in $\mathcal{C}$ and the rearrangement of topology is accounted for by the functional $\mathcal{R}$, defined as
\begin{equation}
  \mathcal{R} \left[\chi \right] := \mathcal{D}_{\chi} \left[\chi \right] + \mathcal{F} \left[\chi \right] \ .
\end{equation}
Here, the flux term
\begin{equation}
  \mathcal{F} \left[\chi \right] := \int_{\Omega} \frac{1}{2} \, \beta \, \norm{\nabla \chi}^2 \ \mathrm{d} V
\end{equation}
accounts for the convective rearrangement with the regularization parameter $\beta \re{> 0}$.
It thus serves as gradient penalization for the density variable and also controls the members size via the parameter $\beta$.
Additionally, the source term
\begin{equation}
  \mathcal{D}_{\chi} \left[\chi \right] := \int_{\Omega} p^\mathrm{diss}_\chi \, \chi \ \mathrm{d} V
\end{equation}
accounts for local rearrangement.
Analogously to \eqref{eq:pdiss}, the non-conservative term for local rearrangement is assumed to be derivable from an associated dissipation function according to
\begin{equation}
p^\mathrm{diss}_\chi := \pf{\Delta^{\mathrm{diss}}_\chi}{\dot{\chi}} \ .
\end{equation}
For the dissipation function, we follow \cite{jantos2019accurate} and chose
\begin{equation}
\Delta^{\mathrm{diss}}_\chi = \frac{1}{2} \, \eta\, \dot{\chi}^2 \ .
\end{equation}
The viscosity parameter $\eta \re{> 0}$ controls the velocity of evolution of topology.
In this manner, the Hamilton functional \eqref{eq:HamiltonFunctinalTTO} is able to couple microstructure evolution and topology optimization.
We propose that an optimal structure can be found if this functional becomes stationary.

The stationary condition with respect to all variables
\begin{equation}
  \delta\mathcal{H} = \delta_{\boldface{u}} \mathcal{H} + \delta_{\boldsymbol{\varepsilon}^{\mathrm{p}}} \mathcal{H} + \delta_{\chi} \mathcal{H} = 0 \qquad \forall\; \delta_{\boldface{u}}, \delta_{\boldsymbol{\varepsilon}^{\mathrm{p}}}, \delta_{\chi}
\end{equation}
yields the following system of governing equations
\begin{equation}
  \label{eq:generalFieldEquations}
  \begin{cases}
    \displaystyle \delta_{\boldface{u}} \mathcal{H} = 0 = \int_{\Omega} \pf{\Psi}{\boldsymbol{\varepsilon}} : \delta\boldsymbol{\varepsilon} \ \mathrm{d} V - \int_{\Omega} \boldface{b}^{\star} \cdot \delta\boldface{u} \ \mathrm{d} V - \int_{\Gamma_{\sigma}} \boldface{t}^{\star} \cdot \delta\boldface{u} \ \mathrm{d} A & \forall\; \delta\boldface{u} \\[4mm]
    \displaystyle \delta_{\boldsymbol{\varepsilon}^{\mathrm{p}}} \mathcal{H} = 0 = \int_{\Omega} \pf{\Psi}{\boldsymbol{\varepsilon}^{\mathrm{p}}} : \delta\boldsymbol{\varepsilon}^{\mathrm{p}} \ \mathrm{d} V + \int_{\Omega} \pf{\Delta^{\mathrm{diss}}}{\dot{\boldsymbol{\varepsilon}}^\mathrm{p}} : \delta\boldsymbol{\varepsilon}^{\mathrm{p}} \ \mathrm{d} V + \delta_{\boldsymbol{\varepsilon}^{\mathrm{p}}} \, \mathcal{C} & \forall\; \delta\boldsymbol{\varepsilon}^{\mathrm{p}} \\[4mm]
    \displaystyle \delta_{\chi} \mathcal{H} = 0 = \int_{\Omega} \pf{\Psi}{\chi} \ \delta\chi \ \mathrm{d} V - \int_{\Omega} \eta \, \dot\chi \, \delta\chi \ \mathrm{d} V - \int_{\Omega} \beta \ \nabla \chi \cdot \nabla\delta\chi \ \mathrm{d} V + \delta_{\chi} \, \mathcal{C} & \forall\; \delta\chi
  \end{cases}
\end{equation}
where each equation belongs to one of the independent system variables, cf. also~\cite{Junker2021hyperelastic} for a general microstructure consideration in case of finite deformations.
Here, the standard notation $\delta\boldsymbol{\varepsilon}:=\mathrm{sym}(\nabla\otimes\delta\boldface{u})$ is used.
The first condition is identified as the weak form of the balance of the linear momentum where the stress is given by $\boldsymbol{\sigma} = \partial\Psi / \partial\boldsymbol{\varepsilon}$.
The second condition constitutes as governing equation for the plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}$ and the last equation is the field equation for topology optimization.

\subsection{Specification of the energetic quantities and \re{the} constraints}
The system of governing equations \eqref{eq:generalFieldEquations} establish\re{es} the \re{general} framework for the optimization process.
However, by specification of the free energy density $\Psi$, the dissipation function $\Delta^{\mathrm{diss}}$, and the constraint functional $\mathcal{C}$ the characteristics of the surrogate material model for plasticity and the density variable are defined.

For the free energy, we follow the classical approach of elasto-plastic materials and combine it with the relative density $\rho(\chi)$ in \eqref{eq:chiInterpolation}.
This gives
\begin{equation}
\Psi(\boldsymbol{\varepsilon}^\mathrm{e},\chi) = \rho(\chi) \, \Psi_0  := \frac{1}{2} \, (\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{\mathrm{p}}): \chi^3 \, \dsE_0 : (\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{\mathrm{p}})
\end{equation}
where the stiffness tensor of the full material is given by $\mathbb{E}_0$ and the energy of the virtually full material is given by
\begin{equation}
  \label{eq:fullHelmholtzEnergy}
  \Psi_0 := \frac{1}{2} \, (\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{\mathrm{p}}): \mathbb{E}_0 : (\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{\mathrm{p}}) \, .
\end{equation}
Consequently, we obtain for the stresses
\begin{equation}
\label{eq:stress}
\boldsymbol{\sigma} = \pf{\Psi}{\boldsymbol{\varepsilon}}= \chi^3 \, \dsE_0 : (\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{\mathrm{p}}) \ .
\end{equation}
The derivative of $\Psi$ with respect to $\boldsymbol{\varepsilon}^{\mathrm{p}}$ thus yields
\begin{equation}
\pf{\Psi}{\boldsymbol{\varepsilon}^{\mathrm{p}}} = - \boldsymbol{\sigma}
\end{equation}
and the derivative of $\Psi$ with respect to $\chi$ yields
\begin{equation}
\label{eq:DFTopOpt}
p:= -\pf{\Psi}{\chi} = - 3 \, \chi^2 \, \Psi_0
\end{equation}
as driving force or sensitivity for the topology optimization, respectively.
The driving force $p$ is non-zero for all conditions with $\Psi_0\not=0$ since $\chi>\chi_\text{min}$.
Furthermore, the evolution of plastic strains influences $\Psi_0$ and, in turn, the driving force $p$ and thus the update condition for optimization, cf.~\eqref{eq:generalFieldEquations}$_3$.

The following equations formulate the surrogate material model for the evolution of plastic strains in the context of thermodynamic topology optimization regarding three assumptions.
The first one concerns the dissipation function.
In a classical elasto-plastic material model, the dissipation function is defined as
\begin{equation}
\label{eq:dissfunc}
\Delta^{\mathrm{diss}} = r \, \norm{\dot{\boldsymbol{\varepsilon}}^\mathrm{p}}
\end{equation}
with the yield limit $r$.
This approach yields a rate-independent formulation.
Details on such an approach can be found, e.\,g.{}, in \cite{junker2013principle,junker2021extended}.
However, this physically motivated choice would contradict our intention to prevent the path-dependence and the related dissipative material behavior.
Therefore, we assume that the dissipation-free evolution of \re{plastic strains within} the surrogate model is acquired by postulating a vanishing dissipation function, i.\,e.{}
\begin{equation}
\Delta^{\mathrm{diss}}=0 \, .
\end{equation}

The second assumption is that the yield condition is included by demanding
\begin{equation}
\norm{\bfsigma^{\mathrm{dev}}} = r
\end{equation}
during plastic evolution where the stress deviator $\bfsigma^{\mathrm{dev}}= \boldsymbol{\sigma} - 1/3 \ \tr{\boldsymbol{\sigma}} \boldface{I}$ is computed by
\begin{equation}
  \bfsigma^{\mathrm{dev}} = \mathbb{P} : \boldsymbol{\sigma}
\end{equation}
with the projection tensor $\mathbb{P}$.
The threshold value $r$ will be defined phenomonologically and needs to be combined with the relative density $\rho\left(\chi\right)$ according to \cite{Duysinx1998Topology} for physical consistency.
Therefore, ideal plasticity is determined by a constant material parameter\re{,} e.\,g.{} the yield stresses $\sigma^\mathrm{Y}$\re{,} which yields
\begin{equation}
 \label{eq:r_ideal}
  r = \chi^3 \, \sigma^\mathrm{Y} \, .
\end{equation}
Hardening can be included by choosing a non-constant $r = r(\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}})$.
To this end, we propose linear hardening by defining
\begin{equation}
  \label{eq:r_linHard}
  r = r(\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}) = \chi^3 \left( \sigma^\mathrm{Y} + h \, \norm{\boldsymbol{\varepsilon}^{\mathrm{p}}} \right)
\end{equation}
with the slope of hardening curve $h$ and exponential hardening according to \cite{junker2017numerical} by
\begin{equation}
  \label{eq:r_expHard}
  r = r(\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}) = \chi^3 \left( \sigma^\mathrm{Y} + h_{1} \, \norm{\boldsymbol{\varepsilon}^{\mathrm{p}}} + \frac{1}{\re{\kappa}} \, \left( h_{1} - h_{0} \right) \left( e^{-\re{\kappa} \, \norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} - 1 \right) \right) \, \re{.}
\end{equation}
\re{Here,} $h_0$ denotes the initial and $h_1$ the end slope of the hardening curve and $\re{\kappa}$ controls the transition from $h_0$ to $h_1$.
Since our approach is equivalent for different definitions of $r$, we always use the general notation $r$ as yield criterion in the following.

The third assumption is that the plastic strains are volume-preserving, i.\,e.{}
\begin{equation}
  \label{eq:constraint_volumePreservation}
  \tr{\boldsymbol{\varepsilon}^{\mathrm{p}}} = 0 \, \Leftrightarrow \, \boldface{I} : \boldsymbol{\varepsilon}^{\mathrm{p}} = 0 \, .
\end{equation}
\re{This assumption combined with the above equations and definitions yield in} a non-linear algebraic equation for the determination of the plastic strains \re{in the following \secref{sec:stationarityConditionPlasticStrains}.}
Regarding \re{the volume-preservation} condition, \re{this} equation \re{would be} ill-posed due to a projection onto the deviator subspace.
Hence, no unique solution exists and a special numerical solution \re{would be} needed for the solution.
\re{However, we found to account for the volume-preservation of the plastic strains in a more efficient way, so that the resulting non-linear algebraic equation is well-posed: we make use of the hydrostatic stress by}
\re{
\begin{equation}
  \tr{\boldsymbol{\sigma}} = \boldface{I} : \boldsymbol{\sigma} = \boldface{I} : \chi^3 \, \dsE_0 : \left( \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{\mathrm{p}} \right)
\end{equation}
and apply the volume-preserving \eqref{eq:constraint_volumePreservation} so that with $\boldface{I} : \chi^3 \, \dsE_0 : \boldsymbol{\varepsilon}^{\mathrm{p}} = 0$ the constraint reads
\begin{equation}
  \boldface{I} : \boldsymbol{\sigma} = \boldface{I} : \chi^3 \, \dsE_0 : \boldsymbol{\varepsilon} \, .
\end{equation}
}

The limitation of \re{the} stress \re{norm by the yield threshold} and the volume preservation are included through the constraint functional by using the Lagrange parameters $\lambda_\mathrm{\sigma}$ and $\lambda_\mathrm{V}$, respectively.

It remains to identify the constraints for the density variable $\chi$ to finally formulate the constraint functional $\mathcal{C}$.
The first constraint is given by the interval in which $\chi$ is defined: values of $\chi$ that are negative are not reasonable; same limitation holds true for values of $\chi$ that are larger than one.
Consequently, we demand $\chi \in [\chi_{\mathrm{min}},1]$ where the lower bound is set to a small value $1 \gg \chi_{\mathrm{min}} > 0$ due to numerical reasons.
These bounds are taken into account by use of a Karush Kuhn Tucker parameter $\gamma$.
Furthermore, the volume of the topology relative to the total design volume is prescribed by the parameter $v_{0}$.
Consequently, it has to hold
\begin{equation}
\int_\Omega \chi \ \mathrm{d} V = v_{0} \, \Omega
\end{equation}
which is included to the constraint functional by use of a Lagrange parameter $\lambda_\mathrm{\chi}$.

Combining these four constraints, i.\,e.{} norm of the stress deviator being equivalent to the yield threshold $r$, volume preservation of the plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}$, bounded interval for $\chi$, and control of the total relative structure volume $v_{0}$, we finally arrive at
\begin{multline}
\label{eq:ConstraintsFunctional}
  \mathcal{C} := \lambda_\mathrm{\sigma} \int_{\Omega} \left(\norm{\bfsigma^{\mathrm{dev}}} - r\right) \ \mathrm{d} V   + \lambda_\mathrm{V} \int_{\Omega} \left( \re{\boldface{I} : \boldsymbol{\sigma} - \boldface{I} : \chi^3 \, \dsE_0 : \boldsymbol{\varepsilon}} \right) \ \mathrm{d} V \\
  + \int_{\Omega} \gamma \, \chi \ \mathrm{d} V + \lambda_\mathrm{\chi} \left(\int_{\Omega} \chi \ \mathrm{d} V - v_{0} \, \Omega \right)  \ .
\end{multline}

\subsection{The stationarity condition with respect to the plastic strains}
\label{sec:stationarityConditionPlasticStrains}
It remains to appropriately analyze the stationarity condition of the Hamilton functional with respect to the plastic strains.
This conditions enables us to compute the plastic strains which, in combination with the total strain, specify the stress state.
To this end, we use the specifications for a vanishing dissipation function $\Delta^{\mathrm{diss}}$ and the constraint functional \eqref{eq:ConstraintsFunctional} to evaluate \eqref{eq:generalFieldEquations}$_2$ as
\begin{eqnarray}
\label{eq:StatPlastLoc}
&& \int_\Omega \left( -\boldsymbol{\sigma} + \lambda_\mathrm{\sigma} \left[ \pf{\norm{\bfsigma^{\mathrm{dev}}}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} - \re{\pf{r}{\bfepsp}} \right] - \lambda_\mathrm{V} \,\re{ \boldface{I} : \chi^3 \, \dsE_0} \right) :\delta\boldsymbol{\varepsilon}^{\mathrm{p}} \ \mathrm{d} V = 0 \qquad \forall \; \delta \boldsymbol{\varepsilon}^{\mathrm{p}} \notag \\
&& \Rightarrow\quad - \boldsymbol{\sigma} + \lambda_\mathrm{\sigma} \left[ \pf{\norm{\bfsigma^{\mathrm{dev}}}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} - \re{\pf{r}{\bfepsp}} \right] - \lambda_\mathrm{V} \, \re{\boldface{I} : \chi^3 \, \dsE_0} \notag\\
&& = - \boldsymbol{\sigma} -\lambda_\mathrm{\sigma} \left[ \frac{\bfsigma^{\mathrm{dev}}}{\norm{\bfsigma^{\mathrm{dev}}}} :\mathbb{P}: \chi^3 \, \dsE_0 + \re{\pf{r}{\bfepsp}} \right] - \lambda_\mathrm{V} \, \re{\boldface{I} : \chi^3 \, \dsE_0} = 0 \, .
\end{eqnarray}
Solving \eqref{eq:StatPlastLoc} for the plastic strains constitutes our surrogate model for the plastic material behavior.
A detailed derivation of the Lagrange multipliers is deferred to App.~\ref{appendix:DerivationSurrogateModel}.
There, we show that the governing equation for the plastic strains is given as
\begin{eqnarray}
  \label{eq:s_epsP}
  \boldface{s} &:=& - \boldsymbol{\sigma} + \frac{r^3}{\displaystyle \bfsigma^{\mathrm{dev}} :  \chi^3 \, \dsE_0 : \bfsigma^{\mathrm{dev}} + \pf{r}{\boldsymbol{\varepsilon}^{\mathrm{p}}} : \bfsigma^{\mathrm{dev}} \, r} \, \left[\frac{\bfsigma^{\mathrm{dev}}:\chi^3 \, \dsE_0}{r} + \pf{r}{\boldsymbol{\varepsilon}^{\mathrm{p}}} \right] \\
  && + \frac{\boldface{I} : \chi^3 \, \dsE_0 : \boldsymbol{\varepsilon}}{\boldface{I} : \mathbb{E}_0 : \boldface{I}} \, \boldface{I} : \mathbb{E}_0 = \boldsymbol{0} \notag
\end{eqnarray}
which is a non-linear algebraic equation.
The derivative of the yield criterion $r$ is defined as
\begin{equation}
  \label{eq:dr}
  \displaystyle \pf{r}{\boldsymbol{\varepsilon}^{\mathrm{p}}} = \begin{cases}
  \displaystyle \boldsymbol{0} & \text{ideal plasticity} \\
  \displaystyle \pf{r \left(\norm{\bfepsp} \right)}{\boldsymbol{\varepsilon}^{\mathrm{p}}} = \pf{r \left(\norm{\bfepsp} \right)}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} \, \pf{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} = r' \, \frac{\boldsymbol{\varepsilon}^{\mathrm{p}}}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} & \text{hardening}
\end{cases}
\end{equation}
where the term $r'$ for the defined types of hardenings reads
\begin{equation}
  \displaystyle r' = \begin{cases}
  \displaystyle h & \text{linear hardening} \\
  \displaystyle  h_1 - \left( h_{1} - h_{0} \right) \, e^{-\re{\kappa} \, \norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} & \text{exponential hardening}
\end{cases} \, .
\end{equation}
In case of ideal plasticity with \re{$r = \text{constant}$} and the derivative from \eqref{eq:dr}\re{,} we can reduce \eqref{eq:s_epsP} to
\begin{equation}
  \label{eq:s_epsP_ideal}
  \boldface{s}^{\mathrm{\re{ip}}} := - \boldsymbol{\sigma} + \frac{r^2}{\bfsigma^{\mathrm{dev}} :  \mathbb{E}_0 : \bfsigma^{\mathrm{dev}}} \, \bfsigma^{\mathrm{dev}}:\mathbb{E}_0 + \frac{\boldface{I} : \chi^3 \, \dsE_0 : \boldsymbol{\varepsilon}}{\boldface{I} : \mathbb{E}_0 : \boldface{I}} \, \boldface{I} : \mathbb{E}_0 = \boldsymbol{0} \, .
\end{equation}

\textit{Remark:} it is worth mentioning that we do not receive a differential equation for the internal variable as it is usually the case.
This routes back to assuming a dissipation-free evolution of the plastic strains which, in turn, are determined by energy minimization.

Components of the plastic strain tensor only evolve to compensate high stresses which are greater than the yield stress $\sigma^\mathrm{Y}$.
Therefore, it is mandatory to identify a suitable criterion for distinguishing whether an elastic or plastic material behavior is present.
Since the purpose of the modified surrogate plasticity model is to display the same material behavior for loading like a classical material model for elasto-plasticity, we make use of the indicator function that would result from the dissipation function in \eqref{eq:dissfunc} via a Legendre transformation, cf. \cite{junker2013principle}.
This indicator function reads
\begin{equation}
\label{eq:YieldFunctionSigma}
\Phi_{\sigma} = \norm{\bfsigma^{\mathrm{dev}}} - r \le 0
\end{equation}
where elastic behavior is present for $\Phi_{\sigma} <0$ and plastic behavior for $\Phi_{\sigma} =0$.

Fitting the characteristic\re{s} of the classical elasto-plastic material model, physical unloading from a plastic state can be detected by this indicator function when the stress decreases once again below the yield threshold $r$.
The elastically stored energy is released first and the residual, plastic strains remains.
In this way, the hysteresis loop in the stress/strain diagram of a physical material evolves.

This behavior must be suppressed by the surrogate material model as discussed above.
Virtual unloading from a plastic state should \re{immediately} result in a decrease of plastic strains.
Thus, the plastic strains are reduced first and only if no plastic strains are present anymore, the elastically stored energy is released.
In this way, the loading branch in the stress/strain curve is followed both for loading and virtual unloading.

Consequently, the stress is not a suitable measure for the indicator function related to the surrogate model.
\re{Hence}, the strains are \re{identified as suitable} measure.
\re{We therefore,} reformulate the indicator function \eqref{eq:YieldFunctionSigma} in terms of strains.
\re{To this end, the yield threshold $r$ is compared to the linear stress $\boldsymbol{\sigma}^{\star} = \chi^3 \, \dsE_0 : \boldsymbol{\varepsilon}$ which occurs depending on the total strain $\boldsymbol{\varepsilon}$.}
Therefore, we can present the yield function as
\begin{equation}
  \label{eq:YieldFunctionPlasticStrains}
  \Phi_{\varepsilon} = \norm{\mathbb{P} : \boldsymbol{\sigma}^{\star}} - r \qquad \text{where } \qquad
  \Phi = \begin{cases}
    < 0    & \text{elastic} \\
    \geq 0 & \text{plastic}
  \end{cases} \, .
\end{equation}

\subsection{The stationarity condition with respect to the density variable}
Finally, the evolution of the density variable needs to be formulated.
Therefore, it remains to investigate the governing equation for the density variable $\chi$ which is given by \eqref{eq:generalFieldEquations}$_3$.
Making use of the constraint functional $\mathcal{C}$ in \eqref{eq:ConstraintsFunctional} and the driving force for topology optimization $p$ in~\eqref{eq:DFTopOpt}, the stationarity with respect to $\chi$ takes the form
\begin{equation}
\label{eq:TopOptWeak}
\int_\Omega \left( -p - \eta \, \dot{\chi} + \gamma + \lambda_\mathrm{\chi} \right) \delta\chi \ \mathrm{d} V - \int_\Omega \beta \, \nabla\chi \cdot\nabla\delta\chi \ \mathrm{d} V = 0 \qquad \forall \; \delta\chi
\end{equation}
which is a parabolic differential equation and shows some similarities to phase field equations, cf.~\cite{bartels2021cahn} for instance.
Analogously to the stationarity with respect to the displacements in \eqref{eq:generalFieldEquations}$_1$, this equation \eqref{eq:TopOptWeak} is the weak form of the associated Euler equation (which is the balance of linear momentum for the displacements).
Therefore, one possibility for numerical evaluation would be given by direct application of the finite element method.
A comparable approach has been presented in \cite{junker2015variational}.
However, it has turned out that this procedure is much more time consuming than applying the numerical method that has been presented in \cite{jantos2019accurate} due to the complex constraints of the bounded interval for $\chi$ and the prescribed total density $v_{0}$.
Therefore, in order to apply the method of the previous work in \cite{jantos2019accurate} which reduces the numerical efforts by approximately one order of magnitude, we transform \eqref{eq:TopOptWeak} to its strong form by integration by parts.
This results in
\begin{equation}
\label{eq:TopOptStrong}
  \begin{cases}
    \displaystyle \eta \, \dot\chi \, \in - p + \beta \, \mathop{}\!\mathbin\bigtriangleup\chi + \lambda_\mathrm{\chi} + \gamma & \forall\; \boldface{x} \in \Omega \\[4mm]
    \displaystyle \boldface{n} \cdot \nabla\chi = 0 & \forall\; \boldface{x} \in \partial\Omega
  \end{cases}
\end{equation}
where \eqref{eq:TopOptStrong}$_2$ is the Neumann boundary condition for the density variable.
It ensures conservation of the prescribed total structure volume.
Meanwhile, the change of the density variable is defined by \eqref{eq:TopOptStrong}$_1$ and accounts for the Laplace operator which is defined as
\begin{equation}
\mathop{}\!\mathbin\bigtriangleup\chi:= \frac{\partial^2 \chi}{\partial x^2} + \frac{\partial^2 \chi}{\partial y^2} + \frac{\partial^2 \chi}{\partial z^2} \, .
\end{equation}
The transient characteristic of this term require the specification of an initial value for $\chi(\boldface{x},t=0)=\chi_\mathrm{ini}\,\forall \;\boldface{x}\in\Omega$, which will be introduced with the numerical treatment in \secref{sec:updateDensity}.

\FloatBarrier
\section{Numerical implementation}
\label{sec:numerics}
In summary, the following system of coupled differential-algebraic equations needs to be solved:
\begin{equation}
\label{eq:numericalFieldEquations}
\begin{cases}
  \displaystyle 0 = \int_\Omega \boldsymbol{\sigma} \cdot \delta\boldsymbol{\varepsilon} \ \mathrm{d} V - \int_\Omega \boldface{b}^* \cdot \delta\boldface{u} \ \mathrm{d} V - \int_{\partial\Omega} \boldface{t}^* \cdot \delta\boldface{u} \ \mathrm{d} A & \forall \; \delta\boldface{u} \\[4mm]
  \displaystyle \boldface{0}  = \boldface{s} \text{, see \eqref{eq:s_epsP}} & \forall \; \boldface{x} \in \Omega\\[4mm]
  \displaystyle \dot{\chi} \, \in \, \frac{1}{\eta} \left( - p + \beta \, \mathop{}\!\mathbin\bigtriangleup\chi + \lambda_\mathrm{\chi} + \gamma \right) & \forall \; \boldface{x} \in \Omega
\end{cases}
\end{equation}
\re{The numerical implementation based on this solution is written in Julia programming language \cite{Julia} and published as open-access file in \cite{Kick2022Implementation}.}
It is worth mentioning that we use for now on the usual Voigt notation for the stresses and strains which reduces, for instance, the double contraction to a scalar product in \eqref{eq:numericalFieldEquations}$_1$ and \eqref{eq:numericalFieldEquations}$_2$.
\re{}

The numerical solution of the system of equations of the displacement field $\boldface{u}$, the microstructural plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}$ and the topology density $\chi$ is a sophisticated task due to the inherent non-linearities, constraints, and strong coupling.
Therefore, instead of applying a monolithic update scheme, cf.~\cite{junker2015variational}, we solve the equations in a staggered manner.
This can be interpreted as operator split which has turned beneficial in our previous works as in~\cite{jantos2019accurate} and also for adaptive finite element usage in~\cite{vogel2020adaptive}.
Here, both the finite element method (FEM) and the finite difference method (FDM) are \re{employed} for the solution.
This combination in the staggered process is referred to as neighbored element method (NEM)\re{,} cf.~\cite{jantos2019accurate}.
According to the staggered process, our method \re{shows similarities to} conventional mathematical optimization methods which are composed of alternating structure computation and optimization algorithm.

During the iterative solution of \eqref{eq:numericalFieldEquations}\re{,} each iteration step corresponds to an update step of the thermodynamic topology optimization.
In this way, an evolutionary update \re{of,} e.\,g.{}\re{,} the density field takes place.
For this purpose, we employ a standard discretization in pseudo-time, given as
\begin{equation}
\Delta t := t_{n+1} - t_n
\end{equation}
where $t_{n+1}$ refers to the current iteration step and $t_n$ to the previous iteration step.

\subsection{\re{U}pdate of the displacements}
Due to the operator split, a standard finite element approach is employed for updating the displacements and the stress in~\eqref{eq:numericalFieldEquations}$_1$ is evaluated as
\begin{equation}
\label{eq:numStress}
\boldsymbol{\sigma}_{n+1} = \chi^3_n \, \mathbb{E}_0 \cdot(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}^{\mathrm{p}}_n)
\end{equation}
so that this current stress is based on the plastic strains of the previous iteration $n$.
Thus, the stress and the resulting displacement field evolve through the optimization process.
To this end, the displacement field is approximated using the Galerkin-Ansatz
\begin{equation}
  \label{eq:GalerkinAnsatz}
  u_{k} = N_{o} \, u_{o}^\comp{k} = \boldface{N} \cdot \hat\boldface{u}^\comp{k}
\end{equation}
with the shape function $\boldface{N}$ and the nodal displacement $\hat\boldface{u}^\comp{k}$ in the spatial direction $k$.
Therefore, the weak form of the balance of linear momentum in~\eqref{eq:numericalFieldEquations}$_1$ transforms to
\begin{equation}
  \label{eq:updateDisplacements}
  \int_\Omega \boldsymbol{\sigma} \cdot \delta\boldsymbol{\varepsilon} \ \mathrm{d} V = \delta\hat\boldface{u} \cdot \int_\Omega \boldface{B}^\mathrm{T} \cdot \boldsymbol{\sigma} \ \mathrm{d} V = 0 =: \delta\hat\boldface{u} \cdot \boldface{r} \qquad \forall\; \delta\hat\boldface{u}
\end{equation}
when body forces are neglected.
Here, $\boldface{B}$ denotes the usual operator matrix including the spatial derivatives of the shape function.
The quantity $\delta\hat\boldface{u}$ is the global column matrix of nodal virtual displacements which also includes the Dirichlet boundary conditions.
Finally, the global residual column matrix is denoted by $\boldface{r}$ and, accordingly, the nodal displacements will be found from $\boldface{r} = \boldsymbol{0}$.
The global residual $\boldface{r}$ is assembled in usual manner by
\begin{equation}
  \boldface{r} := \AssemargPJ{e}{} \boldface{r}_e \quad\text{where}\quad \boldface{r}_e := \int_{\Omega_e} \boldface{B}^\mathrm{T}_e \cdot \boldsymbol{\sigma} \ \mathrm{d} V
\end{equation}
denotes the residual column matrix for each element $e$.
More details on the finite element method can be found in standard textbooks, e.\,g.{}, \cite{wriggers2008nonlinear}.

\re{Since our numerical implementation (cf. \cite{Kick2022Implementation}) of t}he thermodynamic topology optimization including plasticity has been coded in Julia \cite{Julia}\re{,} we made use of the finite element toolbox Ferrite~\cite{Ferrite}.
Ferrite uses a gradient-based equation solver as it is the standard for many finite element programs.
Consequently, the iterative solution  process for  $\boldface{r}=\boldsymbol{0}$ is performed by
\begin{equation}
  \boldface{r}^{i+1} = \boldface{r}^{i} + \pf{\boldface{r}}{\hat\boldface{u}} \cdot \Delta\hat\boldface{u} = \boldsymbol{0}
\end{equation}
where the iteration number is given by $i$.
The increment $\Delta\hat\boldface{u}$ updates the displacement field iteratively for fixated plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}_n$ and density field $\chi_n$.
The required element tangent is computed as
\begin{equation}
\pf{\boldface{r}_{e}}{\hat\boldface{u}_{e}} = \int_{\Omega_{e}} \boldface{B}_{e}^\mathrm{T} \cdot \chi_n^3 \, \mathbb{E}_0 \cdot \boldface{B}_{e} \ \mathrm{d} V
\end{equation}
with the column matrix of displacements for each finite element $e$ denoted as $\hat\boldface{u}_{e}$.
Then, the assembled tangent is constructed by
\begin{equation}
\pf{\boldface{r}}{\hat\boldface{u}} = \AssemargPJ{e}{} \pf{\boldface{r}_{e}}{\hat\boldface{u}_{e}} \ .
\end{equation}

\textit{Remark:} It is worth mentioning that we used the tensors package \cite{JuliaTensors} of Julia in our numerical implementation which is optimized for using tensors of higher order.
Therefore, we did not perform a finite element programming in standard form, i.\,e.{}, by using the Voigt notation, but used the full tensor notation.
This, of course, also effects the dimensions of the other quantities, i.\,e.{}, the $\boldface{B}$ operator is an array with three indices.
For a more usual presentation, we presented the formulas by using the Voigt notation and deferred our array-based programming using the tensors package to \appref{appendix:FerriteFEM}.

\subsection{\re{U}pdate of the plastic strains}
\label{ssec:UpdateEpsP}
The plastic strains are defined, as usual, for each integration point.
According to the discretization we employ for the density variable, all integration points in the same finite element are evaluated with the same value for the density variable $\chi$.
More details are given in \secref{sec:updateDensity} when we discuss the numerical treatment for the density variable.

The plastic strains are determined from solving~\eqref{eq:numericalFieldEquations}$_2$ which is a non-linear algebraic equation.
Within the update scheme of the plastic strains\re{,} we employ the operator split with $\boldsymbol{\sigma}={\boldsymbol{\sigma}}(\boldsymbol{\varepsilon}_{n+1},\boldsymbol{\varepsilon}^{\mathrm{p}}_{n+1},\chi_n)$ accounting for the element-wise density from the last iteration $n$ and the updated value of the plastic strains.
For the numerical implementation we make use of Newton's method to find the roots \re{of $\boldface{s}$} and define the Newton iterator $i$.
The Newton method for~\eqref{eq:s_epsP} reads
\begin{equation}
  \boldface{s} + \pf{\boldface{s}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} \cdot \Delta \boldsymbol{\varepsilon}^{\mathrm{p}} = \boldsymbol{0}
\end{equation}
and the plastic strains are iteratively updated according to
\begin{equation}
  \label{eq:NewtonUpdateESepsp}
  \boldsymbol{\varepsilon}^{\mathrm{p}}_{n+1} \leftarrow \boldsymbol{\varepsilon}^{\mathrm{p}}_{i+1} = \boldsymbol{\varepsilon}^{\mathrm{p}}_{i} - \left[ \pf{\boldface{s}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} \right]^{-1} \cdot \boldface{s} \ .
\end{equation}
The analytical tangent reads
\begin{eqnarray}
  \label{eq:ds_epsP}
  \pf{\boldface{s}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} &=& - \chi^3 \, \dsE_0 + \frac{1}{\left(\displaystyle \bfsigma^{\mathrm{dev}} \cdot \chi^3 \, \dsE_0 \cdot \bfsigma^{\mathrm{dev}} + \pf{r}{\bfepsp} \cdot \bfsigma^{\mathrm{dev}} \, r \right) ^2}  \\
  && \Bigg[ 3 \, r^2 \, \pf{r}{\bfepsp} \left( \bfsigma^{\mathrm{dev}} \cdot \chi^3 \, \dsE_0 \cdot \bfsigma^{\mathrm{dev}} + \pf{r}{\bfepsp} \cdot \bfsigma^{\mathrm{dev}} \, r \right) \notag \\
  && - r^3 \bigg( - 2 \, \chi^3 \, \dsE_0 : \chi^3 \, \dsE_0 : \bfsigma^{\mathrm{dev}} + \re{\frac{\partial^2 r}{\partial \bfepsp \partial \bfepsp}} : \bfsigma^{\mathrm{dev}} \, r \notag \\
  && - \pf{r}{\bfepsp} : \chi^3 \, \dsE_0 \, r + \pf{r}{\bfepsp} : \bfsigma^{\mathrm{dev}} : \pf{r}{\bfepsp} \bigg) \Bigg] \otimes \left[ \frac{1}{r} \, \bfsigma^{\mathrm{dev}} \cdot \chi^3 \, \dsE_0 + \pf{r}{\bfepsp} \right] \notag \\
  && - \frac{r}{\displaystyle \bfsigma^{\mathrm{dev}} \cdot \chi^3 \, \dsE_0 \cdot \bfsigma^{\mathrm{dev}} + \pf{r}{\bfepsp} \cdot \bfsigma^{\mathrm{dev}} \, r} \, \left[ \mathbb{P} : \chi^3 \, \dsE_0 : \chi^3 \, \dsE_0 \, r + \bfsigma^{\mathrm{dev}} : \chi^3 \, \dsE_0 \otimes \pf{r}{\bfepsp} \right] \notag
\end{eqnarray}
where the yield criterions $r$ was defined in \eqref{eq:r_linHard} and \eqref{eq:r_expHard} as well as its first derivatives in \eqref{eq:dr}.
The second derivative of the yield criterion $r$ reads
\begin{equation}
  \re{\frac{\partial^2 r}{\partial \bfepsp \partial \bfepsp}} = \begin{cases}
    \displaystyle \boldsymbol{0} & \text{ideal plasticity} \\
    \begin{split}
      \displaystyle \ppf{r \left(\norm{\bfepsp} \right)}{\re{\boldsymbol{\varepsilon}^{\mathrm{p}}}} &= r' \, \ppf{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}}{\re{\boldsymbol{\varepsilon}^{\mathrm{p}}}} + \pf{r'\left(\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}\right)}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} \, \pf{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} \otimes \frac{\boldsymbol{\varepsilon}^{\mathrm{p}}}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} \\
      &= r' \, \ppf{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}}{\re{\boldsymbol{\varepsilon}^{\mathrm{p}}}} + r'' \, \frac{\boldsymbol{\varepsilon}^{\mathrm{p}} \otimes \boldsymbol{\varepsilon}^{\mathrm{p}}}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}^2}
    \end{split} & \text{hardening}
  \end{cases}
\end{equation}
where we make use of
\begin{equation}
  \ppf{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}}{\re{\boldsymbol{\varepsilon}^{\mathrm{p}}}} = \frac{\mathbb{I}}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} - \frac{\boldsymbol{\varepsilon}^{\mathrm{p}} \otimes \boldsymbol{\varepsilon}^{\mathrm{p}}}{\norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}^3} \, .
\end{equation}
Furthermore, $r''$ is defined in terms of the type of the hardening as
\begin{equation}
  \displaystyle r'' = \begin{cases}
  \displaystyle 0 & \text{linear hardening} \\
  \displaystyle  \re{\kappa} \, \left( h_{1} - h_{0} \right) \, e^{-\re{\kappa} \, \norm{\boldsymbol{\varepsilon}^{\mathrm{p}}}} & \text{exponential hardening}
\end{cases} \, .
\end{equation}
The initial value for the plastic strains is chosen as $\boldsymbol{\varepsilon}^{\mathrm{p}}_{\mathrm{ini}}=\boldsymbol{0}$ \re{at the beginning of each iteration step}.
The convergence is defined such that all components of $\boldface{s}$ must be numerically zero, $\re{\max\{\boldface{s}\}} \leq 10^{-8}$ for instance.

It turns out, that the components of $\boldface{s}$ are small for each integration point located at every element with a small density variable $\chi$.
For this reason, the value of plastic strains computed by the described method are not as accurate as for larger density values.
Therefore, we propose to factorize equation \eqref{eq:s_epsP} with $1 / \chi_{n}$ so that it reads
\begin{equation}
  \label{eq:s_epsP_scaled}
  \tilde\boldface{s} := \frac{1}{\chi_{n}} \, \boldface{s}  = \boldsymbol{0}
\end{equation}
and its tangent \eqref{eq:ds_epsP} can be denoted as
\begin{equation}
  \label{eq:ds_epsP_scaled}
  \pf{\tilde\boldface{s}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} = \frac{1}{\chi_{n}} \, \pf{\boldface{s}}{\boldsymbol{\varepsilon}^{\mathrm{p}}} \, .
\end{equation}
The roots of any equation remain the same during factorization so that the scaling is only a numerical technique which has no influence on the magnitude of the resulting value but on the precision.
An overview of this numerical update algorithm is given in \algoref{alg:updatePlasticStrain}.
\begin{algorithm}[!htb]
  \caption{Compute update of plastic strains by Newton method}
  \label{alg:updatePlasticStrain}
  \input{./03_figures/alg_update_epsP.tex}
\end{algorithm}

However, to numerically stabilize convergence, it is purposeful to compute only plastic strains for stresses that \re{differ significantly} from the current yield criterion $r$.
Therefore, we propose to keep constant plastic strains within the plastic case if the stress is \re{close} the current yield criterion~$r$.
The criterion \re{for update plastic} is the trial stress defined as
\begin{equation}
  \boldsymbol{\sigma}_{\mathrm{trial}} := \chi^3 \, \dsE_0 \cdot \left(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}^{\mathrm{p}}_{n} \right)
\end{equation}
\re{with} the plastic strains from the last iteration $n$.
To this end, the classic indicator function defined in \eqref{eq:YieldFunctionSigma} depending on this trial stress is evaluated.
If the \re{relative} value $\Phi_\sigma(\boldsymbol{\sigma}_{\mathrm{trial}}) / r$ is less than \SI{1}{\percent} then the current plastic strains are set equal to the plastic strains from the last iteration: $\boldsymbol{\varepsilon}^{\mathrm{p}}_{n+1} \leftarrow \boldsymbol{\varepsilon}^{\mathrm{p}}_{n}$.
\re{Otherwise, $\boldface{s} = \boldsymbol{0}$ is solved for the updated values of $\boldsymbol{\varepsilon}^{\mathrm{p}}_{n+1}$.}

\re{In summary, the numerical implementation of the complete update scheme with all cases can be viewed as Julia code in \cite{Kick2022Implementation}.}

\subsection{\re{U}pdate of the density variable}
\label{sec:updateDensity}
Each value of the density field is evaluated for one finite element $e$ as discrete subvolume.
The evolution of the density variable is described by the transient partial differential equation in~\eqref{eq:numericalFieldEquations}$_3$ which needs to be discretized both in time and space for numerical evaluation.
Various strategies can be used for this purpose, e.\,g.{}, a finite element approach would be possible.
However, due to constraint of bounded interval for density $\chi$ and prescribed design volume $v_{0}$, a direct FE approach consumes a remarkable amount of computation time, cf.~\cite{junker2015variational}, where such a procedure has been discussed.
A more advantageous numerical treatment for this equation has therefore been presented in~\cite{jantos2019accurate} which is based on a generalized FDM along with an operator split.
More details on the numerical performance of this method, also regarding important aspects like convergence behavior and robustness, have been investigated in~\cite{vogel2020adaptive}.
In this paper, we make use of the published finite difference strategies and therefore only recall the fundamental update \re{strategy} and refer to the original publications \re{as well as our code (cf. \cite{Kick2022Implementation})} for a detailed explanation.

The transient character of the evolution equation demands the definition of the initial value for the density variable for each element.
As naive guess, we set each discretized density variable to $\chi_\mathrm{ini} = v_{0}$.
Therefore, the constraint of the given prescribed structure volume is identically fulfilled.

The change of density is driven by the driving force $p$ in equation \eqref{eq:TopOptStrong}.
Considering the operator split, the driving force $p$ is based on the Helmholtz free energy $\Psi_{0,n+1} := \Psi_{0,n+1} \left( \boldsymbol{\varepsilon}_{n+1}, \boldsymbol{\varepsilon}^{\mathrm{p}}_{n}, \chi_{n} \right)$.
High values of the driving force $p$ result in increasing densities and low values result in decreasing densities, respectively.
Since the actual value of the driving force is of no significance it is thus suitable to normalize the driving force with the weighted driving force (cf. equation (36) in \cite{jantos2019accurate}) by
\begin{equation}
   p_{\mathrm{w}} := \frac{\displaystyle \sum_e \left(\chi_{e} - \chi_{\mathrm{min}}\right) \, \left(1 - \chi_{e} \right) \, p_{e}}{\displaystyle \sum_e \left(\chi_{e} - \chi_{\mathrm{min}}\right) \, \left(1 - \chi_{e} \right)}
\end{equation}
to define the dimensionless driving force $\bar p := p/p_{\mathrm{w}}$.

Subsequently, the update scheme is employed according to \cite{jantos2019accurate}.
Then, the discretized evolution equation for the density variable for each element is given by
\begin{equation}
  \label{eq:updateDensity}
  \chi_{n+1} = \chi_{n} + \Delta t \, \frac{1}{\eta} \, \left[ - \bar p_{n+1} + \beta \, \mathop{}\!\mathbin\bigtriangleup\chi_{n} + \lambda_{\mathrm{\chi}} + \gamma \right]
\end{equation}
analogously to equation (49) \cite{Junker2021hyperelastic}.
Due to this, we are able to account for the regularization parameter $\beta$ in length unit squared and the viscosity $\eta$ in time unit as general optimization parameters.

To determine the value of the Lagrange parameter $\lambda_{\mathrm{\chi}}$ for the volume constraint, the update equation \eqref{eq:updateDensity} is solved iteratively by a simple bisection algorithm analogously to Alg.~1 in \cite{jantos2019accurate}.
This process also determines $\gamma$.
\re{Both are implemented in \cite{Kick2022Implementation} with the density update scheme.}
\subsection{Optimization process}
The presented update schemes take place in a global optimization process.
As proposed\re{,} we denote this staggered process of FEM and FDM as NEM, cf. \cite{jantos2019accurate}:
first the update of the displacements $\boldface{u}_{n+1}$ is solved by the finite element method for fixed values of the plastic strains at the previous global iteration step $\boldsymbol{\varepsilon}^{\mathrm{p}}_{n}$ and fixed values of the density variable $\chi_{n}$.
After updating the displacements, both the update of the plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}_{n+1}$ and the update of the density variable $\chi_{n+1}$ are performed using the updated displacements $\boldface{u}_{n+1}$.
The updated value for the plastic strains and the density variable are used for updating the displacements in the succeeding global iteration step $n \leftarrow n+1$.
In each iteration\re{,} the current stress and therefore the displacement field lag behind the physical reality because of the staggered process.
This could be counteracted by \re{further updates of the displacement field $\boldface{u}_{n+1}$ and the plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}_{n+1}$ by FEM}
for several times \re{(loops)} before going on with the next update of the topology.
The flowchart of the thermodynamic topology optimization including plasticity is given in \figref{fig:flowchart_TTO}.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/flowchart_TTO.pdf}
  \caption{Flowchart of the proposed numerical implementation of the thermodynamic topology optimization including plasticity.}
  \label{fig:flowchart_TTO}
\end{figure}

\FloatBarrier
\section{Numerical results}
We present several aspects of our novel thermodynamic topology optimization including plasticity by investigation of various numerical experiments.
We begin with the presentation of the general functionality of the proposed surrogate material model for plasticity on the material point level.
Afterwards, we show the impact of the material model on the optimized construction parts by means of analyzing several quasi\re{-}2D and 3D boundary value problems.
All results are based on \re{our} numerical implementation \re{\cite{Kick2022Implementation}} \re{in} Julia \cite{Julia}.
We use the material parameter for steel summarized in \tabref{tab:materialParameters}.
\begin{table}[!htb]
  \centering
  \caption{Material parameters for steel.}
  \label{tab:materialParameters}
  \begin{tabular}{L{3.1cm}C{0.6cm}R{1.3cm}L{1.0cm}C{0.03mm}|C{0.03mm}L{4.5cm}C{0.7cm}R{1.1cm}L{1.0cm}}
    Young's modulus & $E_0$ &  \num{210000} & [\si{\mega\pascal}]
    &&& slope of hardening & $h$/$h_1$ & \num{1000} & [\si{\mega\pascal}] \\
    Poisson's ratio & $\nu$ & \num{0.3} & [-]
    &&& initial slope of hardening & $h_0$ & \num{40000} & [\si{\mega\pascal}] \\
    yield stress & $\sigmaY_\mathrm{exp}$ & \num{300} & [\si{\mega\pascal}]
    &&& transition of hardening & $\re{\kappa}$ & \num{300} & [-]\\
  \end{tabular}
\end{table}
\re{The yield stress for modelling results from the material parameter with $\sigma^\mathrm{Y} = \sqrt{2 / 3} \, \sigmaY_\mathrm{exp}$.}
The hardening parameters are chosen according to \cite{junker2017numerical}.
An overview of the different material models used in the following is given in \figref{fig:plot_models_mm_types} on material point level.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/plot_models_mm_types}
  \caption{\re{Overview of material behavior: elastic and different types of plastic models.}}
  \label{fig:plot_models_mm_types}
\end{figure}

\subsection{Surrogate model for plasticity}
The computation of plastic strains takes place at the microstructural level.
To investigate the results of the proposed surrogate material model for plasticity, we present a first result at the material point and thus without topology optimization.
Consequently, we prescribe the strain as linear function of \num{100} load steps with tension and pressure loading and unloading.
For this, we determine the strain tensor depending on the load step $l$ according to
\begin{equation}
\label{eq:PrescribedStrain}
\boldsymbol{\varepsilon}(l) = \varepsilon_{11}(l) \, \left(
\begin{array}{rrr}
  1   & 0.6   & 0.6   \\
  0.6 & -\nu  & -0.1  \\
  0.6 & -0.1  & -\nu  \\
\end{array} \right) \, .
\end{equation}
To present a result that is representative, the diagonal entries correspond to the material parameters given above (\tabref{tab:materialParameters}), i.\,e.{}, we use the Poisson's ratio of steel, and the shear components have been chosen randomly.
The maximum value of the component in $11$-direction is set to $\varepsilon_{11,\max}(l) = 0.005 \, [-]$.
The numerical results for the surrogate model for plasticity at the material point are given as stress/strain diagram exemplary for ideal plasticity.

Matching the scalar-valued comparison of the indicator function, the von Mises stresses are plotted above the norm of strains in \figref{fig:plot_spmodel_mm_yielding}.
It indicates that the intended material behavior is displayed: first, the stress/strain curve of the proposed material model increases linearly in the elastic region.
The end points of the elastic region are indicated by $\varepsilon_\mathrm{exp}^\mathrm{Y}$ and $\sigmaY_\mathrm{exp}$, respectively.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/plot_spmodel_mm_yielding.pdf}
  \caption{Overview of the decision criteria for microstructural update cases using the example of ideal plasticity. The update of the plastic strains is executed if the plastic case is detected.}
  \label{fig:plot_spmodel_mm_yielding}
\end{figure}%
Then, the stress reaches the yield stress level $r$, here $\sigmaY_\mathrm{exp}$, in the plastic case.
This behavior coincides to classical plasticity models.
However, the remarkable difference is that the unloading case is also included in \figref{fig:plot_spmodel_mm_yielding}.
Here, no hysteresis is observed but with decreasing strains the stress level is maintained until the strains indicate the elastic region.
The result is thus independent of the unloading history.
Correspondingly, the increase or decrease of plastic strains in the surrogate material model directly reacts on the increase or decrease of strains in the plastic case.

An important difference of our novel surrogate model for dissipation-free plasticity to classical elasto-plastic material models is that we do not formulate our model by using an ordinary differential equation.
Consequently, path-dependence, as intended, is excluded in our model.
Of course, there exists no proof that the different formulations, ODE for classical models vs. algebraic equation for our model, give same results even when only the loading case is considered for which we demand a similar material behavior.
To investigate the quality of our novel surrgate model in this regard, we compare the surrogate material model and the hysteresis curve for a classical elasto-plastic model accounting for \re{one} component of the stress/strain state.
Thus, both curves are shown in \figref{fig:plot_models_mm_compare}.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/plot_models_mm_compare.pdf}
  \caption{Material point curve of the surrogate model and the classic elasto-plastic model for ideal plasticity. The new idea about the surrogate model is the ``virtual'' compared to the physical unloading.}
  \label{fig:plot_models_mm_compare}
\end{figure}%
Here, the behavior for loading and unloading can be observed in greater detail.
As a result, the surrogate material model deviates from the purely physical classical elasto-plastic material behavior exactly as intended.
Both models show \re{the} identical physical loading path but differ in unloading.
While the classic model results in the typical hysteresis by dissipation during physical unloading, the virtual unloading follows back the loading path \re{in our surrogate model}.
Therefore, the proposed surrogate model displays a physically \re{reasonable} plastic material behavior but without considering dissipation.

\textit{Remark:} It is worth mentioning that we obtain exactly the behavior as for hyperelastic surrogate model in the 1D case.
\re{However,} this holds true for each individual tensor component which differ in different stress levels in the plastic regime which are determined by the specific strain state.
Consequently, our surrogate material model yields the intended results also for the 3D case in which the calibration of a hyperelastic model is a very challenging task, if possible at all.

Another investigation of the quality of the surrogate model is discussed by the results of a FEM simulation.
To this end, we choose a fix density distribution of the clamped beam (defined in \secref{sec:benchmarkProblems}) given by optimization results (\figref{fig:pic_ClampedBeam_tto_compare_chi}, \num{1} loop).
For this structure and boundary value problem\re{,} both for the surrogate model and classic elasto-plasticity a simulation is applied \re{in which} we ramp the maximum displacement up over \num{20} load steps.
All computations are performed for all plasticity types: ideal, linear hardening and exponential hardening.
The resulting distribution of plastic strains and its relative difference is plotted in \figref{fig:pic_models_fem_compare}.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/pic_models_fem_compare.pdf}
  \caption{Comparing the classic and the surrogate model by a FEM simulation on a given structure which is loaded in \num{20} steps to the maximum load of $u=\SI{-0.05}{\milli\metre}$ accounts for different plastic characteristics show very small differences in the plastic strains.}
  \label{fig:pic_models_fem_compare}
\end{figure}%
The maximum deviation is \re{always} less than \SI{1}{\percent}.
Considering the mathematical difference of the two models, the difference of computed plastic strain is unexpectedly low.

This allow\re{s} us to validate that the surrogate model \re{along with its implementation} address the proposed aspects on \re{the} material point \re{level} and \re{also} confirms \re{accuracy} with\re{in} the FEM.

\FloatBarrier
\subsection{Optimization results with surrogate model for plasticity}
\subsubsection{Benchmark problems and optimization parameters}
\label{sec:benchmarkProblems}
To demonstrate the functionality of the consideration of plasticity in the thermodynamic topology optimization, several boundary value problems are tested.
To this end, we present all considered design spaces with the respective boundary conditions and symmetry planes.

The clamped beam in \figref{fig:pic_ClampedBeam_bvp} is fixated at both sides and chosen in analogy to Maute \textit{et al.}{}~\cite{Maute1998Adaptive}.
The quasi\re{-}2D classical Messerschmitt-Bölkow-Blohm (MBB) beam shown in \figref{fig:pic_MBB_bvp} is simply supported at the lower corner nodes.
Both models are loaded centrally (without symmetr\re{y} plane) on the design space from above.
\begin{figure}[!htb]
  \centering
  \begin{minipage}[b]{0.46\textwidth}
    \centering
    \includegraphics{./03_figures/pic_ClampedBeam_bvp_sym.pdf}
    \caption{Dimensions of design space and boundary value conditions for the quasi\re{-}2D clamped beam.}
    \label{fig:pic_ClampedBeam_bvp}
  \end{minipage}
  \hspace{0.06\textwidth}
  \begin{minipage}[b]{0.46\textwidth}
    \centering
    \includegraphics{./03_figures/pic_MBB_bvp_sym.pdf}
    \caption{Dimensions of design space and boundary value conditions for the quasi\re{-}2D MBB.}
    \label{fig:pic_MBB_bvp}
  \end{minipage}
\end{figure}%
As 3D example, we investigate the boundary value problem given in \figref{fig:pic_3DCantilever_bvp} and denote it as 3D cantilever.
The corners of one side are fixated and the load is exerted at the bottom of the opposite side.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/pic_3DCantilever_bvp_sym.pdf}
  \caption{Dimensions of design space and boundary value conditions for the 3D cantilever.}
  \label{fig:pic_3DCantilever_bvp}
\end{figure}%
All models are discretized by hexahedral finite element meshes with element edge size $e_{\mathrm{size}}$ and linear shape functions.
The thickness of the quasi\re{-}2D models is discretized by one finite element with size $e_{\mathrm{size}}$.

It is worth mentioning that in contrast to topology optimization of linear elastic materials, our results depend in a non-linear way on the amplitude of load (which might be provided either by external forces or prescribed displacements).
Here, the load conditions are applied as prescribed displacements $\boldface{u}^{\star}$ \re{where $\boldface{u}^{\star}_{\max}$ is chosen such that} plasticity \re{evolves} during optimization.

Our novel surrogate model allows to account for a physically reasonable computation of the plastic strains without repeating the entire loading history for each optimization step which is usually necessary to estimate the sensitivity.
Therefore, it is worth mentioning that maximum loading, i.\,e.{}, the loading for which the structure is optimized, can be employed instantaneously.
This is a remarkable difference to other optimization schemes including plasticity.
Since the solution of the finite element problem consumes the highest amount of computation time, our novel approach enables us to save numerical costs \re{by reducing the number of necessary FEM simulations per iteration to even one or a few loops.}

The density variable can be varied in the interval $\left[\chi_{\min}, 1\right]$ where the minimum value is set to $\chi_{\min} = \num{0.001}$.
Therefore, the minimal material stiffness is given by $\num{e-9}\times \mathbb{E}_0$.
The regularization parameter is chosen as $\beta = 2 \, e_{\mathrm{size}}^2 \, \si{\square\milli\metre}$ and the viscosity for all simulations is set to $\eta = \SI{15}{\second}$, corresponding to our previous work \cite{jantos2019accurate}.
All necessary model and optimization parameters for the different boundary value problems are collected in \tabref{tab:modelParameters}.
\begin{table}[!htb]
  \centering
  \caption{Model and optimization parameters.}
  \label{tab:modelParameters}
  \begin{tabular}{l|rlclccc}
    boundary value problem & \#$_\mathrm{elements}$ & $e_{\mathrm{size}} \ [\si{\milli\metre}]$ & $\boldface{u}^{\star} \ [\si{\milli\metre}]$  & $v_{0} \ [-]$ & $\eta \ [\si{\second}]$  & $\beta \ [\si{\square\milli\metre}]$ \\
    \hline \hline
    quasi\re{-}2D clamped beam & \num{5000}\hspace*{2.0mm} & \hspace*{4.8mm}\num{0.02} & \num{0.05} & \hspace*{1.8mm}\num{0.5} & \num{15} & \num{0.80e-3} \\
    quasi\re{-}2D MBB beam & \num{4800}\hspace*{2.0mm} & \hspace*{4.8mm}\num{0.025} & \num{0.02} & \hspace*{1.8mm}\num{0.5} & \num{15} & \num{1.25e-3} \\
    3D cantilever & \num{26364}\hspace*{2.0mm} & \hspace*{4.8mm}\num{1/26} & \num{0.06} & \hspace*{1.8mm}\num{0.15} & \num{15} & \num{2.96e-3} \\
  \end{tabular}
\end{table}

\re{As mentioned,} the stresses lag behind the strains due to the staggered process.
In order to better approximate physics, we compute \numrange{1}{5} FEM simulations within one optimization iteration and before updating the topology for the next time.
This additional simulations are denote\re{d} as loops in the following.

The illustrations of the field data are created with Paraview \cite{Paraview}.
Even if the models make use of symmetry planes, the results are presented as whole (mirrored) in some instances.
The resultant structures are obtained by using the isovolume filter for the density variable $\chi$ with the minimum threshold set to \num{0.5}.
This is the average value of the interval in which $\chi$ has been defined.

\FloatBarrier
\subsection{Optimal structures}
\label{sec:optimizedStructures}
We investigate the impact of inclusion of plasticity on the resultant optimal structure.
To this end, the optimization results are compared with results of thermodynamic topology optimization for a linear elastic material behavior.
This can be achieved while setting the yield stress to an unphysically high value, i.\,e.{} $\sigmaY_\mathrm{exp} = \SI{500000}{\mega\pascal}$.
This ensures that no plastic deformation is active since the von Mises norm of the stress is below this value for all boundary value problems considered.
The results obtained from this elastic optimization are\re{, of course,} consistent with results obtained in our previous publications, cf. \cite{jantos2019accurate}, for instance.
All structures are presented for the converged iteration step.
The structures with shades of green correspond to the thermodynamic topology optimization including plasticity (ideal or hardening) whereas the gray structure is the result for a purely linear elastic topology optimization.

Due to loading, high plastic strains \re{may} occur in the entire design space.
Two result regions with lower stress intensities in topology are possible: i) thicker cross-section areas reduce the maximum value of the averaged stress such that the remaining stress is limited by the yield criterion $r$\re{,} or ii) vanishing \re{sub}structure\re{s} \re{because no stresses occur for void material.}
For an example of the distribution of plastic strains, we refer to \secref{sec:evolutionProcess} \re{in which a} structure is \re{computed where} the highest stresses and thus plastic deformations are present at the constrained boundaries in terms of external loading and supports.
Consequently, we observe thicker member sizes here, cf. the center and supports of the clamped beam in \figref{fig:pic_ClampedBeam_tto_compare_chi} and the center of the MBB in \figref{fig:pic_MBB_tto_compare_chi}.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/pic_ClampedBeam_tto_compare_chi.pdf}
  \caption{Overview of optimized structures \re{and stiffness $S$} for the clamped beam with the surrogate material model accounting for different plastic characteristics and several numbers of inner loops.}
  \label{fig:pic_ClampedBeam_tto_compare_chi}
\end{figure}%
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/pic_MBB_tto_compare_chi.pdf}
  \caption{Overview of optimized structures \re{and stiffness $S$} for the MBB with the surrogate material model accounting for different plastic characteristics and several numbers of loops.}
  \label{fig:pic_MBB_tto_compare_chi}
\end{figure}%
Otherwise, high plastic strains occur in the complete middle part of the design space of the clamped beam.
For this reason, a large void area below the truss corresponding to the loading can be seen in \figref{fig:pic_ClampedBeam_tto_compare_chi}.
This void region is even wider than in the elastic optimization result.
For an increasing number of FEM loops even the lower trusses disappear completely.
This has an impact on the total structure regarding the general topology: due to prescribed total structure volume, the angles and thicknesses of some bars change significantly for the plastic optimization.
Consequently, remarkably different structures are computed when plastic material behavior is considered.

Based on the staggered optimization process, stresses are overestimated.
Therefore, we propose to achieve stresses that are more physical due to the increasing number of additional FEM loops.
Thereby, the displacement field follows the stresses and therefore the plastic strains \re{better reflects the physical reality.}
The optimization results observed for the clamped beam and the MBB confirm the assumed influence, see \figref{fig:pic_ClampedBeam_tto_compare_chi} and \figref{fig:pic_MBB_tto_compare_chi} from left to right.
Both thicker trusses with reduction of structure and rearrangement of thickness in further trusses seem possible.
As a result, plastic strains, corresponding to the displacement field, are crucial for the \re{optimized} structures
Therefore, the number of loops need to be chosen wisely.

\textit{Remark:} It should be mentioned, that the number of \re{\numrange{1}{5}} loops used here is significantly less than the number of load steps required for a path-dependent classic plasticity model\re{, e.\,g.{} \num{20} load increments in \figref{fig:pic_models_fem_compare}}.
Therefore, this proposed method is still fast and efficient.

A special characteristic of our approach is that we can model both ideal plasticity and hardening.
Therefore, we are able to determine that specifications of plasticity result in different structures.
Some structures show small differences but especially with \num{3} loops differences are obvious.
This can be seen for the clamped beam in \figref{fig:pic_ClampedBeam_tto_compare_chi} and for the MBB in \figref{fig:pic_MBB_tto_compare_chi} in a vertical comparison.
For instance, the MBB with \num{3} loops has more truss elements for linear hardening than for ideal plasticity and even more with exponential hardening.
This means, a precisely defined material behavior is important.
Therefore, it is a great advantage that real material behavior can be reproduced in an accurate way by this approach, and can be entered in the optimization process directly in this way.

The optimization results of the 3D cantilever seen in \figref{fig:pic_3DCantilever_tto_compare_chi} proves the functionality and applicability of our approach for fully 3D boundary problems.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/pic_3DCantilever_tto_compare_chi.pdf}
  \caption{Overview of optimized structures \re{and stiffness $S$} for the 3D cantilever with the surrogate material model accounting for different plastic characteristics and several numbers of loops.}
  \label{fig:pic_3DCantilever_tto_compare_chi}
\end{figure}

\FloatBarrier
\subsection{Convergence behavior}
Another aspect of analysis is to discuss the evolution of the optimization objective which is to minimize the compliance of the structure.
Since a compliance minimization analogously causes a stiffness maximization, we use the latter for presentation.
The stiffness is computed in analogy to other works on topology optimization by $S = 1/\boldface{f}\cdot\hat{\boldface{u}}$.
Consequently, we expect a decreasing function for $S$ when the reaction force increases during the evolution of the structure.
The order of magnitude of stiffness is very different for elastic and plastic optimization.
For a convincing representation a logarithmic stiffness axis is chosen.
We define convergence as soon as the relative stiffness changes less than \num{1e-5} for the first time and less than \num{1e-4} for further three succeeding iteration steps.
This rather strict convergence criterion is chosen to exclude a wrong detection of convergence in the plastic case.

The stiffness and iteration step of convergence is plotted for the clamped beam and the MBB with \num{3} loops in \figref{fig:plot_u3_stiffness}, for instance.
\begin{figure}[!htb]
  \centering
  \begin{subfigure}[b]{0.49\textwidth}
    \includegraphics{./03_figures/plot_ClampedBeam_stiffness_u3_log.pdf}
    \caption{Clamped beam (3 loops)}
    \label{fig:plot_ClampedBeam_u3_stiffness}
  \end{subfigure}
  \begin{subfigure}[b]{0.49\textwidth}
    \includegraphics{./03_figures/plot_MBB_stiffness_u3_log.pdf}
    \caption{MBB (3 loops)}
    \label{fig:plot_MBB_u3_stiffness}
  \end{subfigure}
  \caption{Convergence of stiffness evolution during the optimization process.
  The first time that the convergence criterion has been reached is indicated by a vertical line.}
  \label{fig:plot_u3_stiffness}
\end{figure}%
We still see the usual evolution of the stiffness during topology optimization which is that the stiffness increases while a discrete black/white structure evolves.
In the elastic case, the maximum stiffness converges towards a constant value.

The onset of plasticity includes remarkable reduction of stiffness since locally higher strains do not result in higher or stress: the yield stress $\sigma^\mathrm{Y}$ is the fixed upper limit for ideal plasticity, and the increase of stress is slowed down with hardening.
This is a physically reasonable behavior.
Therefore, the stiffness of structures including plasticity is lower than of those which behave purely elastically.
This becomes particularly clear with the clamped beam in \figref{fig:plot_ClampedBeam_u3_stiffness} where \re{larger} values of plastic \re{strains are observed}.
In general, the (absolute value of the) differences in the stiffness plots corresponds to the dissipated energy due to \re{the} plastic \re{formation of} deformations.

Furthermore, the plastic strains are even lower for hardening than for ideal plasticity.
This is caused by the yield criterion $r$ which allows the stresses to increase in a defined manner with hardening.
Therefore, the plots also show a greater stiffness especially for exponential hardening.
\re{Structures with a higher stiffness are thus more similar to elastically optimized structures, cf. the clamped beam in \figref{fig:pic_3DCantilever_tto_compare_chi} with exponential hardening.}

It is remarkable that sometimes plastic optimizations converge in less iteration steps than the elastic optimizations\re{, cf. \figref{fig:plot_u3_stiffness}.}
\re{The number of convergence iterations is a major factor for the difference in computation time in plastic and elastic optimizations.
This can be seen by comparing the runtimes for elastic and plastic optimizations with \num{1} loop in \tabref{tab:convergenceRuntime}.
\begin{table}[!htb]
  \centering
  \caption{\re{Convergence iteration steps and relative runtime according to elastic optimization.}}
  \label{tab:convergenceRuntime}
  \re{\begin{tabular}{l|lrrrr}
    boundary value problem & type of plasticity & \multicolumn{2}{c}{\makecell{convergence\\iteration step}} & \multicolumn{2}{c}{rel. runtime} \\

     & & \textit{1 loop} & \textit{3 loops} & \textit{1 loop} &\textit{3 loops} \\
    \hline \hline
     & elastic & \multicolumn{2}{c}{\num{96}} & \multicolumn{2}{c}{\num{1.0}} \\
     & ideal plastic & \num{167}\hspace*{2.4mm} & \num{87}\hspace*{2.4mm} & \num{3.25}\hspace*{1.8mm} & \num{5.12}\hspace*{1.8mm} \\
     & linear hardening & \num{295}\hspace*{2.4mm} & \num{242}\hspace*{2.4mm} & \num{6.15}\hspace*{1.8mm} & \num{14.75}\hspace*{1.8mm} \\
    \multirow{-4}{*}{quasi\re{-}2D clamped beam} & exp. hardening & \num{237}\hspace*{2.4mm} & \num{157}\hspace*{2.4mm} & \num{4.96}\hspace*{1.8mm} & \num{9.00}\hspace*{1.8mm} \\
    \hline
     & elastic & \multicolumn{2}{c}{\num{273}} & \multicolumn{2}{c}{\num{1.0}} \\
     & ideal plastic & \num{292}\hspace*{2.4mm} & \num{191}\hspace*{2.4mm} & \num{1.10}\hspace*{1.8mm} & \num{2.26}\hspace*{1.8mm} \\
     & linear hardening & \num{201}\hspace*{2.4mm} & \num{227}\hspace*{2.4mm} & \num{0.71}\hspace*{1.8mm} & \num{2.50}\hspace*{1.8mm} \\
    \multirow{-4}{*}{quasi\re{-}2D MBB} & exp. hardening & \num{125}\hspace*{2.4mm} & \num{373}\hspace*{2.4mm} & \num{0.44}\hspace*{1.8mm} & \num{3.64}\hspace*{1.8mm} \\
  \end{tabular}}
\end{table}
It is obvious that computation time also increases with the number of loops (cf. \tabref{tab:convergenceRuntime} for \num{3} loops) which remain less than \num{5} loops in our method.}
Therefore, with the surrogate model the needed computational resources for a plastic optimization \re{is comparable to an elastic optimization which is applicable in engineering practice.}

\FloatBarrier
\subsection{\re{Structure e}volution during \re{the} optimization process}
\label{sec:evolutionProcess}
The evolution of the structure and the plastic strains $\boldsymbol{\varepsilon}^{\mathrm{p}}$ during the optimization process is exemplary presented for the clamped beam with \num{3} loops with ideal plasticity and linear hardening in \figref{fig:pic_ClampedBeam_tto_evolution}.
\begin{figure}[!htb]
  \centering
  \includegraphics{./03_figures/pic_ClampedBeam_tto_evolution.pdf}
  \caption{Evolution of increasing/decreasing plastic strains and structure during the optimization process for the clamped beam (3 loops) with ideal plasticity and linear hardening.}
  \label{fig:pic_ClampedBeam_tto_evolution}
\end{figure}%
Especially for the clamped beam with linear hardening\re{,} the \re{structure and void areas} under the influence of plastic strains can be observed, as explained in \secref{sec:optimizedStructures}.
Furthermore, we see that the value of plastic von Mises strains is lower for hardening than for ideal plasticity for the optimization of the clamped beam.

It is worth mentioning that the amount of plastic strains also reduces during the optimization while stiffness increases and thus strains are locally reduced (again).
This can be seen when comparing the iteration steps $n=42$ and $n=100$ at the area of external displacement and support for ideal plasticity in \figref{fig:pic_ClampedBeam_tto_evolution}.
Therefore, it is a crucial property of the proposed material model to reduce plasticity without dissipation.
This proves that the proposed surrogate material model for plasticity without dissipation operates as expected during the optimization process.
It is thus possible to consider the plastic strain evolution simply by considering the current strain while avoiding the repeated computation of the entire loading path.

\FloatBarrier
\section{Conclusion and outlook}
A novel approach to the thermodynamic topology optimization including plasticity was presented.
To avoid the computation of a complete load path for estimating the plastic strain state in every optimization step, a novel surrogate material model was developed.
To this end, the model was constructed to be dissipation-free such that the plastic strains result from pure energy minimization.
The resultant system of governing equations followed as stationarity conditions from an extended Hamilton functional.
The system comprised the field equations for the displacements and the density variable, and an algebraic equation for the plastic strains.
In the algebraic equation for plastic strains\re{,} \re{arbitrary} types of plasticity can be included by defining the \re{related} yield criterion: \re{exemplary we used} ideal plasticity, linear hardening and exponential hardening.
For the numerical implementation, we employed the neighbored element method for solving the weak form of the balance of linear momentum and the strong form of the evolution equation for the density variable.
Thereby, optimization is \re{solved as} evolutionary process in a staggered manner.
We presented both a material point and FEM investigation to demonstrate the general functionality of the novel material model and various finite boundary value problems for optimization.
Significant deviations between optimized structures for purely elastic materials and the surrogate model for plastic deformations could be detected.
Also differences can be observed for different numbers of FEM loops during one iteration step as well as with ideal plasticity, linear or exponential hardening.
All optimizations result in reliable convergence \re{and with a suitable number of iteration steps}.
During the optimization process, our surrogate material model allows both to predict the microstructural state both for increasing and decreasing strain states due to topology optimization: the plastic strains always correspond to a state evolved during pure loading as is the case for the optimized component during real application.
A remarkable numerical advantage is a computation runtime for the optimization including plasticity is comparable to that for an elastic optimization.

These findings provide the following insights for future research: with the staggered process\re{,} the physical reality is always mapped with a time delay and the optimization is based on these results.
We tried to compensate this delay by additional FEM loops within one optimization iteration.
Therefore, it would be particularly interesting for further research to investigate a monolithic \re{treatment of} thermodynamic topology optimization.

\FloatBarrier
\section*{Acknowledgment}
We highly acknowledge the funding of this research by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the project grant JU 3096/2-1.

redpajama_set_name: RedPajamaArXiv