<figure><img src="image_1.png" /></figure>

HAL
open science

# Stress-constrained topology optimization using approximate reanalysis with on-the-fly reduced order modeling

M. Xiao, J Ma, D. Lu, B. Raghavan, W. Zhang

► To cite this version:

M. Xiao, J Ma, D. Lu, B. Raghavan, W. Zhang. Stress-constrained topology optimization using approximate reanalysis with on-the-fly reduced order modeling. Advanced Modeling and Simulation in Engineering Sciences, 2022, 9 (1), pp.17. 10.1186/s40323-022-00231-x. http://hal-03810665

HAL Id: hal-03810665

https://hal.science/hal-03810665

Submitted on 11 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Stress-constrained topology optimization using approximate reanalysis with on-the-fly reduced order modeling

Manyu Xiao<fnref n="1" />*, Jun Ma<fnref n="1" />, Dongcheng Lu<fnref n="1" />, Balaji Raghavan<fnref n="2" /> and Weihong Zhang<fnref n="3" />

## Abstract

Most of the methods used today for handling local stress constraints in topology optimization, fail to directly address the non-self-adjointness of the stress-constrained topology optimization problem. This in turn could drastically raise the computational cost for an already large-scale problem. These problems involve both the equilibrium equations resulting from finite element analysis (FEA) in each iteration, as well as the adjoint equations from the sensitivity analysis of the stress constraints. In this work, we present a paradigm for large-scale stress-constrained topology optimization problems, where we build a multi-grid approach using an on-the-fly Reduced Order Model (ROM) and the p-norm aggregation function, in which the discrete reduced-order basis functions (modes) are adaptively constructed for adjoint problems. In addition to reducing the computational savings due to the ROM, we also address the computational cost of the ROM learning and updating phases. Both reduced-order bases are enriched according to the residual threshold of the corresponding linear systems, and the grid resolution is adaptively selected based on the relative error in approximating the objective function and constraint values during the iteration. The tests on 2D and 3D benchmark problems demonstrate improved performance with acceptable objective and constraint violation errors. Finally, we thoroughly investigate the influence of relevant stress constraint parameters such as the p norm factor, stress penalty factor, and the allowable stress value.

**Keywords:** Stress constraints, Topology optimization, PCA, ROM, Adjoint equations

## Introduction

Mechanical stress is a critical factor that affects performance, service life, fatigue resistance and safety of structural components, and is inarguably an important design criterion. Traditional topology optimization formulations (such as the minimum compliance problem with a volume constraint [1–3]) typically do not consider stress constraints, which could result in the phenomenon of high stress concentrations, leading to a final "optimized" design that all too often fails to meet real engineering requirements, eventually requiring manual adjustment or shape optimization as a subsequent step.

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence:
manyuxiao@nwpu.edu.cn</footnote>

<footnote n="1">¹Xi’an Key Laboratory of Scientific Computation and Applied Statistics, School of Mathematics and Statistics, Northwestern Polytechnical University, Xi’an 710072, China</footnote>

<footnote n="2">²Laboratoire de Génie Civil et Génie Mécanique EA 3913, Institut National des Sciences Appliquées de Rennes, Rennes, France</footnote>

<footnote n="3">³State IJR Center of Aerospace Design and Additive Manufacturing, Northwestern Polytechnical University, Xi’an 710072, China</footnote>
In the past two decades, many researchers have worked on stress-constrained topology optimization [4,5] in order to ensure that the design results are closer to engineering reality. However, when compared with the stress-constraint-free problem, stress-constrained topology optimization design presents a unique set of difficulties, which may be broadly classified as follows:

1. The stress singularity phenomenon. This was first studied for the truss optimization problem. Cheng and Jiang [6,7] gave a mathematical explanation for the singularity in the topology optimization problem, attributing it to the discontinuity of the stress constraint. Duysinx and Bendsøe [8] linked stress singularity to the phenomenon of non-vanishing stresses as the design variables tended towards zero, in other words, a region with a very small relative density could end up with a non-zero strain which then gives rise to non-zero (or even very high) stress, despite the fact that it represents a hole (void area). Today, the singularity problem is typically resolved by adopting the relaxation approach [7–9], QP-relaxation method [10] and stress penalization [11].

2. Large number of local stress constraints. Since stress is a local parameter, the greater the number of elements for the domain discretization, the more there will be stress evaluation points, resulting in a large number of local non-convex stress constraints, which in turn will lead to a steep rise in computational cost for the sensitivity analysis of local stress constraints [4,12]. To reduce the computational complexity, local stress constraints are usually transformed into a global stress constraint or some form of cluster-based stress constraints by using a condensation function that approximates the value of the maximum local function. Currently, the two most common types of condensation functions are the Kreisselmeier-Steinhauser function (K-S) [4, 12–16] and the P-norm function [9, 16–18], and they can merge a mass of local stresses into a single constraint formula.

3. Highly nonlinear nature of the stress constraint. Any change in density value of the adjacent region will significantly affect the stress level within some key regions [13,17] such as depressions and corners. This requires the optimization procedure to be able to effectively reduce or eliminate the stress concentration phenomenon, and that the solution algorithm must maintain numerical consistency with the optimization procedure in order to ensure stable convergence of the overall procedure [19,20].

In view of the three broad classes of problems associated with stress-based topology optimization listed above, researchers have put forward a variety of effective methods to deal with them.

However, with modern-day sophisticated numerical methods and the recent development of advanced manufacturing technologies, the performance requirements for mechanical parts in the aerospace industry have become increasingly demanding [21]. Modeling with highly refined finite element models typically results in a series of ultra-large-scale problems with degrees of freedom ranging from tens of millions to even billions. Given the presence of non-convex, non-linear stress constraints in large-scale topology optimization, the main hindrances then are the solution of large-scale linear systems during the FEA analysis (equilibrium equation) in each iteration, as well as the sensitivity analysis of stress constraint functions due to the non self-adjoint nature [22]. Stress coag-
ulation alone is then inadequate to improve the design, meaning that we need another approach [23].

One highly effective and proven solution is Approximate Reanalysis (AR) using a Reduced Order Model (ROM) technique [24–27]. Due to the non-self-adjointness of the problem, we need to solve two high-dimensional systems during each iteration of the optimization process, the first being the equilibrium equations from finite element analysis (FEA), the other being the adjoint equations for the sensitivity analysis of the stress constraints. To the authors’ knowledge, virtually no work has been reported on coupling stress-constrained topology optimization with approximate reanalysis using a ROM.

In this work, we first apply two previously proposed ROM approaches to the stress-
constrained topology optimization problem, the first being on-the-fly ROM using Prin-
cipal Components Analysis (PCA) [28] and the other being multi-fidelity ROM using the
multigrid method (MG) [29]. In addition, based on the proposed multi-fidelity ROM, we
develop two improved strategies for further reductions in computational costs, the first
one being the Adap-ROM where the residual threshold for the approximate solution of
the equilibrium equations can be adaptively selected according to the relative error of the
objective function value during the iteration; the second strategy being Bi-ROM where we
simultaneously apply the ROM to both the equilibrium equations and adjoint equations
(during the sensitivity analysis of the stress constraint). Instead of the classical Optimality
Criteria (OC) method, the Method of Moving Asymptotes (MMA) [30] is used as the
optimization algorithm in this work, and the P-norm function in [9] is adopted for coag-
ulating the local stress constraints to reduce the computational cost involved. The stress
singularity phenomenon is avoided by using stress penalization as in [11]. Recent years
there are more papers about stress constraint [31–33].

This paper is organized in the following manner: in “Theoretical formulation”, the
theoretical formulation is presented with a classical minimum compliance problem con-
strained by both stress and volume, followed by adjoint equation-based sensitivity analysis
of both the objective function and the stress constraint function. “Projection-based ROM
approaches” first introduces two projection-based ROM approaches including on-the-fly
ROM construction using PCA and multi-fidelity ROM construction using the multigrid
method, then combines both approaches into a single algorithm. In “Two improved strate-
gies”, we propose two improved strategies (i.e. Adap-ROM and bi-ROM). “Numerical
tests” details the numerical investigations using a benchmark 2D L-shaped beam exam-
ple. In “3D L-shaped beam”, we first apply a variety of ROM methods to the 3D compliance
minimization problem of L-shaped beam, and then analyze the impact of the p norm fac-
tor p, the stress penalty factor q and the allowable maximum stress value σ on the design
results. The paper ends with discussion, concluding comments and recommendations for
future work.

**Theoretical formulation**

Topology optimization model

The essence of topology optimization is mathematically determining the distribution of
material in a geometrically fixed design domain under a set of constraints, such that
the resulting structure achieves some optimal combination of mechanical performance
objectives. Here, we consider the classical minimum compliance problem under stress
and volume constraints, expressed as follows:

$$ \begin{array}{l} \min_{\rho} c(\rho) = f^T u = u^T K u \\ \text{s.t.} \quad \left\{ \begin{array}{l} K(\rho)u = f \\ v_f = \frac{\sum_{e=1}^{N_e} \rho_e V_e}{V_0} \le \hat{v}_f \\ \sigma_{\max}^{vM} = \max(\sigma_e^{vM}) \le \hat{\sigma} \\ 0 \le \rho_e \le 1, e = 1, \dots, N_e \end{array} \right. \end{array} \quad (1) $$

where, $\rho$ is the vector consisting of the design variables (i.e. element (e) density), and $N_e$ is the total number of design variables. $\rho$ denotes the compliance of the structure (objective function), $f$ and $u$ are external load vector and displacement vector of the structure respectively, $K$ is the global stiffness matrix of the structure. $V_e$ is the volume of element e and $V_0$ is the initial volume of the structure, occupying the entire design domain, $v_f$ and $\hat{v}_f$ is the material volume fraction and its allowable value respectively. $\sigma_e^{vM}$ represents the von Mises stress for an arbitrary element, $\sigma_e^{vM}$ and $\hat{\sigma}$ are the maximum von Mises stress of structure and its constraint limit value.

To suppress the generation of gray elements (intermediate values between 0 and 1) and obtain a clear black-and-white pattern, we use a modified solid isotropic material with penalization (SIMP) interpolation [31], where the stiffness is penalized by a power function relation between the Young's modulus $E_e$ and element relative density $\rho_e$, given by:

$$ E_e(\rho_e) = E_{\min} + \rho_e^{\omega}(E_0 - E_{\min}), \quad (2) $$

where, $\omega > 1$ is the SIMP penalty factor. is the Young's modulus of solid material while $E_{\min}$ is a vanishingly small positive value to avoid numerical instabilities such as a singular global stiffness matrix. In this paper, $\omega$ assumes an empirical value of 3 while $E_{\min} = 10^{-9}$.

To alleviate the singular solution problem, in a manner similar to stiffness penalization in SIMP model, we penalize the stress vector in following scheme [12]:

$$ \sigma_e = (\rho_e)^q \sigma_e^0 = (\rho_e)^q D_0 B u_e, \quad (3) $$

where, $q < 1$ is the stress penalty factor, $D_0$ is the constitutive matrix of solid materials and $B$ is strain-displacement matrix. $\sigma_e^0 = D_0 B u_e$ is the stress vector (Voight notation) calculated at the geometric center of the FE element (e), and the $\sigma_e^0$ has a different form in 2D and 3D case, which is given by:

$$ \sigma_e^0 = \begin{cases} (\sigma_{ex}^0, \sigma_{ey}^0, \tau_{exy}^0)^T & (2D \text{ case}) \\ (\sigma_{ex}^0, \sigma_{ey}^0, \sigma_{ez}^0, \tau_{exy}^0, \tau_{eyz}^0, \tau_{ezx}^0)^T & (3D \text{ case}). \end{cases} \quad (4) $$

The corresponding von Mises stress is then calculated as follows [34]:

$$ \sigma_e^{vM} = \begin{cases} (\sigma_{ex}^2 + \sigma_{ey}^2 - \sigma_{ex}\sigma_{ey} + 2\tau_{ey}^2)^{\frac{1}{2}} & (2D \text{ case}) \\ (\sigma_{ex}^2 + \sigma_{ey}^2 + \sigma_{ez}^2 - \sigma_{ex}\sigma_{ey} - \sigma_{ey}\sigma_{ez} - \sigma_{ez}\sigma_{ex} + 3\tau_{exy}^2 + 3\tau_{eyz}^2 + 3\tau_{ezx}^2)^{\frac{1}{2}} & (3D \text{ case}). \end{cases} \quad (5) $$

In order to reduce the calculation cost caused by massive local stress constraints, we adopt the P-norm function [14] to coagulate the $N_e$ local stress constraints in Eq. (1) so that to form a global stress constraint, given by:

$$ \sigma^{PN} = \left( \sum_{e=1}^{N_e} (\sigma_e^{vM})^p \right)^{\frac{1}{p}}, \quad (6) $$
where, $p$ may be interpreted as the P-norm factor. And it determines the difference between the original function and its approximate value. When larger values are used, more weight is given to peak stress. Theoretically, as the value of the parameter $p$ approaches infinity, we have:

$$ \lim_{p \to \infty} \sigma^{PN} = \max(\sigma_e^{vM}) = \sigma_{\max}^{vM}. \qquad (7) $$

Thus, the global stress constraint in Eq. (1) is transformed into,

$$ \sigma^{PN} \le \hat{\sigma}. \qquad (8) $$

However, it's impossible to choose an infinite value of $p$ and global measure function $\sigma^{PN}$ may not approximate $\sigma_{\max}^{vM}$ very well. Therefore, we use following Eq. (8) to modify the constraint in Eq. (7),

$$ \tilde{\sigma}^{PN} = \xi \sigma^{PN} \le \hat{\sigma}, \qquad (9) $$

where, adjusted parameter $\xi = \sigma_{\max}^{vM}/\sigma^{PN}$.

Generally, a high value of $p$ could create problems with numerical accuracy since $(\sigma_e^{vM})^p$ could be very large. Actually when a high value of $p$ is used, the peak stress is weighted more heavily, it often leads to oscillation or even divergence during the optimization iterations. The oscillation results from the high non-linearity of the global stress function. One work-around is to normalize the stress constraint, which is given by:

$$ \tilde{\sigma}_{norm}^{PN} = \xi \left( \sum_{e=1}^{N_e} (\sigma_e^{vM})^p \right)^{\frac{1}{p}} \le 1. \qquad (10) $$

## Sensitivity analysis

Using the adjoint method, the sensitivity of the objective function in Eq. (1) can be elegantly derived as follows:

$$ \frac{\partial c}{\partial \rho_e} = -\mathbf{u}^T \frac{\partial \mathbf{K}}{\partial \rho_e} \mathbf{u}, \qquad (11) $$

where, Eq. (10) is the well-known self-adjoint gradient as seen in [29].

Using the chain rule, we obtain the sensitivity of the stress constraint function $\tilde{\sigma}_{norm}^{PN}$

$$ \frac{\partial \tilde{\sigma}_{norm}^{PN}}{\partial \rho_e} = \frac{\partial \tilde{\sigma}_{norm}^{PN}}{\partial \sigma_e^{vM}} \left( \frac{\partial \sigma_e^{vM}}{\partial \sigma_e} \right)^T \frac{\partial \sigma_e}{\partial \rho_e}, \qquad (12) $$

where, the three terms are respectively calculated as follows:

1. the derivative of normalized P-norm function with respect to von Mises stress.

$$ \begin{aligned} \frac{\partial \tilde{\sigma}_{norm}^{PN}}{\partial \sigma_e^{vM}} &= \frac{\xi}{p \hat{\sigma}} \left( \sum_{e=1}^{N_e} \left( \frac{\sigma_e^{vM}}{\hat{\sigma}} \right)^p \right)^{\left(\frac{1}{p}-1\right)} \times p \left( \frac{\sigma_e^{vM}}{\hat{\sigma}} \right)^{p-1} \\ &= \frac{\xi}{\hat{\sigma}} \left( \sum_{e=1}^{N_e} \left( \frac{\sigma_e^{vM}}{\hat{\sigma}} \right)^p \right)^{\left(\frac{1}{p}-1\right)} \times \left( \frac{\sigma_e^{vM}}{\hat{\sigma}} \right)^{p-1}. \end{aligned} \qquad (13) $$

2. The derivative of the von Mises stress with respect to stress components.

For the 2D case:

$$ \left\{ \begin{aligned} \frac{\partial \sigma_e^{vM}}{\partial \sigma_{ex}} &= \frac{1}{2\sigma_e^{vM}} (2\sigma_{ex} - \sigma_{ey}) \\ \frac{\partial \sigma_e^{vM}}{\partial \sigma_{ey}} &= \frac{1}{2\sigma_e^{vM}} (2\sigma_{ey} - \sigma_{ex}) \\ \frac{\partial \sigma_e^{vM}}{\partial \sigma_{ez}} &= \frac{2\tau_{exy}}{\sigma_e^{vM}}. \end{aligned} \right. \qquad (14) $$
For the 3D case:

$$ 
\left\{ 
\begin{aligned}
\frac{\partial \sigma_e^{vM}}{\partial \sigma_{ex}} &= \frac{1}{2\sigma_e^{vM}} (2\sigma_{ex} - \sigma_{ey} - \sigma_{ez}) \\
\frac{\partial \sigma_e^{vM}}{\partial \sigma_{ey}} &= \frac{1}{2\sigma_e^{vM}} (2\sigma_{ey} - \sigma_{ex} - \sigma_{ez}) \\
\frac{\partial \sigma_e^{vM}}{\partial \sigma_{ez}} &= \frac{1}{2\sigma_e^{vM}} (2\sigma_{ez} - \sigma_{ex} - \sigma_{ey}) \\
\frac{\partial \sigma_e^{vM}}{\partial \tau_{exy}} &= \frac{3\tau_{exy}}{\sigma_e^{vM}}, \quad \frac{\partial \sigma_e^{vM}}{\partial \tau_{eyz}} = \frac{3\tau_{eyz}}{\sigma_e^{vM}}, \quad \frac{\partial \sigma_e^{vM}}{\partial \tau_{ezx}} = \frac{3\tau_{ezx}}{\sigma_e^{vM}}
\end{aligned}
\right. 
\qquad (15) $$

3. The derivative of the stress components with respect to the design variables:

$$ \frac{\partial \sigma_e}{\partial \rho_e} = q \rho_e^{q-1} \mathbf{D}_0 \mathbf{B} \mathbf{u}_e + \rho_e^q \mathbf{D}_0 \mathbf{B} \frac{\partial \mathbf{u}_e}{\partial \rho_e}, \qquad (16) $$

where, the term $\frac{\partial \mathbf{u}_e}{\partial \rho_e}$ is obtained from the equilibrium equations:

$$ \mathbf{K} \mathbf{u} = \mathbf{f}. \qquad (17) $$

By calculating the derivative of both sides of this equation, we have:

$$ \frac{\partial \mathbf{K}}{\partial \rho_e} \mathbf{u}_e + \mathbf{K} \frac{\partial \mathbf{u}_e}{\partial \rho_e} = \mathbf{0}. \qquad (18) $$

Giving,

$$ \frac{\partial \mathbf{u}_e}{\partial \rho_e} = -\mathbf{K}^{-1} \frac{\partial \mathbf{K}}{\partial \rho_e} \mathbf{u}_e. \qquad (19) $$

Substituting Eq. (18) into Eq. (15) and then Eq. (15) into Eq. (11) gives,

$$ \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \rho_e} = \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \sigma_e^{vM}} \left( \frac{\partial \sigma_e^{vM}}{\partial \sigma_e} \right)^T \times \left( q \rho_e^{q-1} \mathbf{D}_0 \mathbf{B} \mathbf{u}_e - \rho_e^q \mathbf{D}_0 \mathbf{B} \mathbf{K}^{-1} \frac{\partial \mathbf{K}}{\partial \rho_e} \mathbf{u}_e \right). \quad (20) $$

Here, an adjoint variable $\lambda$ is defined as:

$$ \lambda^T = \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \sigma_e^{vM}} \left( \frac{\partial \sigma_e^{vM}}{\partial \sigma_e} \right)^T \mathbf{D}_0 \mathbf{B} \mathbf{K}^{-1}. \qquad (21) $$

Thus, we obtain the adjoint equations:

$$ \mathbf{K} \lambda = \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \sigma_e^{vM}} \mathbf{B}^T \mathbf{D}_0^T \frac{\partial \sigma_e^{vM}}{\partial \sigma_e}. \qquad (22) $$

Then, by solving Eq. (22), we obtain the adjoint variable $\lambda$ which is inserted into Eq. (20) so that a final gradient of $\tilde{\sigma}_{\text{norm}}^{PN}$ can be expressed as,

$$ \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \rho_e} = \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \sigma_e^{vM}} \left( \frac{\partial \sigma_e^{vM}}{\partial \sigma_e} \right)^T q \rho_e^{q-1} \mathbf{D}_0 \mathbf{B} \mathbf{u}_e - \rho_e^q \lambda^T \frac{\partial \mathbf{K}}{\partial \rho_e} \mathbf{u}_e, \qquad (23) $$

where, the terms $\partial \tilde{\sigma}_{\text{norm}}^{PN}/\partial \sigma_e^{vM}$ and $\partial \sigma_e^{vM}/\partial \sigma_e$ are determined by Eq. (13) and Eq. (14) (or Eq. (15)) respectively. And as we know that,

$$ \frac{\partial \mathbf{K}}{\partial \rho_e} = \omega \rho_e^{\omega-1} \tilde{E}_e \mathbf{K}^0. \qquad (24) $$

So we can substitute Eq. (24) into equation Eq. (23) to get

$$ \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \rho_e} = \frac{\partial \tilde{\sigma}_{\text{norm}}^{PN}}{\partial \sigma_e^{vM}} \left( \frac{\partial \sigma_e^{vM}}{\partial \sigma_e} \right)^T q \rho_e^{q-1} \mathbf{D}_0 \mathbf{B} \mathbf{u}_e - \omega \rho_e^{q+\omega-1} \lambda^T \tilde{E}_e \mathbf{K}^0 \mathbf{u}_e. \qquad (25) $$
# Projection-based ROM approaches

## On-the-fly ROM construction using PCA

For large-scale topology optimization problems, the main computational workload stems from the numerical solution of the equilibrium equation at each iteration:

$$ Ku = f, \qquad (26) $$

where, the computation of this high-precision full-field solution involves the inversion of a high-dimensional linear systems that can consist of up to millions or billions of degrees of freedom (DOF) [30]. Since the coefficient matrix K is typically ill-conditioned (due to the SIMP interpolation) and also has a large number of DOFs, traditional methods depending on iteration methods like PCG may still be inadequate to significantly reduce the computational cost.

To alleviate this solution process, the authors call on a previously-proposed method called on-the-fly ROM construction using PCA [28]. The gist of this approach is mapping the displacement field solution of the large-scale problem Eq. (26) to a low-dimensional space spanned by an appropriately obtained orthonormal basis $\Phi(\Phi^T\Phi = I)$ that is constructed on-the-fly and updated by implementing PCA for a set of stored solution snapshots obtained during the previous iterations. Here, we briefly present the salient features of the approach.

Given a snapshot matrix consisting of $N_b$ displacement solutions at previous $N_b$ successive iterations,

$$ U = [u_1, u_2, \dots, u_{N_b}] \in \mathbb{R}^{N \times N_b}, \qquad (27) $$

by zero-mean processing for matrix U, we obtain its centered version $\tilde{U}$,

$$ \tilde{U} = [u_1 - \tilde{u}, u_2 - \tilde{u}, \dots, u_{N_b} - \tilde{u}], \quad \tilde{u} = \frac{1}{N_b} \sum_{i=1}^{N_b} u_i, \qquad (28) $$

then performing singular value decomposition (SVD) on the centered matrix $\tilde{U}$, we have,

$$ \bar{U} = \Psi_{N \times N_b} \Sigma_{N_b \times N_b} V_{N_b \times N_b}^T, \qquad (29) $$

the reduced basis $\Phi$ is obtained by truncating the first $M \ll N_b$ columns corresponding to the largest $N_b$ singular values in matrix $\Psi$ which is composed of $N_b$ left singular vectors $\phi_i$, $i = 1, \dots, N_b$,

$$ \Phi \leftarrow [\phi_1, \phi_2, \dots, \phi_M], \qquad (30) $$

usually, the value of M (i.e. the size of reduced basis $\phi$), is determined by following energy rule:

$$ 1 - \frac{\sum_{i=1}^{M} \kappa_i^2}{\sum_{i=1}^{N_b} \kappa_i^2} \le \varepsilon_{\text{trunc}}, \qquad (31) $$

where, $\kappa_i$ denotes the ith largest singular value of centered matrix $\tilde{U}$ and $\varepsilon_{\text{trunc}}$ is a user-defined truncation error (in this paper 1%).

Now, by mapping the displacement field $u$ on to the low-dimensional space spanned by the reduced basis $\Phi$ constructed above, we have the following ROM expression,

$$ u \approx \tilde{u} = \Phi\alpha + \tilde{u}. \qquad (32) $$

Substituting the approximate solution $\tilde{u}$ into Eq. (26) and then projecting the equations onto the basis, we get a reduced system with a much smaller dimensionality,

$$ \Phi^T K \Phi \alpha = \Phi^T (f - K \bar{u}), \qquad (33) $$
the high computational cost of solving full-scale system Eq. (26) may then be greatly reduced through a significantly less expensive solution of Eq. (33). Finally, a linear combination of projection coefficient $\alpha$ as in Eq. (32) is used to recover the full-filed solution.

The approximate accuracy of $\tilde{u}$ is monitored during the optimization via the residual equilibrium error which is given by,

$$e_{rb} = \frac{\|K\tilde{u} - f\|}{\|f\|}, \qquad (34)$$

when the force residual of a certain optimization loop iteration (after at least $N_b$ iterations have been completed so as to be able to obtain the basis) is lower than a pre-specified threshold $\hat{e}_{rb}$, we accept the reduced-basis solution in place of the full-field solution. Once the force residual exceeds the tolerance $\hat{e}_{rb}$, we are forced to recompute a full-field solution of Eq. (26) by the direct method, and then use this new solution to replace the “oldest” stored displacement vector in current snapshot matrix. In order to reduce the storage requirements for large scale problems, recalculating PCA basis for a single snapshot update, the incremental PCA can be considered as shown in the reference [35]. In this way, the original basis is updated and improved using PCA for the next loop iteration, which forms the so-called on-the-fly construction and updating of the ROM.

Here two points have to be noted. Firstly, the construction of the snapshot matrix, requiring the solution of Nb linear systems, may be done by an iterative solver (see [28]), avoiding full-size K construction. Secondly, the projection of the linear problem in Eq. (33) also does not require the assembly of the global stiffness matrix K, as the reduced matrix $\Phi^T K \Phi$ and the RHS term $\Phi^T K$ may be assembled directly by projecting individual element contributions.

## Multi-fidelity ROM construction using multigrid

Topology optimization with ROM-based reanalysis (as seen in “On-the-fly ROM construction using PCA”) still suffers from a significant computational burden, as we need to repeatedly update the reduced basis to obtain a "sufficiently good" prediction of the displacement field when the accuracy of the ROM is insufficient. The new snapshot is then treated as a full-field solution calculated on a fine fixed FE mesh (i.e. single-fidelity), which should be effectively computed in some way. So, based on the above mentioned method, [26] proposed an approach of multi-fidelity ROM construction using the multi-grid method, where the original equilibrium equations are solved on grids of various resolutions, and then quickly obtain a multi-fidelity snapshot solution for constructing and updating the reduced basis. Next, A typical "V" cycle MG iteration of recursive form is given in Algorithm 1 as follows ($\ell = 1, 2, \dots, n-1$):

<table><thead><tr><th>Algorithm 1</th><th>MG (K<sup>(ℓ)</sup>, f<sup>(ℓ)</sup>, u<sup>(ℓ)</sup>, ℓ)</th></tr></thead><tbody><tr><td>1:</td><td>procedure INITIALIZATION: (u<sup>(ℓ)</sup>, {K<sup>(ℓ+1)</sup>}, {P<sup>(ℓ+1)</sup>}, {W<sub>ℓ</sub>})</td></tr><tr><td>2:</td><td>if ℓ + 1 = n then</td></tr><tr><td>3:</td><td>u<sup>(ℓ+1)</sup> ← (K<sup>(ℓ+1)</sup>)<sup>-1</sup>f<sup>(ℓ+1)</sup></td></tr><tr><td>4:</td><td>else</td></tr><tr><td>5:</td><td>K<sup>(ℓ)</sup>u<sup>(ℓ)</sup> = f<sup>(ℓ)</sup>: u<sup>(ℓ)</sup> ← W<sub>ℓ</sub>(K<sup>(ℓ)</sup>, f<sup>(ℓ)</sup>, u<sup>(ℓ)</sup>)</td><td>▷ Pre-smoothing, apply β<sub>1</sub>(β<sub>1</sub> ≤ 3)</td></tr><tr><td>6:</td><td>f<sup>(ℓ+1)</sup> = (P<sup>(ℓ+1)</sup>)<sup>T</sup>(f<sup>(ℓ)</sup> - K<sup>(ℓ)</sup>u<sup>(ℓ)</sup>), u<sup>(ℓ+1)</sup> = 0</td><td></td></tr><tr><td>7:</td><td>u<sup>(ℓ+1)</sup> ← MG((K<sup>(ℓ+1)</sup>, f<sup>(ℓ+1)</sup>, u<sup>(ℓ+1)</sup>, ℓ + 1))</td><td>▷ Recursive iteration</td></tr><tr><td>8:</td><td>u<sup>(ℓ)</sup> ← u<sup>(ℓ)</sup> + P<sup>(ℓ)</sup><sub>l+1</sub>u<sup>(ℓ+1)</sup></td><td>▷ Coarse-grid correction</td></tr><tr><td>9:</td><td>K<sup>(ℓ)</sup>u<sup>(ℓ)</sup> = f<sup>(ℓ)</sup>: u<sup>(ℓ)</sup> ← W<sub>ℓ</sub>(K<sup>(ℓ)</sup>, f<sup>(ℓ)</sup>, u<sup>(ℓ)</sup>)</td><td>▷ Post-smoothing, apply β<sub>2</sub>(β<sub>2</sub> ≤ 3)</td></tr></tbody></table>
In this paper, we set the level number $n = 2$ and $\beta_1 = \beta_2 = 3$, i.e. a simple two-grid “V” cycle. And we choose the FE displacement solution obtained from the previous main optimization loop as the initial guess $\mathbf{u}_\ell$, details of the coarse grid operators $K^{(\ell+1)} = (P_{\ell+1}^\ell)^\mathrm{T}K^\ell P_{\ell+1}^\ell$ and prolongation operators can be seen in [24]. The smooth operators use the PCG method. Algorithm 1 starts from $MG(K^{(1)}, f^{(1)}, \mathbf{u}^{(1)}, 1)$, where $K^{(1)} = K$, $f^{(1)} = f$. For the sake of accuracy, we repeat the above “V” cycle loop twice (initial guess for the next loop is selected as the MG solution of the previous loop).

