Source-linked AI summary

Robust topology optimization of isogeometric Kirchhoff--Love shells by the Speed Method

Philipp Eisenhardt, Fernando Hübner Scherer, Barbara Wohlmuth, Andrei Constantinescu

arXiv:2610.30389v1math.NAcond-mat.mtrl-scics.NA

TL;DR

The paper targets robust topology optimization of Kirchhoff–Love shells across non-matching multi-patch geometries, using a CAD-compatible, volume-based speed method and stochastic uncertainty models. Its multi-patch hemisphere example recovers a literature reinforcement pattern with fewer elements, and uncertain loading results show reduced compliance variability at increased expected compliance.

  • Problem

    The paper addresses robust topology optimization for Kirchhoff–Love shells across non-matching multi-patch geometries, where IGA naturally satisfies the shells’ C1-continuity requirement.

    Source: §1 · p. 2
  • Method

    A NURBS level set and speed-method volume-based derivative optimize CAD-compatible shell shapes, while Monte Carlo samples of Matérn Gaussian random fields model loading and material uncertainty.

    Source: §1 · p. 2, §2.3 · p. 15, §4 · p. 31
  • Results

    40 × 40 cubic NURBS elements per patch recover the multi-patch hemisphere reinforcement pattern of 128 × 128 × 1 solid elements; increasing α from 0 to α = 5 reduces double-parabola compliance standard deviation roughly 73%, with 6% higher expectation.

    Source: §4 · p. 30
  • Takeaways & Limitations

    The framework can also accommodate peak von Mises stress and Conditional Value-at-Risk objectives because the shape derivative enters through a realization-dependent scalar weight.

    Source: §4 · p. 30

    Monte Carlo sampling scales each iteration’s cost with sample size, although independent samples permit parallelization.

    Source: §4 · p. 31

Abstract

from arXiv · show

This paper presents a framework for robust topology optimization of Kirchhoff-Love (KL) multi-patch shells using isogeometric analysis and non-matching meshes and provides a CAD-compatible method. The material distribution is optimized via a NURBS-based level set method, enabling accurate boundary representation and seamless export to standard CAD/CAM workflows. The shape derivative is computed using the speed method, yielding a volume-based formulation that avoids boundary approximations and aligns with the isogeometric principle of coarse, high-order discretizations. This eliminates the need for local refinement, adaptive quadrature, or trimming as the descent direction for the compliance and volume objectives are obtained through a biharmonic equation. The framework is extended to robust optimization under uncertainty, accounting for stochastic loading and material properties via Monte Carlo sampling and Gaussian Random Fields. Numerical examples demonstrate the method's efficiency and robustness across single- and multi-patch shells, reproducing literature results with fewer elements than traditional FEM or refined IGA approaches. The developed methods and tools can be directly extended to problems involving stress-based optimization and risk-averse metrics. The numerical implementations are open-source and available on GitHub.

1. Introduction

The paper introduces CAD-compatible topology optimization for Kirchhoff–Love multi-patch shells, combining NURBS level sets with a volume-based speed-method shape derivative. It situates this approach among existing shell optimization methods and extends it to uncertainty.

1. Introduction

  • IGA naturally supplies the C1 continuity required by KL shells and connects with finite-element frameworks through Bézier extraction.KL shells are used in engineering structures, biological tissues, and models of complex material behavior.Source: §1 · p. 2
  • The framework combines NURBS level sets and a volume-based speed-method derivative for CAD-compatible optimization of KL multi-patch shells with non-matching meshes.The formulation avoids boundary approximations, trimming, local refinement, and adaptive quadrature, while obtaining descent directions through a biharmonic equation.Source: §1 · p. 2
  • Earlier KL shell topology optimization used boundary-based shape derivatives requiring fine refinement, while this paper also addresses loading and material uncertainty.Prior work includes parametric, Free Form Deformation, Adaptive Bubble, and level-set approaches; the paper presents single- and multi-patch geometries in Section 3.Source: §1 · p. 2

2. Derivation of the Shape Derivatives for Kirchhoff–Love shells by the speed method

This section derives speed-method shape sensitivities for Kirchhoff–Love shell optimization, covering compliance and volume objectives and their level-set, isogeometric implementation. It extends the deterministic procedure to robust optimization with stochastic loading and material properties using scaled deterministic derivatives.

2. Derivation of the Shape Derivatives for Kirchhoff–Love shells by the speed method

  • Under the Kirchhoff–Love kinematic assumption, mid-surface displacement characterizes shell deformation, and small-displacement strains split into membrane and bending components.The shell geometry is represented by a parametrized surface mapping between parametric and material domains.Source: §2 · pp. 3–4

