Structural Design Using Laplacian Shells

TL;DR

This paper presents a novel method for designing lightweight, structurally robust shell objects by iteratively adjusting local thickness based on stress analysis, ensuring smooth, intersection-free boundaries through Laplace's equation solutions.

Contribution

It introduces a shape parametrization using Laplace's equation and a gradient-free optimization approach for efficient, intersection-free shell design under complex forces.

Findings

Successfully designed lightweight shells for various 3D models

Validated structural robustness through physical experiments

Ensured smooth, intersection-free boundaries in shell shapes

Abstract

We introduce a method to design lightweight shell objects that are structurally robust under the external forces they may experience during use. Given an input 3D model and a general description of the external forces, our algorithm generates a structurally-sound minimum weight shell object. Our approach works by altering the local shell thickness repeatedly based on the stresses that develop inside the object. A key issue in shell design is that large thickness values might result in self-intersections on the inner boundary creating a significant computational challenge during optimization. To address this, we propose a shape parametrization based on the solution to the Laplace's equation that guarantees smooth and intersection-free shell boundaries. Combined with our gradient-free optimization algorithm, our method provides a practical solution to the structural design of hollow objects with a single inner cavity. We demonstrate our method on a variety of problems with arbitrary 3D models under complex force configurations and validate its performance with physical experiments.

Figures

Tables

Table 1. Table 1: Performance of our algorithm on a variety of models. Columns 2 to 5 are the number of finite elements in the volumetric mesh and the number of vertices on the boundary, skeleton and remaining regions, respectively. Column 6 reports the number of force configurations for which the shell object is optimized. In all examples, the target factor of safety is 90 % percent 90 90\% of that of the fully solid versions.

Model	Elements	Vertices			Force	Iteration	Time [s]	Volume [cm3]
		Boundary	Skeleton	Internal	Configurations			Initial	Optimized	Reduction [%]
Cactus	43727	1856	233	5858	1	29	11.60	68.99	26.51	61.6
Sea Horse	99054	7502	105	11814	1	42	157.23	6.42	1.48	77.0
Deer Head	89521	7500	834	8969	20	47	92.68	27.80	9.14	67.1
Lilium	41962	3389	2085	4513	3	37	24.20	-	91.35	-
Octopus	130896	12502	893	14035	4	53	70.69	75.32	26.00	65.5
Test Beam	87284	6000	976	9730	2	38	76.84	20.38	9.32	54.3
Pitcher	61881	4997	735	6703	820	92	1788.39	175.85	31.09	82.3
Spot	56504	2930	310	7247	2189	50	1037.77	367.75	132.04	64.1
Mug	48604	3000	2679	3575	2510	76	1212.17	104.12	46.26	55.6
Shark	69243	5757	680	7522	4282	33	1100.37	66.45	12.79	80.8

Equations

$$\left[\begin{array}{ccc}{\mathcal{L}}_{b,b}& {\mathcal{L}}_{b,s}& {\mathcal{L}}_{b,in}\\ {\mathcal{L}}_{s,b}& {\mathcal{L}}_{s,s}& {\mathcal{L}}_{s,in}\\ {\mathcal{L}}_{in,b}& {\mathcal{L}}_{in,s}& {\mathcal{L}}_{in,in}\end{array}\right]\left[\begin{array}{c}{\mathbf{T}}_b\\ {\mathbf{T}}_s\\ {\mathbf{T}}_{in}\end{array}\right]=\left[\begin{array}{c}\cdot \\ \cdot \\ 0\end{array}\right].$$

$${\mathcal{L}}_{in,in}{\mathbf{T}}_{in}=-{\mathcal{L}}_{in,b}{\mathbf{T}}_B-{\mathcal{L}}_{in,s}{\mathbf{T}}_S.$$

$$\underset{{\mathbf{T}}_B}{\mathrm{minimize}}$$

$$\mathbf{K}({\mathbf{T}}_B){\mathbf{u}}_l={\mathbf{f}}_l\forall l\in {\mathcal{B}}_L,$$

$${\sigma }_{cr}({\mathbf{T}}_B)\le {\sigma }_y/k,$$

$$T_c\le {\mathbf{T}}_B\le T_u$$

$${\mathbf{K}}_e={\mathbf{K}}_e^{void}+{\rho }_e^{\beta }({\mathbf{K}}_e^{solid}-{\mathbf{K}}_e^{void}).$$

$${\mathbf{\sigma }}_l={\mathbf{C}}_{\mathbf{g}}\mathbf{B}{\mathbf{u}}_l,$$

$${\sigma }_i=\underset{l}{\max }\left({\sigma }_i^{vm}\right)\forall l\in {\mathcal{B}}_L,$$

$${\tau }_j=\sum _i{\sigma }_i\left(\frac{1/{d_{ij}}^q}{{\sum }_{j}1/{d_{ij}}^q}\right),\forall i\in \mathcal{V}\text{ and }\forall j\in \mathcal{B}$$

