# A statistical decomposition of geometric imperfections applied to robust topology optimisation

Noal Florian, Alain Etienne, Marco Montemurro, Jean-Yves Dantan, Julien Gardan

► To cite this version:

Noal Florian, Alain Etienne, Marco Montemurro, Jean-Yves Dantan, Julien Gardan. A statistical decomposition of geometric imperfections applied to robust topology optimisation. Reliability Engineering and System Safety, 2026, 265 (Part A), pp.111442. 10.1016/j.ress.2025.111442. hal-05502746

HAL Id: hal-05502746

https://hal.inrae.fr/hal-05502746v1

Submitted on 10 Feb 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# A statistical decomposition of geometric imperfections applied to robust topology optimisation

Noal Florian a,*, Alain Etienne a,*, Marco Montemurro b,c,†, Jean-Yves Dantan a,†, Julien Gardan d,‡

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

## ARTICLE INFO

### Keywords:

- Topology optimisation

- Robust optimisation

- Monte Carlo method

- Additive manufacturing

- Geometric imperfections

- Density based models

## ABSTRACT

The present paper deals with the integration of geometric imperfections due to the additive manufacturing process in the framework of a density-based topology optimisation method. Specifically, the model presented in this work splits the geometric imperfections into modes using a strategy based on a statistical decomposition of local geometric patterns. The aleatory uncertainty related to the lack of repeatability of the manufacturing process is propagated through the Monte Carlo method and a robust topology optimisation problem formulation is proposed. The method is applied to the classic problem of the minimisation of a cost function formulated as a weighted sum of the mean and the standard deviation of the structural compliance subject to a constraint on the volume. The effectiveness of the approach is tested on 2D benchmark problems taken from the literature and complemented by a sensitivity analysis on the convergence and accuracy of Monte Carlo estimates. Moreover, a modification of cost function weights allows to support decision-making depending on uncertainty fearfulness. The proposed method has proven effective in mitigating uncertainties arising from geometric imperfections during the design activity.

## 1. Introduction

Topology optimisation (TO) is widely used in multiple industrial sectors in order to find innovative solutions of intricate shape fulfilling different design criteria. Since its infancy, topology optimisation has seen the development of different methods, including the combination of a density-based approach with a variety of penalisation schemes, such as the isotropic material with penalisation (SIMP) scheme [1], the rational approximation of material properties (RAMP) [2], component-based methods [3,4], methods wherein the topology and/or anisotropy descriptors are represented through computer-aided design compatible entities [5,6], the level-set method [7], etc.

TO methods are today extensively used to design products of complex shapes that can be manufactured through various additive manufacturing (AM) technologies [8]. Nevertheless, AM technologies are affected by uncertainty, leading to deficiencies in the certification of such manufactured components [9]. Uncertainty can be mitigated using two strategies: robust design optimisation (RDO) and reliability-based design optimisation (RBDO). RDO aims to minimise robust measures of the system performance. Generally, the cost function is built as a

weighted sum of the mean and standard deviation, but it is possible to combine probabilities on the constrains functions too [10]. On the other hand, RBDO aims at obtaining a design that performs without failure for a required duration; to achieve this goal, the used metrics are probability of failure [11]. In this paper, the RDO is coupled with TO, hence the resulting method will be referred to as robust topology optimisation (RTO).

Uncertainty, from a modelling perspective, comes from the geometry, physical parameters as well as initial and boundary conditions [12]. For instance, uncertainty in boundary conditions is studied in [13], where uncertainty on the intensity and orientation of external forces are mitigated thanks to perturbation method. Mitigating uncertainty is also applied to dynamic structures where the intensity and the frequency of harmonic excitation are affected by uncertainty [14]. In [15], the Young's modulus and the Poisson's ratio are affected by aleatory and epistemic uncertainty; the uncertainties are aggregated and mitigated in a two-stage analysis. It is noteworthy that few research works have been conducted on the quantification of the uncertainty factor affecting, simultaneously, material properties, geometry, initial

<footnote>a Arts et Métiers Institute of Technology, Université de Lorraine, LCFC, 4 rue Augustin Fresnel, F-57078, Metz, France</footnote>

<footnote>b Université de Bordeaux, Arts et Métiers Institute of Technology, CNRS, INRA, Bordeaux INP, I2M UMR 5295, F-33405, Talence, France</footnote>

<footnote>c Institut Universitaire de France (IUF), 7500, Paris, France</footnote>

<footnote>d Institut de Recherche, ESTP, 28 Avenue du Président Wilson, F-94230, Cachan, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: florian.noal@ensam.eu (N. Florian).</footnote>
and boundary conditions. In their study, Javed and Rout [16] analysed the influence of four factors, namely Young's modulus, external forces, feature size radius and variations in the geometry. Among them, the variations in the geometry have the greatest influence. Therefore, a particular care must be put to geometry variations and imperfections due to the manufacturing process, which is the main focus of this paper.

That is why our work builds upon and extends several key developments in the field of RTO mitigating geometric imperfections, although it is not a new topic. For instance, filtering techniques inspired by mathematical morphology have been integrated into density-based TO frameworks, as shown in [17], where erosion and dilation filters are applied to the final density field to simulate thin and thick geometrical deviations from the nominal design. This approach is reformulated as a min/max compliance minimisation problem, where the topology is updated based on eroded, intermediate, and dilated designs, and later extended to improve numerical stability [18]. Moreover, in the last decade, this geometrically robust method has been successfully applied to more complex physical problems including various design requirements, such as stress-related requirements [19] or self-contacting structures [20]. They are particularly effective in eliminating undesirable geometric patterns like one-node hinges, even though they do not always ensure strict control over length scale or local smoothness.

This idea has been extended to the “standard” framework of density-based TO method, which is described [21,22]. In this framework, a smooth Heaviside filter is used to reduce intermediate density levels. Two hyperparameters control its behaviour: $\beta$ that handles filter sharpness (the higher $\beta$, the better the approximation) and $\eta$ that is a threshold that discriminates densities and enforces convergences towards a black-and-white design. Introducing uncertainty into $\eta$, as in [23], allows for the modelling of geometric imperfections. Compared to previous methods, the geometric imperfection, which affects the parameter $\eta$, is now modelled and represented with probability distributions, namely a uniform probability distribution. This approach has been extended in [24] where $\eta$ is represented by a random field. This allows spatial variations of geometric defects to be considered. However the chosen uncertainty propagation technique requires a significant computational cost. In order to alleviate this problem, the Monte Carlo simulation (MCS) was replaced with frugal approaches, like stochastic collocation methods, in [25]. Affecting $\eta$ with uncertainty will consider only relatively small variations around part boundaries. Accordingly, the centre of the convolution filter is affected by uncertainty, miming larger variations in the topology [26]. The centre is represented by a conditional random field, which better captures boundary-specific variations. Another way to include geometric variation in TO is through the mesh of the finite element model [27]. In that case, the position of nodes is affected by uncertainty. The formulated perturbation method has a computational cost lower than methods based on simulations. However, in this Lagrangian approach, one can only consider small geometric perturbations. To overcome this issue, Eulerian representations have been developed. Different TO approaches have been tested, including level-set method [28], and density-based method with SIMP scheme [29]. The geometric imperfections are introduced via fictitious advection or Hamilton-Jacobi equations. Additionally, local smoothness has been addressed using clustering techniques, such as K-mean and a tailored metric that quantifies the frequency of elements in both the horizontal and vertical directions [30]. Nevertheless, different geometric patterns may affect the accuracy of the manufactured geometry [31].

From the experimental side, geometric imperfections modelling is generally done thanks to skin model shape (SMS) approaches [32], which is sought to represent the nominal boundaries, where tolerances are needed, by a set of points. Then, to model defects, the discrete SMS is altered by changing the position of each node, resulting in a large set of parameters [31]. Domain expertises was combined in these representations, wherein the geometric variations are split into numerous defect modes [31]. However, SMS is not perfectly compliant with complex

shapes resulting from TO algorithms and made with additive manufacturing. To address this limitation, the utilisation of a three-dimensional grid of volumetric pixels, called voxels, allows comprehensive characterisation and analysis of complex geometric structures. Voxel-based metrology facilitates accurate and detailed reconstructions of objects, capturing intricate features that may be challenging to represent using traditional point clouds or surface-based methods [33,34].