The algorithms for both approaches (introduced in “Projection-based ROM approaches”) are combined in Algorithm 2, where SF-ROM stands for on-the-fly ROM construction using PCA using single-fidelity (SF) snapshots calculated on a fixed FE mesh, and MF-ROM for multi-fidelity ROM construction using multigrids where multi-fidelity (MF) snapshots are obtained by MG iteration. Full-order model (FOM) represents the full-field solution in Eq. (26) obtained using either the direct method (for SF-ROM) or the MG method (for MF-ROM).

<caption>**Algorithm 2 SF-ROM and MF-ROM**</caption>

1: **procedure** INITIALIZATION: (compute global stiffness matrix $K_i$ at $i_{th}$ iteration loop, given button, set $i = 1$)
2:   **if** $1 \le i \le N_b$ **then**
3:     $\mathbf{u}_i \leftarrow \underset{\text{FOM by direct method}}{\text{K}_i \mathbf{u}_i = f}$
4:     **add** $\mathbf{u}_i$ to snapshot matrix $\mathcal{U}$ and compute its centered version $\bar{\mathcal{U}}$
5:     **if** $i = N_b$ **then**
6:       $\Phi \leftarrow \text{SVD}(\bar{\mathcal{U}})$
7:     **else**
8:       $\alpha \leftarrow \Phi^\mathrm{T} K_i \Phi \alpha = \Phi^\mathrm{T}(f - K_i \bar{\mathbf{u}})$
9:       $e_{rb} \leftarrow \|K_i \bar{\mathbf{u}} - f\|$.
10:      **if** $e_{rb} > \hat{e}_{rb}$ **then**
11:        remove the oldest snapshot from $\mathcal{U}$
12:        **if** button=1 **then**
13:          $\mathbf{u}_i \leftarrow \underset{\text{FOM by direct method}}{\text{K}_i \mathbf{u}_i = f}$ ▷ SF-ROM
14:        **else**
15:          $\mathbf{u}_i \leftarrow \underset{\text{FOM by MG method}}{\text{K}_i \mathbf{u}_i = f}$ ▷ MF-ROM
16:          **add** $\mathbf{u}_i$ to snapshot matrix $\mathcal{U}$ and compute its centered version $\bar{\mathcal{U}}$
17:          $\Phi \leftarrow \text{SVD}(\bar{\mathcal{U}})$
18:        **else**
19:          $\bar{\mathbf{u}}_i \leftarrow \underset{\text{ROM}}{\Phi \alpha + \bar{\mathbf{u}}}$
20:        **Update** the cyclic variable : $i = i + 1$ Back to line 1

## Two improved strategies

### Adap-ROM

Whether we use the SF-ROM or the MF-ROM, a pivotal ROM parameter, the force residual threshold $\hat{e}_{rb}$ must always be predefined by users. Sometimes it's hard to give an appropriate value of $\hat{e}_{rb}$ as it clearly involves a trade-off between calling the ROM solution and the FOM solution during the optimization procedure. On the one hand, if we use a smaller value, then the optimal objective may have better accuracy but would require more CPU time as the FOM would ostensibly be called more frequently, on the other hand, if we use a bigger value, then the optimal objective may have poorer accuracy but would require less CPU time as the ROM would tend to be called more frequently. Therefore, in this paper, we propose an adaptive selection strategy for the value of force residual threshold which is then chosen according to the relative error variation of objective function value during the iteration loop.

$$e_c = \frac{\|c_{i-1} - c_{i-2}\|}{\|c_{i-2}\|} (i \ge 3), \quad (35)$$
$c_{i-1}$, $c_{i-2}$ are the objective function values at two consecutive iterations $i-1_{th}$ and $i-2_{th}$, here the relative error is gradually decreasing because of the convergence behavior of the compliance minimization problem. Details of the strategy are summarized in the following Algorithm 3.

<table><caption>**Algorithm 3 Adaptive parameter selection of force residual threshold**</caption><tr><td>1:</td><td>procedure INPUT:($c_{i-1}$, $c_{i-2}$, $\hat{e}_{rb}^0$)</td><td>▷ Two Objective function values $c_{i-1}$, $c_{i-2}$</td></tr><tr><td>2:</td><td></td><td>▷ Initial force residual threshold $\hat{e}_{rb}^0$</td></tr><tr><td>3:</td><td>procedure INPUT:($\theta$, $\hat{e}_c^0$, $\zeta$)</td><td>▷ Scale factor $\theta$, the relative error threshold $\hat{e}_c^0$</td></tr><tr><td>4:</td><td></td><td>▷ Decay factor $\zeta$</td></tr><tr><td>5:</td><td>if $e_c > \hat{e}_c^0$ then</td><td></td></tr><tr><td>6:</td><td>$\hat{e}_{rb} = 0$</td><td>▷ No ROM</td></tr><tr><td>7:</td><td>if $\zeta \hat{e}_c^0 < e_c \le \hat{e}_c^0$ then</td><td></td></tr><tr><td>8:</td><td>$\hat{e}_{rb} = \hat{e}_{rb}^0$</td><td></td></tr><tr><td>9:</td><td>else</td><td></td></tr><tr><td>10:</td><td>set $\hat{e}_{rb} = \theta \hat{e}_{rb}^0$</td><td>▷ Update the $\hat{e}_{rb}$</td></tr></table>

In this paper, we set $\theta > 1$, so that the force residual threshold changes from a relatively smaller value towards a bigger value as the relative error $e_c$ reduces. More specially, the adaptive selection of $\hat{e}_{rb}$ is roughly divided into three stages: In the first stage, $e_c$ is relatively higher and so the FOM rather than ROM is called on to maintain the stability of the optimization procedure. In the second stage, the variation of the objective function is less, so we begin using the ROM with a smaller value of $\hat{e}_{rb}$. Finally in the third stage, $e_c$ is small enough (meaning that objective is sufficiently converged) therefore a larger value $\hat{e}_{rb}$ of can be used, thus the ROM solution is more readily accepted. Note that the code in Algorithm 3 may be directly inserted between lines 7 and 8 in Algorithm 2 for obtaining a new version of ROM using adaptive parameter selection for the residual threshold.

## Bi-ROM