2.1. Shape optimization problem

  • The optimization seeks a shell domain that improves compliance while accounting for volume, with shape sensitivities derived by perturbing the domain and differentiating the pulled-back state and cost.A Lagrangian formulation accounts for changes in displacement, and its optimality conditions provide the direct and adjoint equations used in the compliance derivative.Source: §2.1 · pp. 5, 7, 8
  • Biharmonic regularization supplies the descent direction because bending strains involve second derivatives of velocity; the level set evolves by Hamilton–Jacobi advection and is reinitialized every five iterations.The material domain is represented implicitly by a level set, and the speed-method shape derivative drives its evolution.Source: §2.1 · pp. 10–11

2.2. Numerical Optimization Algorithm

  • IGA encodes geometry and shape functions with NURBS; multi-patch coupling enforces C1 continuity of displacement and descent direction through penalty methods.The implementation uses FEniCS with tIGAr and PENGoLINS, and the associated numerical results and code are publicly available.Source: §2.2 · pp. 11–13
  • The optimization algorithm uses degree-3 NURBS for u and θ, while ϕ is a B-spline on the same discretization; color coding distinguishes deterministic and uncertainty-specific steps.Source: §2.2 · p. 14

2.3. Extension of the Shape Derivative to Loading and Material Uncertainties

  • Uncertainty in elastic moduli or external loads is modeled by Matérn-covariance Gaussian random fields, while robust compliance accounts for expected compliance and its variability using α ≥ 0.Source: §2.3 · p. 15
  • The robust shape derivative is the expectation of realization-scaled deterministic derivatives and replaces the deterministic derivative in the same regularization and advection workflow.Each uncertainty sample defines a mechanical problem with its own compliance.Source: §2.3 · pp. 15–16

3. Results and Discussion

The examples test the approach on deterministic and stochastic shell problems, including a curved multi-patch hemisphere. For the hemisphere, the method produces a symmetric reinforcement pattern and similar shapes to prior work with fewer elements.

3. Results and Discussion

  • The results section includes deterministic and stochastic examples, compares them with available benchmarks, and points to the GitHub repository for further examples.All physical dimensions are given in coherent SI units unless stated otherwise.Source: §3 · p. 17

Example 1. Multi-Patch Hemisphere

  • The hemisphere benchmark uses four shell patches on a radius-10 half-sphere with an 18° central hole, under distributed loading and clamping at the hole.The shell thickness is 0.02; the remaining boundary is free.Source: Example 1. Multi-Patch Hemisphere · p. 17, Figure 2 · p. 17
  • With 40 × 40 cubic NURBS KL elements, the hemisphere design achieved similar performance to ’s 128 × 128 × 1 solid IGA elements.The optimized pattern has rib-like supports and is relatively symmetric despite no symmetry constraints.Source: Figure 3 · p. 18, Example 1. Multi-Patch Hemisphere · p. 18
  • The hemisphere objective stabilizes after approximately 30 iterations.Figure 4 tracks objective, compliance, and volume fraction against iteration number.Source: Figure 4 · p. 18

Deterministic Point Loading

The double parabola examples examine how deterministic point-load direction affects optimized shell topology. The resulting supports differ across vertical, horizontal, and diagonal loading cases.

Deterministic Point Loading

  • The double parabola has side length 5 and height 0.6, with three concentrated loads applied at its center in vertical, horizontal, and diagonal directions.The models use cubic NURBS and the same geometry, though the volume Lagrange multiplier differs for vertical versus horizontal loads to keep volume fractions comparable.Source: Deterministic Point Loading · p. 19
  • Vertical loading yields cross-shaped supports reinforced along the z-axis, horizontal loading yields asymmetric load-aligned struts, and diagonal loading yields trusses connecting the load point to clamped boundaries.The optimized topologies correspond to the three distinct concentrated-load directions.Source: Deterministic Point Loading · p. 19, Figure 6 · p. 20

Uncertainty in Loading Direction

The study varies a robustness parameter for a double parabola subject to uncertain loading direction. Higher parameter values reduce compliance variance but increase expected compliance, yielding a risk–performance trade-off.