Taking inspiration from the RTO approaches developed in the structural optimisation community and those used to represent geometric imperfections in the manufacturing process community, this paper introduces an original method to include geometric imperfections and the associated uncertainty in an RTO framework. Previous studies focus mainly on slight modifications of existing TO methods and algorithms to include geometric uncertainty through classical techniques widely used in the field of RDO. Unlike the approaches available in the literature, in this paper we proposed a methodology consistent with the work done in the manufacturing community and merge it to TO, which is a first step towards validation and uncertainty quantification [35] and ultimately certification [9].

In the light of the above (non-exhaustive) literature review, in this paper, we propose a new methodology together with an original modelling framework incorporating domain expertise directly into the representation of geometric imperfections in the context of TO problems. Specifically, we affect the geometry through a combination of statistical modes that characterise and penalise local geometric patterns in a differentiated manner. This allows for a more nuanced and physically meaningful treatment of imperfections, setting our contribution apart from existing methods. This model is integrated into a classic density-based TO algorithm using the SIMP scheme and applied to a compliance minimisation problem subject to a constraint on the volume fraction of the solid phase. Additionally, the sensitivity of the optimised solution to the model parameters, the weights used in the formulation of the cost function and the sample size of the Monte Carlo simulations, is also investigated.

The paper is organised as follows. Section 2 introduces the fundamentals of density-based TO method using the SIMP scheme for the classical problem of compliance minimisation under a constraint on the volume. The proposed model for geometric imperfections is discussed in detail in Section 3. Section 4 presents numerical results on 2D benchmark problems taken from the literature. Conclusions and prospects are drawn in Section 5 in fine.

## 2. Uncertainty propagation in robust topology optimisation under geometric imperfections

### 2.1. Fundamentals of the density-based topology optimisation method

In this subsection, the classical density-based approach with SIMP interpolation scheme is briefly explained for the problem of compliance minimisation subject to a constraint on the volume. This description is a first step to understand how RTO works. First, a design space, $\Omega_{ds} \subset \mathbb{R}^D$, is defined by the designer. A non-design space can also be considered where the addition or the removal of materials is prohibited, leading to regions always full ($\Omega_{full}$) or free ($\Omega_{free}$) of materials. The total space, $\Omega_{tot}$, is discretised using a mesh of which the element size is defined by the user. Let $\mathcal{L} = \{\Lambda_1, \dots, \Lambda_{n_e}\}$ denote the set containing the $n_e$ elements forming the mesh. The objective of a density-based method is to find the best material arrangement by assigning a density level to each element in $\mathcal{L}$. Each density varies between zero and one; density equal to one indicates the presence of material, whereas density equal to zero indicates the absence of material. The density-based TO algorithm allows intermediate values of the density which have no physical meaning. However, to enforce convergence of the density field towards lower and upper bounds, several penalty schemes have been developed, such as the SIMP scheme [1]. Nevertheless, optimised topologies obtained with the density-based method are prone to two main issues, viz. the
dependence of the solution on mesh refinement and the occurrence of the chequerboard effect [36]. To overcome such issues, convolution and filters are commonly used in the literature and applied in this paper. SIMP scheme is a power law that penalises intermediate densities and it is applied to the element stiffness matrix [21,22]:

$$K_e = [\alpha + (1 - \alpha)\bar{\rho}_e^p] K_e^0, \quad (1)$$

where $\bar{\rho}_e$ is the filtered element density (final density), $K_e$ and $K_e^0$ are the penalised and unpenalised element stiffness matrices, respectively, $p \ge 1$ is a penalisation constant ($p = 3$ in this paper). Finally, $\alpha > 0$ ($\alpha = 10^{-6}$ in this paper) is a ratio between void and full stiffness matrices and used to avoid singularities in the finite element model. In this work we use a convolution filter, which can also be employed to prescribe a minimum length scale through the definition of the radius $r^*$. The convolution filter alter the element densities (decision variables) such that:

$$\hat{\rho}_e = \frac{\sum_{i \in I_e} H(\Lambda_i, \Lambda_e) \rho_i v_i}{\sum_{i \in I_e} H(\Lambda_i, \Lambda_e) v_i}, \quad (2)$$

where $\hat{\rho}_e$ is the density of the generic element after applying this filter. $H(\Lambda_i, \Lambda_e) = \max(0, r^* - \|O_i - O_e\|_2)$ is a weight associated with the elements in the neighbourhood of $\Lambda_e$ within a distance $r^*$ from the point $O_e$; accordingly $I_e = \{i | \|O_i - O_e\|_2 \le r^*\}$ is the set containing the indices of such elements, i.e., it is the local support of the element, and $v_i$ is the volume of the element $\Lambda_i$. Nevertheless, this filter is not effective in addressing the issue of intermediate density levels. To this end, the Heaviside filter is used such that [18]:

$$\bar{\rho}_e = \frac{\tanh(\beta\eta) + \tanh(\beta(\bar{\rho}_e - \eta))}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))}, \quad (3)$$

where $\bar{\rho}_e$ is the final density for the element $\Lambda_e$ used in Eq. (1) to penalise the element stiffness matrix, whereas $\beta$ and $\eta$ are two parameters used to control the filter sharpness and the threshold value to filter intermediate densities, as discussed in the Introduction. The strategy used in this work to set the value of these parameters is presented in Section 4.

To solve the optimisation problem, physics and metrics must be defined. In this article, only static analyses in the framework of linear elasticity are considered and solved by the finite element method. Therefore, a linear system needs to be solved, such that:

$$K(\bar{\rho})u(\bar{\rho}) = f, \quad (4)$$

where $K$ is the global stiffness matrix, which is a function of the final density field, $u$ is the vector of generalised nodal displacements and $f$ is the vector of generalised nodal forces. The metric used here is the generalised compliance, which is related to the total potential energy, which is a self-adjoint function under non-zero Neumann–Dirichlet boundary conditions [37]. Therefore, the optimisation problem reads:

$$\begin{array}{ll} \min_{\rho} & C(\bar{\rho}) = f^T u - u^T r, \\ \text{subject to} & V(\bar{\rho}) - V^* = 0, \\ & 0 \le \rho \le 1, \\ & K(\bar{\rho})u(\bar{\rho}) = f, \end{array} \quad (5)$$

where $r$ and $u_d$ are nodal reactions and imposed displacements, respectively. The volume of the part is constrained to be equal to $V^*$ and is expressed as:

$$V(\bar{\rho}) = \sum_{i=1}^{n_e} \bar{\rho}_i v_i \quad \text{and} \quad \frac{\partial V}{\partial \rho_e} = \sum_{j \in I_e} \frac{\partial \bar{\rho}_j}{\partial \rho_e} \frac{\partial \bar{\rho}_j}{\partial \rho_e} v_j. \quad (6)$$

The sensitivity of the compliance reads as follows [37]:

$$\frac{\partial C}{\partial \rho_e} = -u^T \frac{\partial K}{\partial \rho_e} u = - \sum_{j \in I_e} \frac{\partial \bar{\rho}_j}{\partial \rho_e} \frac{\partial \bar{\rho}_j}{\partial \rho_e} \frac{p \bar{\rho}_e^{p-1}}{\alpha + (1-\alpha) \bar{\rho}_e^p} u_j^T K_j u_j, \quad (7)$$

where $u_j$ is the vector of the generalised displacements of the $j$th element. Partial derivatives $\partial \bar{\rho}_j / \partial \bar{\rho}_i$ and $\partial \bar{\rho}_i / \partial \rho_e$ appearing in Eqs. (6) and (7) read:

$$\frac{\partial \bar{\rho}_j}{\partial \rho_e} = \frac{H(\Lambda_j, \Lambda_e) v_e}{\sum_{i \in I_j} H(\Lambda_i, \Lambda_j) v_i}, \quad (8)$$

$$\frac{\partial \bar{\rho}_e}{\partial \bar{\rho}_e} = \beta \frac{1 - \tanh(\beta(\bar{\rho}_e - \eta))^2}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))}. \quad (9)$$

The above formulae are used to compute gradients. They are employed to feed the optimisers used in this paper, for instance the method of moving asymptotes (MMA) [38] from the open-source library NLOpt [39].

## 2.2. Robust formulation with Monte Carlo simulation