Recall that in “Introduction”, we have mentioned that stress-constrained topology optimization problems are generally non-self-adjoint, and will yield an additional adjoint equation (see (Eq. (22)) derived from sensitivity analysis of the stress constraint function, that needs to be solved. Taking advantage of the fact that the adjoint equation is of the same scale as the FE equilibrium equations, it would seem logical to apply the ROM idea given in Algorithm 2 to the adjoint equations, giving:

$$ \lambda = \Xi\mu + \tilde{\lambda}, \qquad (36) $$

where, $\Xi$ is the reduced basis, $\mu$ represents the projection coefficients and $\tilde{\lambda}$ is the mean snapshot for the set of adjoint solutions. An important feature here is that the construction of the ROM for the adjoint equation and the equilibrium equations in this paper are completely independent, which means the snapshot number and residual threshold of the ROM for the adjoint equations could be completely different from those for the equilibrium equations.

## Numerical tests

In this section, we first apply the two ROM methods that were introduced in “Projection-based ROM approaches”: SF-ROM and MF-ROM to a 2D benchmark stress-constrained compliance minimization problem. Next, using MF-ROM, we demonstrate the effectiveness of the two improved strategies that were presented in “Two improved strategies”. In order to better understand the impact of the parameters associated with
the stress constraints on optimal design, we use different values for the p norm factor $p$, the stress penalty factor $q$ and the allowable maximum stress $\hat{\sigma}$ for numerical experiments using a benchmark 3D L-shaped beam. Unless otherwise stated, the variables and geometric parameters are dimensionless for all test examples in this paper. Pre-defined parameters: Young's modulus $E_0 = 1$ for the solid material with a lower limit $E_{min} = 10^{-9}$ for void material, Poisson's ratio $\nu = 0.3$. The penalty factor of stiffness $\omega = 3$, density filter $r_{min} = 2.5$ and the allowable volume fraction $\dot{v}_f = 0.2$. Since the objective function (here compliance) is conveniently represented using a separable convex approximation by the Method of Moving Asymptotes (MMA), which is easily adapted to different types of topology optimization problems, we retain the MMA as the optimizer for iteration loop. The optimization iterations stop either when the convergence condition is satisfied or once the maximum allowable number of iterations is reached. Here, the convergence criterion is defined as follows:

$$ \frac{\|\rho_{i-1} - \rho_i\|_2}{\|\rho_i\|_2} \le Tol, \qquad (37) $$

where, $\|\cdot\|_2$ denotes the binary norm of a vector. $i$ denotes the current iteration step, and $Tol$ is the allowable convergence error, in this paper, $Tol = 10^{-4}$.

## 2D L-shaped beam

Consider the classical optimization problem of the 2D L-shaped beam [4, 5, 8, 10], the design domain and boundary conditions are illustrated in Fig. 1. The structure's top face on the left side is fixed and an external load $f = 1$ is applied on the top right corner A. We use 16384 bilinear quadrilateral elements to discretize the design domain, each element having unit volume and thickness. To avoid stress concentration at the loading node A, the concentrated load is evenly distributed over the 6 adjacent nodes between node A and B, as shown in Fig. 1. For this 2D case, we set the p norm factor $p = 8$, the stress penalty factor $q = 0.8$ and the allowable maximum stress. The maximum iteration number is set to 200. Note that for all test cases shown hereafter, the iterations were halted due to the maximum allowable number of optimization loop iterations having been attained.

Table 1 shows the optimal topology and corresponding stress distribution obtained using the reference method (without ROM), as well as the SF-ROM method and the MF-ROM method using a fixed snapshot number $N_b$ and various values for the force residual threshold $\hat{e}_{rb}$. Figures 2 and 3 compare the iteration histories for the convergence of the compliance and stress (von Mises stress $\sigma_{max}^{vM}$ and P-norm stress $\sigma^{PN}$) between the SF-ROM and MF-ROM methods.

From Table 1 a macroscopic view leads us to the premature conclusion that the optimal topology using SF-ROM and MF-ROM seem basically consistent with that using the reference method; but on closer inspection more structural branches occur and the arc of reentry of the L-shaped beam is more distinct when using a bigger value of $\hat{e}_{rb}$. And it should be noted that for both SF-ROM and MF-ROM, stress concentration still exists at the location of the applied force. Figure 2 shows that the objective optimization is converged regardless of whether or not we use SF-ROM or MF-ROM. Figure 3 shows that stress constraints using both the SF-ROM and MF-ROM are active and gradually satisfied during the overall optimization procedure. Furthermore, due to the constraint violation shown in Fig. 3, the convergence history of the objective function is non-monotonous as shown in Fig. 2.
<table><caption>**Table 1** Comparison of optimal topology and stress distribution for 2D L-shaped beam obtained using reference method (No ROM) and SF-ROM and MF-ROM under various $\hat{e}_{rb} = 0.1, 0.5$ and fixed $N_b = 4$</caption><thead><tr><th>Method</th><th>Optimal topology</th><th>Stress distribution</th></tr></thead><tbody><tr><td>Reference</td><td></td><td></td></tr><tr><td>SF-ROM N<sub>b</sub> = 4, e<sub>rb</sub> = 0.1</td><td></td><td></td></tr><tr><td>SF-ROM N<sub>b</sub> = 4, e<sub>rb</sub> = 0.5</td><td></td><td></td></tr><tr><td>MF-ROM N<sub>b</sub> = 4, e<sub>rb</sub> = 0.1</td><td></td><td></td></tr><tr><td>MF-ROM N<sub>b</sub> = 4, e<sub>rb</sub> = 0.5</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig. 1 Design domain and boundary conditions of L-shaped beam</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 2 Comparison of compliance iterative history between SF-ROM and MF-ROM usingand</figcaption></figure>

The above design results are summarized in Table 2, where the relative error of the optimal compliance (c*) characterizes the approximate accuracy of the objective function value using ROM,

$$ error = \frac{|c_{refe}^* - c_{rom}^*|}{|c_{refe}^*|}, \qquad (38) $$

where, $c_{refe}^*$ and $c_{rom}^*$ represent the optimal compliance obtained by using the reference method and the ROM method respectively.

From Table 2, as we increase the value of $\hat{e}_{rb}$, the CPU clock time using both SF-ROM and MF-ROM show a significant decrease compared to the reference method, which is easy to explain: the cheaper reduced-basis solution (ROM) is called more frequently, thus the computational efficiency of the overall optimization routine is improved. However, the optimal compliance using both ROM methods has a very poor accuracy if we use a larger value of the force residual threshold ($\hat{e}_{rb} = 0.5$). This indicates that we must make
<figure><img src="image_7.png" /><figcaption>Fig. 3 Comparison of stress iterative history between SF-ROM (a) and MF-ROM (b)</figcaption></figure>

<table><caption>**Table 2** Comparison of design results between SF-ROM and MF-ROM (after 200 iterations)</caption><thead><tr><th>Method</th><th>N<sub>b</sub></th><th>ê<sub>rb</sub></th><th>Calls to ROM</th><th>CPU running time (/s)</th><th>Optimal compliance (c*)</th><th>Relative error of c* (%)</th></tr></thead><tbody><tr><td>Reference</td><td>–</td><td>–</td><td>0/200</td><td>1884.26</td><td>291.48</td><td>–</td></tr><tr><td rowspan="2">SF-ROM</td><td rowspan="2">4</td><td>0.1</td><td>135/200</td><td>1624.22</td><td>292.3</td><td>0.28</td></tr><tr><td>0.5</td><td>164/200</td><td>1355.04</td><td>301.5</td><td>3.44</td></tr><tr><td rowspan="2">MF-ROM</td><td rowspan="2">4</td><td>0.1</td><td>129/200</td><td>1541.80</td><td>294.19</td><td>0.93</td></tr><tr><td>0.5</td><td>164/200</td><td>1296.32</td><td>300.11</td><td>2.96</td></tr></tbody></table>

<table><caption>**Table 3** Comparison of design results between MF-ROM and Adap-ROM (after 200 iterations)</caption><thead><tr><th>Method</th><th>N<sub>b</sub></th><th>ê<sub>rb</sub></th><th>Calls to ROM</th><th>CPU running time (/s)</th><th>Optimal compliance (c*)</th><th>Relative error of c* (%)</th></tr></thead><tbody><tr><td>Reference</td><td>–</td><td>–</td><td>0/200</td><td>1884.26</td><td>291.48</td><td>–</td></tr><tr><td rowspan="2">MF-ROM</td><td>4</td><td>0.05</td><td>105/200</td><td>1677.06</td><td>291.10</td><td>0.13</td></tr><tr><td>4</td><td>0.5</td><td>164/200</td><td>1296.32</td><td>300.11</td><td>2.96</td></tr><tr><td>Adap-ROM</td><td>4</td><td>0 → 0.05 → 0.5</td><td>116/200</td><td>1652.98</td><td>291.68</td><td>0.07</td></tr></tbody></table>

a trade-off between efficiency and accuracy. Furthermore, by comparing the performance of the MF-ROM with the SF-ROM, we find that MF-ROM can further drive down the computational cost with a minor loss of objective accuracy when using the same set of values of $N_b$ and $ê_{rb}$.

Next, we use the first improved strategy: adaptive parameter selection of the force residual threshold $ê_{rb}$ to modify the MF-ROM method in “Projection-based ROM approaches”. Some parameters need to be predefined: initial force residual threshold $ê_{rb}^0 = 0.05$ and its scale factor $\theta = 10$, initial limit of objective relative error (see Eq. (32)) $ê_{rb}^0 = 2 \times 10^{-3}$ and its decay factor $\zeta = 10^{-1}$. The optimized design results are summarized in Table 3 and Fig. 4 shows the optimal topology using the MF-ROM and Adap-ROM as well as the reference method.

From Table 3, we see that topology optimization using Adap-ROM has a clear advantage in terms of accuracy compared to the MF-ROM method using a bigger force residual ($ê_{rb} = 0.5$), while requiring less computational time compared to the reference method, as well as the MF-ROM method using a smaller $ê_{rb} = 0.05$ (as expected), which shows that this improved method allows a good compromise between accuracy and speed. From Fig. 4 we see that the optimal topology obtained using Adap-ROM is closer to that obtained using the reference method, compared to the optimal topology obtained using MF-ROM,
<figure><img src="image_8.png" /><figcaption>Fig. 4 Comparison of optimal topology</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5 Change history of using Adap-ROM</figcaption></figure>

<table><caption>Table 4 Comparison of design results between MF-ROM and Bi-ROM (after 200 iterations)</caption><thead><tr><th>Method</th><th>N<sub>b</sub></th><th>ê<sub>rb</sub></th><th>Calls to ROM</th><th>CPU running time (/s)</th><th>Optimal compliance (c*)</th><th>Relative error of c* (%)</th></tr></thead><tbody><tr><td>Reference</td><td>–</td><td>–</td><td>0/200</td><td>1884.26</td><td>291.48</td><td>–</td></tr><tr><td>MF-ROM</td><td>4</td><td>0.1</td><td>129/200</td><td>1541.80</td><td>294.19</td><td>0.93</td></tr><tr><td rowspan="2">Bi-ROM</td><td>4</td><td>0.1</td><td>132/200</td><td>1338.92</td><td>293.88</td><td>0.82</td></tr><tr><td>6</td><td>0.05</td><td>148/200</td><td>1338.92</td><td>293.88</td><td>0.82</td></tr></tbody></table>

which is again expected. Next we give the change history of relative error $e_c$ generated by using Adap-ROM in Fig. 5a, which allows us to decide when to use different force residual thresholds. And Fig. 5b gives the evolutionary curve of the force residual $e_{rb}$.

Then, for the same 2D benchmark stress-constrained topology optimization problem, we use another improved strategy Bi-ROM to modify the MF-ROM method in “Projection-based ROM approaches”. For the equilibrium equations (Eq. (23)), we set $N_b = 4$ and $e_{rb} = 0.1$, for adjoint equations (Eq. (21)), we set $N_b = 6$ and $e_{rb} = 0.05$. The design results are summarized in Table 4 and Fig. 6 shows the optimal topology obtained using MF-ROM and Bi-ROM as well as reference method.

From Table 4 we can see that Bi-ROM is less time consuming compared to the MF-ROM and the accuracy of the optimized compliance is better. Figure 6 shows that the Bi-ROM appears to yield a clearer optimal topology than the MF-ROM method.

Figure 7 shows the iteration history for the force residual $e_{rb}$ of the equilibrium equations as well as the adjoint equations, using the Bi-ROM method.
<figure><img src="image_10.png" /><figcaption>Fig. 6 Comparison of optimal topology</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7 Evolutionary curve of force residual $e_{rb}$ using Bi-ROM method for equilibrium equations(left) and adjoint equations(right)</figcaption></figure>

A conclusion may be directly drawn from Fig. 8 that for the Bi-ROM, the ROM solution of the adjoint equations is more readily “acceptable” than that of the equilibrium equations, even for smaller values of $e_{rb}$. This is because the residual evolutionary curve of equilibrium equations frequently oscillated above and below the threshold line (red) while the residual evolutionary curve of equilibrium equations basically stays below the threshold line as the iteration progresses.

Finally, we combine the Adap-ROM and Bi-ROM methods together (denoted by AdapBi-ROM) to further improve the MF-ROM method. Some parameters are pre-defined like before: initial limit of objective function relative error $\hat{e}_c^0 = 2 \times 10^{-3}$ decay factor $\zeta = 10^{-1}$, for the equilibrium equations: $N_b = 4$, initial force residual threshold $\hat{e}_c^0 = 0.05$ with scale factor $\theta = 10$, for adjoint equations $N_b = 6$, initial force residual threshold $\hat{e}_c^0 = 0.01$ and its scale factor $\theta = 10$. The optimized design results are summarized in Table 5 as follows.

From Table 5 we can see, the combined strategy AdapBi-ROM shows better performance than any one of the individual strategies(Adap-ROM or Bi-ROM) from a holistic perspective since this combined method can not only significantly reduce the computation time but also ensure a high accuracy of the value of objective function.
<table><caption>**Table 5** Comparison of design results among AdapBi-ROM, Adap-ROM and Bi-ROM (after 200 iterations)</caption><thead><tr><th>Method</th><th>N<sub>b</sub></th><th>ê<sub>rb</sub></th><th>Calls to ROM</th><th>CPU running time (/s)</th><th>Optimal compliance (c*)</th><th>Relative error of c* (%)</th></tr></thead><tbody><tr><td>reference</td><td>–</td><td>–</td><td>0/200</td><td>1884.26</td><td>291.48</td><td>–</td></tr><tr><td>Adap-ROM</td><td>4</td><td>0 → 0.05 → 0.5</td><td>116/200</td><td>1652.98</td><td>291.68</td><td>0.07</td></tr><tr><td>Bi-ROM</td><td>4</td><td>0.5</td><td>167/200</td><td>1294.97</td><td>299.09</td><td>2.61</td></tr><tr><td>AdapBi-ROM</td><td>4</td><td>0 → 0.05 → 0.5</td><td>112/200</td><td>1294.97</td><td>299.09</td><td>2.61</td></tr><tr><td></td><td>6</td><td>0.1</td><td>153/200</td><td>1352.96</td><td>291.52</td><td>0.01</td></tr><tr><td></td><td>6</td><td>0 → 0.01 → 0.1</td><td>96/200</td><td>1352.96</td><td>291.52</td><td>0.01</td></tr></tbody></table>

## Discussion

The multi-fidelity reduced-order-model (MF-ROM) approach clearly enhances the performance of the topology optimization procedure with a superior reduction of computational effort and CPU time compared with single-fidelity reduced-order-model (SF-ROM) for all the test cases studied, noting that a small loss of accuracy of the objective function is expected when using MF-ROM. By using the first proposed improved strategy (Adap-ROM) of adaptively selecting the force residual threshold ($ê_{rb}$) in MF-ROM, we found a trade-off can be made between optimization efficiency and objective accuracy. We also note that the second proposed improved strategy (Bi-ROM) of applying MF-ROM method to both adjoint equations derived from sensitivity analysis of stress constraint function and equilibrium equations in FE analysis can reduce the CPU running time in a more significant manner. In addition, even better accuracy of objective function can be obtained compared with the original MF-ROM method (for this 2D case at least), this may be because using a ROM for the adjoint equations corrects the convergence of the objective function. Finally, we combined both strategies and this can obtain better optimized results because there is a significant improvement in efficiency with a higher accuracy of approximation for the objective function.

## 3D L-shaped beam

In the previous section, we investigated the application of reanalysis using ROM in stress-constrained topology optimization through a 2D numerical example. In this section, we apply the previously described ROM methods to the 3D compliance minimization problem of the L-shaped beam. To better understand the influence of the stress constraint on the optimization results, we then consider making an in depth analysis of some of the parameters associated with the stress constraint based on the AdapBi-ROM method. The L-shaped beam is shown in Fig. 8, where the structure's top surface is fixed and a vertically downward load $f = 1$ is applied along the top right side.

We use standard hexahedral elements to discretize the design domain with each element having unit volume and thickness. The maximum number of iterations is set to 300.

## ROM methods on 3D L-shaped beam

Here, we set $L = 80$ and $W = 50$, thus the design domain is discretized by a total of 32768 elements. 5 ROM methods were tested for this 3D case, The corresponding design results are given in Table 6.
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Design domain and boundary conditions of 3D L-shaped beam</figcaption></figure>

<table><caption>**Table 6** Comparison of design results using several ROM methods (after 200 iterations)</caption><thead><tr><th>Method</th><th>N<sub>b</sub></th><th>ê<sub>rb</sub></th><th>Calls to ROM</th><th>CPU running time (/s)</th><th>Optimal compliance (c*)</th><th>Relative error of c*(%)</th></tr></thead><tbody><tr><td>Reference</td><td>-</td><td>-</td><td>0/300</td><td>23406.28</td><td>45.02</td><td>-</td></tr><tr><td>SF-ROM</td><td>4</td><td>0.05</td><td>224/300</td><td>21388.52</td><td>45.02</td><td>0</td></tr><tr><td>MF-ROM</td><td>4</td><td>0.05</td><td>225/300</td><td>19682.72</td><td>45.02</td><td>0</td></tr><tr><td>Adap-ROM</td><td>4</td><td>0.5</td><td>260/300</td><td>19578.27</td><td>45.07</td><td>0.11</td></tr><tr><td>Bi-ROM</td><td>4</td><td>0 → 0.05 → 0.5</td><td>233/300</td><td>19655.47</td><td>45.01</td><td>0.02</td></tr><tr><td>AdapBi-ROM</td><td>4</td><td>0.5</td><td>260/300</td><td>18663.07</td><td>45.09</td><td>0.16</td></tr><tr><td>AdapBi-ROM</td><td>6</td><td>0.1</td><td>281/300</td><td>18663.07</td><td>45.09</td><td>0.16</td></tr><tr><td>AdapBi-ROM</td><td>4</td><td>0 → 0.05 → 0.5</td><td>233/300</td><td>18896.49</td><td>45.00</td><td>0.04</td></tr><tr><td>AdapBi-ROM</td><td>6</td><td>0 → 0.01 → 0.1</td><td>254/300</td><td>18896.49</td><td>45.00</td><td>0.04</td></tr></tbody></table>

From Table 6 we see that the combined AdapBi-ROM method behaves well both as far as computational time saving as well as final objective accuracy are concerned, as expected, in the 3D large-scale test case. Next, we will use this method as the ROM approach of choice to conduct a parameter study.

## Parametric study

For this, we set $L = 50$ and $W = 6$ in Fig. 8, thus the design domain is discretized by a total of 9600 elements. Some parameters are predefined as follows: initial limit of objective relative error $\hat{e}_c^0 = 5 \times 10^{-3}$ and its decay factor $\zeta = 10^{-2}$, for equilibrium equations $N_b = 4$; initial force residual threshold $\hat{e}_c^0 = 0.04$ and its scale factor $\theta = 20$, for adjoint equations $N_b = 6$; initial force residual threshold $\hat{e}_c^0 = 0.01$ and its scale factor $\theta = 50$.
<table><caption>**Table 7** Optimization results of AdapBi-ROM and reference method using varying *p* (after 300 iterations)</caption><thead><tr><th rowspan="2">p</th><th colspan="2">Calls to ROM</th><th colspan="2">CPU running time (/s)</th><th colspan="2">Optimal compliance (c*)</th></tr><tr><th>Reference</th><th>AdapBi-ROM</th><th>Reference</th><th>AdapBi-ROM</th><th>Reference</th><th>AdapBi-ROM</th></tr></thead><tbody><tr><td rowspan="2">2</td><td>0</td><td>120/300</td><td rowspan="2">2269.06</td><td rowspan="2">2097.43</td><td rowspan="2">80.04</td><td rowspan="2">79.92</td></tr><tr><td>0</td><td>34/300</td></tr><tr><td rowspan="2">4</td><td>0</td><td>192/300</td><td rowspan="2">2401.96</td><td rowspan="2">2020.23</td><td rowspan="2">78.73</td><td rowspan="2">78.79</td></tr><tr><td>0</td><td>85/300</td></tr><tr><td rowspan="2">6</td><td>0</td><td>207/300</td><td rowspan="2">2354.69</td><td rowspan="2">1768.36</td><td rowspan="2">78.81</td><td rowspan="2">78.86</td></tr><tr><td>0</td><td>141/300</td></tr><tr><td rowspan="2">8</td><td>0</td><td>202/300</td><td rowspan="2">2277.97</td><td rowspan="2">1953.88</td><td rowspan="2">78.38</td><td rowspan="2">78.63</td></tr><tr><td>0</td><td>130/300</td></tr><tr><td rowspan="2">12</td><td>0</td><td>201/300</td><td rowspan="2">2209.22</td><td rowspan="2">1871.71</td><td rowspan="2">79.59</td><td rowspan="2">78.45</td></tr><tr><td>0</td><td>147/300</td></tr><tr><td rowspan="2">20</td><td>0</td><td>227/300</td><td rowspan="2">2159.00</td><td rowspan="2">1560.70</td><td rowspan="2">80.33</td><td rowspan="2">80.30</td></tr><tr><td>0</td><td>255/300</td></tr><tr><td rowspan="2">30</td><td>0</td><td>231/300</td><td rowspan="2">2180.82</td><td rowspan="2">1558.75</td><td rowspan="2">80.30</td><td rowspan="2">80.29</td></tr><tr><td>0</td><td>252/300</td></tr></tbody></table>

Influence of *p* norm factor *p*

Here, we vary the *p* norm factor from 2 to 30, and while fixing the stress penalty factor at *q* = 0.8 and the allowable maximum stress at $\hat{\sigma}$. Table 7 gives the design results.

It can be seen from Table 7 that the CPU runtime trend for both the reference method as well as the combined AdapBi-ROM method is not monotonic. More computation time is spent when using a smaller *p* value (like 2 and 4) while less time is spent when using a larger *p* value (like 20 and 30) for these two methods. However, both methods still have some differences in terms of the changing trend of the CPU runtime, where the CPU time of the reference method goes down while that of the AdapBi-ROM method goes up, when we vary the *p* value from 6 to 8.

Interestingly, we find that the optimal compliance obtained using both methods has the same trend of first dropping and then increasing as the value of *p* increases, which means that we can likely obtain the smallest optimal compliance (meaning the “stiffest” possible structure) by using a mid-range value for *p*. In general, to better balance CPU time and structural performance, the *p* value should neither be too large nor too small.

Figures 9, 10 and 11 show the optimal topology and stress iterative history as well as the volume fraction iteration history for different values of the P-norm *p* norm factor with the AdapBi-ROM method.

We see from Fig. 9 that the optimal topologies obtained using different values of *p* based on AdapBi-ROM method are visually different when using *p* = 4 and *p* = 6 (or *p* = 12 and *p* = 20).

From Fig. 10, we can see that the stress value during the entire procedure satisfy its constraint limit progressively faster with a larger value of *p*. Meanwhile, it seems that the larger the value of *p*, the smaller the value of the final maximum von Mises stress $\sigma_{max}^{VM}$. For *p* = 20 and *p* = 30, the curves are below the stress limit within 50 iterations.

Fig. 11 shows that the volume fraction during the iterations can rapidly, and in a stable manner, satisfy its respective constraint limiting value, provided *p* > 2.
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Comparison of optimal topology using different value of *p* based on AdapBi-ROM method</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Comparison of stress iterative history using different value of *p* based on AdapBi-ROM method</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 11 Comparison of volume fraction iterative history using different value of p based on AdapBi-ROM method</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12 Comparison of compliance iterative history using different values of p based on AdapBi-ROM method</figcaption></figure>

Next, we give a comparison of the compliance (objective function) iteration history using different values of p based on AdapBi-ROM method in Fig. 12 where it appears that the larger the p value, the faster the rate of convergence of the optimization.

### Investigation of varying stress penalty factor *q*

We vary the stress penalty factor *q* while fixing the p norm factor *p* = 8 and the allowable maximum stress at. Table 8 gives the design results.

It can be seen from Table 8 that the trends for CPU time and optimal compliance for both the reference method and the AdapBi-ROM method are exactly the same and not monotonic as the *q* value decreases. Therefore, the ROM retains the same impact of the stress penalty factor parameter as for the reference case. Interestingly, the CPU time has a general tendency to first increase and then decrease while the optimal compliance has the opposite tendency, implying a potential conflict between CPU time and optimal
<table><caption>**Table 8** Optimization results of AdapBi-ROM method and reference method using varying *q* (after 300 iterations)</caption><thead><tr><th rowspan="2">q</th><th colspan="2">Calls to ROM</th><th colspan="2">CPU running time (/s)</th><th colspan="2">Optimal compliance (c*)</th></tr><tr><th>Reference</th><th>AdapBi-ROM</th><th>Reference</th><th>AdapBi-ROM</th><th>Reference</th><th>AdapBi-ROM</th></tr></thead><tbody><tr><td rowspan="2">0.9</td><td>0</td><td>228/300</td><td rowspan="2">2150.16</td><td rowspan="2">1584.64</td><td rowspan="2">80.00</td><td rowspan="2">79.98</td></tr><tr><td>0</td><td>236/300</td></tr><tr><td rowspan="2">0.8</td><td>0</td><td>202/300</td><td rowspan="2">2277.97</td><td rowspan="2">1953.88</td><td rowspan="2">78.38</td><td rowspan="2">78.63</td></tr><tr><td>0</td><td>130/300</td></tr><tr><td rowspan="2">0.5</td><td>0</td><td>146/300</td><td rowspan="2">2180.45</td><td rowspan="2">1825.66</td><td rowspan="2">79.90</td><td rowspan="2">79.74</td></tr><tr><td>0</td><td>151/300</td></tr><tr><td rowspan="2">0.3</td><td>0</td><td>198/300</td><td rowspan="2">2232.93</td><td rowspan="2">1884.87</td><td rowspan="2">80.76</td><td rowspan="2">80.76</td></tr><tr><td>0</td><td>128/300</td></tr><tr><td rowspan="2">0.1</td><td>0</td><td>226/300</td><td rowspan="2">2197.40</td><td rowspan="2">1698.84</td><td rowspan="2">80.26</td><td rowspan="2">80.27</td></tr><tr><td>0</td><td>172/300</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 13** Comparison of optimal topology using different value of *q* based on AdapBi-ROM method</figcaption></figure>

compliance for several *q* values. This means that *q* needs to be properly determined to better balance the computational efforts and structural performance.

Just like in “Influence of p norm factor p”, we give the optimal topology and stress iteration history as well as the volume fraction iteration history for different values of the stress penalty factor, using the AdapBi-ROM method, as shown in Figs. 13, 14 and 15, respectively.

From Fig. 13 we can see that the optimal topologies using different values of *q* based on AdapBi-ROM method are visually distinguishable from each other, which can be seen from the test case with *q* = 0.5 and *q* = 0.8.

From Fig. 14 we note that the stress constraint is not satisfied when using relatively smaller values of *q* (0.1 – 0.3) which means the stress penalization is not sufficient. In
<figure><img src="image_18.png" /><figcaption>Fig. 14 Comparison of stress iterative history using different value of *q* based on AdapBi-ROM method</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15 Comparison of volume fraction iterative history using different value of *q* based on AdapBi-ROM method</figcaption></figure>

addition, using a larger value of *q* can cause the maximum von Mises stress to satisfy its constraint limiting value as soon as possible during the procedure. In the case of *q* = 0.9, the curves end up below the stress limit after about 25 iterations.

Figure 15 shows that the volume fraction during the iterations can rapidly and in a stable manner satisfy the constraint for all 5 different values of *q* tested here.

Next, we compare the compliance iteration histories using different values of *q* based on the AdapBi-ROM method in Fig. 16 where different convergence rates of optimization can be obtained by varying *q*, but no monotonicity is observed as *q* increases. It seems plausible that a midrange value of *q* would lead to a relative lower convergence rate as shown below.

### Varying the allowable maximum stress

In this case, we vary the allowable maximum stress $\hat{\sigma}$ while simultaneously fixing the *p* norm factor *p* at 8 and the stress penalty factor *q* at 0.8. Table 9 gives the design results.
<figure><img src="image_20.png" /><figcaption>**Fig. 16** Comparison of compliance iterative history using different values of *q* based on AdapBi-ROM method</figcaption></figure>

<table><caption>**Table 9** Optimization results of AdapBi-ROM method and reference method using varying $\hat{\sigma}$ (after 300 iterations)</caption><thead><tr><th rowspan="2">$\hat{\sigma}$</th><th colspan="2">Calls to ROM</th><th colspan="2">CPU running time (/s)</th><th colspan="2">Optimal compliance (c*)</th></tr><tr><th>Reference</th><th>AdapBi-ROM</th><th>Reference</th><th>AdapBi-ROM</th><th>Reference</th><th>AdapBi-ROM</th></tr></thead><tbody><tr><td rowspan="2">0.15</td><td>0</td><td>112/300</td><td rowspan="2">2413.22</td><td rowspan="2">2237.27</td><td rowspan="2">65.20</td><td rowspan="2">65.09</td></tr><tr><td>0</td><td>35/300</td></tr><tr><td rowspan="2">0.3</td><td>0</td><td>96/300</td><td rowspan="2">2408.48</td><td rowspan="2">2250.57</td><td rowspan="2">96.30</td><td rowspan="2">98.64</td></tr><tr><td>0</td><td>26/300</td></tr><tr><td rowspan="2">0.45</td><td>0</td><td>136/300</td><td rowspan="2">2664.67</td><td rowspan="2">2332.58</td><td rowspan="2">78.58</td><td rowspan="2">78.88</td></tr><tr><td>0</td><td>103/300</td></tr><tr><td rowspan="2">0.65</td><td>0</td><td>229/300</td><td rowspan="2">2240.31</td><td rowspan="2">1566.90</td><td rowspan="2">80.18</td><td rowspan="2">80.18</td></tr><tr><td>0</td><td>251/300</td></tr><tr><td rowspan="2">0.75</td><td>0</td><td>227/300</td><td rowspan="2">2171.39</td><td rowspan="2">1573.41</td><td rowspan="2">80.00</td><td rowspan="2">80.00</td></tr><tr><td>0</td><td>248/300</td></tr></tbody></table>

From Table 9 we see that the trends for the CPU time are the same for both the reference as well as the AdapBi-ROM methods, first increasing and then decreasing as the $\hat{\sigma}$ value increases. A similar trend is seen for the optimal compliance for both methods.

Just like in the previous investigations, we give the optimal topology and stress iteration history as well as the volume fraction iteration history for different values of the allowable maximum stress based on the AdapBi-ROM method, as shown in Figs. 17, 18 and 19, respectively.

From Fig. 17 we can see that the optimal topologies using different values of $\hat{\sigma}$ based on the AdapBi-ROM method are clearly very different from each other. When using smaller values of $\hat{\sigma}$ ($\hat{\sigma}$ = 0.15, $\hat{\sigma}$ = 0.3), the final design structural branches are discontinuous and in a “fractured” state which means that the optimized design has not yet converged sufficiently despite having completed maximum iterations. Interestingly, we find a clear curvature forming at the corner of the L-shaped beam by using these two smaller values
<figure><img src="image_21.png" /><figcaption>**Fig. 17** Comparison of optimal topology using different value of $\hat{\sigma}$ based on AdapBi-ROM method</figcaption></figure>

of $\hat{\sigma}$. As for other values of $\hat{\sigma}$, the optimal topologies are also different from each other (see the case of $\hat{\sigma} = 0.45$ and $\hat{\sigma} = 0.65$).

From Fig. 18 we can see that the stress constraint is not satisfied in the case of the smallest value of $\hat{\sigma}$ ($\hat{\sigma} = 0.15$). Also the maximum von Mises stress constraint is increasingly easier to satisfy during the procedure as $\hat{\sigma}$ increases, but we note that an overtly large value of $\hat{\sigma}$ could potentially render the stress constraint inactive (as seen in the case of $\hat{\sigma} = 0.65$ and $\hat{\sigma} = 0.75$).

From Fig. 19 we can see that the volume constraint is not satisfied for smaller values of $\hat{\sigma} = 0.15$ or $\hat{\sigma} = 0.3$ despite having been satisfied for larger values of $\hat{\sigma}$.

Finally, we compare the compliance iteration history using different values of $\hat{\sigma}$ and the AdapBi-ROM method in Fig. 20, where the compliance clearly has different convergence rate depending on the values of $\hat{\sigma}$ and a relatively large value of $\hat{\sigma}$ could potentially accelerate convergence.

## Conclusions and perspectives

In this paper, we first presented two novel ROM-based Approximate Reanalysis approaches: on-the-fly ROM using Principal Components Analysis (SF-ROM) and multi-fidelity ROM using multigrid method (MF-ROM) and applied them to the non-self-adjoint stress-constrained topology optimization problem, and assessed their performance in terms of efficiency and reduction in computation time, using a 2D benchmark problem. Furthermore, we improved the second approach with two improved strategies: a so-called Adap-ROM method where the residual threshold of the equilibrium equations is adaptively selected according to the relative error of the objective function value during the iteration; and the so-called Bi-ROM method where the ROM is simultaneously applied
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 18** Comparison of stress iterative history using different value of $\hat{\sigma}$ based on AdapBi-ROM method</figcaption></figure>

to both the equilibrium equations as well as the adjoint equations obtained from the sensitivity analysis of the stress constraint. The first improved strategy achieves a balance between efficiency and accuracy, while the second improved strategy can further reduce the computational effort needed. To amplify the effect of the ROM, we combined both strategies for this 2D numerical example showing a good balance between efficiency and accuracy.

We then tested various versions of the ROM methods in this paper through a 3D L-shaped beam. Then using the combined version of Adap-ROM and Bi-ROM, we studied various parameters associated with the stress constraint, including the P-norm factor, the
<figure><img src="image_27.png" /><figcaption>**Fig. 19** Comparison of volume fraction iterative history using different value of $\hat{\sigma}$ based on AdapBi-ROM method</figcaption></figure>

stress penalty factor as well as the allowable maximum stress value, in order to better understand the influence of the stress constraint on the ROM-based optimization results. We showed that the aforementioned parameters have a significant effect on the final optimized results (efficiency, structural performance etc), and therefore a judicious selection of these key parameters is required.

Furthermore, from the point of view of storage, the amount of storage required for the proposed on-the-fly ROM approach is much smaller compared to that needed for the full-field problem with finite element simulation calculation during the optimization
<figure><img src="image_28.png" /><figcaption>Fig. 20 Comparison of compliance iterative history using different values of $\hat{\sigma}$ based on AdapBi-ROM method</figcaption></figure>

process. The proposed method only needs to store a small dimension snapshot matrix of small dimensionality to obtain the reduced basis. The only thing that needs to be done sequentially is the on-the-fly updating of the reduced basis using the residual errors, in order to ensure the quality of approximation.

A logical extension of this work would focus on efficient large-scale dynamic topology optimization. In addition, other non-self-adjoint problems besides the numerical examples presented in this paper need to be addressed.

**Acknowledgements**

The first author expresses her sincere gratitude to the associate professor Kai Long at North China Electric Power University for his help in testing the codes in its initial stages.

**Author contributions**

MX, WZ developed the main idea of the two novel approaches and introduced JM, and DL to implement benchmark test case, analyze numerical results and prepare the concept manuscript. BR provided his feedback on the manuscript. All authors read and approved the final manuscript.

**Authors' information**

Dr. Manyu Xiao, Associate professor, School of Mathematics and Statistics of Northwestern Polytechnical University, Xi'an, Shaanxi, People's Republic of China. Mr. Jun Ma, Dongcheng Lu, Master students, School of Mathematics and Statistics of Northwestern Polytechnical University, Xi'an, Shaanxi, People's Republic of China. Dr. Balaji Raghavan, Associate professor, Laboratoire de Génie Civil et Génie Mécanique EA 3913, Institut National des Sciences Appliquées de Rennes, Rennes, France. Prof. Weihong Zhang, Academician of the Chinese Academy of Sciences, State IJR Center of Aerospace Design and Additive Manufacturing, Northwestern Polytechnical University, Xi'an, 710072, People's Republic of China.

**Funding**

This research study was supported by the National Natural Science Foundation of China (Grant No. 12032018, Grant No. 11620101002 and Grant No. 11972166) and the Natural Science Foundation of Shaanxi Province, China (Grant No. 2021JM-043).

**Availability of data materials**

The data sets and models used and/or analysed during the current study are available from the corresponding author on reasonable request.

**Declarations**

**Competing interests**

The authors declare that they have no competing interests.
References

1. Sigmund O. A 99 line topology optimization code written in matlab. Structural and Multidisciplinary Optimization. 2001;21(2):120-7.

2. Andreassen E, Clausen A, Schevenels M, Sigmund O. Efficient topology optimization in matlab using 88 lines of code. Structural and Multidisciplinary Optimization. 2011;43(1):1-16.

3. Rietz A. Sufficiency of a finite exponent in simp (power law) methods. Structural and Multidisciplinary Optimization. 2001;21(2):159-63.

4. E. Holmberg, T. Bo, and A. Klarbring. Stress constrained topology optimization. Springer-Verlag, 48(1), 2013.

5. Le C, Norato J, Bruns T. Stress-based topology optimization for continua. Structural and Multidisciplinary Optimization. 2010;41(4):605-20.

6. Cheng GD, Guo X. ε-relaxed approach in structural topology optimization. Structural optimization. 1997;13(4):258-66.

7. Cheng GD, Jiang Z. Study on topology optimization with stress constraints. Engineering Optimization. 1992;20(2):129-48.

8. P. Duysinx and M. P. Bendsøe. Topology optimization of continuum structures with local stress constraints. *International Journal for Numerical Methods in Engineering*, 1998.

9. P. Duysinx and O. Sigmund. *New developments in handling stress constraints in optimal material distribution*, pages 1501–1509. AIAA, 2008.

10. M. Bruggi and P. Venini. A mixed fem approach to stress constrained topology optimization. *International Journal for Numerical Methods in Engineering*, 73:1693–1714, 03 2008.

11. Bruggi M. On an alternative approach to stress constraints relaxation in topology optimization. Structural and Multidisciplinary Optimization. 2008;36(2):125-41.

12. Verbart A, Langelaar M, Keulen FV. A unified aggregation and relaxation approach for stress-constrained topology optimization. Structural & Multidisciplinary Optimization. 2017;55(2):1-17.

13. París J, Navarrina F, Colominas I, Casteleiro M. Topology optimization of continuum structures with local and global stress constraints. Structural and Multidisciplinary Optimization. 2009;39(4):419-37.

14. Yang RJ, Chen C. Stress-based topology optimization Structural optimization. 1996;12(2):98-105.

15. J. París, F. Navarrina, I. Colominas, and M. Casteleiro. Block aggregation of stress constraints in topology optimization of structures. *Advances in Engineering Software*, 41(3):433-441, 2010. Advances in optimum engineering design.

16. Yang DX, Liu HL, Zhang WS, et al. Stress-constrained topology optimization based on maximum stress measures. *Computers & Structures*. 2018;198:23-39.

17. J. H. Rong, T. T. Xiao, and L. H. Yu et al. Continuum structural topological optimizations with stress constraints based on an active constraint technique. *International Journal for Numerical Methods in Engineering*, 108(4), 2016.

18. P. Duysinx, L. Van Miegroet, E. Lemaire, and O. Brüls et al. Topology and generalized shape optimization: Why stress constraints are so important? https://doi.org/10.1051/ijsmdo/2008034, 2, 12 2008.

19. Wang MY, Li L. Shape equilibrium constraint: a strategy for stress-constrained structural topology optimization. Structural and Multidisciplinary Optimization. 2013;47(3):335-52.

20. Wang S, Sturler ED, Paulino GH. Large-scale topology optimization using preconditioned krylov subspace methods with recycling. International Journal for Numerical Methods in Engineering. 2007;69(12):2441-68.

21. Mitjana F, Cafieri S, Bugarin F, Segonds S, Castanie F, Duysinx P. Topological gradient in structural optimization under stress and buckling constraints. Applied Mathematics and Computation. 2021;409: 126032.

22. Collet M, Noël L, Bruggi M, Duysinx P. Topology optimization for microstructural design under stress constraints. Structural and Multidisciplinary Optimization. 2018;58(6):2677-95.

23. Noël L, Duysinx P. Shape optimization of microstructural designs subject to local stress constraints within an xfem-level set framework. Structural and Multidisciplinary Optimization. 2017;55(6):2323-38.

24. Amir O, Stolpe M, Sigmund O. Efficient use of iterative solvers in nested topology optimization. Structural and Multidisciplinary Optimization. 2010;42(1):55-72.

25. C. Gogu. Improving the efficiency of large scale topology optimization through on-the-fly reduced order model construction. *International Journal for Numerical Methods in Engineering*, 2015.

26. Y. Choi, G. Oxberry, D. White, T. Kirchdoerfer. Accelerating topology optimization using reduced order models. In The World Congress of Structural and Multidisciplinary Optimization, May 20-24. Beijing, China. 2019;04:2019.

27. Amir O, Bendsøe MP, Sigmund O. Approximate reanalysis in topology optimization. *International Journal for Numerical Methods in Engineering*. 2010;78(12):1474-91.

28. Xiao M, Lu D, Breitkopf P, Raghavan B, Dutta S, Zhang W. On-the-fly model reduction for large-scale structural topology optimization using principal components analysis. Structural and Multidisciplinary Optimization. 2020;62(1):209-30.

29. M. Xiao, D. Lu, P. Breitkopf, B. Raghavan, and S. Dutta. Multi-grid reduced-order topology optimization. *Structural and Multidisciplinary Optimization*, 2020.

30. Svanberg K. The method of moving asymptotes—a new method for structural optimization. *International Journal for Numerical Methods in Engineering*. 1987;24(2):359-73.

31. Bendsøe MP. Optimal shape design as a material distribution problem. Structural optimization. 1989;1(4):193-202.

32. Andreassen E, Clausen A, Schevenels M, et al. Efficient topology optimization in matlab using 88 lines of code. Structural and Multidisciplinary Optimization. 2011;43(1):1-16.

33. Aage N, Andreassen E, Lazarov B, Sigmund O. Giga-voxel computational morphogenesis for structural design. Nature. 2017;550(7674):84-6.
34. Wang D, Lee J, Holland K, Bibby T, Beaudoin S, Cale T. Von mises stress in chemical-mechanical polishing processes. Journal of The Electrochemical Society. 1997;144(3):1121-7.

35. Phalippou P, Bouabdallah S, Breitkopf P, Villon P, Zarroug M. 'on-the-fly' snapshots selection for proper orthogonal decomposition with application to nonlinear dynamics. Computer Methods in Applied Mechanics and Engineering. 2020;367: 113120.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.