Uncertainty in Loading Direction

  • The optimization targets a 50% volume fraction for loading directions normally distributed around the diagonal load, with standard deviation 10°.The compliance-standard-deviation penalty parameter α varies from 0 to 5.Source: Uncertainty in Loading Direction · p. 20
  • Increasing α reduces compliance variance but increases expected compliance, showing the trade-off between these outcomes.Figure 8 plots expected compliance against compliance standard deviation.Source: Figure 8 · p. 21, Uncertainty in Loading Direction · p. 21
  • Higher α values produce thicker diagonal reinforcements and improve resistance to less probable force orientations normal to the expected loading diagonals.The load orientation varies with a Gaussian standard deviation of 10°; α ranges from 0 to 5.Source: Uncertainty in Loading Direction · p. 21, Figure 7 · p. 21

Deterministic Surface Loading

The surface-loading example applies a distributed load to the double-parabola geometry. Its optimized material distribution differs from the point-load design, favoring rigidity across the surface.

Deterministic Surface Loading

  • The example applies a downward distributed load of magnitude 12500 and uses the geometry and initialization of the preceding examples.The weak-phase ersatz parameter is set to 0.05 to ensure well-posedness.Source: Deterministic Surface Loading · p. 22
  • Distributed surface loading produces a more even material distribution than vertical point loading, which concentrates material near the center.The distributed-load design addresses rigidity across the whole surface.Source: Deterministic Surface Loading · p. 22, Figure 10 · p. 23

Distributed Surface Loading by GRF

The double-parabola example applies a spatially varying GRF surface load and compares optimized geometries across robustness levels. Higher robustness produces wavy diagonal contours, while regularization and coarse geometries limit fine details.

Distributed Surface Loading by GRF

  • The applied surface loading is a Matérn-covariance GRF with correlation length ρ = 0.25, f0 = 12500, and σ = 0.4.Figure 11 shows three samples of the applied load amplitudes.Source: Distributed Surface Loading by GRF · p. 23, Figure 11 · p. 23
  • At higher robustness levels α, diagonal geometries develop wavy contours attributed to reinforcements.The reported fine details seen in prior distributed-load studies do not develop with the current regularized, relatively coarse geometries.Source: Distributed Surface Loading by GRF · p. 23
  • Very small localized connections can pose manufacturing problems, so geometries without these details may be preferable depending on the application.Source: Distributed Surface Loading by GRF · p. 24
  • The α = 0 cost function is significantly lower than at higher robustness levels because of the variance penalty.Figure 13 tracks compliance, volume fraction, standard deviation, and cost function through the iterations.Source: Distributed Surface Loading by GRF · p. 24

Double Parabola with Lognormal Young’s Modulus

This double-parabola example studies how spatial variation in Young’s modulus affects optimized geometry, modeling the modulus as a lognormal random field so it remains positive. The resulting topology is qualitatively similar to the earlier result.

Double Parabola with Lognormal Young’s Modulus

Double Parabola with Geometric Perturbation by GRF

The geometric-uncertainty example perturbs an optimized level set with a Gaussian random field and fine-tunes the design for robust compliance. The paper also demonstrates optimization on a multi-patch cylinder under torsion.

Double Parabola with Geometric Perturbation by GRF

  • The geometric uncertainty model perturbs the ideal level set with a GRF, with noise level γ tunable to represent fabrication uncertainty.The examples fine-tune deterministic optimization results for an uncertainty level adapted to the fabrication process and problem geometry.Source: Double Parabola with Geometric… · pp. 25–26
  • The double-parabola fine-tuning uses diagonal loading F = 4.5 · 10^4 and target volume fraction φ* = 0.45.Figure 16 visualizes the unperturbed level set, perturbation field, and their resulting sum.Source: Double Parabola with Geometric… · p. 27, Figure 16 · p. 27
  • Perturbations visibly change the fine arms, which are more sensitive because they are thinner than the other reinforcements and are consequently removed.The comparison is between the deterministic result and the uncertain geometries; Figure 17 distinguishes them with red and blue contours/surfaces.Source: Double Parabola with Geometric… · p. 27, Figure 17 · p. 27

Example 3. Multi-Patch Cylinder

  • The multi-patch cylinder example applies a torsion-type line load of magnitude 255 and enforces target volume fraction ϕ* = 0.5.The initialization ensures the loaded edge is fully in the material phase; Figure 18 shows the mesh and initial shape.Source: Example 3. Multi-Patch Cylinder · p. 28, Figure 18 · p. 28
  • By iteration 100, the cylinder’s geometry is visually similar to its form after 30 iterations, when the cost function had its largest decrease.Figure 19 shows the iteration-100 design, while Figure 20 plots cost-functional, compliance, and volume-fraction histories.Source: Example 3. Multi-Patch Cylinder · p. 28, Figure 19 · p. 28, Figure 20 · p. 29