This subsection extends the optimisation problem of Eq. (5) to include, propagate [40] and mitigate aleatory uncertainties. First of all, let $\Xi = [\Xi_1, \dots, \Xi_{n_{rv}}]$ denote a random vector containing random variables mapping uncertain physical parameters, geometry, etc. RTO objective is to propagate $\Xi$ through a physical model, compute measures and optimise the topology to satisfy the requirements of the problem at hand. In this paper, the random variable $Y$ quantifies the compliance of the part and $\mathbb{P}_Y$ is its probability distribution. Therefore, the problem (5) is reformulated as:

$$\begin{array}{ll} \min_{\rho} & \bar{C} = k m_1(Y) + (2-k)\sqrt{\mu_2(Y)}, \\ \text{subject to} & V(\bar{\rho}) - V^* = 0, \\ & 0 \le \rho \le 1, \end{array} \quad (10)$$

where $m_1(Y)$ and $\mu_2(Y)$ are, respectively, the mean and the variance of $Y$ and $k \in [0, 2]$ is a weight used to balance the contributions due to the mean and the standard deviation. $\bar{C}$ is a robust measure for compliance. Of course, other measures could be used, e.g., probabilistic threshold [41]. Numerous techniques can be found in the literature and impact the accuracy of the estimate. For instance, MCS is replaced with a stochastic reduced order model in [42]: this approach allows reducing the size of the sample, while being a direct substitute for MCS. Nevertheless a full MCS is the most reliable and robust approach even if it requires high computational costs. The estimate of both statistical moments, when MCS is performed, read:

$$\hat{m}_1(Y) = \frac{1}{n_{\text{mcs}}} \sum_{i=1}^{n_{\text{mcs}}} C(w_1^i, \dots, w_{n_{\text{rv}}}^i), \quad (11)$$

and:

$$\hat{\mu}_2(Y) = \frac{1}{n_{\text{mcs}} - 1} \sum_{i=1}^{n_{\text{mcs}}} \left[ C(w_1^i, \dots, w_{n_{\text{rv}}}^i) - \hat{m}_1(Y) \right]^2, \quad (12)$$

where $w_i^j$ is the $i$th realisation (out of $n_{\text{mcs}}$) of the $j$th random variable and $w^i = [w_1^i, \dots, w_{n_{\text{rv}}}^i]$ is a vector that contains the $i$th realisation for the $n_{\text{rv}}$ random variables. By using the expansion $(\mathbb{V}(Y) = \mathbb{E}[Y^2] - \mathbb{E}[Y]^2)$ and invariant $(\mathbb{V}(Y - K) = \mathbb{V}(Y))$ properties of the variance, round-off errors of Eq. (12) are limited [43]. Accordingly, Eq. (12) can be rewritten as follows:

$$\hat{\mu}_2(Y) = \frac{n_{\text{mcs}}}{n_{\text{mcs}} - 1} \left[ \frac{1}{n_{\text{mcs}}} \sum_{i=1}^{n_{\text{mcs}}} [C(w^i) - C_{\text{ref}}]^2 - \left( \frac{1}{n_{\text{mcs}}} \sum_{i=1}^{n_{\text{mcs}}} C(w^i) - C_{\text{ref}} \right)^2 \right], \quad (13)$$

where $C_{\text{ref}}$ is a constant that should be as close as possible to the mean. Since the mean is not known, $C_{\text{ref}}$ is equal to $C(\bar{\rho})$, which is the compliance of the design without uncertainties. In other words, it is the compliance of the design characterised by the field $\bar{\rho}$ and calculated without any perturbation due to $\Xi$. Throughout the rest of this article, the optimised solution obtained under these conditions is referred to as reference design.

The constraint in Eq. (10) is formulated in terms of the volume of the part without uncertainties. In RTO, the constraint function on the volume can be replaced with an estimation of mean volume, or with a probabilistic constraint, where the occurrence of an undesired event
is bounded [11]. To solve the optimisation problem, sensitivities of Eqs. (11) and (13) must be calculated. Their expressions read:

$$ \frac{\partial \hat{m}_1}{\partial \rho_e} = \frac{1}{n_{mcs}} \sum_{i=1}^{n_{mcs}} \frac{\partial C(\mathbf{w}^i)}{\partial \rho_e}, \qquad (14) $$

$$ \frac{\partial \hat{\mu}_2}{\partial \rho_e} = \frac{2}{n_{mcs}(n_{mcs} - 1)} \left[ n_{mcs} \sum_{i=1}^{n_{mcs}} [C(\mathbf{w}^i) - C_{\text{ref}}] \frac{\partial C(\mathbf{w}^i)}{\partial \rho_e} - \left( \sum_{i=1}^{n_{mcs}} [C(\mathbf{w}^i) - C_{\text{ref}}] \right) \left( \sum_{i=1}^{n_{mcs}} \frac{\partial C(\mathbf{w}^i)}{\partial \rho_e} \right) \right]. \qquad (15) $$

To avoid memory issues, increment operators is applied and care was taken as multiprocessing is used to reduce computation time due to the non-atomicity of these operations [44]. It is noteworthy that, for each iteration of the optimisation process, a new sample of the random vector $\Xi$ is drawn. In the literature, some authors draw an initial sample and keep the same during the optimisation [26,45]. We decided to “resample” at each iteration to avoid that the optimiser overfit the initial sample.

### 3. Geometric imperfections modelling

In this section, the geometric imperfections model is detailed together with its integration into the RTO framework presented in the above section.

Most of the geometric imperfections models applied to TO do not include information about local geometric features. However, manufacturing processes, like additive manufacturing, are characterised by defects that can be decomposed as a combination of elementary modes [31]. The objective of the proposed model is to consider local geometric features in the robust optimisation. As mentioned in the Introduction, the proposed model is compliant with density-based TO methods because of the need to integrate a full verification, validation and uncertainty quantification framework. Indeed, today new metrics are defined for tolerance analysis based on three-dimensional X-ray computed tomography data and voxelisation of geometrical models [34]. Experimental data and nominal geometry are both voxelised, resulting in the same topological descriptor as the density-based method. Accordingly, they can be compared for tolerance analysis purposes. Voxel modelling of uncertain geometric imperfections has already been considered for verification, validation and uncertainty quantification of multiscale models [35], strengthening the certification of additive manufactured parts. The above works have encouraged us to adopt a voxel-based, and, hence, a density-based modelling approach.

For the density-based method, geometric imperfections are modelled mainly by affecting the parameters of the penalisation schemes involved in the TO process with aleatory uncertainty. Nevertheless, from a physical perspective, it is preferable to introduce uncertainty in the final density field $\bar{\rho}$ rather than in the parameters of the penalisation scheme, in order to be as close as possible to the effective geometric imperfections due to the AM process. Accordingly, the perturbed density field $\bar{\rho}$ can be expressed as follows:

$$ \bar{\rho}_e = \bar{\rho}_e + \Delta(\bar{\rho}), \qquad (16) $$

where $\bar{\rho}$ is a vector collecting all the element densities affected by uncertainty via the operator $\Delta : \mathbb{R}^{n_e} \to \mathbb{R}$. Indeed, this operator models the geometric deviation (i.e. the imperfections) with respect to the nominal design due to the manufacturing process. This model could be characterised by either partial differential equations, expert knowledge, data or any combinations. In Eq. (16), $\bar{\rho}_e$ falls in the interval $[0, 1]$, meaning that the element densities values obtained after applying the perturbation are adjusted when they are outside this interval. It is noteworthy that, in the proposed framework, the nominal value of the final density of the non-design regions is set equal to zero. Nevertheless, geometric perturbations can also occur in these regions. Accordingly, for elements belonging to the non-design regions the perturbed density is $\bar{\rho}_e = \Delta(\bar{\rho})$.

In our implementation, only the elements falling in the neighbourhood of the element e are inputs of the model, which signifies that the model considers local pattern deviations from the manufacturing process. $S_e$ denotes the set collecting the indices of the elements falling in the neighbourhood of element e as follows:

$$ S_e = \{ i \mid \|O_i - O_e\|_{p_{\text{loc}}} \le r_{\text{loc}} \}, \qquad (17) $$

where $r_{\text{loc}}$ is the radius of the ball centred at $O_e$ equipped with the p-norm, the parameter $p_{\text{loc}}$ must be properly tuned. The use of p-norm allows sketching “square” windows regardless of the value of $r_{\text{loc}}$, provided that $p_{\text{loc}}$ is sufficiently large, e.g., $p_{\text{loc}} > 10$. The densities belonging to $S_e$ are collected in the vector $\bar{\rho}_e$, although the density of the element e does not belong to this set and, hence, it is not stored in the vector $\rho_e$, i.e.