$$d_{ij}=\{\begin{array}{ll}dist(i,j),& \text{if }0\le dist(i,j)<R,\\ \infty ,& \text{otherwise}.\end{array}$$

$${{}^{t+1}{\mathbf{T}}_B}^{\prime }=max\left(\frac{T_{\Sigma }}{{\sum }_j{{\tau }_j}^{\kappa }}{\mathbf{\tau }}^{\kappa },1\right).$$

$${\mathbf{T}}_B=\alpha {}^{t+1}{\mathbf{T}}_B+(1-\alpha ){}^t{\mathbf{T}}_B,$$

$$\mathcal{T}=\mathcal{F}\mathcal{W}$$

$$\mathcal{W}=({\mathcal{F}}^T\mathcal{F}+r\mathbf{I}{)}^{-1}({\mathcal{F}}^T\mathcal{T})$$

Full text

{CCSXML} <ccs2012> <concept> <concept_id>10010147.10010371.10010396.10010402</concept_id> <concept_desc>Computing methodologies Shape analysis</concept_desc> <concept_significance>500</concept_significance> </concept> <concept> <concept_id>10010147.10010371.10010396.10010397</concept_id> <concept_desc>Computing methodologies Mesh models</concept_desc> <concept_significance>300</concept_significance> </concept> <concept> <concept_id>10010405.10010432.10010439.10010440</concept_id> <concept_desc>Applied computing Computer-aided design</concept_desc> <concept_significance>300</concept_significance> </concept> </ccs2012>

\ccsdesc

[500]Computing methodologies Shape analysis \ccsdesc[300]Computing methodologies Mesh models \ccsdesc[300]Applied computing Computer-aided design

\printccsdesc

††volume: 38††issue: 5

1 Introduction

As additively manufactured (AM) parts find their way into industrial applications as functional components, lightweighting and structural optimization methods have become prevalent in shape design [WWY*∗*13, LSZ*∗*14, CBNJ*∗*15, LSD*∗*16, ZKWG16]. In many such methods, the resulting shape usually contains complex internal structures, even for very simple loading configurations. These structures create multiple isolated volumes inside the object and thus, result in manufacturability as well as functionality problems (Figure 1). They prevent easy removal of material encapsulated inside these cavities such as internal supports in fused-deposition modeling (FDM) or excess material in selective laser sintering (SLS) and stereolithography (SLA). The object needs to be cut into pieces or perforated with multiple holes to access each cavity during post-processing to clear the encapsulated material [MAB*∗*15]. However, when extensive, these destructive processes may significantly degrade the mechanical performance and invalidate the structure completely. Furthermore, complex internal structures disrupt the inner cavity that could be used for functional purposes such as housing electronics, mechanisms or wiring [GZN*∗*15, SWT*∗*17]. Shell structures provide a compelling alternative where the object is interpreted as the solid enclosed between two 2-manifold surfaces creating a single connected cavity inside.

We propose a new structural optimization approach for designing minimum weight shell objects. Our approach takes as input (1) a 3D shape represented by its boundary surface mesh and (2) description of the external forces that the object may experience during its use, and produces a minimum weight 3D shell structure that withstands any of the prescribed force configurations (Figure Structural Design Using Laplacian Shells).

For shell structure design, a reasonable approach would be to compute the optimal structure by offsetting the original boundary surface inwards in the normal direction to create the inner boundary. However, this strategy fails to guarantee that the resulting offset surface is free of self-intersections. Even advanced approaches that aim to minimize such self-intersections through local alterations of the offset directions [MAB*∗*15, ZXZ*∗*17] do not guarantee self-intersection free offset surfaces for large thickness variations in high-curvature regions. Additionally, such an approach requires each candidate shell structure to be remeshed to perform finite element analysis (FEA) for evaluating the structural performance under the prescribed force configurations.

Our approach overcomes these challenges using a temperature field–a 3D scalar field obtained by solving the steady-state heat conduction–which serves as a proxy to the varying shell thickness within the shape boundary. This temperature field is the solution to the Laplace equation and the resulting offset surface is by design intersection-free even for large thickness variations. Additionally, the temperature field defined on a constant volumetric mesh allows easy and accurate transition to material distribution within the current shape hypothesis. This capability enables each step of the shape optimization to perform FEAs without requiring costly resmeshing operations. Our approach addresses classical structural design problems with fixed and known external forces [ACS*∗*11] as well as the more general class of problems ranging from multiple external force cases [JHM09] to force location uncertainties [UMK17]. Driven by a gradient-free optimization algorithm, our method provides a practical solution for designing robust hollow structures with a single inner cavity while preserving the appearance of the input model.

The main contributions of the proposed work are:

• •

a novel formulation for shell structure design involving structural mechanics,

• •

a heat-based shape parametrization method that allows large variations in thickness while guaranteeing self-intersection-free boundaries in the resulting structure,

• •

a gradient-free shape optimization approach to arbitrary 3D problems with complex force configurations including multiple loads as well as uncertainties in force locations.

2 Related Work

Our review is comprised of the studies that focus on design for fabrication, structural analysis, and shell object synthesis. We emphasize approaches involving structural optimization for additive fabrication.

Design for fabrication

A large body of work has investigated optimization techniques offering design aids to create shapes that meet the prescribed structural objectives and fabrication constraints [LEM*∗*17]. Recent examples include designing for deformation behavior [BBO*∗*10, STC*∗*13, PZM*∗*15, MZL*∗*17], meta-materials [ZKBT17, MHSL18, IFW*∗*16, ZSCM17], and lightweighting [Ben89, LSZ*∗*14, WLQ*∗*17]. Broader methods that can handle variety of requirements have also been investigated in [CLD*∗*13, CBNJ*∗*15, MHR*∗*16, SXZ*∗*17]. Our problem falls under the category of optimum weight structure design subject to external forces, however our approach aims to generate a specific family of geometries – shell objects.

Other design approaches has been explored to generate models that do not require internal support structures to fabricate them. Topology optimization methods [Lan16, ADE*∗*17], slice-based hollowing methods [WLW*∗*18] and specialized infill structures [WWZW16, LFC*∗*18] are investigated. Although these methods are well-suited for FDM or SLA printing processes, resulting complex structures with large number of disconnected cavities are not desirable for processes such as SLS or polyjet where excess material is typically unavoidable. Our shape parametrization is complementary to these approaches in that the presented support constraints may be facilitated in our optimization to generate shell models requiring minimum or possibly no internal supports to fabricate, thereby making it practical for all of the above mentioned fabrication processes.

Structural analysis

In design optimization problems involving structural mechanics, structural soundness of a candidate design is determined by stress and deformation analysis through FEA [Din86]. Simple, low-cost elemental structures such as trusses [SHOW02, Ros07, KIL*∗*18] and beams [WWY*∗*13, JTSW17, WYY*∗*18] as well as higher degree of freedom elements such as tetrahedral [LSZ*∗*14, SXZ*∗*17] and hexahedral [BS03, HTK13] elements are commonly employed based on the families of geometries in consideration. For thin-walled structures, more specialized elemental structures such as plate and shell elements are widely used [BBWR04]. They allow modeling the bending behavior of thin features accurately with fewer finite elements. This makes them suitable for applications including realistic cloth simulation [TWS06, LCBD*∗*18] and design for sheet metal fabrication [PP02]. However, they are limited to surface-like structures with one dimension (i.e., thickness) being significantly smaller than the other two [BBWR04]. In our approach, we use tetrahedral elements in FEA to analyze candidate shell designs. This allows our optimizer to explore in a larger design space by providing flexibility to generate shell models with wide range of thickness values.

Zhou et al.[ZPZ13] introduce worst-case structural analysis by extending modal analysis used in dynamic systems to static problems in order to identify weak regions of the structure that may fail under arbitrary force configurations. Langlois et al.[LSD*∗*16] present a stochastic finite element model to predict the failure probabilities of objects under scenarios where the loading is stochastic in nature (such as dropping and collisions). For problems with uncertainties in load direction and location spanning a small subspace, Schumacher et al.[SZB18] uses parametric descriptions of the uncertainties to determine the worst-case load.

For scenarios where uncertainty in force contact location is large and cannot be described parametrically, Ulu et al.[UMK17] present a data-driven technique to predict the stress distribution for each possible force configuration and determine the one that creates the highest stress within the shape. Wang et al.[WUSK18] improve this approach by incorporating a computationally tractable experimental design method to select data samples. Our analysis uses a similar approach to deal with contact location uncertainty.

Shell object synthesis

Support structures [HL18], moldable objects [NAI*∗*18] and models with desired deformation behavior [ZLW*∗*16] are among the recent applications of shell design in the context of 3D fabrication. A common approach in generating such shell structures is to offset the object’s surface to create the boundary of the inner cavity. For constant thickness offsets, Minkowski operations [CK10, MHCL15], boolean operations of volumetric primitives [PK08], and use of particle sets [MCS*∗*18] have been investigated.

When the input is a polygonal surface, an offset surface can be generated by shifting the original vertices in the normal direction [QS03]. Building upon this approach, Musialski et al.[MAB*∗*15] present a method to create varying thickness offset surfaces. Their optimizer alters the local shell thickness to minimize a set of objective functions. However, their method is not streamlined for design problems involving structural mechanics; (1) self-intersections occur in areas of high curvature as well as regions where offsets get too large and (2) each candidate shell needs to be remeshed for FEA. The latter problem is addressed in [MHR*∗*16] and [ZXZ*∗*17] by discretizing the shape using shell elements and adjusting their thicknesses without altering the volumetric mesh. While these methods are well-suited for optimization of thin-shells, the analysis accuracy suffers for large thickness values due to their selection of particular finite element type. Additionally, self-intersection problems may still persist for complex geometries with large thicknesses. Driven by similar motivations, we undertake both of these challenges in shell structure design.

Our approach is similar to traditional topology optimization methods [BS03, WWG03] in that we define the material distribution on a fixed volumetric mesh. However, in order to enforce the resulting geometry to be a shell and address the self-intersection problems while doing so, we use a smooth scalar field defined on this volumetric mesh as a proxy to the varying shell thickness values.

3 Algorithm

A shell structure can be interpreted as a solid enclosed between two 2-manifold surfaces–an outer boundary $\mathcal{B}$ (enclosing volume $\mathcal{R}$) and an inner boundary $\mathcal{B}_{i}$ (enclosing volume $\mathcal{R}_{i}$) . For a given $\mathcal{B}$, our design problem aims to find an optimal inner boundary $\mathcal{B}_{i}$ such that the resulting shell structure $\mathcal{M}=\mathcal{R}\setminus\mathcal{R}_{i}$ has as low a mass as possible while remaining robust under the forces it experiences. We next describe our approach to addressing this problem.

3.1 Overview

Figure 2 illustrates our approach. Given an input 3D shape and prescribed boundary and loading configurations (Figure 2(a)), our system optimizes the shell thickness by manipulating a scalar field, i.e., temperature field, defined on the boundary of the object. At each step, governed by the boundary temperature distribution, our system computes the resulting steady-state temperature field inside the object. In this temperature field, the isosurface at a certain cut-off temperature forms the inner boundary of the shell structure where the region encapsulated between the outer and the inner boundaries is solid and the region inside the inner boundary is treated as void. Then, for the current material distribution, our method computes the maximum stress encountered across the entire structure. Based on the stresses, optimization updates the boundary temperature distribution to minimize mass (Figure 2(b-c)). At the end, a minimum weight shell structure satisfying the imposed constraints is obtained (Figure 2(d)). Algorithm 1 summarizes our approach.

3.2 Temperature Field

We obtain the temperature field inside the object by creating a temperature gradient between the boundary surface $\mathcal{B}$ and the shape skeleton $\mathcal{S}$ (Figure 3). The skeleton is assigned a low temperature value and kept constant at all times while the boundary takes larger temperature values and they are adjusted to manipulate the steady state temperature field inside the object. The isosurface at a pre-defined cut-off temperature $T_{c}$ constitutes the inner boundary of the shell object. Here, the skeleton at the medial axis of the object serves as the inner bound for the isosurface and therefore, upper bound for the resulting shell thicknesses. We adopt the approach presented in [TAOZ12] to generate the skeleton. Due to its Laplacian-smoothing-based contraction process, this method provides smooth skeleton approximations that are less sensitive to surface details compared to other 3D skeleton generation algorithms [MAB*∗*15]. Please refer to [ATM*∗*16] for a detailed comparison of skeleton generations methods available in the literature.

The steady-state heat distribution within a volume can be computed by solving the Laplace equation $\nabla^{2}T=0$ subject to Dirichlet boundary conditions $T|_{\mathcal{B}}=T_{B}$ and $T|_{\mathcal{S}}=T_{S}$ where $T$ is the temperature function. In the discrete setting, the linear system can be written as $\boldsymbol{\mathcal{L}}\boldsymbol{T}=\boldsymbol{0}$ where $\boldsymbol{\mathcal{L}}\in\mathbb{R}^{n_{v}\times n_{v}}$ is the discrete laplacian of the volumetric mesh $\mathcal{V}$ with $n_{v}$ vertices and $\boldsymbol{T}$ is the vector of per-vertex temperature values. Note that $\mathcal{V}$ is composed of the boundary vertices $b\in\mathcal{B}$, the skeleton vertices $s\in\mathcal{S}$ and the internal vertices (i.e.,  Steiner points) $in\in\mathcal{I}$, i.e., $\mathcal{V}=\mathcal{B}\cup\mathcal{S}\cup\mathcal{I}$.

Reordering the vertices as $\boldsymbol{T}=[\boldsymbol{T}_{b},\boldsymbol{T}_{s},\boldsymbol{T}_{in}]^{\prime}$, the Laplace equation with Dirichlet boundary conditions can be reformulated as

[TABLE]

Temperature field for the internal vertices can then be computed by solving the bottom block in Equation (1) for $\boldsymbol{T}_{b}=\boldsymbol{T}_{B}$ and $\boldsymbol{T}_{s}=\boldsymbol{T}_{S}$,

[TABLE]

As the volumetric mesh, the boundary vertices and the skeleton vertices are constant during the process, it is sufficient to factorize the $\boldsymbol{\mathcal{L}}_{in,in}$ once as a preprocessing step and use it to solve for new $\boldsymbol{T}_{B}$.

One might argue that the skeleton temperatures could also be altered to provide more flexibility in shape parametrization as the cost in computing the steady state temperature distribution in Equation (2) does not change. It is technically correct that each unique Dirichlet boundary conditions defined by the combination of $\boldsymbol{T}_{B}$ and $\boldsymbol{T}_{S}$ result in a unique solution to the Laplace equation and therefore, a unique temperature distribution inside the object. However, the iso-surface at $T_{c}$ defining the inner boundary $\mathcal{B}_{i}$ of the resulting shell can be the same for different $\boldsymbol{T}_{B}$ and $\boldsymbol{T}_{S}$ combinations. To avoid such redundancy, in our approach, we keep $\boldsymbol{T}_{S}$ constant and modify only $\boldsymbol{T}_{B}$ to create candidate shell structures.

3.3 Shape Optimization

We tackle the following stress constrained mass minimization problem

[TABLE]

where $M$ is the total mass of the solid enclosed between $\mathcal{B}$ and $\mathcal{B}_{i}$ created at the isosurface $T_{c}$ and $T_{u}$ is a large number defining the upper bound for the boundary temperature. $\boldsymbol{K}$ is the global stiffness matrix, $\boldsymbol{f}_{l}$ and $\boldsymbol{u}_{l}$ represent the nodal force and displacement vectors when the external force is applied to surface node $l$ in a user-defined contact region $\mathcal{B}_{L}$. The object fails if the critical stress $\sigma_{cr}$; maximum stress observed among all loading configurations, ever exceeds the allowable stress $\sigma_{y}/k$, where $k$ is the safety factor and $\sigma_{y}$ is the yield strength of the material.

Finding the critical stress requires a set of FEAs–one for each distinct loading configuration. As the candidate shell structure generated at each step of the optimization does not necessarily comply with the volumetric mesh $\mathcal{V}$, each candidate shell needs to be remeshed to perform FEA in the conventional form. However, discontinuous nature of remeshing often results in computational challenges in the optimization problem. Moreover, in addition to its direct costs, remeshing introduces a computational overhead of repeated stiffness matrix computation. Namely, for each unique volumetric mesh, stiffness matrix needs to be recomputed to solve the linear elasticity problem. We address these problems by approaching it similar to topology optimization [BS03], material design [STC*∗*13, XLCB15] and microstructure design [SBR*∗*15]. Each element in the discretized domain $\mathcal{V}$ is associated with a density variable $\rho_{e}$ representing whether the element $e$ is full ($\rho_{e}=1$) or void ($\rho_{e}=0$). Elements that reside between $\mathcal{B}$ and $\mathcal{B}_{i}$, i.e., elements with all their vertices having temperatures larger than $T_{c}$, are deemed to be solid while elements that are completely inside $\mathcal{B}_{i}$, i.e., elements with all their vertices having temperatures smaller than $T_{c}$, are considered to be void. For the elements that lie on the isosurface where at least one of their vertices is inside $\mathcal{B}_{i}$, we adopt the common approach of allowing intermediate densities $\rho_{e}\in[0,1]$. We determine the density of an element based on the fraction of its volume residing between $\mathcal{B}$ and $\mathcal{B}_{i}$. Figure 4 illustrates the possible configurations that an element is assigned an intermediate density value. For an element with volume $V_{e}$, the density is computed as $\rho_{e}=V_{e}^{in}/V_{e}$ where $V_{e}^{in}$ is the portion between $\mathcal{B}$ and $\mathcal{B}_{i}$. Assuming linear isotropic materials and small deformations, the elemental stiffness matrix $\boldsymbol{K}_{e}$ can be related to $\rho_{e}$ and the stiffness matrix for base material $\boldsymbol{K}_{e}^{solid}$ as

[TABLE]

Here, $\boldsymbol{K}_{e}^{void}=\epsilon\boldsymbol{K}_{e}^{solid}$ is the stiffness matrix assigned to the void regions to avoid singularities in FEA. We penalize the intermediate density values using a penalization factor $\beta$ [Ben89]. We use $\epsilon=10^{-8}$ and $\beta=3$. For the volumetric mesh $\mathcal{V}$ with $m$ elements, one can assemble $\rho_{e}(\boldsymbol{T}_{B})$ into vector $\boldsymbol{\rho(\boldsymbol{T}_{B})}\in\mathbb{R}^{m}$ and construct the global stiffness matrix $\boldsymbol{K}(\boldsymbol{\rho})$ in order to determine the displacements $\boldsymbol{u}_{l}$ from $\boldsymbol{K}\boldsymbol{u}_{l}=\boldsymbol{f}_{l}$. Then, the stress-displacement relationship can be written as

[TABLE]

where $\boldsymbol{\sigma}_{l}\in\mathbb{R}^{6m}$ captures the unique six elements of the elemental stress tensor and $\boldsymbol{B}$ is the strain-displacement matrix that depends only on the elements’ rest shapes. Block-diagonal matrix $\boldsymbol{C}_{g}\in\mathbb{R}^{6m\times 6m}$ is constructed with elemental elasticity tensors $\boldsymbol{C}_{e}(\boldsymbol{\rho})$ on the diagonal. For each element, $\boldsymbol{C}_{e}$ can be computed analogous to $\boldsymbol{K}_{e}$ in Equation (4). In our formulation, we use 10-node quadratic tetrahedral elements. Note that vertices of $\mathcal{V}$ constitute the corner nodes of the quadratic elements and additional middle nodes are generated for FEA purposes only.

The approach formulated in Equation (4)-(5) is useful because it preserves the same discretization throughout the optimization. Additionally, it complies well with our heat-based shell object parametrization in Equation (1)-(2) by sharing the same discretization.

We solve the optimization problem in Equation (3) using a gradient-free iterative approach inspired by [BT15]. Figure 5 illustrates a single iteration in our method. Given the material distribution at step $t$ of the optimization, we start by computing distribution of maximum stresses across $\mathcal{V}$, $\boldsymbol{\sigma}$ as

[TABLE]

where $\sigma_{i}$ is the maximum stress that node $i$ experiences for all possible force configurations. Here, for each loading configuration $l$, we compute the von Mises stress at node $i$ $\sigma_{i}^{vm}$ by extrapolating the corresponding components of $\boldsymbol{\sigma}_{l}$ using the element shape functions. Note that $\max(\boldsymbol{\sigma})$ constitutes $\sigma_{cr}$ in Equation (3). We then calculate the effective boundary stress (EBS) by projecting the stresses at the internal nodes onto the boundary nodes only. The effective boundary stress distribution allows us to estimate the boundary temperature, thereby the shell thickness for the next step $t+1$ of the optimization. Thus, the higher EBS at a vertex, the more critical it is deemed and therefore the larger temperature it is assigned.

Effective Boundary Stress

We compute the effective boundary stress by distributing the stress at internal vertices to their closest boundary vertices proportional to the distance between them. The idea behind such a heuristic is based on our observation that stresses created by external forces at a certain region of a shell object can be manipulated effectively by adjusting (i.e., thickening or thinning) the shell around that particular region. Therefore, the distance-based projection allows us to adjust the boundary node temperatures and thus, the local thicknesses accordingly. EBS at boundary vertex $j$, can be defined as

[TABLE]

where

[TABLE]

Here, $dist(i,j)$ is the graph distance between vertices $i$ and $j$, $R$ is the influence depth and $q$ is the influence exponent. The main reason behind using graph distance as opposed to Euclidean distance here is that the graph distance takes the shape boundaries into account by constraining the paths between vertices to be within the shape. Figure 6 illustrates the contributions of internal vertices on boundary vertices for two example cases. Note that $R$ should be selected large enough that each internal node should be inside of the influence region of at least one boundary node. We use $R=10$ and $q=3$ for our examples.

Boundary Temperature

Algorithm 2 describes how the boundary temperature $\boldsymbol{T}_{B}$ is updated at each step of the optimization. As the total mass $M$ of the shell object is linearly proportional to $\boldsymbol{T}_{B}$, the change in mass can be controlled by the change in a temperature budget $T_{\Sigma}$; a maximum cap on the total sum of the boundary temperatures at any step of the optimization. Hence, $T_{\Sigma}$ is adjusted by a step size $h$ based on the critical stress value. $T_{\Sigma}$ is increased for $\sigma_{cr}$ larger than the allowable stress and it is reduced otherwise. The temperature budget is then distributed among the boundary vertices proportional to their corresponding EBS as

[TABLE]

Here, $\boldsymbol{\tau}\in\mathbb{R}^{n_{b}}$ is a vector storing EBS values where $n_{b}$ denotes the number of vertices in $\mathcal{B}$ and $\kappa$ is the proportion exponent. We use $\kappa=5.0$ in our algorithm.

As a final step, we linearly blend boundary temperatures from the previous step with the current step as

[TABLE]

where $\alpha$ is the blending factor. This helps attenuate drastic jumps between two consecutive iterations that could be created due to local stress concentrations and result in smooth transitions throughout the optimization. We found out that $\alpha=0.5$ works well for all our examples.

As the optimization approaches toward the final result, we observe frequent direction changes in $T_{\Sigma}$ across consecutive iterations. This usually results in back and forth jumps between two states. We address this problem by halving the step size in every direction change. Hence, the step size provides an effective indicator for convergence. We initialize the optimization with $h=0.1$ and terminate it when $h$ is smaller than $10^{-8}$.

4 Force Location Uncertainty

In the case of a small number of force configurations, the maximum stress distribution in Equation (6) and $\sigma_{cr}$ can be calculated by performing an FEA per force configuration. However, for cases where a large number of force configurations need to be considered, a brute force approach in this way becomes restrictively impractical. A vast majority of such problems can be generalized under force location uncertainties where the external forces’ contact locations exhibit significant variations during the use of the object.

For problems with force location uncertainties, we adopt a similar approach to [UMK17] in estimating the maximum stress distribution $\boldsymbol{\sigma}\in\mathbb{R}^{n_{v}}$ and $\sigma_{cr}$. We perform FEAs for only a small number of force configurations and use them to construct a mapping between the nodal forces and the resulting stress distributions and estimate the stress distributions for the remaining force configurations.

Suppose uncertainty in the external force locations is defined such that the normal forces are allowed to make contact within a user-specified union of contact regions $\mathcal{B}_{L}\subseteq\mathcal{B}$. We start by uniformly sampling a number of force instants on $\mathcal{B}_{L}$ such that the geodesic distance between the closest samples is maximized. Here, we use approximate geodesic distances [CWW13] for computational efficiency. We assume that each force is distributed to a small circular area around the contact point to avoid stress singularities. Therefore, it can be represented as a sparse vector $\boldsymbol{f}_{l}^{\prime}$ of size $n_{b}$ where each element corresponding to a boundary surface node that lies inside this circular area takes a non-zero value equal to the magnitude of the nodal force component. Von Mises stress corresponding to each force instant forms a vector $\boldsymbol{\sigma}^{vm}$ of size $n_{v}$.

Using a small number of training force samples $p$ (i.e., $p\ll n_{b}\ll n_{v}$), it is not possible to represent the relationship between these two high dimensional data using a simple mapping function. Thus, we project both the force and von Mises stress data onto lower dimensional spaces; we use the Laplacian basis for forces and the principal components for stresses. We compute the Laplacian basis functions as the first $s$ eigenvectors of the surface graph Laplacian $\boldsymbol{\mathcal{L}}_{\mathcal{B}}\in\mathbb{R}^{n_{b}\times n_{b}}$. In all our examples, we use $s=15$ eigenvectors to construct our lower dimensional Laplacian basis. For stresses, we use principal component analysis (PCA) to obtain the lower dimensional basis. A PCA on the stress data results in $(p-1)$ principal vectors.

Lower dimensional representations of the force instants and the corresponding stresses allow us to construct a simple mapping between these two spaces. We build the following quadratic regression model with L2 regularization

[TABLE]

where $\boldsymbol{\mathcal{T}}$ is the lower dimensional representation of nodal von Mises stresses, $\boldsymbol{\mathcal{F}}$ is a matrix storing lower dimensional force representation including the quadratic terms and $\boldsymbol{\mathcal{W}}$ is the coefficient matrix that models the quadratic map. The coefficient matrix can then be computed as

[TABLE]

where $\boldsymbol{I}$ is identity matrix and $r$ controls regularization. Using this map, we estimate the nodal von Mises stresses $\boldsymbol{\sigma}^{vm}$ for each force instant. This allows us to compute an estimation for $\boldsymbol{\sigma}$ using Equation (6). Figure 7 illustrates $\boldsymbol{\sigma}$ obtained through our estimation in comparison to the ground truth for an example model.

Note that computed $\boldsymbol{\sigma}$ here is only an approximation of the actual values, thus cannot be used directly to obtain $\sigma_{cr}$. We adopt the hierarchical search approach described in [UMK17] to compute $\sigma_{cr}$ accurately and guarantee the structural soundness of the resulting shell under prescribed force location uncertainties.

5 Results and Discussion

5.1 Shape Parametrization

In our heat based shape parametrization, the shell structure is represented by a temperature field created on the boundary surface of the input mesh, i.e., the number of design parameters are identical to the number of boundary vertices. This provides our optimizer enough flexibility to perform localized changes with large variations on the shell thickness while preserving smoothness in the resulting internal surfaces. Figure 8 illustrates an example set of shells that can be achieved by our shape parametrization on a simple sphere model. Although the number of design parameters could be large for models with intricate surface details such as the octopus (Figure Structural Design Using Laplacian Shells) or the sea horse (Figure 11), our parametrization preserves a certain level of smoothness inherently due to the smooth temperature field created inside the object.

Another advantage of our shape parametrization method is that the resulting temperature field is easily transitioned into the density based representation from which the structural analysis can be constructed directly. This allows us to use the same volumetric mesh throughout the entire optimization without requiring costly remeshing operations. Moreover, during this transition, only a small number of intermediate density elements are created (Figure 5(a)), which is highly desirable in structural optimization as the analysis accuracy is only minimally affected.

Skeleton and Guarantee for a Single Hole

Skeletons used in our examples are shown in green in Figure 11 and Figure 12. Note that the skeletons are essentially a set of connected vertices in our formulation. Therefore, the variety of geometric representations including polylines, open or closed surfaces, meso-skeletons [TAOZ12] or combinations of these can be used as skeletons in our approach. In our examples, we employ curvilinear skeletons for the deer head, sea horse, cactus, spot, shark and octopus, a meso-skeleton for the mug, and an open surface mesh for the beam model.

As the skeleton defines the inner bound for the resulting shell, users may constrain the design space by manually manipulating the skeleton geometry. For example, in the pitcher model, we add a spherical surface to the curvilinear skeleton inside the inner chamber to set the minimum capacity in the final result. As the isosurface defining the inner boundary is restricted to reside between the skeleton and the outer boundary surface, the resulting structure is guaranteed to have the internal part of the skeleton hollow. This can be useful in designing dedicated housings to incorporate electronics or other instrumentation.

Importantly, because the solution to the heat equation is a harmonic function, the maximum principle [ABW01] guarantees that all voids are connected to the skeleton and all solid areas are connected to the outer boundary. This means that an optimized object will never have more voids than the number of connected components in the skeleton. To ensure the creation of a single hole, we restrict the number of skeletons to one in our approach.

Although the internal skeletons can be obtained easily for closed meshes, they are not well-defined for open surfaces. Figure 9 demonstrates an example case where the input to our algorithm is an open surface mesh and a solid shell structure is to be created by thickening the input surface. The surface is fixed at its boundary and three different loads are applied to top of the surface. In this case, we create a surface to serve as the skeleton (shown in green) in our algorithm by projecting the mesh to a lower plane. Then, the volumetric mesh is created in the space enclosed between these two surfaces, on which the temperature field is solved. The resulting shell boundary is shown in yellow. The combination of this boundary and the input surface constitutes the resulting solid shell.

Comparison

Figure 10 compares the performance of our shape parametrization with the reduced order offset surfaces approach of Musialski et al. [MAB*∗*15]. The offset directions in this method are computed such that the self-intersections are avoided as much as possible. However, high curvature regions are still susceptible to such problems for large thickness values. Therefore, there is no guarantee that the candidate shell structures generated during the optimization will be free of self-intersections. Although self-intersections might be tolerated for certain optimization problems, design optimization involving structural mechanics are highly sensitive to such complications as each candidate design needs to be evaluated by a set of FEAs. In our method, because the inner shell boundary is a level-set of a smooth and continuous temperature field, the resulting shell is guaranteed to be self-intersection free. Note that the inner boundaries shown in Figure 10 are results of random sets of feasible design variables. For the offset surfaces, we use $36$ manifold harmonics to generate the results. In both cases, we use the same skeleton structure generated using [TAOZ12]. Parts of the skeleton that are outside of the boundary are manually corrected.

5.2 Structural Optimization

5.2.1 Multiple Problem Configurations

Figure 11 illustrates the results of our method on various 3D problems where the shell structures are required to withstand a number of fixed and known force configurations. Displacement boundary conditions are shown in blue and the force contact regions are shown in red. Our shell structure optimization algorithm detects the failure-prone parts of the objects and adjusts the local shell thickness accordingly. For example, in the deer head model, the branching regions of the antler where the stresses are high due to large bending moments are thickened while the rest of the head is left rather hollow. Similarly, in the cactus model, since the trunk is already structurally sound under a load applied to the long arm, only the region where the arm branches out from the trunk is optimized to have a large shell thickness.

Our formulation in Equation (6) also allows optimization of structures for different displacement boundary conditions as well as for any arbitrary loading configuration at a cost of increased computational complexity. In such cases, each set of Dirichlet boundary conditions needs to be incorporated into the stiffness matrix and a linear solve needs to be performed per problem configuration at each step of the optimization. Figure 14 illustrates an example case where the beam is optimized to withstand two sets of boundary conditions as well as loading configurations.

In all of the examples, given the boundary and loading conditions, the target factor of safety for optimization in Equation (3) is set to be $90\%$ of the factor of safety of the fully solid, original models. This effectively means the optimization tries to preserve at least $90\%$ of the structural strength of the fully solid models. Under this relatively stringent condition, our algorithm achieves $54\%$ to $77\%$ mass reduction. Table 1 reports the reduction in volume together with various other metrics relevant to these models.

5.2.2 Force Location Uncertainty

In Figure 12, we demonstrate the results of our method on a set of example problems where there are uncertainties in the force configurations. In these examples, the displacement boundary conditions are assumed to be known and fixed and the forces are compressive and normal to the boundary surface.

Similar to the previous cases, a thicker shell is obtained around the parts of the objects where high stresses may develop under the given force location uncertainties. Notice the handles of the pitcher and mug models or the tail connecting the shark body to the base plate.

Although the computational cost is higher and the optimization takes longer to converge compared to the simpler loading configurations mentioned earlier, a similar volume reduction performance is obtained. For the problems with force location uncertainties, we achieved $55\%$ to $82\%$ reduction in volume (Table 1).

Comparison

In Figure 13, we compare our approach with the build-to-last method [LSZ*∗*14] and the lightweighting method in [UMK17]. As the build-to-last method creates honeycomb-like structures, it tends to generate large number of unconnected internal voids. For the shark model, the number of such components turns out to be $48$. Therefore, unless a special technique is employed to avoid internal supports in 3D printing, post-processing could be very cumbersome and costly. Our approach, however, generates a shell structure with a single connected hole.

In comparison to [UMK17], our optimum result weighs $14.7\%$ less when it is optimized for $20N$ load that is allowed to be applied anywhere on the surface of the shark in normal direction. As their approach is a reduced order method, the performance is affected highly by the number of material basis being used. For $15$ material modes, the optimization can perform more global alterations than local changes. On the other hand, our method can perform more localized alterations using all the vertex temperatures and thereby hollowing out both the fins as well as the tail.

5.3 Validation and Performance

Fabrication

We 3D printed our optimum shells using polyjet (OBJET Connex 350), SLS (ProMaker P2000 HT) and FDM (Strarasys F170) technologies. As our method creates a single connected void inside the object, internal supports (or excess material) can be removed by piercing a single access hole to empty the material trapped inside the shell. For models printed using SLS and polyjet, we used water jet for this purpose. However, for cases where a large amount of material can be accessed through only a small passage (such us the shark or deer head models), we observed that such an approach might not be effective. In those cases, a 3D printing technology allowing soluble support material is a more practical solution. We used soluble support in our FDM prints as an example (see accompanying video).

Physical Tests

Figure 14 illustrates a beam model we designed to perform tests to physically evaluate the strength of the 3D printed model. The model is required to withstand $115N$ force in three-point bending and tensile loading configurations. As the stress is concentrated in the middle section under three-point bending compared to the more uniformly distributed stresses in tensile loading case, our algorithm thickens the middle while keeping the ends as thin as possible.

To validate the optimization result, we performed three-point bending and tensile tests on our optimized beam model. We used an INSTRON universal testing machine for our tests. For comparison, we performed the same two tests on an identically weighing uniform thickness beam. Figure 15 shows the test setups as well as the resulting force-displacement curves. For the uniform thickness model, we measured $96.3N$ and $780.8N$ of failure (yielding) forces for three-point bending and tensile tests, respectively. Although its tensile performance is better, the uniform thickness model fails to satisfy the design constraint of $115N$ for the three-point bending case. Our optimum result, however, can withstand $153.4N$ and $433.7N$ forces for these cases, thereby satisfying the design constraint for both. Note that, for the same mass, our method shifts the material towards the center to improve the performance for the three-point bending test and satisfy the constraint while sacrificing its tensile performance (while satisfying both imposed constraints).

Performance

Table 1 shows the performance of our algorithm. We tested our method on a computer with an 3.7GHz Intel Xeon W-2145 CPU and 32GB memory. Our tests include the optimization of various 3D models for problem configurations of different complexities. As in most of the structural optimization approaches, FEAs constitute the majority of the computational cost in our approach. Therefore, the performance is heavily affected by the size of the volumetric mesh as well as the number of problem configurations. Sea horse (large number of elements) vs. cactus (small number of elements) highlights the impact of the mesh size on the optimization time. For the same number of force configurations, computation time per iteration increases $\sim 9.5\times$ on average when the number of elements is increased by $\sim 2\times$. Similarly, problems with force location uncertainty (Pitcher, Spot, Mug and Shark) take longer to optimize compared to problems with a single (Cactus, Sea Horse) or a small number of problem configurations (Deer Head, Lilium, Octopus and Test beam). However, our estimation based approach in determining $\boldsymbol{\sigma}$ and $\sigma_{cr}$ allows us to achieve $\sim 5\times$ acceleration on average over a brute force approach.

In all of our examples, the optimization converges under $100$ iterations. Figure 16 illustrates the convergence profiles for Cactus and Shark as examples of classical structural design problems with a single force configuration and problems with force location uncertainties, respectively. Note that larger step sizes used at the initial stages allow the optimizer to reach reasonable solutions quickly by removing large amount of material. In later stages, however, smaller step sizes result in more local and intricate alterations requiring larger number of iterations. Although, faster convergence can be achieved by relaxing the stopping criteria, this comes at the cost of limited flexibility in optimization.

5.4 Limitations and Future Work

Our algorithm can create very thin shells when required. However, the density values associated with the elements on such thin regions could be very low resulting in inaccurate stress calculations. We address this problem by enforcing a single layer of boundary elements to be solid at all times. The thickness of this layer, however, is largely dictated by the quality of the volumetric mesh. Linear elasticity model may fail to predict the nonlinear buckling failure modes of shell structures. In the future, our analysis could be extended or complemented with nonlinear corotational methods to capture such failure modes more accurately.

The temperature distribution in our algorithm depends on the geometry of the skeleton created. We found mean curvature skeletons to work well for organic shapes. However, for man-made shapes such as voluminous mechanical objects, the generated skeleton may limit the design space and restrict the quality of the resulting shell structure. For such cases, a manual adjustment tool such as a sketching interface would be beneficial to the user. For open meshes, we use a simple heuristic of projecting the input surface to a plane in creating the skeleton. Although this approach works well for smooth surfaces with low curvature, computation of the skeleton for complex open surfaces remains an open problem.

For very large boundary temperatures, our formulation is bound to create a small internal void around the skeleton (such as the spout of the pitcher model in Figure 12). Such small voids might result in manufacturability issues in 3D printing and more importantly can create stress concentration problems. A natural extension to our approach would be to limit the minimum hole size by thickening the skeleton using morphological operations.

6 Conclusion

We present a lightweight shell structure optimization method for 3D objects. We propose a heat based shape parametrization method to create shell structures with large thickness variations. With this method, we show that smooth internal surfaces can be created without self-intersections. A rapid transition between our heat based shape parametrization and the density based representation provides a practical solution to the computationally demanding design problem involving structural mechanics. Combined with the data-driven critical stress estimation approach, we demonstrate that our method can be applied to a generalized set of problems where there is uncertainty in the force locations. We evaluate the performance of our algorithm on a variety of 3D models. Our results show that significant mass reductions can be achieved by optimizing the shell thickness while ensuring that the object is structurally sound against a wide range of force configurations.

Acknowledgements

The authors would like to thank Nurcan Gecer Ulu, Walter Hsiao and Saigopal Nelaturi for insightful discussions and David Johnson for his help in 3D printing. We are grateful to the designers whose 3D models we used: kyoday (Thingiverse) for sea horse, yeg3d (Thingiverse) for deer head, xiaoxunyue2016 (Thingiverse) for octopus, HappyMoon (Thingiverse) for mug, Maninder Singh Chanay (Grabcad) for pitcher, Keenan Crane for spot, Sorkine-Hornung et al. for the cactus and Lu et al. for the shark. This work is partly funded by NCDMM America Makes Project #4058.