Torsion with axial loading by GRF

The cylinder torsion example extends the loading model with a stochastic axial component represented by a Matérn GRF. The figure compares optimized shapes for deterministic loading and robustness levels α = 0 and 10.

Torsion with axial loading by GRF

  • A Matérn GRF is superimposed as a random axial loading component to represent possible axial loads alongside the prior pure axial loading examples.Source: Torsion with axial loading by GRF · p. 29
  • The loading field shown uses correlation length ρ = 0.2 and smoothness γ = [1, 3], corresponding to harmonic and biharmonic SPDE-generated GRFs.Source: Torsion with axial loading by GRF · p. 29
  • Figure 21 compares optimized shapes for the deterministic case and robustness levels α = 0 and 10 at fixed correlation length ρ and smoothness γ.Source: Figure 21 · p. 29

4. Conclusion

The framework combines CAD-compatible, volume-based shape optimization with a biharmonic descent direction for Kirchhoff–Love multi-patch shells. It reproduces a literature reinforcement pattern with fewer elements and examines robust designs under loading, material, and geometric uncertainty.

4. Conclusion

  • The volume-based speed-method derivative avoids boundary approximations, while biharmonic regularization supplies the descent direction on the same multi-patch discretization.This avoids trimming, local refinement, and selective integration; the numerical implementation is open-source.Source: §4 · pp. 29–31
  • 40 × 40 cubic NURBS elements per patch recover the reinforcement pattern in obtained with 128 × 128 × 1 solid elements.The framework is validated on single- and multi-patch geometries under point, line, and surface loading, and under several sources of uncertainty.Source: §4 · p. 30
  • Increasing α from 0 to 5 reduces compliance standard deviation by roughly 73% while increasing its expectation by 6%.For uncertain loading direction in the double parabola, the robust design produces thicker diagonal reinforcements aligned with the less probable load orientations.Source: §4 · p. 30
  • Under material uncertainty, moderate dispersion and long correlation length leave the topology largely unchanged, with boundaries adapting only locally as α increases.The sampled fields vary slowly relative to the structural members, producing nearly uniform effective stiffness across reinforcements.Source: §4 · p. 30
  • Topology changes appreciably when uncertainty alters load character, whereas member-scale geometric perturbations remove the thinnest diagonal reinforcements.The geometric perturbations affect features most sensitive to boundary changes.Source: §4 · pp. 30–31
  • Coarse discretization and biharmonic regularization limit finer geometric detail; resolving it would require finer meshes, weaker regularization, or trimming-based integration.The authors note that small localized connections can be difficult to fabricate.Source: §4 · p. 31
  • Monte Carlo sampling makes each iteration's cost scale with sample size, while Python-based Bézier extraction dominates per-iteration overhead in the current implementation.The samples are independent and amenable to parallelization; variance reduction or multilevel sampling could reduce cost.Source: §4 · p. 31

Appendix A Derivation of result (20)

The appendices derive shape and strain sensitivities using transformations to a fixed parametric domain, then describe Gaussian random-field generation through SPDEs. They note boundary-related covariance error for bounded-domain approximations and restrict the treatment to integer smoothness parameters.

Appendix A Derivation of result (20)

  • The shape-derivative derivation maps the expression to parametric coordinates, applies product and chain rules, and pushes the result to the physical domain.The final step uses surface divergence.Source: Appendix A · p. 32

Appendix B Derivation of results (24), (25) and (27)

  • The strain-derivative derivations rewrite membrane and bending strains over the original parametric domain using transformed parametric gradients and Hessians.The appendix differentiates the membrane strain over pseudo-time and then develops the bending-strain expression.Source: Appendix B · pp. 32–34

Appendix C Generation of GRFs by SPDE Sampling

  • On bounded domains, boundary conditions are required, and the resulting approximation can introduce errors relative to Matérn covariance near the boundary.The stated Matérn covariance result applies when Equation 65 is solved on an unbounded domain.Source: Appendix C · p. 36
  • The treatment limits random-field smoothness to integer k, while pointing to a method for computing fields with general smoothness.The text motivates the restriction by noting that general smoothness is not required for most engineering applications.Source: Appendix C · p. 36
  • The harmonic and biharmonic SPDEs generate fields with smoothness parameters γ = 1 and γ = 3, respectively, and are compared on a non-matching two-patch domain.The accompanying figure shows sample fields and averaged covariance against the ground-truth Matérn covariance for a sample size of 4 · 104.Source: Appendix C · p. 37, Figure 22 · p. 37
Loading 2610.30389v1…