$$ \bar{\rho}_e = \{ \bar{\rho}_i \mid i \in S_e \wedge i \neq e \}. \qquad (18) $$

Let $n_\#$ be the cardinality of the set $S_e$, i.e. the number of elements belonging to $S_e$. Therefore the operator $\Delta$ is a function of local elements, thus only neighbouring geometric features are ipso facto considered and Eq. (16) is rewritten as:

$$ \bar{\rho}_e = \bar{\rho}_e + \Delta(\bar{\rho}_e). \qquad (19) $$

The aim of the proposed geometric imperfections model is twofold:

* It allows detecting topology boundaries that cannot be detected when using density-based TO algorithms, which is a shortcoming when geometric imperfections are considered [28]. This aspect is handled by introducing a threshold on the local standard deviation (STD).

* It allows distinguishing local features, such as boundary curvatures; and assign randomness to these imperfection modes. In our model, modes are discriminated through a statistical analysis of local elements. Since manufacturing processes, for instance additive manufacturing [31], exhibit imperfection mode, this aspect is crucial in our method.

Fig. 1 illustrates values of skewness and STD for some representative local geometric features and boundary layout in two dimensions. Cases A and B exhibit straight boundaries, leading to null skewness and high STD. In cases C and D, the presence of a fillet radius in the vicinity of the element $A_e$ leads to a reduced standard deviation but an increase in skewness. Invariance in rotation is shown in cases A and B. Cases C and D are also invariant in rotation. However, it is possible to distinguish between cases C and D. If the fillet is rounding an interior corner of the topology (case C), then the skewness is negative; conversely, when the fillet is rounding an exterior corner (case D) the skewness is the same in intensity yet opposite sign. Finally, case E represents a connection between two topological branches, leading to a reduction in standard deviation and skewness. Consequently, the proposed geometric imperfection model is well suited for density-based TO algorithms and can account for the local boundary curvature. Given the aforementioned observations and analysis, we infer that the deviation model in Eq. (19) can be stated as follows:

$$ \Delta(\bar{\rho}_e) = \Delta_e = f(\text{STD}(\bar{\rho}_e), a) \left( k_{\text{STD}} \text{STD}(\bar{\rho}_e) + k_{\text{Sk}} \left| \text{Sk}(\bar{\rho}_e) \right| \right) \Xi_e^{\text{AM}}, \qquad (20) $$

with:

$$ f(\bar{\rho}_e, a) = \begin{cases} 1 & \text{if } \text{STD}(\bar{\rho}_e) \ge a, \\ 0 & \text{otherwise.} \end{cases} \qquad (21) $$

In the above formulae, STD and Sk are the functions used to compute the standard deviation and the skewness of the vector $\bar{\rho}_e$, respectively. To discriminate imperfection modes, two weight factors are used: $k_{\text{STD}}$ that is associated with the value of the standard deviation; and $k_{\text{Sk}}$ for the skewness. Finally, $\Xi_e^{\text{AM}}$ is a random variable allocated to the density $\bar{\rho}_e$ and is related to the manufacturing process. Since experimental data are not available to correctly characterise the
<figure><img src="image_4.png" /><figcaption>Fig. 1. Examples of local geometric features and their associated standard deviation and skewness.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Matching of an imperfection and the element size.</figcaption></figure>

parameters of the present model, they have been selected in order that most deviation values $\Delta_e$ fall between -1 and +1. Here the random variables are uncorrelated in position and can be considered as a discretised version of random fields where correlation lengths are sufficiently small. The parameter $a$ in Eq. (21) is set as $a = 10^{-4}$.

To integrate the model of geometric imperfections of Eqs. (19) and (20) in the density-based TO algorithm, it suffices to replace the final density $\bar{\rho}_e$ with the perturbed density $\bar{\rho}_e$ in the expression of the penalised stiffness matrix of the generic element of Eq. (1) and to apply the chain rule for derivative of the generic response function $z$ as follows:

$$ \frac{\partial z}{\partial \bar{\rho}_e} = \sum_{r \in S_e} \sum_{j \in I_e} \frac{\partial z}{\partial \bar{\rho}_r} \frac{\partial \bar{\rho}_r}{\partial \bar{\rho}_j} \frac{\partial \bar{\rho}_j}{\partial \bar{\rho}_e} . \quad (22) $$

In Eq. (22) the partial derivative $\partial \bar{\rho}_r / \partial \bar{\rho}_j$ is equal to:

$$ \frac{\partial \bar{\rho}_r}{\partial \bar{\rho}_j} = \begin{cases} 1 & \text{if } r = j, \\ \frac{\partial \Delta(\bar{\rho}_r)}{\partial \bar{\rho}_j} & \text{otherwise.} \end{cases} \quad (23) $$

The term $\partial \Delta(\bar{\rho}_r) / \partial \bar{\rho}_j$ of Eq. (23) can be calculated according to the following steps. Firstly, the STD partial derivative of Eq. (23) is computed as follows:

$$ \frac{\partial \mathrm{STD}(\bar{\rho}_r)}{\partial \bar{\rho}_j} = \frac{\partial \mathrm{Var}(\bar{\rho}_r)}{\partial \bar{\rho}_j} \frac{1}{2\sqrt{\mathrm{Var}(\bar{\rho}_r)}}, \quad (24) $$

with:

$$ \begin{aligned} \frac{n_{\#} - 1}{n_{\#}} \mathrm{Var}(\bar{\rho}_r) &= \frac{1}{n_{\#}} \sum_{i \in S_r} \bar{\rho}_i^2 - \left( \frac{1}{n_{\#}} \sum_{i \in S_r} \bar{\rho}_i \right)^2, \\ &= \frac{1}{n_{\#}} \sum_{i \in S_r} \bar{\rho}_i^2 - (\hat{m}_1(\bar{\rho}_r))^2, \end{aligned} \quad (25) $$

and

$$ \frac{\partial \mathrm{Var}(\bar{\rho}_r)}{\partial \bar{\rho}_j} = \frac{2}{n_{\#} - 1} \left( \bar{\rho}_j - \frac{2}{n_{\#}} \sum_{i \in S_r} \bar{\rho}_i \right). \quad (26) $$

Secondly, the expressions of the skewness and of its partial derivative are written as follows:

$$ \mathrm{Sk}(\bar{\rho}_r) = \frac{\frac{1}{n_{\#}} \sum_{i \in S_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r))^3}{\left( \frac{1}{n_{\#}-1} \sum_{i \in S_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r))^2 \right)^{3/2}}, \quad (27) $$

$$ \begin{aligned} \frac{\partial \mathrm{Sk}(\bar{\rho}_r)}{\partial \bar{\rho}_j} = \frac{3}{2} \left( \sum_{i \in S_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r))^3 \right) & \frac{\left( \frac{2}{n_{\#}} \sum_{i \in \mathcal{T}_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r)) + 2 \left( \frac{1}{n_{\#}} - 1 \right) (\bar{\rho}_j - \hat{m}_1(\bar{\rho}_r)) \right)}{n_{\#}(n_{\#} - 1) \left( \frac{1}{n_{\#} - 1} \sum_{i \in S_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r))^2 \right)^{3/2+1}} \\ & - \frac{3}{n_{\#}} \frac{\sum_{i \in \mathcal{T}_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r))^2 + 3 \left( \frac{1}{n_{\#}} - 1 \right) (\rho_j - \hat{m}_1(\bar{\rho}_r))^2}{n_{\#} \left( \frac{1}{n_{\#} - 1} \sum_{i \in S_r} (\bar{\rho}_i - \hat{m}_1(\bar{\rho}_r))^2 \right)^{3/2}}, \end{aligned} \quad (28) $$

where $\mathcal{T}_r = \{i \in S_r | i \neq j\}$. However, in Eqs. (27) and (28) the skewness is biased and $n_{\#}^2/(n_{\#}-1)(n_{\#}-2)$ is added to the right-hand side of these equations. In Eq. (20) two more terms still need to be derived. Firstly, the absolute value is approximated with an expression that is a
<figure><img src="image_6.png" /><figcaption>Fig. 3. Geometry and boundary conditions of the 2D benchmark problem.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Optimised topologies obtained without uncertainty with MMA 4(a) and OC 4(c) and the juxtaposition 4(b).</figcaption></figure>

function of the error function (erf) [46], such as

$$|Sk(\bar{\rho}_r)| \approx Sk(\bar{\rho}_r) \erf(10Sk(\bar{\rho}_r)). \qquad (29)$$

Subsequently, the partial derivative of the above function is calculated as follows:

$$\frac{\partial |Sk(\bar{\rho}_r)|}{\partial \bar{\rho}_j} = \frac{\partial Sk(\bar{\rho}_r)}{\partial \bar{\rho}_j} \left( \erf(10Sk(\bar{\rho}_r)) + \frac{2Sk(\bar{\rho}_r) \exp(-100Sk(\bar{\rho}_r))}{10\sqrt{\pi}} \right). \qquad (30)$$

Finally, the indicator function in Eq. (21) is replaced with the following smooth function:

$$f(\bar{\rho}_r, a) = \frac{1}{\exp(10(a - \text{STD}(\bar{\rho}_r))) + 1}, \qquad (31)$$

and the partial derivative of the aforementioned equation is written as follows:

$$\frac{\partial f(\bar{\rho}_r, a)}{\partial \bar{\rho}_j} = \frac{\partial \text{STD}(\bar{\rho}_r)}{\partial \bar{\rho}_j} \frac{10 \exp(a - \text{STD}(\bar{\rho}_r))}{(\exp(10(a - \text{STD}(\bar{\rho}_r))) + 1)^2}. \qquad (32)$$

Initially, to simplify the problem, the partial derivative $\partial\Delta(\bar{\rho}_r)/\partial\bar{\rho}_j$ is set to zero.

Since the TO problem is solved using a fixed mesh, the imperfection size is limited by the element size. When the imperfections are larger

than the element size, they can be adjusted to the actual mesh size by using bigger imperfection patterns, as explained in [16] and illustrated in Fig. 2. If the imperfections are smaller than the element size, then multi-scale topology optimisation might be a solution [47].

Finally, as the goal of the proposed strategy is to include geometric imperfections in the TO process, one might inquire as to what occurs outside the boundary of the design space. It is not possible to prevent imperfections resulting from the manufacturing process from extending beyond the boundaries of the prescribed design space. This is simply an input provided by the designer. In the present paper, we add a non-design space that surrounds the prescribed one, yet this is limited to the region of the frontier where no boundary conditions are applied. In the literature only a few studies consider the spread of geometrical imperfections. To the best of the authors' knowledge, only [26] has addressed this issue, but using a different model for geometric imperfections.

## 4. Numerical experiments and results

In this section a 2D benchmark problem is introduced to analyse the influence of the optimised solution on the following factors: repeatability due to the MCS, sensitivity of the solution to the robust objective cost function, by modifying the parameter k of Eq. (10), and the behaviour of our model by altering the two weight factors in Eq. (20).

Fig. 3 shows the geometry, the boundary conditions and the definition of the design and non-design regions for the benchmark problem. The design space is a rectangle having dimension L = 160 mm and H = 40 mm. The nodes located at x = 0 are clamped, whereas a point force F = (0, -1) N is applied to the node located at (x, y) = (L, H/2). As stated above, a non-design space encompasses the design space except on the faces where a boundary condition is applied. The mesh of the finite element model is made of 160 × 40 elements with four nodes and two degrees of freedom per node for the design region, plus 2 × 2 × 160 elements for both non-design regions (top and underside). The mesh is generated with Gmsh [48]. Finite element analysis is executed with the version 0.6 of Fenicsx [49]. Young's modulus and Poisson's ratio are equal to E = 1 MPa and ν = 0.3, respectively. The value of radius for the convolution filter is equal to r* = 1.3 mm, which is the minimal value to avoid chequerboard patterns. It is noteworthy that the plane stress hypothesis is tacitly assumed for all the analyses presented in the following of the paper. To ease the numeric implementation, the code is based on the object-oriented programming paradigm of Python language [50].

The solution search for problem (10) is carried out by considering two different approaches. In the first one, we employ optimality criteria (OC) [22] and the value of the parameter η (see Eq. (3)) is kept constant and set to 0.5 unless otherwise stated. In the second one, the method of moving asymptotes (MMA) algorithm [38] from NLOpt Python package [39] is used. For OC the “move limit” is set to 0.2. Regarding
<figure><img src="image_8.png" /><figcaption>Fig. 5. Optimised topologies obtained with MMA 5(c) and OC 5(d) for robust topology optimisation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Bubble charts of standardised mean and standard deviation compliance for various values of *k* using the OC and the MMA algorithm.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Optimised topology obtained with $k = 0.05$ and the OC 7(a) for robust topology optimisation.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. $s_d$ versus *k* for the OC and the MMA algorithm.</figcaption></figure>

the MMA algorithm, the default values of the main parameters have been considered. The constraint on the volume in Eq. (10) has been set by considering a volume fraction equal to 40% of the volume of $\Omega_{ds}$. It is noteworthy that for the MMA algorithm the equality constraint on the volume has been transformed into two inequality constraints by considering a tolerance of 5% on the predefined volume fraction. The stopping criterion of the algorithms is based on the number of iterations, because the cost function could be noisy. This is inconsequential as the cost function shows similar behaviour in convergence rates. Therefore the number of iterations is chosen accordingly. Finally, a continuation scheme is used to increase the value of $\beta$ (see Eq. (3)).

incrementally. In the first step, the parameter $\beta$ is set to $\beta = 1$ and remains constant for 250 iterations. Subsequently, $\beta$ is increased to $\beta = 4$. After an additional 20 iterations, $\beta$ is updated to $\beta = 8$. Finally, the value of $\beta$ is doubled every 5 iterations until it reaches $\beta = 64$. Moreover, no continuation scheme is used to update the value of the penalty parameter of the SIMP scheme in Eq. (1), which is kept constant and equal to $p = 3$. The results of the optimisation without uncertainty (deterministic case) are shown in Fig. 4. The compliance is equal to 393 N mm (when $\beta = 64$ the compliance oscillates between 391 N mm and 402 N mm) and 372 N mm when using OC and the MMA algorithm, respectively. The compliance with MMA optimiser is lower because the volume fraction is about 43% ($\le 45\%$) of the design space, but for OC
<figure><img src="image_12.png" /><figcaption>Fig. 9. Box charts of different metrics used in the cost function.</figcaption></figure>

optimiser the fraction is equal to 40%. Beam-like topological features are present in both topologies, which are symmetric with respect to the plane $y = H/2$.

## 4.1. First results

Fig. 5 illustrates the results of the RTO with model nominal values, namely $k_{STD} = 1.6$, $k_{Sk} = 0.2$, $k = 1$ and $n_{mcs} = 4000$. As far as the other parameters involved in the formulation of the RTO problem are concerned, unless otherwise specified, in this section the following values are used: $\Xi_e^{AM} \sim N(0, 0.4)$, $a = 10^{-4}$, $r_{loc} = 1.5$ mm and $p_{loc} = 2$. Both cases show similar trends: the topological branches are thicker, especially where one-node connected hinges occur. Moreover, the top and bottom topological branches exhibit a small slope, which is missing when the optimisation is performed without uncertainty or without a non-design space. Lastly, both designs lose their symmetry, even if it is not easily distinguishable due to the substantial size of Monte Carlo (MC) sample. Due to the size of the MC sample, which is made of 4000 individuals (corresponding to 4001 finite element simulations per iteration), the computation cost is significant. For instance, when using a single central processing unit (CPU) of a Core i9-14900K with 64 GB of RAM, the computation requires approximately 270 h in sequential mode. When multiprocessing is employed with 32 CPUs cores, the computation time is reduced to about 16 h.

The optimised solutions obtained without uncertainty (Fig. 4(c)) and by mitigating uncertainty (Fig. 5(d)) show differences. In fact, the designs obtained using the OC are similar except in the middle region of the domain. Specifically, the design obtained through the mitigation of uncertainty exhibits a greater number of small topological branches, which reduce the compliance of the structure. This often occurs when the continuation scheme is used on the SIMP penalty parameter. When uncertainty on the geometry is considered, these small features arise at the beginning of the optimisation and this is due to the high variability coming from uncertainty estimation (more details on this aspect are given in Section 4.3). Regarding the solutions obtained with the MMA algorithm, the differences between the design obtained with (Fig. 5(c)) and without (Fig. 4(a)) uncertainty are more noticeable, especially

<table><caption>Table 1
Compliance [N mm] of the solutions shown in Figs. 4 and 5.</caption><thead><tr><td>Case</td><td>Reference</td><td>Mean</td><td>Standard deviation</td><td>Maximum</td><td>Minimum</td></tr></thead><tbody><tr><td>Fig. 4(a)</td><td>372</td><td>485</td><td>11.7</td><td>685</td><td>459</td></tr><tr><td>Fig. 5(c)</td><td>385</td><td>490</td><td>6.60</td><td>524</td><td>468</td></tr><tr><td>Fig. 4(c)</td><td>393</td><td>510</td><td>10.3</td><td>666</td><td>480</td></tr><tr><td>Fig. 5(d)</td><td>403</td><td>520</td><td>7.50</td><td>561</td><td>499</td></tr></tbody></table>

in the left part of the design domains where a greater number of small topological branches occur. Table 1 reports the results of the designs shown in Figs. 4 and 5. This table includes statistical measures of designs compliance, namely mean, standard deviation, maximum and minimum values. The statistical measures for the deterministic cases were estimated by propagating uncertainty through our geometric imperfections model, using the same parameters values given at the beginning of this section and with a sample of 4000 individuals. These calculations are useful to understand if aleatory uncertainties are effectively mitigated when applied to the proposed approach. As expected, the compliance is lower if geometric uncertainty is not included in the problem formulation. However, their estimated statistical measures are worse for the range and the STD compliance.

## 4.2. Cost function insights

In this subsection, the behaviour of the cost function is analysed by varying the parameter $k$ of Eq. (10) and by assessing its influence on the optimised solution. Fig. 6 illustrates a bubble chart obtained for different values of $k$ when using both the OC and the MMA algorithm to carry out the solution search for the RTO problem. In this figure, data are scaled with min-max normalisation of the compliance values to ease the reading. For each analysis, the following values of the parameter $k$ have been considered: 0.01, 0.05, 0.1, 0.5, 1.0, 1.5, 2.0. It is noteworthy that, for the first three values, MMA optimiser does not produce a feasible topology. Nevertheless, for the other values of $k$, the same trend is observed for both the OC and the MMA algorithm: the lower $k$, the lower compliance STD; and the higher $k$ the lower mean compliance.
<figure><img src="image_13.png" /><figcaption>Fig. 10. Mean (left) and entropy (right) of the reference designs with the MMA algorithm.</figcaption></figure>

It is noteworthy that the discrepancies between the geometries with the best standard deviation compliance (also the worst mean compliance) and the best mean compliance are marginal. This behaviour is attributable to the fact that for $k \ge 0.1$, the influence of the mean value of the compliance in Eq. (10) is more pronounced than the influence of the standard deviation of the compliance on the optimised solution. However, even if the geometry is not feasible, when lowering $k$ to 0.05, the resulting geometry shows thicker branches everywhere and is more slender. The poorly defined part is in its first quarter, as shown in Fig. 7, and this may be due to:

* Ill-posed problem because minimising only the standard deviation of the compliance leads to a meaningless design. For instance, a solution exhibiting a uniform density field equal to the volume fraction is already variance-free. Particularly, the convergence curves for $k = 0.01$ reveal that the mean and the reference compliance after about twenty iterations increase gradually, but the STD remains almost constant.

* The optimiser is unable to accommodate the volume constraint, although in such cases the convergence curves demonstrate stable metrics prior to the implementation of the continuation scheme on $\beta$. This is due to the occurrence of a considerable number of elements characterised by intermediate values of density (which are eliminated by the Heaviside filter).

It is noteworthy that, the MCS method allows calculating other metrics without further efforts. To check how our model affects the

geometry, we computed the integral of mean deviation $\hat{m}_1(\Delta_e)$ between the reference design and the sampled geometries. If the mesh is regular, then the integral of the mean deviation takes the following form:

$$s_{\Delta} = \int_{\Omega_{tot}} \hat{m}_1(\Delta) d\Omega_{tot} \approx \sum_{e=1}^{n_e} \hat{m}_1(\Delta_e) v_e. \quad (33)$$

Fig. 8 shows values of $s_{\Delta}$ for different $k$ values of the optimised solutions obtained with both the OC and the MMA algorithm. It can be observed that all values are negative, which indicates that our model displays a tendency to overstretch the geometries. This may be a contributing factor to the observed discrepancy between the mean and the reference compliance, as reported in Table 1. In addition, if $k$ is increased (the mean is favoured), then less geometric deviation is observed.

### 4.3. Repeatability

When an MCS is performed, two parameters of this method influence simulation results, namely the sample size and the random seed. If the estimation $P_Y$ is good enough, the optimised results should be insensitive to these factors. To this end, we analyse the sensitivity of the reference, the mean and the STD compliance, by considering five sample size values, namely 4000, 2000, 1000, 500, 100, and five random seeds. Fig. 9 illustrates the box charts representing the variation of the aforementioned metrics. Variability in these metrics decreases when the number of samples increases. Unfortunately, even with large samples, variability is persistent and is representative of geometric changes in
<figure><img src="image_14.png" /><figcaption>Fig. 11. Mean (left) and entropy (right) of the reference designs with the OC.</figcaption></figure>

reference designs. The MMA optimiser seems to be more sensitive to the statistical moments accuracy. Regarding the computing time, using a sample of 100 individuals requires approximately 5 h with a single CPU.

Figs. 10 and 11 illustrate the mean of the 25 reference designs for the OC and the MMA algorithm, respectively. Moreover, these figures also illustrate the entropy: a high entropy means that the different seeds give similar density levels, and vice versa. In all cases, the topologies are nearly consistent amongst themselves. However, a larger sample size culminates in a convergence in the optimised solution shown in Fig. 5, for both the OC and the MMA algorithm. Nevertheless, it is difficult to establish a general guideline to select the appropriate sample size. A large sample size would result in convergent and robust design, but it would require a high computational effort. A small sample size would result in a solution characterised by high variability and, in some cases, an inability to converge (this is evidenced by the solution obtained with the MMA algorithm when considering a sample size of 100). The simulation method and its parameterisation are a source of epistemic uncertainty that should not be ignored. However, the mitigation of this uncertainty comes at the cost of a significant increase in the required sample size [51]. To the authors' best knowledge, it is the first time that the robustness assessment of the MCS is tested by changing simulation parameters at the beginning of the optimisation process for RTO problems.

<table><caption>**Table 2**

Compliance [N mm] of the solutions shown in Fig. 12.</caption><thead><tr><th>Case</th><th>Reference</th><th>Mean</th><th>Standard deviation</th><th>Maximum</th><th>Minimum</th></tr></thead><tbody><tr><td>Fig. 12(a)</td><td>403</td><td>520</td><td>7.50</td><td>561</td><td>499</td></tr><tr><td>Fig. 12(b)</td><td>403</td><td>514</td><td>6.82</td><td>547</td><td>492</td></tr><tr><td>Fig. 12(c)</td><td>396</td><td>437</td><td>2.44</td><td>447</td><td>429</td></tr></tbody></table>

## 4.4. Sensitivity of the optimised solution to the parameters of the geometric imperfections model

In this subsection the behaviour of the geometric imperfections model is investigated, by analysing the influence of the parameters $k_{STD}$ and $k_{Sk}$ of Eq. (20) on the optimised solution. To this end, three different combinations have been considered: $k_{STD} = 1.6$ and $k_{Sk} = 0.2$ (Fig. 12(a)), $k_{STD} = 2$ and $k_{Sk} = 0$ (Fig. 12(b)), and $k_{STD} = 0$ and $k_{Sk} = 0.4$ (Fig. 12(c)). Table 2 references the results obtained with the OC and Fig. 12 illustrates the distinct reference design configurations obtained for each combination of $k_{STD}$ and $k_{Sk}$.

The results shown in this subsection used only the OC, because the MMA algorithm did not converge when $k_{STD} = 0$ and $k_{Sk} = 0.4$. Furthermore, the convergence curves exhibit spiky behaviour with the OC. Moreover, in such a case a decrease of the standard deviation and mean compliance is observed, which indicates that this configuration has negligible influence on the cost function. This is corroborated by a visual analysis of the juxtaposition with the optimised solution that has been obtained without uncertainty. Both designs are close,
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Optimised topologies obtained by varying $k_{STD}$ and $k_{Sk}$ with the OC for robust topology optimisation.</figcaption></figure>

e.g. one-node connected hinges is not resolved. The skewness properties seem to influence part topology. In fact, as the skewness is null for straight boundary configuration (see Fig. 1), the topological branches are straighter compared to the nominal model design (where $k_{STD} = 1.6$ and $k_{Sk} = 0.4$).

Finally, if the parameter $k_{Sk}$ is set to zero, then the optimised design is much closer to the nominal one, which is consistent with previous remarks, i.e. $k_{Sk}$ is the least influential parameter. Particularly, when $k_{STD} \neq 0$, one-node connected hinges do not occur in the final design. Therefore, $k_{STD}$ is the main parameter that has a positive influence on the structure in terms of global thickness.

To see how the geometric imperfections model is altering the topology, Fig. 13 illustrates the integral of mean deviation between the reference designs and the sampled geometries for Figs. 12(b) and 12(c). When $k_{Sk} = 0.4$ and $k_{Sk} = 0$ (see Fig. 13(b)), the values of $\hat{m}_1(\Delta_e)$ at the bottom and top topological branches are smaller than the counterpart obtained when $k_{STD} \neq 0$ and highest values are located at the intersection of the topological branches, which is consistent with our model. The topology obtained in the case $k_{Sk} = 0.4$ and $k_{Sk} = 0$ presents high value of $\hat{m}_1(\Delta_e)$ in the outer and inner parts of the boundaries where topological branches are inclined. This is due to the small value of the parameter $a$ in Eq. (20).

## 5. Conclusions

In this article, we presented a geometric imperfection models, which decompose imperfections into imperfection modes in order to include domain expertise and knowledge. This model has been integrated into a density-based topology optimisation algorithm in order to propose a general framework for robust topology optimisation including geometric uncertainty due to the manufacturing process. To mitigate uncertainties in topology optimisation, robust design optimisation is used and the moments are estimated through the Monte Carlo simulation approach.

To verify that the geometric imperfections model and the optimisation problem perform as expected, a design of experiments is

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** $\hat{m}_1(\Delta_e)$ field for the cases in Figs. 12(b) and 12(c) at the iteration 249.</figcaption></figure>

conducted. Firstly, modifying the weights of the objective function demonstrated that it is possible to optimise the topology based on the designer's aversion to uncertainties. Increasing the weight assigned to the standard deviation results in a globally thicker structure but reduces its average performance (the mean compliance). Thus, compliance serves as a metric to account for geometric imperfections in topological branches. Additionally, a repeatability study revealed that the topology is significantly influenced by the accuracy of the statistic moments estimation and, consequently, on the partial derivatives. Establishing a general rule for selecting an appropriate sample size remains challenging. Finally, the geometric imperfections model was tested by its associated weights. The model proved consistent with the description provided in Section 3. Specifically, assigning a higher value to $k_{STD}$ leads to a uniform increase in topological branches thickness, while increasing $k_{Sk}$ results in a greater number of horizontal topological branches.

As far as the prospects of this work are concerned, the proposed approach should be extended and tested on 3D problems by employing a more efficient method to calculate statistical moments. The model of imperfections should also be generalised by considering different imperfections shapes and mesh sizes and the associated parameters should be determined experimentally based on tomography data or roughness measurements. Additionally, it would be of great interest to modify the indicator function defined in Eq. (20) in order to consider uncertainty also inside the boundary and not only along the boundary. Moreover, the influence of a non-design region, located in the zones where boundary conditions are applied, on the optimal solution should also be considered as well as the integration of other uncertainties sources effects, e.g., material properties, or other constraints and cost functions, such as mass, failure criteria and buckling factors. Finally, in order to better take into account the designer's aversion to uncertainties and his confidence in the model and its parameters, epistemic uncertainty should also be mitigated.

## CRediT authorship contribution statement

**Noal Florian:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Formal analysis, Data curation, Conceptualization. **Alain Etienne:** Writing – review & editing, Validation, Supervision, Software, Resources, Conceptualization. **Marco Montemurro:** Writing – review & editing, Supervision, Methodology, Formal analysis. **Jean-Yves Dantan:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Julien Gardan:** Writing – review & editing, Supervision.
# Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

# Acknowledgements

This paper has been supported by the Agence National de la Recherche (ANR) for the project Robust Design For Additive Manufacturing (RobustAM) under the following reference ANR-21-CE10-0005, which is thankfully acknowledged by the authors. The authors also wish to express their gratitude to the open-source community.

# Data availability

Data will be made available on request.

# References

[1] Bendsøe MP, Sigmund O. Topology optimization: theory, methods, and applications. Springer Science & Business Media; 2003, http://dx.doi.org/10.1007/978-3-662-05086-6.

[2] Stolpe M, Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization. Struct Multidiscip Optim 2001;22(2):116-24. http://dx.doi.org/10.1007/s001580100129.

[3] Norato J, Bell B, Tortorelli D. A geometry projection method for continuum-based topology optimization with discrete elements. Comput Methods Appl Mech Engrg 2015;293:306-27. http://dx.doi.org/10.1016/j.cma.2015.05.005.

[4] Zhang W, Chen J, Zhu X, Zhou J, Xue D, Lei X, Guo X. Explicit three dimensional topology optimization via moving morphable void (mmv) approach. Comput Methods Appl Mech Engrg 2017;322:590-614. http://dx.doi.org/10.1016/j.cma.2017.05.002.

[5] Giulio Costa MM, Pailhès J. Nurbs hyper-surfaces for 3d topology optimization problems. Mech Adv Mater Struct 2021;28(7):665-84. http://dx.doi.org/10.1080/15376494.2019.1582826.

[6] Montemurro M, Mas A, eddine Zerrouq S. Topology and anisotropy optimization of continua using non-uniform rational basis spline entities. Comput Methods Appl Mech Engrg 2024;420:116714. http://dx.doi.org/10.1016/j.cma.2023.116714.

[7] Allaire G, Jouve F, Toader A-M. Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004;194(1):363-93. http://dx.doi.org/10.1016/j.jcp.2003.09.032.

[8] Reddy K SN, Ferguson I, Frecker M, Simpson TW, Dickman CJ. Topology optimization software for additive manufacturing: A review of current capabilities and a real-world example. In: Volume 2A: 42nd design automation conference of international design engineering technical conferences and computers and information in engineering conference. American Society of Mechanical Engineers; 2016, http://dx.doi.org/10.1115/DETC2016-59718.

[9] Venturi F, Taylor R. Additive manufacturing in the context of repeatability and reliability. J Mater Eng Perform 2023;32(15):6589-609. http://dx.doi.org/10.1007/s11665-023-07897-3.

[10] Ma Y-Z, Li C-X, Wang Y-Y, Zhang Z-Y, Li H-S, Ding A-N, Rui X-T. Robust design optimization of a multi-body system with aleatory and epistemic uncertainty. Reliab Eng Syst Saf 2024;245:110029. http://dx.doi.org/10.1016/j.ress.2024.110029.

[11] Schuëller G, Jensen H. Computational methods in optimization considering uncertainties - an overview. Comput Methods Appl Mech Engrg 2008;198(1):2-13. http://dx.doi.org/10.1016/j.cma.2008.05.004, Computational Methods in Optimization Considering Uncertainties.

[12] Roy CJ, Oberkampf WL. A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing. Comput Methods Appl Mech Engrg 2011;200(25):2131-44. http://dx.doi.org/10.1016/j.cma.2011.03.016.

[13] Peng X, Li J, Jiang S, Liu Z. Robust topology optimization of continuum structures with loading uncertainty using a perturbation method. Eng Optim 2018;50(4):584-98. http://dx.doi.org/10.1080/0305215X.2017.1337754.

[14] Wang X, Shi Y, Meng Z, Yang B, Long K. Uncertainty-oriented topology optimization of dynamic structures considering hybrid uncertainty of probability and random field. Reliab Eng Syst Saf 2024;110744. http://dx.doi.org/10.1016/j.ress.2024.110744.

[15] Meng Z, Wu Y, Wang X, Ren S, Yu B. Robust topology optimization methodology for continuum structures under probabilistic and fuzzy uncertainties. Internat J Numer Methods Engrg 2021;122(8):2095-111. http://dx.doi.org/10.1002/nme.6616.

[16] Javed A, Rout B. Tolerance range selection of topologically optimized structures with the effects of uncertainties of manufacturing process. Proc Inst Mech Eng C: J Mech Eng Sci 2014;228(18):3441-61. http://dx.doi.org/10.1177/0954406214528484.

[17] Sigmund O. Manufacturing tolerant topology optimization. Acta Mech Sin 2009;25(2):227-39. http://dx.doi.org/10.1007/s10409-009-0240-z.

[18] Wang F, Lazarov BS, Sigmund O, methods Onprojection. Convergence and robust formulations in topology optimization. Struct Multidiscip Optim 2011;43(6):767-84. http://dx.doi.org/10.1007/s00158-010-0602-y.

[19] da Silva GA, Beck AT, Sigmund O. Topology optimization of compliant mechanisms with stress constraints and manufacturing error robustness. Comput Methods Appl Mech Engrg 2019;354:397-421. http://dx.doi.org/10.1016/j.cma.2019.05.046.

[20] Frederiksen AH, Sigmund O, Poulios K. Topology optimization of self-contacting structures. Comput Mech 2024;73(4):967-81. http://dx.doi.org/10.1007/s00466-023-02396-7.

[21] Andreassen E, Clausen A, Schevenels M, Lazarov BS, Sigmund O. Efficient topology optimization in matlab using 88 lines of code. Struct Multidiscip Optim 2011;43(1):1-16. http://dx.doi.org/10.1007/s00158-010-0594-7.

[22] Ferrari F, Sigmund O. A new generation 99 line matlab code for compliance topology optimization and its extension to 3d. Struct Multidiscip Optim 2020;62(4):2211-28. http://dx.doi.org/10.1007/s00158-020-02629-w.

[23] Lazarov BS, Schevenels M, Sigmund O. Robust design of large-displacement compliant mechanisms. Mech Sci 2011;2(2):175-82. http://dx.doi.org/10.5194/ms-2-175-2011.

[24] Schevenels M, Lazarov B, Sigmund O. Robust topology optimization accounting for spatially varying manufacturing errors. Comput Methods Appl Mech Engrg 2011;200(49):3613-27. http://dx.doi.org/10.1016/j.cma.2011.08.006.

[25] Lazarov BS, Schevenels M, Sigmund O. Topology optimization considering material and geometric uncertainties using stochastic collocation methods. Struct Multidiscip Optim 2012;46(4):597-612. http://dx.doi.org/10.1007/s00158-012-0791-7.

[26] Jansen M, Lombaert G, Diehl M, Lazarov BS, Sigmund O, Schevenels M. Robust topology optimization accounting for misplacement of material. Struct Multidiscip Optim 2013;47(3):317-33. http://dx.doi.org/10.1007/s00158-012-0835-z.

[27] Jansen M, Lombaert G, Schevenels M. Robust topology optimization of structures with imperfect geometry based on geometric nonlinear analysis. Comput Methods Appl Mech Engrg 2015;285:452-67. http://dx.doi.org/10.1016/j.cma.2014.11.028.

[28] Chen S, Chen W. A new level-set based approach to shape and topology optimization under geometric uncertainty. Struct Multidiscip Optim 2011;44(1):1-18. http://dx.doi.org/10.1007/s00158-011-0660-9.

[29] Sato Y, Izui K, Yamada T, Nishiwaki S, Ito M, Kogiso N. Reliability-based topology optimization under shape uncertainty modeled in eulerian description. Struct Multidiscip Optim 2019;59(1):75-91. http://dx.doi.org/10.1007/s00158-018-2051-y.

[30] Xia M, Sun D, Yang S. A novel methodology for robust topology optimization considering manufacturing errors and topology deviations. IEEE Trans Magn 2022;58(9):1-4. http://dx.doi.org/10.1109/TMAG.2022.3172854.

[31] Dantan J-Y, Huang Z, Goka E, Homri L, Etienne A, Bonnet N, Rivette M. Geometrical variations management for additive manufactured product. CIRP Ann 2017;66(1):161-4. http://dx.doi.org/10.1016/j.cirp.2017.04.034.

[32] Zhu Z, Qiao L, Anwer N. An improved tolerance analysis method based on skin model shapes of planar parts. Procedia CIRP 2016;56:237-42. http://dx.doi.org/10.1016/j.procir.2016.10.072, the 9th International Conference on Digital Enterprise Technology - Intelligent Manufacturing in the Knowledge Economy Era.

[33] Moroni G, Petrò S, Polini W. Geometrical product specification and verification in additive manufacturing. CIRP Ann 2017;66(1):157-60. http://dx.doi.org/10.1016/j.cirp.2017.04.043.

[34] Moroni G, Petrò S. Segmentation-free geometrical verification of additively manufactured components by x-ray computed tomography. CIRP Ann 2018;67(1):519-22. http://dx.doi.org/10.1016/j.cirp.2018.04.011.

[35] Gorgulianslan RM, Grandhi RV, Choi H-J, Choi S-K. Prediction assessment and validation of multiscale models for additively manufactured lattice structures under uncertainty. J Mech Sci Technol 2019;33(3):1365-79. http://dx.doi.org/10.1007/s12206-019-0238-9.

[36] Rozvany GIN. A critical review of established methods of structural topology optimization. Struct Multidiscip Optim 2009;37(3):217-37. http://dx.doi.org/10.1007/s00158-007-0217-0.

[37] Montemurro M. On the structural stiffness maximisation of anisotropic continua under inhomogeneous neumann-dirichlet boundary conditions. Compos Struct 2022;287:115289. http://dx.doi.org/10.1016/j.compstruct.2022.115289.

[38] Svanberg K. The method of moving asymptotes—a new method for structural optimization. Internat J Numer Methods Engrg 1987;24(2):359-73. http://dx.doi.org/10.1002/nme.1620240207.

[39] Johnson SG. The NLopt nonlinear-optimization package. 2007, https://github.com/stevengj/nlopt.
[40] Helton J, Johnson J, Oberkampf W. An exploration of alternative approaches to the representation of uncertainty in model predictions. Reliab Eng Syst Saf 2004;85(1):39-71. http://dx.doi.org/10.1016/j.ress.2004.03.025, Alternative Representations of Epistemic Uncertainty.

[41] Beyer H-G, Sendhoff B. Robust optimization – a comprehensive survey. Comput Methods Appl Mech Engrg 2007;196(33):3190-218. http://dx.doi.org/10.1016/j.cma.2007.03.003.

[42] Torres AP, Warner JE, Aguiló MA, Guest JK. Robust topology optimization under loading uncertainties via stochastic reduced order models. Internat J Numer Methods Engrg 2021;122(20):5718-43. http://dx.doi.org/10.1002/nme.6770.

[43] Tony F. Chan GHG, Leveque RJ. Algorithms for computing the sample variance: Analysis and recommendations. Amer Statist 1983;37(3):242-7. http://dx.doi.org/10.1080/00031305.1983.10483115.

[44] Gorelick M, Ozsváltd. High performance python: practical performant programming for humans. O'Reilly Media; 2020, URL https://www.oreilly.com/library/view/high-performance-python/9781492055013/.

[45] De S, Hampton J, Maute K, Doostan A. Topology optimization under uncertainty using a stochastic gradient-based approach. Struct Multidiscip Optim 2020;62(5):2255-78. http://dx.doi.org/10.1007/s00158-020-02599-z.

[46] Bagul YJ, Chesneau C. Sigmoid functions for the smooth approximation to the absolute value function. Moroc J Pure Appl Anal 2021;7(1):12-9. http://dx.doi.org/10.2478/mjpaa-2021-0002.

[47] Bertolino G, Montemurro M. Two-scale topology optimisation of cellular materials under mixed boundary conditions. Int J Mech Sci 2022;216:106961. http://dx.doi.org/10.1016/j.ijmecsci.2021.106961.

[48] Geuzaine C, Remacle J-F. Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities. Internat J Numer Methods Engrg 2009;79(11):1309-31. http://dx.doi.org/10.1002/nme.2579.

[49] Baratta IA, Dean JP, Dokken JS, Habera M, Hale JS, Richardson CN, Rognes ME, Scroggs MW, Sime N, Wells GN. DOLFINx: the next generation FEniCS problem solving environment. 2023, http://dx.doi.org/10.5281/zenodo.10447666, preprint.

[50] Python Software Foundation. Python language reference, version: 3.10.12. 2023, http://www.python.org.

[51] Karanki D, Rahman S, Dang V, Zerkak O. Epistemic and aleatory uncertainties in integrated deterministic and probabilistic safety assessment: Tradeoff between accuracy and accident simulations. Reliab Eng Syst Saf 2017;162:91-102. http://dx.doi.org/10.1016/j.ress.2017.01.015.