GENERATIVE DESIGN SHAPE OPTIMIZATION WITH SINGULARITIES AND DISCONNECTION PREVENTION FOR COMPUTER AIDED DESIGN AND MANUFACTURING

Abstract

Methods, systems, and apparatus, including medium-encoded computer program products, for computer aided design of physical structures using generative design processes. A method includes obtaining one or more load cases and one or more design criteria for a modeled object; iteratively modifying a three dimensional shape of the modeled object in accordance with the one or more design criteria and the one or more load cases, the iteratively modifying comprising regulating shape change velocities for an implicit surface representation of the three dimensional shape that exceed a reference velocity, where the reference velocity is set based on a mean and a standard deviation of a shape derivative on the implicit surface; and providing the three dimensional shape of the modeled object for use in manufacturing a physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

Description

BACKGROUND

This specification relates to computer aided design of physical structures, which can be manufactured using additive manufacturing, subtractive manufacturing and/or other manufacturing systems and techniques.

Computer Aided Design (CAD) software has been developed and used to generate three-dimensional (3D) representations of objects, and Computer Aided Manufacturing (CAM) software has been developed and used to evaluate, plan and control the manufacture of the physical structures of those objects, e.g., using Computer Numerical Control (CNC) manufacturing techniques. Typically, CAD software stores the 3D representations of the geometry of the objects being modeled using a boundary representation (B-Rep) format. A B-Rep model is a set of connected surface elements specifying boundaries between a solid portion and a non-solid portion of the modelled 3D object. In a B-Rep model (often referred to as a B-Rep), geometry is stored in the computer using smooth and precise mathematical surfaces, in contrast to the discrete and approximate surfaces of a mesh model, which can be difficult to work with in a CAD program.

CAD programs have been used in conjunction with subtractive manufacturing systems and techniques. Subtractive manufacturing refers to any manufacturing process where 3D objects are created from stock material (generally a “blank” or “workpiece” that is larger than the 3D object) by cutting away portions of the stock material. Such manufacturing processes typically involve the use of multiple CNC machine cutting tools in a series of operations, starting with a roughing operation, an optional semi-finishing operation, and a finishing operation. In addition to CNC machining, other subtractive manufacturing techniques include electrode discharge machining, chemical machining, waterjet machining, etc. In contrast, additive manufacturing, also known as solid free form fabrication or 3D printing, refers to any manufacturing process where 3D objects are built up from raw material (generally powders, liquids, suspensions, or molten solids) in a series of layers or cross-sections. Examples of additive manufacturing include Fused Filament Fabrication (FFF) and Selective Laser Sintering (SLS). Other manufacturing techniques for building 3D objects from raw materials include casting and forging (both hot and cold).

In addition, CAD software has been designed so as to perform automatic generation of 3D geometry (generative design) for a part or one or more parts in a larger system of parts to be manufactured. This automated generation of 3D geometry is often limited to a design space specified by a user of the CAD software, and the 3D geometry generation is typically governed by design objectives and constraints, which can be defined by the user of the CAD software or by another party and imported into the CAD software. The design objectives (such as minimizing the waste material or weight of the designed part) can be used to drive the geometry generation process toward better designs. The design constraints can include both structural integrity constraints for individual parts (i.e., a requirement that a part should not fail under the expected structural loading during use of the part) and physical constraints imposed by a larger system (i.e., a requirement that a part not interfere with another part in a system during use). Further, examples of design constraints include maximum mass, maximum deflection under load, maximum stress, etc.

Some CAD software has included tools that facilitate 3D geometry enhancements using lattices and skins of various sizes, thicknesses and densities, where lattices are composed of beams or struts that are connected to each other or directly to solid parts at junctions, and skins are shell structures that overlay or encapsulate the lattices. Such tools allow redesign of a 3D part to be lighter in weight, while still maintaining desired performance characteristics (e.g., stiffness and flexibility). Such software tools have used lattice topologies of various types that can be used to generate lattice structures that can be manufactured.

Moreover, the inputs to a generative design process can include a set of input solids (B-Rep input) that specify boundary conditions for the generative design process, but many modern generative design solvers do not operate directly on the exact surface boundary representation of their input solids. Instead, B-Reps are sampled and replaced with volumetric representations such as level sets or tetrahedral or hexahedral meshes, which are significantly more convenient and efficient for the physical simulations and material synthesis computed by the solver. The set of input solids can include “preserve bodies”, which should always be present in the design and which represent interfaces to other parts of the systems or locations on which boundary conditions should be applied (for example mechanical loads and constraints). Other regions in which geometry should or should not be generated can also be provided in a similar manner, such as input solids that define “obstacle bodies”, which represent regions where new geometry should not be generated.

SUMMARY

This specification describes technologies relating to computer aided design of physical structures using generative design processes, where the three dimensional (3D) models of the physical structures can be produced with singularities and disconnection prevention, and where the physical structures thus produced can then be manufactured using additive manufacturing, subtractive manufacturing and/or other manufacturing systems and techniques.

In general, one or more aspects of the subject matter described in this specification can be embodied in one or more methods, including: obtaining, by a computer aided design program, a design space for a modeled object, for which a corresponding physical structure will be manufactured, one or more design criteria for the modeled object, and one or more in-use load cases for the physical structure; iteratively modifying, by the computer aided design program, a generatively designed three dimensional shape of the modeled object in the design space in accordance with the one or more design criteria and the one or more in-use load cases for the physical structure, wherein the iteratively modifying includes performing numerical simulation of the modeled object in accordance with a current version of the three dimensional shape and the one or more in-use load cases to produce a current numerical assessment of a physical response of the modeled object, computing shape change velocities for an implicit surface in a level-set representation of the three dimensional shape, changing the shape change velocities in accordance with a polynomial function that has been fit to at least a portion of the shape change velocities above a reference velocity, updating the level-set representation using the shape change velocities to produce an updated version of the three dimensional shape of the modeled object, and repeating at least the performing, the computing, the changing and the updating until a predefined number of shape modification iterations have been performed or until the generatively designed three dimensional shape of the modeled object in the design space converges to a stable solution for the one or more design criteria and the one or more in-use load cases; and providing, by the computer aided design program, the generatively designed three dimensional shape of the modeled object for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

The one or more methods can include setting the reference velocity based on a mean and a standard deviation of a shape derivative on the implicit surface.

The one or more methods can include fitting the polynomial function being a cubic polynomial to velocity values above the reference velocity.

The iteratively modifying can include, after the updating and before the repeating: checking for an excessive change having been made during the updating; making the current version of the three dimensional shape be the updated version of the three dimensional shape for a next iteration to undue the excessive change; and slowing shape changes for the next iteration of the iteratively modifying.

Checking for the excessive change can include comparing a change in the one or more design criteria resulting from the updating with a predefined limit on an amount of change allowed for the one or more design criteria in a single iteration of the iteratively modifying.

Slowing the shape changes for the next iteration can include reducing a target volume change for the generatively designed three dimensional shape of the modeled object for the next iteration of the iteratively modifying.

The iteratively modifying can include, before performing the numerical simulation: identifying any elements, generated from the current version of the three dimensional shape for the numerical simulation, that are partially but not entirely within the implicit surface; and removing any of the identified elements that are not connected to any nodes of the current version of the three dimensional shape with Dirichlet boundary conditions, before performing the numerical simulation.

The iteratively modifying can include, before performing the numerical simulation: identifying any elements, generated from the current version of the three dimensional shape for the numerical simulation, that are partially but not entirely within the implicit surface; and setting densities of the identified elements equal to respective volume fractions of the identified elements, wherein a volume fraction of an identified element is a fractional amount of the identified element that falls within the implicit surface; wherein performing the numerical simulation includes penalizing stiffness of the identified elements according to the densities.

These and other methods described herein can be implemented using a non-transitory computer-readable medium encoding a computer aided design program operable to cause one or more data processing apparatus to perform the method(s). In some implementations, a system includes: a non-transitory storage medium having instructions of a computer aided design program stored thereon; and one or more data processing apparatus configured to run the instructions of the computer aided design program to perform the method(s). Further, such systems can include an additive manufacturing machine, or other manufacturing machines, and the one or more data processing apparatus can be configured to run the instructions of the computer aided design program to generate instructions for such machines (e.g., toolpath specifications for the additive manufacturing machine) from the three dimensional model, and manufacture the physical structure corresponding to the object with the machines using the instructions (e.g., the additive manufacturing machine using the toolpath specifications).

Particular embodiments of the subject matter described in this specification can be implemented to realize one or more of the following advantages.

Singularities, e.g., caused by sharp re-entrant corners or bad meshing, and geometric disconnections, e.g., caused by the use of Ersatz materials, can be mitigated or outright avoided. The techniques described can be implemented for automatic detection and correction, allowing for shape and topology optimization to continue where otherwise it may halt as a result of detecting a disconnection or singularity. Objects can be optimized even in limited design spaces. Identified regions of improvement, i.e., regions containing singularities or disconnections, can be corrected for to prevent affecting the stress constraint defined for the object during design optimization.

The process can be performed automatically and in conjunction with other techniques, such as the level-set method, conventional body-fitted solvers, backup-and-recovery, load case-specific advection to prevent disconnected geometry, Proportional-Integral-Derivative (PID) controllers, other adaptive controllers and related technologies, including adaptive PID tuning, PID auto-tuning, the latter of which can be used to meet arbitrary design constraints.

The details of one or more embodiments of the subject matter described in this specification are set forth in the accompanying drawings and the description below. Other features, aspects, and advantages of the invention will become apparent from the description, the drawings, and the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A shows an example of a system usable to design and manufacture physical structures.

FIG. 1B shows an example of a process of designing and manufacturing physical structures.

FIG. 2A shows a graphical representation of an example of geometry mapping between an initial design configuration and a current design configuration.

FIG. 2B shows a graphical representation of an example of a data mapping from solid mesh to level-set grid.

FIG. 2C shows an example of a problem set up for a level-set based topology optimization process; FIG. 2D shows a narrow-band velocity example for the problem set up from FIG. 2C; and FIG. 2E shows the velocity from FIG. 2C, after extension.

FIG. 3A shows an example of a process that generates one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes.

FIG. 3B shows a graphical representation of examples of initial seeding geometries applied to a design space and resulting designs for a Michell type arch problem after a generative design process is executed on the design space.

FIG. 3C shows a graphical representation of an example of a bubble insertion history for intermediate designs for a Michell type arch problem during execution of a generative design process.

FIG. 4A shows an example of a process that generates one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes implementing arbitrary constraint handling with controlled convergence.

FIG. 4B shows examples of graphs tracking progressive reduction and increase of a target constraint value using B-splines, respectively.

FIG. 4C shows a graphical representation of example variations of a value μ with a target volume change given different voxel sizes.

FIG. 4D shows examples of graphs tracking target volume versus actual volume fraction during an iterative optimization process with approximate volume control and without adaptive control.

FIG. 4E shows examples of graphs tracking target volume versus actual volume fraction during an iterative optimization process with PID control (but without adaptive PID control).

FIG. 4F shows an example of a process of adaptively modifying parameter values for a PID controller.

FIG. 4G shows an example of a graph that tracks different measures used in constraint normalization.

FIG. 4H shows examples of graphs tracking convergence history of constraints with and without line search.

FIG. 5A shows an example of a process that generates one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes that solve for safe life fatigue constraints on the body.

FIG. 5B shows an example of an SN-curve tracking fatigue strength (stress) on a body over a number of cycles.

FIG. 6A shows an example of a process of iteratively modifying a generatively designed three dimensional shape of a modeled object in the design space in accordance with the one or more design criteria, including at least one damage tolerant fatigue constraint.

FIG. 6B shows graphical representations of geometries in which thickness of the geometries are computed according to different measuring techniques.

FIG. 6C shows an example of a graph tracking curves of sigmoid functions at different rates.

FIG. 6D shows an example of a table describing fatigue crack properties.

FIG. 6E shows an example of a process of generating one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes that solve for damage tolerant fatigue constraints on the body.

FIG. 7A is an example of a process of iteratively modifying a generatively designed three dimensional shape of a modeled object in the design space in accordance with the one or more design criteria, including a stress constraint.

FIG. 7B shows a graphical representation of an example of a build angle measured on an object.

FIG. 7C shows an example of a process of generating one or more portions of a three dimensional model of an object to be manufactured, using a stress constraint.

FIG. 8A shows an example of a process of iteratively modifying a generatively designed three dimensional shape of a modeled object in the design space in accordance with the one or more design criteria, while avoiding excessive sudden changes and minimizing the likelihood of disconnections.

FIG. 8B is a graphical representation of an example of geometry disconnection during optimization.

FIG. 8C is a graphical representation of an example of a geometry with simulation elements classified based on intersection of the elements with the geometry.

FIG. 9 shows a schematic diagram of a data processing system including a data processing apparatus, which can be programmed as a client or as a server.

Like reference numbers and designations in the various drawings indicate like elements.

DETAILED DESCRIPTION

FIG. 1A shows an example of a system 100 usable to design and manufacture physical structures. A computer 110 includes a processor 112 and a memory 114, and the computer 110 can be connected to a network 140, which can be a private network, a public network, a virtual private network, etc. The processor 112 can be one or more hardware processors, which can each include multiple processor cores. The memory 114 can include both volatile and non-volatile memory, such as Random Access Memory (RAM) and Flash RAM. The computer 110 can include various types of computer storage media and devices, which can include the memory 114, to store instructions of programs that run on the processor 112, including Computer Aided Design (CAD) program(s) 116, which implement three-dimensional (3D) modeling functions and includes one or more generative design processes for topology optimization using numerical simulation, and including material or microstructure shape optimization techniques, geometrical or macrostructure shape optimization techniques, or both (e.g., using one or more level-set based topology optimization processes).

The numerical simulation performed by the CAD program(s) 116 can simulate one or more physical properties and can use one or more types of simulation to produce a numerical assessment of physical response (e.g., structural response) of the modelled object. For example, finite element analysis (FEA), including linear static FEA, finite difference method(s), and material point method(s) can be used. Further, the simulation of physical properties performed by the CAD program(s) 116 can include Computational Fluid Dynamics (CFD), Acoustics/Noise Control, thermal conduction, computational injection molding, electric or electro-magnetic flux, and/or material solidification (which is useful for phase changes in molding processes) simulations. Moreover, the CAD program(s) 116 can potentially implement hole and/or fixture generation techniques to support clamping during manufacturing and/or manufacturing control functions.

As used herein, CAD refers to any suitable program used to design physical structures that meet design requirements, regardless of whether or not the CAD program is capable of interfacing with and/or controlling manufacturing equipment. Thus, CAD program(s) 116 can include Computer Aided Engineering (CAE) program(s), Computer Aided Manufacturing (CAM) program(s), etc. The CAD program(s) 116 can run locally on computer 110, remotely on a computer of one or more remote computer systems 150 (e.g., one or more third party providers' one or more server systems accessible by the computer 110 via the network 140) or both locally and remotely. Thus, a CAD program 116 can be two or more programs that operate cooperatively on two or more separate computer processors in that one or more programs 116 operating locally at computer 110 can offload processing operations (e.g., generative design and/or numerical simulation operations) “to the cloud” by having one or more programs 116 on one or more computers 150 perform the offloaded processing operations.

The CAD program(s) 116 present a user interface (UI) 122 on a display device 120 of the computer 110, which can be operated using one or more input devices 118 of the computer 110 (e.g., keyboard and mouse). Note that while shown as separate devices in FIG. 1A, the display device 120 and/or input devices 118 can also be integrated with each other and/or with the computer 110, such as in a tablet computer (e.g., a touch screen can be an input/output device 118, 120). Moreover, the computer 110 can include or be part of a virtual reality (VR) and/or augmented reality (AR) system. For example, the input/output devices 118, 120 can include a VR/AR input glove 118a and/or a VR/AR headset 120a. In any case, a user 160 interacts with the CAD program(s) 116 to create and modify 3D model(s), which can be stored in 3D model document(s) 130.

In the example shown, an initial 3D model 132 is a seed model for input to a generative design process. In this example, the user 160 has defined a mechanical problem, for a generative design process to operate on, to produce a new 3D model from a starting 3D model 132. In this case, the defined problem is the Michell type arch problem, where the user 160 has specified a domain 134 and loading cases 136. However, this is but one of many possible examples.

In some implementations, the user 160 (or other person or program) can specify a design space for an object to be manufactured, a numerical simulation setup (e.g., load(s) and material(s)) for numerical simulation (e.g., FEA, CFD, Acoustics/Noise Control, thermal conduction, computational injection molding simulations, electric or electro-magnetic flux, material solidification, etc.) of the object, at least one design objective (e.g., minimize material usage) for the object, and at least one design constraint (e.g., a volume constraint) for the object. In some implementations, the inputs for use in numerical simulation and generative design processes can include one or more regions of a current 3D model in which to generate new 3D geometry, loading case(s) defining one or more loads in one or more different directions to be borne by a physical structure being designed, one or more materials (e.g., one or more isotropic solid materials identified as a baseline material model for the design space), one or more seed model types to use as input to a generative design process, one or more generative design processes to use, and/or one or more lattice topologies to use in one or more regions of the design space. Inputs to the generative design and numerical simulation processes can include non-design spaces, different types of components (e.g., rods, bearings, shells), one or more target manufacturing processes and associated parameters, obstacle geometries that should be avoided, preserve geometries that should be included in the final design, and parameters related to various aspects, such as resolution of the design, type of synthesis, etc.

Moreover, the CAD program(s) 116 provide user interface elements in the UI 122 to enable the user 160 to specify the various types of inputs noted above, and all (or various subsets) of these inputs can be used in the generative design and numerical simulation processes described in this document. Further, the user 160 can be enabled by the UI 122 of the CAD program(s) 116 to design a part using traditional 3D modelling functions (to build precise geometric descriptions of the 3D design model) and then use generative design and simulation processes in a design space specified within one or more portions of the 3D design model. Thus, as will be appreciated, many possible types of physical structures can be designed using the systems and techniques described in this document, the UI 122 can be used to create a full mechanical problem definition for a part to be manufactured, and the generative design and numerical simulation processes can accelerate new product development by enabling increased performance without time consuming physical testing.

Further, as described herein, the CAD program(s) 116 implement at least one generative design process, which enables the CAD program(s) 116 to generate one or more portions of the 3D model(s) automatically (or the entirety of a 3D model) based on design objective(s) and constraint(s), i.e., design criteria, where the geometric design is iteratively optimized based on simulation feedback. Note that, as used herein, “optimization” (or “optimum”) does not mean that the best of all possible designs is achieved in all cases, but rather, that a best (or near to best) design is selected from a finite set of possible designs that can be generated within an allotted time (e.g., as specified by a predefined number of shape modification iterations) given the available processing resources. The design criteria can be defined by the user 160, or by another party and imported into the CAD program(s) 116. The design criteria can include structural integrity constraints for individual parts (e.g., a requirement that a part should not fail under the expected structural loading during use of the part) and physical constraints imposed by a larger system (e.g., a requirement that a part be contained within a specified volume so as not to interfere with other part(s) in a system during use).

Various generative design processes can be used, which can optimize the shape and topology of at least a portion of the 3D model. The iterative optimization of the geometric design of the 3D model(s) by the CAD program(s) 116 involves topology optimization, which is a method of light-weighting where the optimum distribution of material is determined by minimizing an objective function subject to design constraints (e.g., structural compliance with volume as a constraint). Topology optimization can be addressed using a variety of numerical methods, which can be broadly classified into two groups: (1) material or microstructure techniques, and (2) geometrical or macrostructure techniques. Microstructure techniques are based on determining the optimum distribution of material density and include the Solid Isotropic Material with Penalization (SIMP) method and the homogenization method. In the SIMP method, intermediate material densities are penalized to favor either having p=0 or p=1, denoting a void or a solid, respectively. Intermediate material densities are treated as composites in the homogenization method.

In contrast, macrostructure techniques treat the material as being homogeneous, and the three dimensional topology of the modeled object being produced is represented as one or more boundaries between one or more solid regions (having the homogenous material therein) and one or more void regions (having no material therein) within the design space (also referred to as the domain or a sub-space of the domain for topology optimization). The shape(s) of the one or more boundaries are optimized during the generative design process, while the topology is changed in the domain as a result of the shape optimization in combination with adding/removing and shrinking/growing/merging the void region(s). Thus, the types of final optimized topologies that can result from a generative design process using a macrostructure technique can depend significantly on the number and sizes of voids within the seed geometry along with the addition and removal of voids during the optimization process.

Note that, while only one seed model 132 is shown in FIG. 1A (where this model 132 includes a complex solid region 132A surrounding many holes 132B of the void region) it should be appreciated that the generative design processes described in this document can employ two or more seed geometries/models for any given generative design process iteration, so as to improve the final result of shape and topology optimization. Further, during the shape and topology optimization process, one or more voids can be introduced into the solid domain and/or one or more solids can be introduced into the void domain, so as to improve the final result of the shape and topology optimization. Thus, the CAD program(s) 116 can include various types of available seed geometries and mid-process geometry introductions, along with a user interface element allowing the user 160 to design their own seed geometries and mid-process geometry introductions. Likewise, the user 160 can run two or more generative design process iterations (saving the results from each) until a preferred generative design is produced.

In various implementations, the CAD program(s) 116 provide generative design shape optimization processes (1) controlled convergence, (2) with damage prevention over loading cycles, (3) with size limited fatigue damage, (4) using a build material strength model, and/or (5) with singularities and disconnection prevention, as described herein. In some implementations the CAD program(s) 116 implement all of the above-listed processes, while in other implementations the CAD program(s) 116 implement a subset of the above-listed processes.

Once the user 160 is satisfied with a generatively designed 3D model, the 3D model can be stored as a 3D model document 130 and/or used to generate another representation of the model (e.g., an .STL file for additive manufacturing). This can be done upon request by the user 160, or in light of the user's request for another action, such as sending the 3D model 132 to an additive manufacturing (AM) machine 170, or other manufacturing machinery, which can be directly connected to the computer 110, or connected via a network 140, as shown. This can involve a post-process carried out on the local computer 110 or a cloud service to export the 3D model 132 to an electronic document from which to manufacture. Note that an electronic document (which for brevity will simply be referred to as a document) can be a file, but does not necessarily correspond to a file. A document may be stored in a portion of a file that holds other documents, in a single file dedicated to the document in question, or in multiple coordinated files.

In any case, the CAD program(s) 116 can provide a document 135 (having toolpath specifications of an appropriate format) to the AM machine 170 to produce a complete structure 138, which includes the optimized topology and shape (in this example, an arch design generated for the Michell type arch problem). The AM machine 170 can employ one or more additive manufacturing techniques, such as granular techniques (e.g., Powder Bed Fusion (PBF), Selective Laser Sintering (SLS) and Direct Metal Laser Sintering (DMLS)), extrusion techniques (e.g., Fused Deposition Modelling (FDM), which can include metals deposition AM). In addition, the user 160 can save or transmit the 3D model for later use. For example, the CAD program(s) 116 can store the document 130 that includes the generated 3D model.

In some implementations, subtractive manufacturing (SM) machine(s) 174 (e.g., a Computer Numerical Control (CNC) milling machine, such as a multi-axis, multi-tool milling machine) can also be used in the manufacturing process. Such SM machine(s) 174 can be used to prepare initial work-pieces on which AM machine(s) 170 will operate. In some implementations, a partially complete structure 138 is generated by the AM machine(s) 170 and/or using casting methods (e.g., investment casting (IC) using ceramic shell or sand casting (SC) using sand cores), and this partially complete structure 138 then has one or more portions removed (e.g., finishing) by the CNC machine 174 in order to form the completed structure. In some implementations, the CAD program(s) 116 can provide a corresponding document 135 (having toolpath specifications of an appropriate format, e.g., a CNC numerical control (NC) program) to the SM machine 174 for use in manufacturing the part using various cutting tools, etc. Moreover, in some implementations, the complete structure 138 is produced in its entirely using SM machine(s) 174.

In various implementations, the CAD program(s) 116 of the system 100 can implement one or more generative design processes as described in this document.

Generative design processes seek an optimal geometric shape, topology, or both. For example, generative design processes seek an optimal geometric shape among alternative designs by minimizing a performance-related objective function subject to constraints:

	minimize	J (s,u(s)) s ϵ	(1)
	such that	gi(s, u (s)) = 0 i = 1, . . . , ng	(2)

where s is a vector of design variables related to a geometric shape of the domain, and u is a vector of state variables (e.g., displacement) that depend on s. Additional constraints (e.g., equilibrium) are denoted by a set g. For simplicity, equality constraints are assumed here. Mathematical programming methods used to minimize (1) can be gradient-based or non-gradient-based. Gradient-based methods (versus non-gradient-based methods) generally use more information associated with design sensitivity, for example:

$\begin{array}{cc}\frac{dJ}{ds}\left(s,u\left(s\right)\right)=\frac{\partial J}{\partial s}+\frac{\partial J}{\partial u}\frac{du}{ds}& \left(3\right)\end{array}$

which is a derivative of the performance-related objective function with respect to the design variables. In lattice-based methods, s represents a lattice thickness. In level-set based topology optimization methods, s represents a boundary of a solid region.

FIG. 1B shows an example of a process of designing and manufacturing physical structures. Design variables are obtained 180, e.g., by CAD program(s) 116, for use in producing a generative 3D model. Different generative design processes can be formulated by using different combinations of design variables, which can include lattices, density fields, and level-sets. In some implementations, the design variables can include various types of inputs, e.g., received through UI 122, such as selection among different generative design synthesis methods made available by CAD program(s) 116 in the system 100. In some implementations, the available generative design synthesis methods include (1) a level-set-based topology optimization that provides a basic level-set method for topology optimization, (2) a lattice and skin optimization that provides a thickness optimization of lattice and skin, (3) a hybrid topology optimization that provides a topology optimization with lattice infill, (4) an inside-out hybrid topology optimization in which the lattice infill is present in a negative space between the topology-optimized design and the original design space, (5) a hollow topology optimization that provides a method for topology optimization with internal hollow regions, and/or (6) a hybrid-hollow topology optimization that provides a method for topology optimization with lattice infill and internal hollow regions. For further details regarding such generative design synthesis methods, see PCT/US2019/060089, filed Nov. 6, 2019, published as PCT Publication No. WO 2020/097216 A1 on May 14, 2020, and titled, “MACROSTRUCTURE TOPOLOGY GENERATION WITH DISPARATE PHYSICAL SIMULATION FOR COMPUTER AIDED DESIGN AND MANUFACTURING,” which is hereby incorporated by reference. Moreover, the available generative design synthesis methods can employ one or more of the systems and techniques described below in connection with FIGS. 3A-8C.

Additional design variables are possible, such as (1) a design space for generative design geometry production, e.g., a boundary representation (B-Rep) 3D model designed or loaded into CAD program(s) 116 that serves as a sub-space of an optimization domain of a described generative design process, and/or (2) a set of input solids that specify boundary conditions for generative design geometry production, e.g., B-Reps selected using UI 122 to specify sub-space(s) that are preserved for use as connection point(s) with other component(s) in a larger 3D model or separate 3D model(s). Different combinations of design variables can be used, e.g., by CAD program(s) 116 in response to input from the user 160. For example, a user 160 may select different generative design synthesis methods to use within respective different design spaces within a single 3D model.

Other design variables can include a setup for numerical simulation, e.g., densities of elements in an FEA model or a homogenized lattice material representation for a selected lattice topology to be used with a topology optimized 3D shape of the part being generatively designed. The design variables can include various design objectives and constraints, such as described in this document. Furthermore, functions can be provided, e.g., by CAD program(s) 116, that assist the user in specifying design variables. For example, a lattice recommender can provide predictions for suitable lattice settings for a given problem using a single solid simulation. In some implementations, the lattice recommender described in PCT Publication No. WO 2017/186786 A1, filed 26 Apr. 2017, titled “METHOD AND SYSTEM FOR GENERATING LATTICE RECOMMENDATIONS IN COMPUTER AIDED DESIGN APPLICATIONS,” which is hereby incorporated by reference, is used.

With the generative design variables specified, one or more 3D model(s) are produced 185, e.g., by CAD program(s) 116, using one or more generative design processes (e.g., using one or more selected generative design synthesis methods). In some implementations, the one or more generative design processes can use described level-set methods, where s, from Equations 1, 2 & 3, represents a boundary of a solid region that is implicitly represented using one or more level-sets, which are signed distance values computed on a Cartesian background grid. In a level-set-based topology optimization method, the outer shape of a structure is represented by a one-dimensional high-level level set function, and a change in shape and configuration is replaced by a change in the level set function value, so as to obtain an optimum structure. The level set function refers to a function that indicates whether each part of the design domain where the initial structure is set corresponds to a material domain (material phase) that forms the structure and is occupied by a material, a void domain (void phase) where a void is formed, or a boundary between these two domains, wherein a predetermined value between a value representing the material domain and a value representing the void domain represents the boundary between the material domain and the void domain.

In some implementations that use level-set based topology optimization methods, one or more octree data structures are used for resolving geometry accurately. Level-set based topology optimization generally involves optimizing the shape of a design domain using a shape derivative, which is the derivative of a constrained minimization problem with respect to the shape. The shape changes are applied on a level-set, which allows topology changes during shape modifications. The outcome of this type of generative design process is the partitioning of the design space into solid and void regions, resulting in an optimized shape, often with topology changes. For this type of level-set-based topology optimization, as well as the variations on this type of level-set-based topology optimization described in this document, one or more of the following approaches can be used.

Linear Elastic Topology Optimization

Consider the linear elastic boundary value problem for a solid body with the domain Ω:

−∇·D∈(u)=ƒ in Ω  (4)

u=0 on Γ_D  (5)

D∈(u)n=t on Γ_N  (6)

where ∈(u) is the linear strain tensor, D is the fourth order constitutive tensor, u is the displacement vector, ƒ is the external load vector and t is the prescribed traction on the Neumann boundary Γ_N with the outward normal n. For simplicity, homogeneous Dirichlet boundary conditions can be assumed on Γ_D. The constrained topology optimization problem can then be

	minimize	J (Ω, u)		(7)
	subject to	−∇ · D∈(u) = f	in Ω	(8)
		u = 0	on ΓN	(9)
		D∈(u)n = t	on ΓN	(10)

where compliance minimization can be used as the objective function

$\begin{array}{cc}J\left(\Omega ,u\right)={\int }_{\Omega }f·\mathrm{ud}\Omega +{\int }_{{\Gamma }_N}\overline{t}·\mathrm{ud}\Gamma =\frac{1}{2}{\int }_{\Omega }D\in \left(u\right):\in \left(u\right)d\Omega & \left(11\right)\end{array}$

FIG. 2A shows a graphical representation of an example of geometry mapping 240 between an initial design configuration and a current design configuration. The solution space in topology optimization can be defined by different perturbations of the geometric shape within the design space. In this context, a linear mapping 242 which maps a given domain 246 Ω, into a perturbed domain 248 Ωt can be defined. With this mapping, a material point with the coordinate x∈Ω can be mapped 244 onto

x _t =x+tδv,t≥0  (12)

where δv is a prescribed constant vector field, and t is a scalar parameter (see FIG. 2A). Note that solving equations using gradient based mathematical programming methods involves using the directional derivative of the objective function in the direction of the velocity field δv

$\begin{array}{cc}\frac{\mathrm{dJ}}{d\Omega }\left(\Omega ,u\left(x\right)\right)\delta v=\underset{t\to 0}{\lim }\frac{J\left({\Omega }_{t,}u\left(x_t\right)\right)-J\left(\Omega ,,u\left(x\right)\right)}{t}& \left(13\right)\end{array}$

More than one approach can be used to obtain the directional derivative of the objective function for use in gradient based optimization methods. Approaches that are suitable for use with gradient based optimization methods include direct differentiation, semi-analytical derivatives, adjoint method, and finite difference. Moreover, further approaches to obtaining values for a directional derivative of an objective function, including approximation techniques, are described in detail below in connection with FIGS. 4A-4F.

Adjoint Method

Evaluating the shape derivative (Equation 13) can require the directional derivative of the state variable u in the direction of the velocity vector δv. This is can be seen by using the chain rule

$\begin{array}{cc}\frac{\mathrm{dJ}}{d\Omega }\delta v=\frac{\partial J}{\partial \Omega }\delta v+\frac{\partial J}{\partial u}\frac{\partial u}{\partial \Omega }\delta v& \left(14\right)\end{array}$

But in some implementations, an adjoint method can be used which involves the formation of a Lagrangian L(δ, u, λ) which depends on domain shape Ω, displacement field u, and Lagrange parameters λ

L(Ω,u,Δ)=J(Ω,u)+λ[∫_Ωƒ+∇·D∈(u)dΩ−∫_ΓDudΓ+∫_ΓNt−D∈(u)ndΓ]  (15)

The stationary condition for the Lagrangian, i.e., δL(Ω, u, λ)=0, can yield a complete set of shape optimization equations. For example, the adjoint problem for compliance minimization (Equation 11) can be given by considering the variation of the Lagrangian with respect to the displacements u. After introducing the cost function (Equation 11) and reformulating the domain term with the divergence theorem

∫_Ωƒ·δudΩ+∫_ΓNt·δudΓ−∫_Ωδu·(∇σ(λ))dΩ−∫_ΓDλ·(D:∇(δu))ndΓ+∫_ΓNδu·σ(λ)ndΓ=0  (16)

the corresponding boundary value problem, referred to as the adjoint problem, can become

−∇·σ(λ)=−ƒ in Ω  (17)

λ=0 in Γ_D  (18)

σ(λ)·n=t in Γ_N  (19)

This can lead to determining that λ=−u is the solution of the adjoint problem. This means that the adjoint problem (Equations 17-19) does not need to be solved explicitly for the compliance minimization problem (Equation 11). Such problems are called self-adjoint, where the solution of the direct problem also yields the adjoint solution. However, this is not often the case and different adjoint problems may have to be solved depending on the nature of the direct problem and the objective function. An advantage of the use of the Lagrangian as is includes the identity:

$\begin{array}{cc}\frac{\mathrm{dJ}}{d\Omega }\left(\Omega ,u\left(x\right)\right)\delta v=\frac{\partial L}{\partial \Omega }\left(\Omega ,u\left(x\right),\lambda \right)\delta v& \left(20\right)\end{array}$

This equation can enable the shape derivative (Equation 13) to be expressed as a boundary integral of the following form:

$\begin{array}{cc}\begin{array}{c}\mathrm{DJ}\left(x,u\left(x\right)\right)\left[\delta v\right]={\int }_{\Omega }f\left(u,\lambda \right)\left(\nabla ·\delta v\right)d\Omega \\ ={\int }_{\Gamma }f\left(u,\lambda \right)\left(\delta v·n\right)d\Gamma \end{array}& \left(21\right)\end{array}$

Without loss of generality, it can be assumed that some boundary variations are not relevant in practical shape optimization. In solid mechanics, the boundary variations can usually be of the form:

δv = 0	on ΓD		(22)
δv = 0	on ΓN	with σn = t,
δv ≠ 0	on ΓN	with σn = 0.

This means that only parts of the boundary TN with no traction are free to move during the shape optimization. In this context, the variation of the Lagrangian (Equation 15) in the direction OP with structural compliance (Equation 11) as the cost function can become:

$\begin{array}{cc}\frac{\partial L}{\partial \Omega }\delta v={\int }_{{\Gamma }_N}\left(2u·f-\nabla u:\sigma \left(u\right)\right)\left(\delta v·n\right)d\Gamma & \left(23\right)\end{array}$

Without restricting δv as stated in Equation 22, the variation of the Lagrangian can contain several more terms. During iterative optimization of the shape, the shape derivative (Equation 23) can be used as gradient information. In order to achieve maximum decrease in the objective function, the boundary perturbation can be chosen as follows

δv=−(2u·ƒ−D∈(u):∈(u)).  (24)

This boundary perturbation can be applied along the direction of the normal δ{circumflex over (v)}=vn, where v is the shape change velocity and is given by

$\begin{array}{cc}v=\frac{\mathrm{dJ}}{d\Omega }\left(\Omega ,u\left(x\right)\right)\delta v={\int }_{\Gamma }\left(2u·f-D\in \left(u\right):\in \left(u\right)\right)\left(\delta v·n\right)d\Gamma & \left(25\right)\end{array}$

Volume Control

Topology optimization using only a compliance minimization objective (Equation 11) can result in the optimum topology covering the full design space. Thus, some form of volume constraint is often required. Moreover, in some implementations, control over volume changes during topology optimization can be important for several reasons: 1) to enforce volume constraints; 2) to provide user control of topology optimization progress, e.g., more volume changes during initial iterations and fewer volume changes during later iterations; and 3) to ensure that arbitrary constraints not having shape derivatives are satisfied.

Note that the presence of shape derivatives for constraints can require modifying the shape change velocity in Equation 25. A modified objective function can be considered where the volume is penalized by a penalty parameter μ in:

$\begin{array}{cc}J\left(\Omega ,u\right)=\frac{1}{2}{\int }_{\Gamma }\left(D\in \left(u\right):\in \left(u\right)d\Omega +\mu V\left(\Omega \right)\right)d\Gamma & \left(26\right)\end{array}$

The corresponding shape derivative (Equation 25) can then be given by:

$\begin{array}{cc}\frac{\mathrm{dJ}}{d\Omega }\left(\Omega ,u\left(x\right)\right)\delta v={\int }_{\Gamma }\left(2u·f-D\in \left(u\right)+\mu \right)\left(\delta v·n\right)d\Gamma & \left(27\right)\end{array}$

wherein μ is constant along the boundary. The velocity term in the shape derivative (Equation 25) can now have an additional term as follows:

v=−(2u·ƒ−D∈(u):∈(u)+μ)  (28)

Augmented Lagrangian Method

In some implementations, an Augmented Lagrangian method is used. Some approaches (see Volume Control above) can have limitations such as creating difficulty in meeting prescribed volume targets. Essentially, the final volume of the design can depend on the value of μ prescribed in Equation 26. In such cases, meeting prescribed design constraint targets can be achieved by using the Augmented Lagrangian method. Consider the following Lagrangian for compliance minimization with a final volume target of V_ƒ(Ω):

$\begin{array}{cc}L\left(\Omega ,u\right)=\frac{1}{2}{\int }_{\Omega }D\in \left(u\right):\in \left(u\right)d\Omega +\mu \left(V\left(\Omega \right)-V_f\left(\Omega \right)\right)+\frac{\lambda }{2}{\left(V\left(\Omega \right)-V_f\left(\Omega \right)\right)}^2& \left(29\right)\end{array}$

The shape derivative can then be given by:

$\begin{array}{cc}\frac{dJ}{d\Omega }\left(\Omega ,u\left(x\right)\right)\delta v={\int }_{\Gamma }\left(2u·f-D\in \left(u\right):\in \left(u\right)+\mu +\frac{\lambda }{2}\left(V\left(\Omega \right)-V_f\left(\Omega \right)\right)\right)\left(\delta v·n\right)d\Gamma & \left(30\right)\end{array}$

where the penalty parameters λ, μ can be updated in an increasing sequence such that they converge to the optimal Lagrange multipliers. In some implementations, one or more heuristic methods are used for updating penalty parameters.

Body-Fitted Solvers

In some implementations, one or more body-fitted mesh based solvers are used. Using such body-fitted mesh based solvers with the level-set method involves mapping data from the solid mesh to a Cartesian grid (note that the inverse mapping is trivial due to the structured nature of the Cartesian grid). This involves two mappings as shown in FIG. 2B, which shows a graphical representation of an example of a data mapping 260 from solid mesh 262 to level-set grid 264.

Data in solid mesh elements (e.g., strain energy, Von Mises stress) can first be mapped to solid mesh nodes. This mapping can be achieved by data averaging. For example, by averaging solid mesh element e_i data in the solid mesh 262 at solid node n_j. Further, data in the solid mesh nodes can be mapped to voxel grid points using linear shape functions. Data at solid mesh nodes n_j can be linearly interpolated to level-set grid point g_i in the level-set grid 264. Having such a mapping can allow the level-set method to be used with complex FEA models solved with body-fitted solvers.

A detailed example of a macrostructure topology optimization process is now described. For simplicity of presentation, the compliance minimization problem is used with a penalized volume (Equation 26). Moreover, for all the detailed examples below, it is presumed that FEA is used for numerical simulation for ease of presentation, but the other numerical simulation types noted above can also be used.

The starting shape can simply be the design space or the intersection of design space with a suitable seed geometry. The initial level-set ψ₀ can be created by converting the starting shape to a signed distance field (SDF). The conversion can be done using an implementation of an OpenVDB toolkit. Other methods are also possible. Next, the shape can be iteratively optimized until the objective function has converged. The FEA model used for simulation includes solid elements throughout the solid region(s) of the generative design space being optimized. In each iteration, the constitutive model D of each element e in the FEA model can be updated by changing the modulus of elasticity according to the relative location of the element with respect to the current level-set. For example, elements outside the level-set (called void elements) are given very low stiffness D_v, while the constitutive model of those inside are set to that of the original solid material D_s. This can be achieved by checking the average level-set of the element nodes:

$\begin{array}{cc}\psi \left(e\right)=\frac{1}{n_e}{\Sigma }_j^{n_e}\psi \left(n_j\right)& \left(31\right)\end{array}$

where n_j denotes the coordinates of element nodes. Once the FEA model is up-to-date with the current geometry represented by the level-set, the boundary value problem (Equations 4-6) can be solved in order to compute an advection velocity. Shape changes can be applied by advecting the level-set ψ using, for example, the Hamilton-Jacobi equation:

$\begin{array}{cc}\frac{\partial \psi }{\partial t}-v❘\nabla \psi ❘=0& \left(32\right)\end{array}$

where v is the shape derivative (Equation 28). Note that one or more heuristic shape update methods can be used, and one or more methods to move the 0^th iso-contour of a level-set can be used given a movement direction specified by the shape derivative. A linear mapping between the FEA nodes and the level-set grid points (see Body-fitted Solvers above) can allow the FEA results, e.g., D∈(u):∈(u), to be transferred to the level-set where the advection velocity v is computed.

Once the objective function has converged, the surface of the final level-set can be extracted using a contouring method which extracts the 0^th iso-contour of ψ. An example of the algorithm is as follows:

Level-Set Algorithm for Level-set Based Topology Optimization
    Input: (Ω, Dv, Ds, μ)
    Output: (Ωs)
    // initiate level-set ψ from design space
1:  ψ0 = fSDF(Ω)
    // iterate until convergence tolerance c is met
2:  i = 0
3:  while i = 0 or |Ji − Ji−1| > c do
    // set constitutive model of FEA elements as void Dv or solid Ds
4:  K                      e                                        ⁢                                          {                                                                                                                                                                  B                                T                                                            ⁢                                                              D                                s                                                            ⁢                              B                                                                                                                                                                          if                                ⁢                                                                                                                                         ψ                                    i                                                                    (                                  e                                  )                                                                                            ≤                              0                                                                                                                                                                                                                          B                                T                                                            ⁢                                                              D                                v                                                            ⁢                              B                                                                                                                                                                          if                                ⁢                                                                                                                                         ψ                                    i                                                                    (                                  e                                  )                                                                                            >                              0
    // formulate and solve FEA problem
5.  K = Σen Ke
6.  u − K−1f
    // compute advection velocity
7:  v                    =                                                                                            1                          2                                                ⁢                        D                                            ∈                                                                        (                          u                          )                                                :                                            ∈                                                                        (                          u                          )                                                +                        μ
    // solve Hamilton- Jacobi and obtain new level-set ψi+1
8:  ∂                                                      ψ                            i                                                                                                    ∂                          t                                                                    -                                              v                        ⁢                                                                              ❘                            "\[LeftBracketingBar]"                                                                                ∇                                                                                      ψ                              i                                                                                                            ❘                            "\[RightBracketingBar]"                                                                                                                =                    0
9:  ψi+1←ψi
    // compute objective
10: J                                              i                        +                        1                                                              =                                                                                            1                          2                                                ⁢                                                                              ∫                            Ω                                                    D                                                                    ∈                                                                        (                          u                          )                                                :                                            ∈                                                                        (                          u                          )                                                ⁢                        d                        ⁢                        Ω
    // increment iteration
11: i ← i + 1
12: end while
    // obtain final level-set as solid region
13: Ωs←∫SDFT(ψi+1)

Note that generalizing this algorithm for arbitrary objectives and constraints can require modifications, such as the following. In line 6, any adjoint problems (see Adjoint Method above) needed for computing shape derivatives of each objective and constraint should be addressed. In line 7, the augmented Lagrangian method (see Augmented Lagrangian Method above) should be used for combining the different shape derivatives to yield a single advection velocity. For further details regarding implementing these modifications to generative design processes generally, see Augmented Lagrangian Algorithm for Constrained Shape Optimization below.

Geometry

Let Σ be a smooth, watertight, oriented surface in Euclidean space, with normal vector field N_Σ pointing in the direction of “free space” outside of Σ. A solid object can be created out of Σ by “thickening” it in the negative normal direction. That is: a small h∈ can be chosen to define

Ω_h:={x−sN_Σ(x): for all x∈[0,h]}  (33)

In plain English, Ω_h consists of all points sandwiched between and an offset of by a distance h in the negative normal direction, which is denoted

Σ_h:={x−hN_Σ(x): for all x∈Σ}  (34)

The boundary ∂χ_h of Ω_h consists of two disjoint surfaces: Σ itself and the offset surface Σ_h.

Note that what is commonly referred to as the “outward” unit normal vector field of ∂Ω_h (which points from inside the solid material defined by Ω_h to outside the solid material) is equal to N_Σ on Σ but points in the opposite direction to N_Σ on Σ_h. Its exact formula at y∈Σ_h is −N_Σ(proj_Σ(y)) where proj_Σ is the mapping that takes a point y to its closest point on Σ—so for example, if it is known that y=x−hN_Σ(x) and h is sufficiently small, then proj_Σ(y)=x.

Deformations of the Geometry

Σ can be allowed to deform. The deformation can be generated by a surface normal speed function which can be denoted as Θ^⊥:Σ→. Such a deformation can be achieved in a number of ways: for example, representing Σ as the zero level set of a function on the background Euclidean space and using advection with respect to (an extension of) Θ^⊥. The methods can be equivalent to first order in the deformation. What matters is the infinitesimal variation of Σ itself, which is precisely the surface normal speed function. Any deformation of Σ will have a notion of the “magnitude” of the deformation, which is a positive scalar ε. For example, if the deformation is generated by advection, then the magnitude of the deformation corresponds to the advection time. The deformed surface can be denoted by Σ_ε.

Once Σ deforms, then Σ_h deforms as well. It is simply “dragged along” in such a way that the offset distance h is maintained between Σ_ε and the deformed version of Σ_h. Consequently the deformation of the thickened object Ω_h is completely determined by the deformation of Σ. Note that it can be shown that the infinitesimal variation of Σ_h at a point y∈Σ_h is just −Θ^⊥ (proj_Σ(y)) where proj_Σ is the mapping that takes a point y to its closest point on Σ.

Optimization by Steepest Descent

An optimal Σ can be found for some objective function: Surfaces→ using a Steepest Descent method. This is based on the shape Taylor formula which says that the variation of Σ by a magnitude E with respect to variations generated by some Θ^⊥:Σ→ as above satisfies

Φ(Σ_ε)≈Φ(Σ)+εDΦ_Σ(Θ^⊥)  (35)

for sufficiently small ε, where DΦ_ε(Φ^⊥) symbolizes the shape derivative of Φ at Σ.

Using a formula for the shape derivative in terms of Θ^⊥, a particular Θ^⊥ can be chosen such that the shape derivative term in the shape Taylor formula is guaranteed to be negative. As a result, the objective function will decrease if Σ is varied with respect to this chosen Θ^⊥, at least for a sufficiently small ε. After Σ is updated by performing the variation (e.g., by advection of the level set function representing ε), the updated shape represents an improvement with respect to the objective function. For further improvement, the procedure can be iteratively repeated until convergence occurs.

Shape Derivative for a Class of Objective Functions

Consider objective functions of the following form. Let Φ₀: Domains→ be some “volumetric” shape function (e.g., one which can be evaluated on a volumetric domain) such as a convex linear combination of the averaged structural compliance with respect to a set of load cases and the total mass. This allows definition of a “surface” shape function by

Φ(Σ)=Φ₀(Ω_h)  (36)

It can be assumed that it is known how to calculate the shape derivative of Φ₀ at an arbitrary domain Ω and with respect to an arbitrary variation Ω_∈ of Ω generated by a boundary normal speed function V^⊥:∂Ω→ with respect to the outward unit normal vector field of Ω. That is, Φ₀ can be standard, meaning that the shape derivative satisfies the Hadamard-Zolésio structure theorem, and providing the formula

Φ₀(Ω_ε)≈Φ₀(Ω)+εDΦ_0,Ω(Θ^⊥) where DΦ_0,Ω(Θ^⊥):=∫_∂ΩG_Ω(x)V^⊥(x)dσ(x)  (37)

for sufficiently small ε. The function G_Ω of Ω, is known as the shape gradient of Φ₀ at the shape Ω. The precise form of G_Ω depends on Φ₀ and can be worked out by calculating the derivative

$\frac{d}{d\epsilon }{\Phi }_0\left({\Omega }_{\epsilon }\right),$

e.g., for averaged structural compliance and volume.

The shape derivative of Φ can be expressed in terms of the shape derivative of Φ₀. This can be done by leveraging the above formula for the shape derivative of Φ₀, applied to Ω_h, almost exactly, and taking into account the nature of the boundary normal speed on the two disjoint surfaces making up ∂Ω_h. That is: Σ can have the boundary normal speed V^⊥:=Θ^⊥, whereas Σ_h can have the boundary normal speed V^⊥:=−Θ^⊥∘proj_Σ. Then, the first-order change in Φ under a variation Σ_ε of Σ generated by the normal speed function Θ^⊥ is given by

Φ(Σ^ε)≈Φ(Σ)+ε(∫_ΣG_Ωh(x)Θ^⊥(x)dσ(x)−∫_ΣhG_Ωh(y)Θ^⊥(proj_Σ(y))dσ(y))  (38)

Thus, the desired shape derivative formula is

DΦ _Σ(Θ^⊥):=∫_ΣG_Ωh(x)Θ^⊥(x)dσ(x)−∫_ΣhG_Ωh(y)Θ^⊥(proj_Σ(y))dσ(y)  (39)

Note that y can be used as the dummy integration variable on Σ_h simply to emphasize that the two integrals are over different spaces and cannot a priori be combined in any simple way.

Extracting a Descent Direction

Recall that the utility of the exact formula for the shape derivative is that it should allow Θ^⊥ to be chosen in such a way that the shape derivative becomes negative, resulting in a decrease of the objective function to first order according to the shape Taylor formula. It is, however, not immediately clear how to do this for the shape derivative of Φ calculated above. This is due to the fact that there are two competing terms (e.g., the integral over Σ and the integral over Σ_h) and it is not clear how these balance out given Θ^⊥.

There are two ways to proceed. The first way is to apply the Hilbert space method of extracting a descent direction. The second way is to apply a change of variables to the integral over Σ_h, allowing it to be expressed as an integral over Σ. This is simple: recall that any point y∈Σ_h can be written in the form y=x−hN_Σ(x); for our purpose now, this should be viewed as a mapping from Σ to Σ_h. Call this mapping n:Σ→Σ_h, with n(x):=x−hN_Σ(x). So by the change of variables formula for surface integrals

$\begin{array}{cc}{\int }_{{\Sigma }_h}{G}_{{\Omega }_h}\left(y\right){\theta }^{\perp }\left(pro{j}_{\Sigma }\left(y\right)\right)d\sigma \left(y\right)={\int }_{n\left(\Sigma \right)}{G}_{{\Omega }_h}\left(y\right){\theta }^{\perp }\left(pro{j}_{\Sigma }\left(y\right)\right)d\sigma \left(y\right)& \left(40\right)\end{array}$ $\begin{array}{cc}={\int }_{\Sigma }{G}_{{\Omega }_h}\left(x-{\mathrm{hN}}_{\Sigma }\left(x\right)\right){\theta }^{\perp }\left(pro{j}_{\Sigma }\left(x-h{N}_{\Sigma }\left(x\right)\right)\right)Jac\left(x\right)d\sigma \left(x\right)& \left(41\right)\end{array}$ $\begin{array}{cc}={\int }_{\Sigma }{G}_{{\Omega }_h}\left(x-h{N}_{\Sigma }\left(x\right)\right){\theta }^{\perp }\left(x\right)Jac\left(x\right)d\sigma \left(x\right)& \left(42\right)\end{array}$

where Jac is the Jacobian of n and is based on the fact that proj_Σ(x−hN_Σ(x))=x for all x∈Σ.

It turns out that the Jacobian of n can be determined. With some help from differential geometry, it can be shown that Jac(x)=1+hH_Σ(x)+h²K_Σ(x) where H_Σ is the mean curvature of Σ and K_Σ is the Gauss curvature of Σ, providing

DΦ _Σ(Θ^⊥)=∫_Σ(G_Ωh(x)−G_Ωh(x−hN_Σ(x))(1+hH_Σ(x)+h²K_Σ(x)))Θ^⊥(x)dσ(x)  (43)

This manipulation of DΦ_Σ(Θ^⊥) provides a nice formula for the shape gradient of Φ at Σ, denoted as

G _Σ(x):=G_Ωh(x)−G_Ωh(x−hN_Σ(x))(1+hH_Σ(x)+h²K_Σ(x)).  (44)

Consequently, the choice of Θ^⊥=−G_Σ makes DΦ_Σ(Θ^⊥) negative, as required for our steepest descent procedure updating Σ towards optimality of Φ.

Augmented Lagrangian Algorithm for Constrained Shape Optimization

Introduction

An algorithm can be developed that solves constrained shape optimization problems where a shape function is to be minimized (e.g., volume or compliance) subject to a mix of equality and inequality constraints. These constraints can use shape functions (e.g., a target volume, an aggregated stress measure, or an aggregated displacement measure). A different algorithm can be used for dealing with non-aggregated (a.k.a. pointwise) constraints, such as the norm of the stress at each point of the shape. Generically, the optimization problems can be of the form

$\begin{array}{cc}\underset{{\Omega }_{\in Adm}}{\min }F\left(\Omega ,u_{\Omega }\right)& \left(45\right)\end{array}$ $\begin{array}{cc}s.t.G_E\left(\Omega ,u_{\Omega }\right)=0& \left(46\right)\end{array}$ $\begin{array}{cc}{G}_I\left(\Omega ,u_{\Omega }\right)=0& \left(47\right)\end{array}$

for u_Ω satisfying a linear elastic PDE in Ω

where Adm is a set of admissible shapes, F is a scalar-valued shape-differentiable shape function, and G_E and G_I are vector-valued shape-differentiable shape functions that can depend explicitly on the shape, or implicitly on the shape through its elastic response under loads.

The Classical Augmented Lagrangian Algorithm

Equality Constraints

The augmented Lagrangian method in the classical setting of optimizing a scalar function of a vector variable can apply most simply to equality-constrained optimization problems of the form

$\begin{array}{cc}\underset{x\in {ℝ}^n}{\min }f\left(x\right)& \left(48\right)\end{array}$ $s.t.g_i\left(x\right)=0\mathrm{for}i=1,\dots ,k$

The augmented Lagrangian method is an enhancement of a so-called penalty method. The penalty method considers a penalized objective function

$L_c\left(x\right):=f\left(x\right)+\frac{c}{2}{{\Sigma }_i\left[g_i\left(x\right)\right]}^2$

and a sequence c:=c_k tending to infinity. Then, one can define x_k:=argmin_x L_ck(x). It is by no means the case that g_i(x_k)=0 for any i. But since the increasing sequence of c-values penalizes the constraints more and more heavily, it can be expected that x_k eventually satisfies the constraints, in the sense that if x_k converges to x* then g_i(x*)=0 for all i. Additionally, it can be expected that x* is a solution of Equation 48.

A problem with the penalty method is that the increasing sequence of c-values can make it harder and harder to minimize L_c(x) numerically due to poor conditioning. The augmented Lagrangian method can remedy this problem. This algorithm minimizes an “augmented Lagrangian”

$\begin{array}{cc}{L}_c\left(x,\mu \right):=f\left(x\right){\sum }_i{\mu }_k^i{g}_i\left(x\right)+\frac{c}{2}{{\sum }_i\left[g_i\left(x\right)\right]}^2& \left(49\right)\end{array}$

and also maintains a sequence of Lagrange multipliers in and increasing penalty parameters c_k. Again, one can define x_k:=argmin_x L_ck(x,μ_k). These sequences are updated in such a way that c_k stabilizes on some large but finite value (hence avoiding the ill-conditioning inherent to general penalty methods). This occurs while x_k converges to a solution of Equation 48 (it is thus optimal and feasible) and μ_k converges to the Lagrange multiplier associated to in the KKT conditions. Since the augmented Lagrangian method iterates on both the primal variable x and the “dual” variable μ, this algorithm is an example of a primal-dual optimization method.

The general framework for the augmented Lagrangian algorithm can require a method to increase the penalty parameters and a method for tightening tolerances

1. Choose a final tolerance Tfinal.
2. Set k = 0. Start with an initial tolerance c0, μ0. Start with an initial tolerance T0
   less restrained than Tfinal.
3. While not converged to within tolerance Tfinal:
   a. Apply an unconstrained optimization algorithm to the problem
      minxLck (x, μk). Stop the algorithm when a convergence to within
      tolerance Tk is reached. Output Xk+1 to be the x-value at this level of
      convergence.
   b. Check constraint satisfaction
      (i)  If the constraints are not satisfied to within tolerance Tk, let Ck+1 :=
           Increase (ck) and μk+1 := μk.
      (ii) If the constraints are satisfied to within tolerance Tk, update the
           Lagrange multipliers by μk+1i := μki + ckgi(xk+1).
   c. Set the tolerance Tk+1 := Tighten (Tk).
   d. Increment k.

Note the feature of the above algorithm that is the update of the Lagrange multiplier. It can be understood as follows. The unconstrained minimum of the augmented Lagrangian at iteration k satisfies ∇L_ck(x_k,μ_k)=0 or

∇ƒ(x_k)+Σ_i(μ_kⁱ+c_kg_i(x_k))∇g_i(x_k)=0  (50)

Of course, the solution x* and the Lagrange multipliers μⁱ*, for the equality-constrained optimization Equation 48 satisfy

∇ƒ(x*)+Σ_iμⁱ*∇g_i(x*)=0  (51)

g ⁱ(x*)=0 for i=1, . . . ,k  (52)

In the augmented Lagrangian algorithm, since it is desirable to ensure that the iterates x_k converge to x, it is desirable to have μ_kⁱ→μ* and g_i(x_k)→0 are consistent with this aim. Therefore, the Lagrange multiplier update can be construed as a fixed-point iteration scheme in the form μ_k+1ⁱ:=μ_kⁱ+c_kg_i(x_k+1) to achieve this.

Inequality Constraints

An elegant extension of the classical equality-constrained augmented Lagrangian algorithm to the classical inequality-constrained case can hinge on two facts. For fact 1, the problem

$\begin{array}{cc}\frac{\min }{x\in {ℝ}^n}f\left(x\right)& \left(53\right)\end{array}$ $s.t.g_i\left(x\right)\le 0\mathrm{for}i=1,\dots ,k$

is equivalent to the problem

$\begin{array}{cc}\underset{\left(x,z\right)\in {ℝ}^n\times {ℝ}^k}{\min }f\left(x\right)& \left(54\right)\end{array}$ $s.t.g_i\left(x\right)+z_i^2=0\mathrm{for}i=1,\dots ,k$

The new z-variables in Equation 54 are known as slack variables and satisfy z_i=√{square root over (−g_i(x*))} at the feasible optimum, where x* is the feasible optimum for Equation 53.

For fact 2, the augmented Lagrangian for Equation 54 is

$\begin{array}{cc}{L}_c\left(x,z,\mu \right):=f\left(x\right)+{\sum }_i{\mu }^i\left[g_i\left(x\right)+z_i^2\right]+\frac{c}{2}{{\sum }_i\left[g_i\left(x\right)+z_i^2\right]}^2& \left(55\right)\end{array}$

and the unconstrained optimization arising in the augmented Lagrangian algorithm can be decomposed and partially solved as follows:

$\begin{array}{cc}\underset{\left(x,z\right)\in {ℝ}^n\times {ℝ}^k}{\min }{L}_c\left(x,z,\mu \right)=\underset{x\in {ℝ}^n}{\min }\underset{z\in {ℝ}^k}{\min }{L}_c\left(x,z,\mu \right)=\underset{x\in {ℝ}^n}{\min }f\left(x\right)+\frac{\min }{x\in {ℝ}^k}{\sum }_i{\mu }^i\left[g_i\left(x\right)+z_i^2\right]+\frac{c}{2}{{\sum }_i\left[g_i\left(x\right)+z_i^2\right]}^2& \left(56\right)\end{array}$ $\begin{array}{cc}=\frac{\min }{x\in {ℝ}^n}f\left(x\right)+{\sum }_i{\mu }^i{g}_i^{+}\left(x,{\mu }^i,c\right)+\frac{c}{2}{{\sum }_i\left[g_i^{+}\left(x,{\mu }^i,c\right)\right]}^2& \left(57\right)\end{array}$

which has introduced

$g_i^{+}\left(x,\mu ,c\right):=\max \left\{g_i\left(x\right),-\frac{{\mu }^i}{c}\right\}.$

This is because the minimization can be performed over z∈ explicitly using elementary calculus techniques and phrase the result in terms of the g_i⁺ function.

The upshot is that one can solve the inequality constrained Equation 53 by applying the augmented Lagrangian algorithm (see Equality Constraints above) to the modified Lagrangian function

$\begin{array}{cc}{\tilde{L}}_c\left(x,\mu \right):=f\left(x\right)+{\sum }_i{\mu }^i{g}_i^{+}\left(x,{\mu }^i,c\right)+\frac{c}{2}{{\sum }_i\left[g_i^{+}\left(x,{\mu }^i,c\right)\right]}^2& \left(58\right)\end{array}$

whose gradient (after some manipulation) is

∇{tilde over (L)}_c(x,μ)=∇ƒ(x)+Σ_i max{0,μⁱ+cg_i(x)}∇g_i(x).  (59)

Note that for completeness, the solution of the minimization over the z-variables appears in Equation 56. Substituting for this solution yields Equation 57. For simplicity, the minimization over the z-variables is rephrased in the form

$\begin{array}{cc}\begin{array}{c}\min \\ Z\in ℝ\end{array}a\left(s+z^2\right)+{b\left(s+z^2\right)}^2& \left(60\right)\end{array}$

where a, b, s are fixed parameters, which can be solved as follows.

The first step is to substitute y:=z² and replace the above minimization problem by a(s+y)+b(s+y)². This new problem is now quite simple, since the function Φ(y):=a(s+y)+b(s+y)² is a simple parabola, and the problem Φ(y) finds the minimum of this parabola constrained to the region y≥0. Therefore the global constrained minimum is either at the global unconstrained minimum of Φ or else at the y=0 boundary of the constrained region, whichever is smaller. This results in either y*=−a/2b−s, or y*=0; from which is selected the result for which Φ*:=Φ(y*) is least. Consequently

$\begin{array}{cc}\begin{array}{cc}{\Phi }_{*}& =\min \left\{\Phi \left(0\right),\Phi \left(-a/2b-s\right)\right\}\\ & =\min \left\{\mathrm{as}+b{s}^2,-a^2/4b\right\}\\ & =a\max \left\{s,-a/2b\right\}+{b\left[\max \left\{s,-a/2b\right\}\right]}^2\end{array}& \left(61\right)\end{array}$

after algebraic re-arrangement. With s=g_i(x) and a=μⁱ and =c/2, Equation 57 is obtained.

Application to Constrained Shape Optimization Problems

The augmented Lagrangian algorithm can be applied to solve the Equations 45-47 with one scalar equality constraint and one scalar inequality constraint for simplicity. Then, Ω∈Adm can be interpreted to mean that Ω satisfies admissibility constraints of the following form. The surface of every Ω∈Adm contains pre-defined ports; every Ω∈Adm contains pre-defined keep-in regions; and every Ω∈Adm avoids pre-defined keep-out regions.

The classical augmented Lagrangian algorithm (see Equality Constraints above) can be mimicked for both equality and inequality constraints but can be adapted for shapes satisfying admissibility constraints. The augmented Lagrangian is defined as

$\begin{array}{cc}{L_c\left(\Omega ,{\mu }^E,{\mu }^I\right):=F\left(\Omega \right)+{\mu }^E{G}_E\left(\Omega \right)+{\frac{c}{2}\left[G_E\left(\Omega \right)\right]}^2+{\mu }^I{G}_I^{+}\left(\Omega ,\mu ,c\right)]}^2+{\frac{c}{2}\left[G_I^{+}\left(\Omega ,\mu ,c\right)\right]}^2& \left(62\right)\end{array}$

with G_I⁺(Ω, μ, c):=max{G_I(Ω), −μ^I/c}. The dependence of F, G_E, G_I on the solution of μ_Ω of the linear elastic PDE has not been indicated for simplicity. To consider only shape functions that are surface or volume integrals of shape-dependent quantities, it is known from general principles that the shape gradient of a function L of this type, evaluated at a given shape Ω, is a scalar-valued function of ∂Ω that is denoted as dL(Ω):∂Ω→. The shape gradient of the augmented Lagrangian at the shape Ω is

dL _c(Ω,μ_E,μ^I)=dF(Ω)+(μ^E+cG_E(Ω))dG_E(Ω)+max{0,μ^I+cG_I(Ω)}dG_I(Ω)  (63)

where dF(Ω), dG_E(Ω), dG_I(Ω) are the shape gradients of the objective and constraint function, respectively, at the shape Ω.

It is also known from general principles that updating the shape Ω to reduce the value of L_c while maintaining the admissibility constraint can be achieved by advecting the level set function representing Ω for a certain time (determined via a line-search or related procedure) with respect to a speed function constructed from an extension of a projection of the shape gradient of L_c (e.g., the shape gradient of L_c can be zeroed out appropriately when a violation of admissibility is detected and use the Adalsteinsson-Sethian velocity extension algorithm to extend the values of dL_c(Ω) from ∂Ω to a narrow band of ∂Ω). This is the basis for an unconstrained, gradient-based shape optimization algorithm for solving the problem min_Ω∈AdmL_c(Ω, μ^E, μ^I) at any fixed values of c, μ^E, μ^I.

The augmented Lagrangian algorithm for solving Equations 45-47 can now be given as follows. It requires a method to increase the penalty parameters and a method for tightening tolerances.

1. Choose a final tolerance Tfinal
2. Set k = 0. Start with initial c0, μ0E, μ0I. Start with an initial tolerance
T0 less restrictive than Tfinal.
3. Initialize the shape Ω0.
4. While not converged to within tolerance Tfinal :
   a. Apply an unconstrained, gradient-based shape optimization
   algorithm                    ⁢                                        to                    ⁢                                        the                    ⁢                                        problem                    ⁢                                                              min                                              Ω                        ∈                        Adm                                                              ⁢                                                                  L                                                  c                          k                                                                    (                                              Ω                        ,                                                  μ                          k                          E                                                ,                                                  μ                          k                          I                                                                    )                                        ⁢                                        until
   convergence occurs to within tolerance Tk. Output the shape
   Ωk+1 at this level of convergence
   b. Check constraint satisfaction.
   (i) If the constraints are not satisfied to within tolerance Tk, let
   ck+1 := Increase (ck) multipliers by μk+1E := μkE and μk+1I :=
   μkI.
   (ii) If the constraints are satisfied to within tolerance Tk,
   update the Lagrange multipliers by μk+1E := μkE + ckGE(Ωk)
   and μk+1I := max{0, μkI + ckGI(Ωk+1)}.
   c. Set the tolerance Tk+1 := Tighten (Tk).
   d. Increment k.

Pre-Advection Operations

In some implementations, one or more operations are performed on the velocity field before the shape is advected, including (1) narrow-band velocity restriction, (2) advection prevention using an advection mask, and/or (3) velocity extension. For the first of these, narrow-band velocity restriction involves restricting the velocity to a narrow-band around the 0^th iso-contour of the level-set. For the second, the advection mask has values of 0 inside the ports (geometric interfaces containing Neumann and Dirichlet boundary conditions) and values of 1 elsewhere in the domain. Shape derivatives and advection velocity is multiplied by the advection mask to prevent any advection inside the ports. With regard to the third of these, the velocity field should be continuous on both negative and positive sides of the 0^th iso-contour of the level-set. However, the objective function (typically strain energy) is often only available inside the negative narrow-band. Velocity extension projects all velocities inside the domain to the positive narrow-band. This can be done by sampling the velocity at a point inside the domain found by moving along the negative the normal a distance equal to the level-set:

v(x)=v(x−ψ(x)n(x))x∉Ω.

An example of the first and third operations is shown in FIGS. 2C-2E.

FIG. 2C an example of a problem set up 200 for a level-set based topology optimization process. The problem set up 200 is intended to generate a pylon for a vehicle, where a design space 202 is specified with respect to ports 204 (also referred to as preserve geometries or preserve sub-spaces) for the pylon to be generated. A cross-section X, which is shown in perspective in FIG. 2C, is also shown in FIGS. 2D and 2E. FIG. 2D shows a narrow-band velocity 210. The number of voxels in the narrow band is denoted by w_nb, and Δs denotes the size of a voxel. FIG. 2E shows the effects 212 of velocity extension with respect to the 0^th iso-contour 214 of the level-set.

The above approaches can be used with various types of level-set-based topology optimizations, which are described further below in connection with FIGS. 3A-8C. Returning to FIG. 1B, all of the generative design processes described in this document can be implemented, e.g., in CAD program(s) 116, to provide both (1) substantial user control over the generative design process, as well as post-processing of the generative design output to form a final, acceptable 3D model of an object, and (2) control functions for providing the generatively designed 3D model for use in manufacturing a physical structure corresponding to the object. Thus, the results of generative design processing can be presented to the user, e.g., in UI 122 on display device 120, along with an option 190 to accept or reject the design.

If the design is rejected, the process of FIG. 1B can return to obtain 180, e.g., by CAD program(s) 116, new design variable(s) for use in producing a new generative 3D model. For example, a new design variable can include selecting 180 a different void addition and/or removal technique for changing the topology of the three dimensional shape during shape optimization. For example, in implementations that use level-set based topology optimization methods, one or more void(s) can be introduced into the material domain based on topological derivatives of an objective function during the process of the structural optimization to allow a change in topology (a change in configuration) such as the introduction of a hole in the material domain. Further details regarding one or more examples of void introduction techniques are described below in connection with FIG. 3A. Moreover, the design variable(s) obtained 180 can include all of the various user inputs described in this document that can affect the generative design process.

Once a design is not rejected 190, the process of FIG. 1B can provide 195, e.g., by CAD program(s) 116, the 3D model for use in manufacturing a physical structure corresponding to the object using one or more computer-controlled manufacturing systems, e.g., AM machine 170, SM machine 174, and/or other manufacturing machines. The providing 195 can involve sending or saving the 3D model to a permanent storage device for use in manufacturing the physical structure corresponding to the object using the one or more computer-controlled manufacturing systems. In some implementations, the providing 195 involves generating 195A, e.g., by CAD program(s) 116, toolpath specifications for the computer-controlled manufacturing system(s) using the 3D model, and manufacturing 195B, e.g., by CAD program(s) 116, at least a portion of the physical structure corresponding to the object with the computer-controlled manufacturing system(s) using the toolpath specifications generated for the additive manufacturing machine.

Note that the 3D model that is provided 195 can be the 3D model produced 185 by a generative design synthesis method or a post-processed version of the generative design output. For example, in some implementations, a 3D mesh model produced by a generative design synthesis method can be converted into a watertight B-Rep 3D model before being provided 195. In some implementations, a generative design output can be post-processed using the systems and techniques described in U.S. Patent Application No. 62/758,053, filed Nov. 9, 2018, and titled “CONVERSION OF GENERATIVE DESIGN GEOMETRY TO EDITABLE AND WATERTIGHT BOUNDARY REPRESENTATION IN COMPUTER AIDED DESIGN,” which systems and techniques are included in U.S. Pat. No. 11,016,470, issued on May 25, 2021, which claims priority to U.S. Patent Application No. 62/758,053 and is titled “CONVERSION OF MESH GEOMETRY TO WATERTIGHT BOUNDARY REPRESENTATION.” Further, in some implementations, the post-processed generative design output can be edited using the systems and techniques described in U.S. application Ser. No. 16/186,136, filed on Nov. 9, 2018, issued as U.S. Pat. No. 10,467,807 on Nov. 5, 2019, and titled, “FACILITATED EDITING OF GENERATIVE DESIGN GEOMETRY IN COMPUTER AIDED DESIGN USER INTERFACE.” Moreover, the generative design processes described below in connection with FIGS. 3A-8C can also be implemented to use the post-processing, editing, and/or provisioning 195 systems and techniques described above. Finally, while described in the context of a CAD program providing multiple options for generative design, each of the generative design processes described in this document can be implemented as a standalone generative design process in a CAD program. Thus, not all of the generative design processes described below in connection with FIGS. 3A-8C need be implemented together in any given implementation.

FIG. 3A shows an example of a process (e.g., as performed by the CAD program(s) 116 of FIG. 1A) that generates one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes, as described in further detail below in connection with FIGS. 3A-8C. The process of FIG. 3A is an example of the defined process 185 from FIG. 1B. A design space, a setup for numerical simulation, and other inputs for a selected generative design process are identified for use in initiating 300 a macrostructure (or geometrical) type of generative model (e.g., a level-set representation of the 3D model) that can be used with one or more selected generative design processes.

When the generative process to be used employs a basic level-set method for topology optimization, the level-set is initiated 300 for the design space. Note that a level-set method is an example of a macrostructure generative modeling technique, where the generative model represents the object being designed as one or more boundaries between one or more solid regions (having the homogenous material therein) and one or more void regions (having no material therein) within the design space. Further, identifying the inputs can involve receiving user input, e.g., via UI 122 on display device 120, importing information from another program or a third party source, and/or one or more of these inputs can be predefined in a given implementation.

The setup for numerical simulation can include one or more physical properties to be simulated and one or more types of simulation to be performed, as discussed above, as well as potentially surrogate modelling or other methods of approximation. In some implementations, the type of numerical simulation is predefined, either for all uses of the program or given a particular context in the program from which the generative design process has been launched. Further, the setup for numerical simulation includes at least one set of loading conditions and/or other physical environment information associated with the type of numerical simulation to be performed.

The design space can be an initial 3D model or one or more portions of an initial 3D model to be used as a starting geometry. In some cases, the design space can be determined as the bounding volume of all initial geometry used to specify loading or other physical environment information associated with the type of numerical simulation to be performed. In some cases, the design space can be unbounded. In some implementations, the portion of the design space to be used as starting geometry can be automatically set by a genetic algorithm or other process. For example, bubble-like holes (e.g., holes 132B) can be put in the domain and a genetic algorithm can be used to vary the bubble size and spacing.

Seeding and the Bubble Method

Initial Seeding

The design space can be initialized using a seeding process, in which the design space is defined by a Boolean intersection between the design space D. and a seed geometry Ω_s, shown as follows:

Ω₀=Ω∩Ω_s  (64)

where Ω₀ is the resulting initialized domain after the seeding geometry Ω_s is applied to the initial design space Ω. The seed geometry can be of a variety of different shapes, e.g., an array of bubbles or a mesh, with parametrized characteristics, e.g., bubble diameter and spacing. The parameters can be user-defined or defined automatically by some seeding process, e.g., bubble seeding, described below.

Initial seeding can result in more efficient optimization and flexibility in design variations. For example, initial seeding can be defined so as to avoid local minima and the need to restart optimizations. Further, initial seeding can aide in creating design variations. In some implementations, a seed geometry is user-defined, e.g., according to a final geometry generated from a previous execution of design process, a random process to randomly initialize the seeding geometry, according to other factors of interest to the user, or some combination of the preceding.

FIG. 3B shows a graphical representation of examples of initial seeding geometries applied to a design space and resulting designs for a Michell type arch problem after a generative design process is executed on the design space. Seeding geometries 314 and 316 represent different geometries for initializing a design space, and can be user-defined or according to some seeding process, e.g., bubble seeding, described below. The seeding geometry 316 is the same as the seed model 132 of FIG. 1A. Final designs 318 and 320 are designs generated by applying the same generative design process after a design space is seeded according to geometries 314 and 316, respectively. The final design 320 is the same as the complete structure 138 of FIG. 1A. The resultant final designs 318 and 320 have different geometries and different physical properties, e.g., different strain energies. Therefore, the seeding geometry can have a significant impact in a final design for different design spaces, even if the generative design process is the same.

Moreover, the other inputs can depend on the type of numerical simulation to be performed and/or the type of generative design process to be used. For example, when a lattice will be used in the generative design process, the other inputs can include, lattice topology, volume fraction, unit size, and thickness. Various types of generative design processes can be used, either alone or in combination, but FIG. 3A shows a representative process for all these various types of generative design processes. After initiation 300 of the macrostructure type generative model, an iterative process of modifying this generative model is begun in order to satisfy the design criteria for the physical structure, e.g., one or more design constraints to be met and one or more design objectives to be maximized.

Numerical simulation of the physical response of the current model (e.g., a level-set representation of the implicit surface of the 3D shape) is performed 302 under the one or more defined loading conditions. In general, the numerical simulation 302 treats each solid region in the current model as a homogenous solid material, and each void region in the current model as having no material therein. However, in some implementations, this treatment can be altered. For example, in hybrid topology optimization, the numerical simulation of the modelled object is performed 302 with at least a portion of the solid region being treated as having at least one void in it (e.g., what the macrostructure generative modeling technique treats as a solid, the numerical simulation treats as partially containing a void in the form of a lattice) or at least a portion of the void region being treated as having at least one solid in it (e.g., what the macrostructure generative modeling technique treats as a void, the numerical simulation treats as partially containing a solid material in the form of a lattice). As another example, in the case of hollow topology optimization, the numerical simulation of the modelled object is performed 302 with at least a portion of the solid region being treated as having a void in it (i.e., what the macrostructure generative modeling technique treats as a solid, the numerical simulation treats as containing a hollow region). Finally, in a combination of these two, in hybrid hollow topology optimization, the numerical simulation of the modelled object is performed 302 with at least a portion of the solid region being treated as having both a partial void region and a complete void in it (i.e., what the macrostructure generative modeling technique treats as a solid, the numerical simulation treats as partially containing a lattice structure around a hollow region). For further details regarding hybrid, hollow and hybrid-hollow topology optimization, see PCT/US2019/060089, filed Nov. 6, 2019, and published as WO 2020/097216 A1 on May 14, 2020.

Results from the simulation are used to update 304 the current model (e.g., the level-set representation) in accordance with a current numerical assessment of the simulated physical response of the current model. For example, shape change velocities can be calculated for the implicit surface in the level-set representation of the 3D shape of the object being modeled, and this level-set representation can be updated 304 using the calculated shape change velocities to produce an updated version of the 3D shape of the modeled object. Note that in various implementations, other operations can be performed before, during or after the updating 304, as part of the iterative modification of the 3D shape of the modeled object, such as described below in connection with FIGS. 4A-8C.

In addition, to produce topology changes in the current model of the object, one or more voids can be inserted 308 into the current model at each or selected ones of the iterative modifications of the current model. For example, one or more bubbles having positions, sizes and shapes determined from the current model can be inserted 308. Moreover, in some implementations, the insertion 308 of void(s) is done only during an earlier portion of the iterative modification and/or only periodically (e.g., at a regular void insertion interval) during the iterative modification. In the example of FIG. 3A, the insertion of void(s) is done only when the current iteration is less than a predefined void insertion cutoff, and then only when the current iteration is equal to the void insertion interval.

Bubble Method

In general, an optional determining 306 and the inserting 308 can be performed as part of a bubble method. A bubble method allows for the topology of the design space to be changed from the inside, as default level-set methods only allow for changes from the boundary. A location for the bubble is identified using the topology derivative, for example computed using a topology-shape sensitivity method, which relates the shape derivative to the topology derivative.

In some implementations, a bubble method is applied with the following features:

1. Position: The shape derivative of the Lagrangian is used as a proxy for the topology derivative. The shape derivative is described in more detail below, with respect to Equation 105.

2. Frequency: Bubbles are not inserted at every iteration. Instead, bubbles are inserted according to user-defined or automatically determined intervals Bubbles are not inserted after the volume reduction iterations have been completed.

3. Size: The size of an inserted bubble at a given iteration is computed according to a user-defined or automatically determined ratio Pb of the current model volume:

V(B_t)=β_bV(Ω_t).  (65)

where V(B_t) is the volume of the bubble inserted at an iteration t, and V(Ω_t) is the volume of the model at the iteration t.

4. Shape: The shape of the inserted bubble is determined based on the distribution of the topology derivative for elements of the current model, as shown below:

$\begin{array}{cc}\frac{dℒ}{d\Omega }\left(e_1\right)<\frac{dℒ}{d\Omega }\left(e_2\right)<\cdots \frac{dℒ}{d\Omega }\left(e_k\right)\cdots <\frac{dℒ}{d\Omega }\left(e_n\right)& \left(66\right)\end{array}$

Where e_k is the k-th element of the current model at a given iteration t. In some implementations, the shape of the inserted bubble is optimized while the overall size of the inserted bubbles increases over multiple iterations.

Then, the volume of elements with the lowest shape derivative are added up, until the necessary bubble volume, i.e., according to Equation 64, is reached. The resultant shape, position, and volume of the bubble B_t generated at iteration t therefore satisfies the following:

B _t =e ₁ ∪e ₂ ∪ . . . ∪e _k s.t. V(B_t)=V(e₁)+V(e₂)+ . . . +V(e_k)  (67)

Returning to FIG. 3A, after the optional determining 306 that the current iteration is less than a void insertion cutoff and is equal to a void insertion interval, one or more voids are inserted 308 into the current model. The void insertion cutoff and the void insertion interval can be predetermined, e.g., by user input or automatically determined, prior to performing 302 the numerical simulation, or prior to initiating 300 the model. The one or more voids can be inserted according to a bubble method, e.g., the bubble method described above with reference to Equations 65-67.

FIG. 3C shows a graphical representation of an example of a bubble insertion history 321 for intermediate models for a Michell type arch problem during execution of a generative design process. Beginning with model Ω₀ 322A at iteration 0, the bubble iteration history 321 includes model Ω₁ 322B at iteration 1; model Ω₁₀ 322C at iteration 10; model Ω₁₁ 322D at iteration 11 (two bubbles inserted); model Ω₂₀ 322E at iteration 20; model Ω₂₁ 322F at iteration 21 (one bubble inserted); model Ω₃₀ 322G at iteration 30; and model Ω₃₁ 322H at iteration 31 (two bubbles inserted), with bubble insertions at t={0, 10, 20, 30}. Also in this example, the bubble volume ratio Pb is set to 0.05, meaning that the volume of the inserted bubbles are each equal to 5 percent of the volume of the current domain, i.e., consistent with Equation 65.

Although the bubble insertion history 321 shows bubbles inserted after a first iteration, bubbles can also be inserted during the first iteration of a generative design process. In some implementations, multiple instances of a bubble method can be executed, with different initial conditions and parameters that result in different design variations. A system implementing the described techniques can provide the results of the design variations for the user to select a preferred variation. In some implementations, the system can automatically select a design variation accordingly to predetermined, e.g., user-provided, criteria.

Returning to FIG. 3A, a check 310 can be made for convergence during each iteration. To determine whether the generative design has converged to a stable solution, the check 310 can identify the condition in which all the design constraints are met and no design objectives have improved significantly since the last one or more iterations. The numerical simulation 302, update 304, any void insertion 308, and check 310 process iterates until convergence. Further, in some implementations, the iterative process ends once a check 312 shows that a predetermined number of shape modification iterations have been completed. Note that the predetermined number of shape modification iterations can be set high enough that it is essentially guaranteed that all design constraints will be met, and using the controlled convergence techniques described below with reference to FIGS. 4A-4H, it can also be essentially guaranteed that convergence to generally optimal values for any design objectives will be achieved within the predefined number of iterations, if not before.

Controlled Convergence

In a controlled convergence technique, a pre-determined time period or number of iterations is defined as the shape and topology optimization begins. For each design constraint, a target much closer to the current value for the constraint is specified for each iteration.

This can result in at least two improvements. First, a user may specify the time period or number of iterations. The convergence rates are then controlled to solve a given problem within the user-specified time period or number of iterations. Therefore, a design process using controlled convergence can arrive at a suitable solution given resource restrictions imposed by the user in how much time is allotted to generating a solution. Second, in specifying a target value for each constraint closer to the current value for the constraint, the likelihood of oscillations and radical changes can be reduced.

FIG. 4A shows an example of a process that generates one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes implementing arbitrary constraint handling with controlled convergence. Controlled convergence is described first, and then arbitrary constraint handling is introduced. A computer-aided design program obtains 416 a design space for a modeled object, for which a corresponding physical structure will be manufactured, one or more design criteria for the modeled object, and one or more in-use load cases for the physical structure, wherein the one or more design criteria include at least one design constraint. For example, the obtaining 416 can be done as part of the initiating 300, as described above with reference to FIG. 3A. As described above with reference to FIGS. 3B-3C, seeding techniques can be applied to the design space which can further improve the generative design process.

An iteration amount is identified 418. The iteration amount can be a time or a count of iterations from a user or calculated otherwise. A series of target values for the at least one design constraint is calculated 420 from an initial target value to a final target value, in accordance with the iteration amount and a function.

Let g^j(s,u(s))−g_T,n^j=0 j=1, . . . , n_g denote a sequence of constraints where g_T,n^j is the final target value at iteration n. Define a target value for each constraint to meet in iteration t as

$\begin{array}{cc}{g}_{T,t}^j=g_{T,n}^j+\frac{N^d\left(\xi \right)}{N^d\left(0\right)}\left(g_0^j-g_{T,n}^j\right)& \left(68\right)\end{array}$

where

$\xi =\frac{t}{n}$

and N^d(ξ) is a B-spline of degree d computed using the recurring relationship

$\begin{array}{cc}\mathrm{Nd}\left(\xi \right)=\frac{1}{d}\left(1-\xi \right)N^{d-1}\left(\xi \right)& \left(69\right)\end{array}$

where N⁰(ξ)=0.

Although this description assumes the use of B-splines as the function for the calculating 420, any suitable smooth function can be used to achieve the same result. Different smooth functions include different orders of functions of the same class, e.g., B-Splines of different orders. The choice in smooth function and order of the function changes how quickly an ending design constraint target value is achieved over the identified 418 iteration amount. Other examples of smooth functions that can be used include polynomials of any order, Lagrange polynomials, and subdivision curves. In some implementations, a system implementing the described techniques can prompt a user for a plurality of reference points, and in response to receiving user-selected reference points, the system can generate a curve, e.g., using interpolation or any appropriate technique, that passes through the points.

FIG. 4B shows examples of graphs 401A, 401B tracking progressive reduction and increase of a target constraint value using B-splines, respectively. B-spline based transition from a starting value of constraint g₀^j to final target g_T,t^j is done during t≤n_v iterations. The final target value is maintained during iterations n_v<t≤n. In the graphs 401A, 401B, the target value of the constraint is given by:

$\begin{array}{cc}{g}_{T,t}^j=\{\begin{array}{cc}{g}_{T,n}^j+\frac{N^d\left(\xi \right)}{N^d\left(0\right)}\left(g_0^j-g_{T,n}^j\right)& \mathrm{if}t\le n_v\\ g_{T,n}^j& \mathrm{otherwise}\end{array}& \left(70\right)\end{array}$

Note particularly that iterations are divided into two parts according to Equation 70. A first part, while the iterations are less than or equal to n_v, includes different target values for different iterations in that first part. In the second part, i.e., after n_v, the target value of the constraint at each iteration is the final target value g_T,n^j.

Graph 401A shows a constraint target g_T,n^j=0.35 from a starting constraint value of 0.9, n_v=30, n=50, and B-Splines of order d=4. Graph 401B shows a constraint target g_T,n^j=0.85 from a starting constraint value of 0.25, n_v=80, n=80, and B-Splines of order d=3. Therefore, optimization methods can be implemented with a variety of different starting and ending constraint values, as well as with different smooth functions and iteration amounts.

Approximate Volume Control

In some implementations, calculating 420 the series of target values for a design constraint includes calculating target changes in a volume fraction of the generatively designed three dimensional shape. Some constraints may have a shape gradient or shape derivative that is not well defined, approximate, or not defined at all. In these cases, a proxy shape derivative can be calculated to improve accuracy and exert more control.

Using a constant value it for all optimization iterations would result in converging to a final volume which depends on the relative magnitude it to velocity v, the latter of which depends on the boundary value problem, given by:

−∇·D∈(u)=ƒ in Ω  (71)

u=0 on Γ_D  (72)

D∈(u)=t on Γ_N  (73)

where Ω is the domain for a solid body, D is the fourth order constitutive tensor of the solid body, u is the displacement vector, ƒ is the external load vector, t is the prescribed traction on the Neumann boundary Γ_N with the outward normal n. For simplicity, only homogeneous Dirichlet boundary conditions are assumed on Γ_D. Instead, a variable value μ_t can be used such that a volume target V_T,t is achieved during iteration t of the topology optimization. Let V_t−1 denote the volume after the (t−1)th iteration. The desired volume change during iteration t reads:

ΔV_t=V_T,t−V_t−1.  (74)

This volume change can be approximated by

ΔV_t≈T∫_Γa(v_se+μ)dΓ  (75)

where a∈{0,1} is the advection mask and T is the time step used in solving the Hamilton-Jacobi equation. Note that the approximation error reduces to zero when smaller time steps are used:

$\begin{array}{cc}\underset{T\to 0}{\lim }\Delta V_t-T{\int }_{\Gamma }a\left(v_{se}+\mu \right)d\Gamma =0.& \left(76\right)\end{array}$

The maximum time step is bounded by the Courant-Friedrichs-Lewy (CFL) condition

$\begin{array}{cc}T\le \frac{C\Delta s}{❘v❘\max }& \left(77\right)\end{array}$ $T=\frac{m_{t}C\Delta s}{❘v❘\max }\mathrm{with}0\le m_t\le 1$

where C is a constant, Δ s is the voxel size and |v|_max is the maximum magnitude of advection velocity, given by

v=—(2u·ƒ−D∈(u):∈(u)+μ  (78)

where variables are defined as they are in Equations 71-73.

Let v_u, v_l denote the bounds of v_se such that v_l≤v_se≤v_u. The maximum velocity magnitude is now given by

|v|_max=v_u+μ≥|u_l+μ|  (79)

|v|_max=−(v_l+μ)≥|u_u+μ|  (80)

|v|_max=—(v_u+μ)≥|u_l+μ|  (81)

|v|_max=v_l+μ|u_u+μ|  (82)

Substituting the values for |v|_max and T in Equation 76 yields:

$\begin{array}{cc}\mu \approx \frac{\int {\Gamma }_{a{v}_{se}}d\Gamma -A{v}_u}{A-\int \Gamma ad\Gamma }\mathrm{if}{v}_u+\mu \ge ❘u_l+\mu ❘& \left(83\right)\end{array}$ $\begin{array}{cc}\mu \approx \frac{\int {\Gamma }_{a{v}_{se}}d\Gamma +A{v}_l}{A+\int \Gamma ad\Gamma }\mathrm{if}-\left(v_l+\mu \right)\ge ❘u_u+\mu ❘& \left(84\right)\end{array}$ $\begin{array}{cc}\mu \approx \frac{\int {\Gamma }_{a{v}_{se}}d\Gamma +A{v}_u}{A+\int \Gamma ad\Gamma }\mathrm{if}-\left(v_u+\mu \right)\ge ❘u_l+\mu ❘& \left(85\right)\end{array}$ $\begin{array}{cc}\mu \approx \frac{\int {\Gamma }_{a{v}_{se}}d\Gamma -A{v}_l}{A-\int \Gamma ad\Gamma }\mathrm{if}{v}_l+\mu \ge ❘u_u+\mu ❘& \left(86\right)\end{array}$ $A=\frac{\Delta V_t}{m_{t}C\Delta s}.$

Note that the number of different cases for μ simplifies to Equations 83 and 84 when the body force term is zero in Equation 78 (ƒ=0). This results in the strain energy component of the velocity being positive, i.e., v_u, v_l∈R⁺, resulting in the maximum velocity magnitude being confined to either Equation 83 or Equation 84. The upper bound on m_t in Equation 77 implies that a value cannot be found for an arbitrary high volume change ΔV_t.

FIG. 4C shows a graphical representation of example variations of a value μ with a target volume change Δ V_t given different voxel sizes A s. Curve 401C corresponds to a voxel size of 2, curve 401D corresponds to a voxel size of 3, and curve 401E corresponds to a voxel size of 4. In these example variations, the upper/lower bounds of μ are found using a bisection algorithm. Essentially, the upper/lower bound is initiated at (0+ΔV_t)/2 and successively iterated until a valid μ is found from Equations 83-86 that satisfies the corresponding assumption in Equations 79-82. Alternatively, an advection time that violates the CFL condition may be used while using multiple advection steps, i.e., set T=1/|v|_max and compute μ from

$\begin{array}{cc}\mu =\frac{\frac{\Delta V_t}{T}-\int \Gamma a{v}_{se}d\Gamma }{\int \Gamma ad\Gamma }& \left(87\right)\end{array}$

In this way, the volume derivative can be used as a proxy shape derivative for many optimization constraints, such as stress, fatigue safety factor, buckling safety factor, and displacement. As described below, accurate volume control can be achieved using adaptive controllers, including PID controllers.

Returning to FIG. 4A, the computer-aided design program iteratively modifies a generatively designed three-dimensional shape of the modeled object in the design space in accordance with the one or more design criteria and the one or more in-use load cases for the physical structure. The iteratively modifying includes performing 422 numerical simulation, computing 424 shape change velocities, and updating 426 level-set representation(s) using the shape change velocities. The modification of the three-dimensional shape can include modification of both a geometry of the 3D shape, as well as a topology of the 3D shape. The numerical simulation is performed 422, for example, as described above with reference to FIG. 3A. The performing 422 includes performing numerical simulation of the modeled object in accordance with a current version of the three dimensional shape and the one or more in-use load cases to produce a current numerical assessment of a physical response, e.g., a structural response, of the modeled object.

The computing 424 includes computing shape change velocities for an implicit surface in a level-set representation of the three dimensional shape in accordance with the respective next (e.g., normalized) target value in a respective series of target values calculated 420 for each design constraint. In some implementations, the techniques described can be applied to density-based methods, like SIMP. The target values can be normalized, the techniques for which are described below. The series of target values begins with the initial target value and ends with a final target value for the design constraint. As described above with reference to FIG. 4B and Equation 70, during a first part of the iterations (up to n_v), the target values are calculated 420 from a smooth function, e.g., a B-Spline. When the iteratively modifying passes nm the next target value used in the computing 424 is the final target value for the respective design constraint.

Volume Control Using Adaptive Controllers, Including Proportional-Integral-Derivative (PID) Controllers

Next, accurate volume control using adaptive controllers is discussed.

In general, an adaptive controller is a technique for providing feedback to continuously compute an error value between a desired value and a measured value, and then applying a correction to some control parameter. A Proportional-Integral-Derivative (PID) controller is a type of adaptive controller that continuously computes an error value between a desired value and a measured value, based on proportional, integral, and derivative terms of the value to some control parameter. Different parameters that effect shape changes (e.g., the relative contribution of a particular constraint gradient) can be controlled using PID controllers. During the iteratively modifying, the proportional, integral and derivative components of the PID controller are adjusted to slow or speedup shape changes implicitly by applying varying amount of control on the controlled parameters in response to an oscillation in the generatively designed three dimensional shape.

The PID components of the PID controller are also adjusted to effect a first level of increase in the measured value in response to repetition of success or failure to satisfy the normalized next target value of the measured value. Increases and reductions to the components can be made using a multiplier value that is based on the average deviation of the measured value of the component from the target value. For example, the multiplier can be 1+abs(deviation). As described below, a PID controller can be used to accurately correct for the target volume of the model for each iteration of an optimization process, including a process implementing controlled convergence.

For convenience, the description below will reference volume fractions denoted by V_T,t^ƒ, V_t^ƒ to indicate target volume fraction and actual volume fraction at the end of iteration t, respectively.

$\begin{array}{cc}{V}_{T,t}^f=\frac{V_{T,t}}{V_o}& \left(88\right)\end{array}$ $V_r^f=\frac{V_t}{V_o}.$

where V₀ denotes the volume of the design space.

Given a maximum number of iterations n and a target final volume fraction V_T,n^ƒ, the target volume fraction for each iteration can be computed using a controlled convergence process described above with reference to Equation 68:

$\begin{array}{cc}{V}_{T,t}^f=V_{T,n}^f+\frac{N^d\left(\xi \right)}{N^d\left(0\right)}\left(V_{T,0}^f-V_{T,n}^f\right)& \left(89\right)\end{array}$

Next, the value for it is computed (Equations 83-86) to achieve this volume target for every iteration. Notice that an error between the desired and actual volume fractions V_T,t^ƒ−V_t^ƒ remains at the end of each process. This is caused by the approximation in Equation 75 which leads to an error of magnitude V_T,t−V_t after iteration t. Note that this error is already accounted for in Equation 74 when computing the volume change for iteration t+1 which is the sum of volume change for the next iteration and the error from the previous iteration

$\begin{array}{cc}\Delta V_{t+1}=\left(V_{T,t+1}-V_{T,t}\right)+\left(V_{T,t}-V_t\right)=V_{T,t+1}-V_t.& \left(90\right)\end{array}$

However, in some cases, this is not sufficient to achieve a prescribed volume target. Therefore, PID controllers are used.

FIG. 4D shows examples of graphs 401F, 401G tracking target volume versus actual volume fraction during an iterative optimization process with approximate volume control and without adaptive control. Graph 401F shows discrepancies between a target and actual volume for a modeled object as the volume is increased over a number of iterations. Graph 401G shows discrepancies between a target and actual volume for a modeled object as the volume is reduced over a number of iterations. In both cases, the target volume and the actual volume do not coincide, and the error can grow larger over longer iterations, with potentially disastrous results.

To address this issue, an adaptive, e.g., PID, controller can be used to adjust the volume target for each iteration, to maintain better control over volume changes. Define the error in volume fraction (between actual and target volume) at iteration t as follows:

e _t ^V =V _T,t ^ƒ −V _t−1 ^ƒ.

The target volume for iteration t is now computed by

$\begin{array}{cc}\begin{array}{c}\Delta V_t=V_0\left(K_p{e}_t^V+K_i{\int }_0^t{e}_{\tau }^{V}d\tau +K_d\frac{{\mathrm{de}}_t^V}{\mathrm{dt}}\right)\\ =V_0\left(K_p{e}_t^V+K_i{\int }_0^t{e}_{\tau }^{V}d\tau +K_d\left(e_t^V-e_{t-1}^V\right)\right)\end{array}.& \left(92\right)\end{array}$

where K_p,K_i,K_d∈R⁺ are the PID parameters. Note that the PID controller is applied on target volume fraction change instead of target volume, to ensure the PID parameters are not dependent on initial domain volume. Unless otherwise stated, K_p=1, K_i=0.1, K_d=0.1 are used hereafter. For simplicity, such PID controllers will be denoted as follows hereafter where the subscript t is dropped for clarity:

ΔV_t=∫PID(e_t^V)  (93)

FIG. 4E shows examples of graphs 401H, 401I tracking target volume versus actual volume fraction during an iterative optimization process with PID control (but without adaptive PID control, as described below). Graph 401H shows the positive volume case (in which volume is added to the model with subsequent iterations), while graph 401I shows the negative volume case (in which volume is subtracted from the model with subsequent iterations). In both cases, the target volume and the actual volume track closely, with the error being smaller than the same optimization performed without a PID controller, as described above with reference to FIG. 4D.

The end result of this process is applying control over a response parameter e_t^V, based on variations of a control parameter ΔV_t. In addition, volume control using PID controllers, and adaptive controller techniques in general, can be applied in combination with all of the other systems and techniques described in this document, including for handling arbitrary equality and inequality constraints, described below, and in combination with the various other shape and topology optimization techniques described in connection with FIGS. 5A-8C.

Adaptive PID Tuning

In some implementations, to provide good control for a broad range of design problems, the K_p, K_i, K_d parameters in an adaptive, e.g., a PID controller are modified in response to controller behavior, by applying a multiplier equal to the average deviation of the iteration result from target. Three key controller statuses monitored for can be:

Oscillations: Oscillations are defined as occurring when, during the observed time period, there exist more than 2 pairs of consecutive success and failure.

Repeat failure or success: A success or failure is considered repeated if it occurs over a number of iterations equal to some predetermined threshold, e.g., 10%, of the specified maximum iterations.

Repeat excessive failure or success: A success or failure is deemed excessive if the relative error between it and the target exceeds the predetermined threshold, e.g., 10%.

PID control theory indicates that the integral term should be used to reduce systematic error, while the derivative term should be used to dampen oscillations. The proportional term guides the speed of convergence. In practice, it is often found the derivative term can be hard to adjust and in fact contribute to oscillatory behavior (a phenomenon known as “derivative kick”).

FIG. 4F shows an example of a process of adaptively modifying parameter values for a PID controller. If oscillations are determined 402, then the proportional and integral terms are reduced. The derivative term is set to a percentage of the proportional term if it is currently at zero, otherwise it is also reduced 404. An aim of these changes is to slow volume change. If repeated success or failure is determined 406, then the proportional term is increased, the integral term is reduced, and the derivative term is set to zero 408. An aim of these changes is to slightly increase volume change to remove small systematic error. If repeated excessive failure or success is determined 410, then the proportional term is not affected, the integral term is increased, and the derivative term is set to zero 412. The aim in this situation is to increase volume change proportionally to the significant systematic error.

For the newly adjusted parameters to take effect, a normal status can be applied for a set number of iterations after the adjustment. Similarly, periods accounted for in the determination of controller statuses can be modified to avoid overcompensation. In some implementations, the design problem is allowed to converge rapidly in the initial optimization stages. Hence, during that interval, only the oscillatory status is adapted for and all others are effectively ignored by the adaptive process.

Therefore, returning to FIG. 4A, as part of the computing 424, the proportional, integral and derivative components of an adaptive controller, e.g., a PID controller can also be adjusted to effect a second level of increase in shape changes in response to the repetition of success or failure having a relative error with respect to the next (e.g., normalized) target value that exceeds a threshold. In this example, the second level of increase in shape changes is greater than the first level of increase in shape changes.

An adaptive controller can also be implemented using any appropriate machine learning technique for predicting changes to a control parameter in response to one or more inputs. For example, a controller can be implemented as a neural network having a plurality of neural network layers. The neural network can include an input layer that receives one or more inputs, and an output layer that outputs a change to a control parameter in response to the inputs. The neural network can include one or more hidden layers between the input and output layer, each hidden layer applying one or more non-linear functions to a received input at the layer, where the functions are weighted according to learned parameter values.

An adaptive controller implemented as a neural network can be trained according to any appropriate supervised learning technique. In the case of supervised learning, a controller can be trained on a training data set that includes inputs to a PID controller paired with output changes to the Proportional, Integral, and Derivative changes of the constraint value in response to the input. The training data can include a subset of all input-output pairs generated from a PID controller, e.g., the previous N iterations where Nis some pre-defined value. The training data can be further modified, e.g., with random oscillations or other variations added.

As described below with reference to Arbitrary Equality and Inequality Constraints, a separate adaptive controller can be implemented for each constraint. In some implementations, the adaptive controller can be implemented as a neural network or other machine learning model that receives, as input, the generated constraint change from each PID controller, and generates, as output, a final volume change considering all constraints. In implementations in which a single machine learning model is used, the training data for training the model can include a subset of input-output pairs, as described above, over a number of iterations, with or without random perturbations.

The adaptive controller can be trained to generate constraint values or a final volume change according to additional inputs besides the proportional, integral, and derivative terms. For example, the additional inputs can be from the topology optimization process, such as features extracted from the current shape of the part being designed for, or features of current stress, strain, and displacement results for the part. An adaptive controller can receive these additional inputs to learn parameter values to generate more accurate controlled values or final shape changes and further stabilize the optimization process.

The adaptive controller can be further trained in response to new test cases and design requirements. Additional training can be done offline, in real-time, or a combination of the two. For example, if a specific optimization task does not converge or is not stable, multiple instances of the optimization process can be run for all test cases but for different controller settings, e.g., different model parameter values weighting the machine learning model, different hyperparameter values, e.g., learning rate, batch size, etc., or both. The settings of the instance with the best performance, e.g., according to some performance metric, can be used in an adaptive controller and corresponding optimization process for producing the final design. The multiple instances can be run offline, e.g., when the number of test cases is small, or adapted online during the optimization process.

Although neural networks are given as an example machine learning technique for implementing adaptive controllers, any suitable technique can be used, e.g., fuzzy logic or any appropriate type of regression model. Moreover, adaptive controllers, as described in this document, can also have one or more extra inputs (other than just a measure of error) that provide more information about the trend and current state of the optimization process

Constraint Normalization During Topology Optimization.

Controlled convergence can be combined with normalization of the constraint values, as described above in reference to the computing 424 of FIG. 4A. The constraint values can have different orders of magnitude, resulting in ill-conditioned optimization problems and gradients that do not exist or are very difficult to implement for complex constraints. For example, the starting value of the fatigue safety factor may be about 10,000, while the target value may be 1. In such situations, the sensitivity of the gradient diminishes as the solution nears the target and a final value of 10 (which has an error of 9/10000 with respect to the starting value) may be deemed acceptable due to the relatively small error. A moving reference value, as described below in the Normalization Algorithm, ameliorates the problem.

Normalization Algorithm

Input: ti, t, n, gT,n, Tk, Ti, Kth, Kcur, Kg
Output: gref,t
// initiate or update reference in every iteration t
1:  if t = 0 then
    // initiate reference value
2:  gref,0={❘"\[LeftBracketingBar]"g0-gT,n❘"\[RightBracketingBar]"⁢if⁢❘"\[LeftBracketingBar]"g0-gT,n❘"\[RightBracketingBar]"≠01⁢otherwise
3:  else
    // check if end of inner loop iteration
4:  if ti ≤ Ti and n − t ≤ Ti then
    // compute moving average for some fixed history length Tk
5:  gavg.t=∑k=1Tkgt-kTk
    // define threshold from moving average
6:  Δ gth = Kth |gavg,t − gT,n|
    // define the threshold from current difference
7:  Δ gcur = Kcur|gi − gT,n|
    // define the min/max allowable change in the reference value
8:  Δ gref = Kggref,t−1
    // count the number of oscillations (n0) during the past t i
    iterations
9:  if n0 > 0 then
    // increase reference value
10: gref,t = max{gref,t−1, min{Δ gth, gref,t−1 +Δ gref}}
    // reset inner loop
11: ti = 0
12: else if gref,t−1 >Δ gth then
    // check if converged and decrease reference value
13: gref,t = min {grefmt−1, max{Δ gth,Δ gcur, gref,t−1 −Δ gref}}
    // reset inner loop
14: ti = 0
15: else if Δ gcur > gref,t−1 then
    // check if excessive violation and increase reference value
16: gref,t = max {gref,t−1, min{Δ gth,Δ gcur}}
    // reset inner loop
17. ti = 0
18: else
    // increment inner loop iteration
19: ti ← ti + 1
20: end if
21: else
    // increment inner loop iteration
22: ti ← ti + 1
23: end if
24: end if

Some features of the Normalization Algorithm are noted, as follows.

The reference value g_ref,t should not be updated in every iteration. Sufficient time should be allowed for the constraint to stabilize after setting a new reference. This is achieved using the notion of an inner loop iteration (lines 3-22 of the Normalization Algorithm). The maximum length of inner loop iterations is given by T_i∈ which can be dynamically adjusted for every design problem, for example, adjusted such that there are at least 6 inner loop iterations with a lower bound of 5 on the inner loop length.

The threshold Δg_th used to determine if the existing reference is too far from the current value is computed using a multiplier K_th∈, e.g., 20, of the moving average of the constraint value. The length of the moving average computation is given by T_k∈, e.g., 10. The threshold is increased if oscillations are detected during the last inner loop, otherwise it is reduced when the current reference value is above the threshold.

Any change to the reference value is bounded by a maximum change Δg_th=K_gg_ref,t−1 computed using a ratio 0<K_g≤1 of the existing value, e.g., 0.6. Any reduced reference value should also be larger than a multiplier of the current violation |g_i−g_T,n|. An example of the multiplier is given by K_cur=1.2.

FIG. 4G shows an example of a graph 401J that tracks different measures used in constraint normalization, as detailed in the Normalization Algorithm and as described immediately above. Specifically, the graph 401J shows how the reference value g_ref,t is updated over time as the constraint value approaches the next target value. The final target value is g_T,n, while g_t represents the current value of the constraint as it approaches the final target value g_T,n oscillating over time. The reference value g_ref,t is larger than both the current and final value for the constraint, and does not change with each iteration. The reference value is adjusted over time, approaching the final target value without crossing the current target value.

In some implementations, a value computed from the Normalization Algorithm can be assigned to a target reference value {tilde over (g)}_ref,t and an adaptive controller (e.g., a PID controller as described above with reference to FIGS. 4D-4F) can be used to enable a smooth transition of the applied reference g_ref,t, in situations in which sudden changes of the reference value might otherwise occur at the start of each inner loop iteration of the Normalization Algorithm. One example application is as follows, for a given iteration t:

Δ{tilde over (g)}_ref,t×{tilde over (g)}_ref,t−g_ref,t−1.

g _ref,t ←g _ref,t+ƒ_PID(Δ{tilde over (g)}_ref,t).  (94)

Arbitrary Constraint Handling

Arbitrary Constraint Handling is a general method that can be applied to any type of optimization constraint, even when exact shape derivatives are not available. Whether proxy or actual shape derivatives are implemented, more accuracy and control can be achieved over the generative design process.

Arbitrary Equality Constraints

Consider an equality constraint of type g(Ω, u(Ω))=0 in an optimization problem where compliance minimization is used as the objective function. In the absence of a shape derivative

$\frac{\mathrm{dg}}{d\Omega },$

the conventional practice in industry is to optimize while monitoring g and to terminate optimization when g≈0. This implies that the end result is an un-converged solution with a lack of control as the solution approaches zero as g→0.

However, adaptive controllers can be used to enforce equality constraints, without shape derivatives. Let gⁱ(Ω,u(Ω))j=1, . . . , n_g, denote a sequence of equality constraints that have been normalized as described above (e.g., using the Normalization Algorithm, above) such that g^j(t)∈[0,1]∀t. Consider the normalized error of each constraint, which can be used to define the error in the constraint for each iteration in a manner similar to Equation 91:

$\begin{array}{cc}{e}_t^j=\frac{g_{t-1}^j-g_{T,t}^j}{g_{\mathrm{ref}}^j}.& \left(95\right)\end{array}$

This can be substituted in Equation 93 to approximate the volume change needed to satisfy constraint g_j in iteration t:

ΔV_t^j≈I^jƒ_PID(e_t^j).  (96)

where I^j=−1 when the constraint is inversely related to the volume change and I_gi=1 otherwise. Control is effectively applied over Δ_t^j to achieve a desirable outcome of the response parameter e_t^j. When multiple constraints are present, each constraint would recommend a different volume change needed to satisfy the constraint. The constraints are grouped into positive and negative values:

V _t ⁺ :={x|x=I ^jƒ_PID(e_t^j),x>0,j=1, . . . ,n_g}.

V _t ⁻ :={x|x=I ^jƒ_PID(e_t^j),x>0,j=1, . . . ,n_g}.  (97)

Then, the single volume change applied in this iteration can be computed with

$\begin{array}{cc}\Delta V_t=\{\begin{array}{c}\max \left(V_t^{+}\right)\mathrm{if}{V}_t^{-}=\varnothing \\ \min \left(V_t^{-}\right)\mathrm{if}{V}_t^{+}=\varnothing \\ \frac{\mathrm{avg}\left(V_t^{+}\bigcup V_t^{-}\right)}{2}\mathrm{if}{V}_t^{-}\ne \varnothing \mathrm{and}{V}_t^{+}\ne \varnothing \\ 0\mathrm{otherwise}\end{array}.& \left(98\right)\end{array}$

The μ value in the shape derivative for volume control for any non-zero ΔV_t can be applied following the method described above for approximate volume control, with reference to Equations 71-87.

Arbitrary Inequality Constraints

The techniques described above with reference to arbitrary equality constraints and Equations 95-98 result in the volume change converging to zero as the constraint becomes satisfied, i.e., ΔV_t^j→0 as e_t^j→0. When inequality constraints are used, this causes a potential problem as inequality constraints can be satisfied with ΔV_t^j≠0. Therefore, when inequality constraints are used, a slack variable called an importance factor can be used to manage relative contribution of each violated constraint. An importance factor can mitigate uncontrolled changes in convergence and constraints interfering with the minimization objective when all constraints are given equal importance. The importance factor regulates the relative importance of different constraints during different iterations of the generative design process.

Let p_t¹, . . . , p_t^ng∈ denote the importance factor of each constraint, Equation 96 can be modified as follows:

ΔV_t^j≈p_t^jI^jƒ_PID(e_t^j)  (99)

where the applied importance factor is computed using predetermined, e.g., user-provided, importance factors {tilde over (p)}^j as follows:

$\begin{array}{cc}{p}_t^j=\frac{{\tilde{p}}^j}{{\max }_{k=0}^{n_g}{s}_t^k{\tilde{p}}^k}& \left(100\right)\end{array}$

where s_t^k=1 if the inequality constraint g^k is violated at iteration t, otherwise s_t^k=0. At each iteration, the applied importance factor is updated according to the sign of s_t^k of all constraints. This allows constraint violations to impact the applied importance factor and hence, the volume change ΔV_t.

To prevent sudden changes to i_t^j due to the change of status of constraint from violated to non-violated, and vice versa, a PID controller can be used, e.g., as described above with reference to Equations 96 and 97 and FIGS. 4D-4F, to stabilize the change in importance factor, as follows:

Δp_t^j=p_t^j−p_t−1^j

p _t ^j ←p _t ^j+ƒ_PID(Δp_t^j).  (101)

Complex generative design problems can be solved according to a combination of the techniques described above, e.g., with controlled convergence, a combination of inequality and equality constraints, and using PID controllers.

Modified Augmented Lagrangian Method for Constraint Handling

The augmented Lagrangian algorithm introduced with reference to Equations 28-30 can be modified as follows to accommodate constraint handling for arbitrary equality and inequality constraints, as described above.

Controlled convergence: The classical augmented Lagrangian method computes the constraint violation term e^j as the difference between the current value of the constraint and the final value, e.g., as follows:

$\begin{array}{cc}{e}^j\left(\Omega ,u\left(\Omega \right)\right)=\frac{g^j\left(\Omega ,u\left(\Omega \right)\right)-g_T^j}{g_{\mathrm{ref}}^j}& \left(102\right)\end{array}$

Controlled constraint convergence is applied as described above, e.g., with reference to FIGS. 4B-4C. Some constraints may converge earlier than others, therefore the amount of changes being made to the design can be gradually reduced, i.e., large changes made during the initial phase with progressively smaller changes later.

Constraints without shape derivatives: Although the shape derivatives for many objectives and constraints can be mathematically computed using the adjoint method (Equations 14-25), implementation of the adjoint method in a commercial finite element solver can be a difficult task. Additionally, sometimes optimization constraints can be assessed by black-box evaluators provided by users. Therefore, proxy shape derivatives can be used following the techniques described above, e.g., with reference to Arbitrary Constraint Handling and FIGS. 4B-4E.

Precise Volume Control: Precise control of volume is often quite important for obtaining a good design output from complex engineering examples. The μ calculation methods introduced above with reference to Approximate Volume Control, Volume Control using Adaptive Controllers, and Equations 71-93 can be integrated into the augmented Lagrangian method, while further enhancing the accuracy using a line search algorithm.

First, constraints can be placed into two groups _v, _nv where the former contains all constraints that are affected by volume changes, and the latter contains constraints that are not, e.g., min/max thickness or centroid constraints. The shape derivative from the augmented Lagrangian method is then modified as follows:

$\begin{array}{cc}\frac{dℒ}{d\Omega }=\frac{d}{d\Omega }+{\mu }^{*}\left({\mu }^j+{\lambda }^j{e}^j\right)\frac{{\mathrm{dg}}^j}{d\Omega }+{\sum }_{g^j\in {𝒢}_{\mathrm{nv}}}\left({\mu }^j+{\lambda }^j{e}^j\right)\frac{{\mathrm{dg}}^j}{d\Omega }& \left(103\right)\end{array}$

where the constraint error e^j is computed using a PID stabilized version of Equation 95 with the importance factor term from Arbitrary Inequality Constraints, described above with reference to Equations 99-101:

$\begin{array}{cc}{e}_t^j=p_t^j{f}_{\mathrm{pid}}\left(\frac{g_{t-1}^j-g_{T,t}^j}{g_{\mathrm{ref},t}^j}\right)& \left(104\right)\end{array}$

When the shape derivative dg^j/dΩ is unavailable, it may be approximated using a suitable proxy shape derivative. For example, any such constraint from the set _v may be approximated using the volume shape derivative as follows:

$\begin{array}{cc}\frac{{\mathrm{dg}}^j}{d\Omega }\approx I^j\frac{\mathrm{dV}}{d\Omega }=I^{j}1& \left(105\right)\end{array}$

where I^j=−1 when the constraint is inversely related to the volume change and I_gi=1 otherwise. Note that the derivative of the volume in the normal direction is unity.

Next, μ* is computed using the concepts described above with reference to Approximate Volume Control, Volume Control using Adaptive Controllers, and Equations 71-93. A target volume change Δv_t is computed for each iteration t using Arbitrary Constraint Handling, using all constraints in _v. Next, μ* is computed using the method described above with reference to Approximate Volume Control and Equations 71-87. Note that Equation 87 needs to be modified as follows to account for arbitrary constraints and objectives:

$\begin{array}{cc}{\mu }^{*}=\frac{\frac{\Delta V_t}{T}-{\int }_{\Gamma }a\left[\frac{d𝒥}{d\Omega }+{\sum }_{g^j\in {𝒢}_{\mathrm{nv}}}\left({\mu }^j+{\lambda }^j{e}^j\right)\frac{{\mathrm{dg}}^j}{d\Omega }\right]d\Gamma }{{\int }_{\Gamma }a\left[\left({\mu }^j+{\lambda }^j{e}^j\right)\frac{{\mathrm{dg}}^j}{d\Omega }\right]d\Gamma }.& \left(106\right)\end{array}$

In some cases, a pitfall should be accounted for: μ*→∞ when

$\left({\mu }^j+{\lambda }^j{e}^j\right)\frac{{\mathrm{dg}}^j}{d\Omega }\to \infty .$

When this is detected, μ* should be set to 1 in such cases. This usually happens as constraints converge towards the target value, i.e., e_j→0.

Updating augmented Lagrangian parameters μ^j, λ^j for iteration t can be stabilized using Adaptive, e.g., PID controllers, as described above with respect to FIG. 4F. For example, see below for the update rule:

Δμ_t^j=μ_t^j−μ_t−1^j

μ_t^j←μ_t^j+ƒ_PID(Δμ_t^j).  (107)

As earlier mentioned, a line search algorithm, e.g., gradient descent or Newton's Method, can be applied to ensure the volume change achieved after advecting the level-set is within some acceptable tolerance with respect to the target volume change Δv_t. Line search is done to find the optimal multiplier l_t such that the advection velocity is given by /dΩ. Note that disabling line search is equivalent to setting l_t=1.

FIG. 4H shows examples of graphs 414A-C tracking convergence history of constraints with and without line search. Graph 414A shows applied versus target volume change Δv_t with line search. Graph 414B shows applied versus target volume change Δv_t without line search. Graph 414C shows the history of velocity multiplier l_t during convergence.

Returning to FIG. 4A, in some implementations, the one or more design criteria include multiple design constraints, which include the first inequality constraint and the second inequality constraint. In these implementations, the design constraints are grouped into a first group and a second group, the first group containing all of the multiple design constraints that are effected by volume changes, and the second group containing one or more remaining ones of the multiple design constraints that are not effected by volume changes. Computing 424 the shape change velocities includes using an augmented Lagrangian method, e.g., as described above, that applies an adjustment factor to a sum of shape change contributions from the first group but does not apply the adjustment factor to a shape change contribution from the second group.

In some implementations, the at least one design constraint without a defined shape gradient includes a first inequality constraint and a second inequality constraint. The first inequality constraint has a first input control parameter to a first proportional-integral-derivative controller and has a first importance factor that is multiplied with a shape change amount provided by the first proportional-integral-derivative controller, e.g., as described above with reference to Arbitrary Inequality Constraints and Equations 99-101. The second inequality constraint has a second input control parameter to a second proportional-integral-derivative controller and has a second importance factor that is multiplied with a shape change amount provided by the second proportional-integral-derivative controller.

Computing 424 the shape change velocities for the implicit surface in accordance with the proxy shape gradient can therefore include adjusting both the first and second importance factors based on whether or not one or more other constraints were violated in a prior iteration of the iteratively modifying. For example, the importance factor can be modified by multiplying the factor with a violation multiplier with respect to the other constraints in prior iterations of the iteratively modifying. In some implementations and as described above, e.g., generally with reference to FIGS. 4D-4F and with reference to Equations 99-101, proportional-integral-derivative control can be used to stabilize adjustments of the first and second importance factors.

Also as described above with reference to FIGS. 4A, 4F, 4G, the oscillation can be more than two pairs of consecutive successes or failures to satisfy the normalized next target value. Further, the repetition of success or failure to satisfy the normalized next target value can be repetition that occurs over a number of iterations equal to or greater than ten percent of the iteration amount and where the threshold is ten percent. In some implementations, the adjusting in response to the oscillation includes reducing the proportional component, reducing the integral component, and reducing or resetting the derivative component. In some implementations, the adjusting in response to the repetition of success or failure includes increasing the proportional component, reducing the integral component, and setting the derivative component to zero. In some implementations, the adjusting in response to the relative error exceeding the threshold includes increasing the integral component and setting the derivative component to zero.

As described above with reference to Adaptive Controllers and FIGS. 4D-4F, the at least one design constraint may not have a defined shape gradient. Therefore, in some implementations, the computing 424 includes computing the shape change velocities for the implicit surface in accordance with a proxy shape gradient output that is adjusted by adaptive control. Specifically, an input control parameter is used that is a measure of error between a normalized current value for the at least one design constraint without the defined shape gradient and a normalized next target value from the respective ones of the target values in the series of target values. The measure of error changes with a reference value that is changed at least once during the iteratively modifying, which can allow for controlled convergence while meeting the critical constraint smoothly.

In some implementations, the at least one design constraint includes a first equality constraint and a second equality constraint, the first equality constraint having a first input control parameter to a first adaptive controller, and the second equality constraint having a second input control parameter to a second adaptive controller, e.g., as described above with reference to Arbitrary Equality Constraints and Equations 95-98. Then, the computing 424 the shape change velocities for the implicit surface in accordance with the proxy shape gradient includes using a maximum shape change amount provided by the first and second adaptive controllers when none of the first and second equality constraints are inversely proportional to shape change and using a minimum shape change amount provided by the first and second adaptive controllers when none of the first and second equality constraints are proportional to shape change. When one constraint is inversely proportional to shape change and at least one constraint is proportional to shape change, an average shape change amount is used.

After the computing 424, the level-set representation(s) are updated 426 using the shape change velocities to produce an updated version of the three dimensional shape of the modeled object. The performing 422, the computing 424, and the updating 426 are repeated until a check 428 determines that a predefined number of shape modification iterations have been performed or that the generatively designed three dimensional shape of the modeled object in the design space converges to a stable solution for the one or more design criteria and the one or more in-use load cases. The updating 426 and the check 428 can be performed, e.g., as described above with reference to FIG. 3A. Then, the computer-aided design program can provide the generatively designed three dimensional shape of the modeled object for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

Fatigue Constraint

Next, a description of optimizing a design based on provided fatigue constraints on the designed body is provided. In addition to the techniques described above, e.g., FIGS. 3A-4H and accompanying description for Controlled Convergence, Seeding, and Arbitrary Constraint Handling, a designed body can be optimized to satisfy fatigue constraints specifying a total life expectancy of the body, or a number of cycles per load case imposed on the body. The process can be performed automatically in addition to an optimization process, e.g., a generative design process as described above.

Two main methods of evaluating this constraint are described: safe life calculation and the damage tolerant approach. The former is about preventing fatigue damage by keeping stresses below an allowable threshold. The latter accepts that there is fatigue damage and aims to ensure a fatigue crack does not cause catastrophic failure before a specific inspection point in time.

Safe Life Fatigue

FIG. 5A shows an example of a process that generates one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes that solve for safe life fatigue constraint(s) on the body. A computer-aided design program obtains 504 a design space for a modeled object, for which a corresponding physical structure will be manufactured, one or more design criteria for the modeled object, one or more in-use load cases for the physical structure, and one or more specifications of material from which the physical structure will be manufactured. The obtaining 504 can be done, for example, as described above with reference to FIG. 3A.

Safe life optimization is based on a material's fatigue properties. Samples are tested to obtain SN-curves available in material databases. The designer can specify the expected amount of loading cycles over the lifetime of the part for each relevant load case. Their cumulative effect guides the allowable design stress. Usually, stress-based safe life fatigue is used. However, in low cycle loading problems, strain life fatigue may be preferred.

Treatment of SN-Curves.

FIG. 5B shows an example of an SN-curve 520 tracking fatigue strength (stress) on a body over a number of cycles. For bodies that include multiple materials, multiple SN-curves can be used, as described in more detail below.

Both damage tolerant and safe life fatigue approaches described herein can make use of material SN-curves. Due to their potential flat and infinite regions, these pose computational challenges and should be interpreted with caution during the optimization process.

SN-curves normally place stress on the y-axis and cycles on the x-axis, as the typical fatigue workflow involves going from a known number of cycles limit to finding the appropriate maximum stress constraint. This makes it natural for SN-curve user-input to be [cycles, stress] points. Such point definition fits the damage tolerant optimization method well and thus only a sorting of points from lowest to highest cycles need be used. However, in the case of safe-life fatigue, the optimization approach goes from stress to cycle data. Simulation of the object being designed provides stress data, for which the predicted number of cycles should be found from the SN-curve. Thus, an important aspect of the method is to invert the user-provided data points into [stress, cycles] pairs and sort those from lowest to highest stresses. Additionally, when an endurance section is detected, i.e., there exists a flat region of the curve at an endurance stress (below which the number of cycles is effectively infinite), the set of points corresponding to that section of the curve is re-sorted from highest stress to lowest stress. The curve reading algorithm then readily returns the highest reported cycle value instead of interpolating to infinity.

When the input data includes data for different materials, a separate number of expected loading cycles is provided for each of the different materials. This includes returning a number of loading cycles from one or more curves fit to the set of data points in a plastic region and an elastic region of the data relating fatigue strength to loading cycles; and then, returning a highest load cycles number from the set of data points in an endurance region of the data relating fatigue strength to loading cycles.

Stress-Based Safe Life Fatigue

Returning to FIG. 5A, the one or more design criteria include a required number of loading cycles for the modeled object for each of the one or more in-use load cases for the physical structure. The input can be provided as cycles for each load case, or as a period of time, e.g., in seconds, minutes, hours, days, weeks, months, years, etc. that can be optionally converted into cycles. The one or more specifications include data relating fatigue strength to loading cycles. The data can be provided as SN-curves for the material. The data can measure stress as an indicator of fatigue stress, but in some implementations strain is used, e.g., when the in-use load cases have a low number of cycles as a requirement.

The program iteratively modifies a generatively designed three dimensional shape of the modeled object in the design space in accordance with the one or more design criteria, the one or more in-use load cases for the physical structure, and the one or more specifications. The iteratively modifying can include modifying both a geometry and a topology for the three dimensional shape of the object.

The iteratively modifying includes performing 506 numerical simulation of the modeled object in accordance with a current version of the three dimensional shape and the one or more in-use load cases to produce a current numerical assessment of a physical response of the modeled object, e.g., as described above with reference to FIGS. 3A & 4A.

A maximized stress or strain element is found 508 for each of the one or more in-use load cases for the physical structure, from the current numerical assessment of the physical response of the modeled object. The element can be a value at a point, location or region of the physical structure. In some implementations and as described below with reference to FIGS. 8A-8C, the finding 508 includes calculating a maximum stress value for an in-use load case based on at least a standard deviation of a stress distribution in the current numerical assessment of the physical response of the modeled object.

Specifically, during topology optimization, each element's stress can be computed at every iteration. After optionally adjusting for singularities (described below in reference to FIGS. 8A-8C, the maximum principal stress across all elements for a given load case can be found and the corresponding number of cycles obtained from the SN curve(s). The ratio of cycles obtained to cycles required represents the damage fraction. Using Miner's rule, the overall safety factor is computed from the damage fractions of all load cases. The fatigue constraint has to work with multi-material problems. Each material's fatigue safety factor is calculated separately and the minimum guides optimization, as shown below:

$\begin{array}{cc}{S}^f=\begin{array}{c}{n}_m\\ \min \\ i\end{array}{S}_{m_i}^f.& \left(108\right)\end{array}$

where S_mi^ƒ is the fatigue safety factor of a material computed using:

$\begin{array}{cc}{S}_{m_i}^f={\left[{\sum }_{\mathrm{lc}=1}^n\frac{C_{\sigma }}{C_{\mathrm{lc}}^{\mathrm{req}}}\right]}^{-1}& \left(109\right)\end{array}$

where the fatigue safety factor S^ƒ is given by the minimum fatigue safety factor of all materials used, C_σ is the number of cycles obtained from the SN-curve for the max principal stress (max_Ωσ^p), and C_lc^reg is the required number of cycles for the given load case lc.

During optimization, the fatigue constraint is defined as S^ƒ−S_T≥0 where S_T^ƒ is the target fatigue safety factor. This can be enforced using methods described above with reference to Controlled Convergence, Arbitrary Constraint Handling, and FIGS. 4A-4H with the shape derivative dS^ƒ/dΩ derived below. In some implementations, the volume shape derivative can be used as the proxy shape derivative, as shown in Equation 105.

Shape Derivative of Aggregated Fatigue Metric

Let Ω be a domain in . In this specification, a stress-based fatigue metric for Ω with respect to a number of load cases that are applied repetitively over time is assessed at each point x∈Ω independently. For purposes of topology optimization, these pointwise measurements are aggregated into a global measurement called an aggregated fatigue metric that approximates the maximum value of the pointwise measurements.

The aggregated fatigue metric takes the following form. Let C:→ be the inverse of the SN curve of the material including Ω and let C_l∈ be a reference value for the l^th load case. This aggregated fatigue metric enables the summation of damage caused by cyclic loading. A target value of 1/(safety factor) is optimized for, where the damage the part can withstand and the actual damage from the load (stress or strain and number of cycles) are summed.

The aggregated fatigue measurement of Ω is defined as:

$\begin{array}{cc}S\left(\Omega \right):={\left({\sum }_{l=1}^L\frac{C\left({\tilde{\sigma }}_l\left(\Omega \right)\right)}{C_l}\right)}^{-1}& \left(110\right)\end{array}$

where σ_l(Ω) is an approximation of the maximum of the magnitude of the stress tensor in Ω under the l^th load case. In other words,

$\begin{array}{cc}\widetilde{{\sigma }_l}\left(\Omega \right):={\left({\int }_{\Omega }{\sigma \left(u_l\right)}^{2p}\right)}^{1/2p}& \left(111\right)\end{array}$

where u_l:Ω→ is the displacement function that satisfies the equations of linear elasticity in Ω with respect to the l^th load case, and σ(u_l) is the associated stress tensor.

The shape derivative of the aggregated fatigue metric S is described in more detail, below.

Proposition 1. The shape derivative of the shape function S with respect to the variation of a shape Ω generated by the normal speed Θ:∂Ω→ is given by:

$\begin{array}{cc}{\mathrm{DS}}_{\Omega }·\Theta =& \left(112\right)\end{array}$ $-{\left[S\left(\Omega \right)\right]}^2{\sum }_{l=1}^L{\frac{C^{\prime }\left({\tilde{\sigma }}_l\left(\Omega \right)\right)}{C_l}\left[{\tilde{\sigma }}_l\left(\Omega \right)\right]}^{1-2p}\left({\int }_{\partial \Omega }\left[{\sigma \left(u_l\right)}^{2p}+\sigma \left(u_l\right):e\left({\lambda }_l\right)\right]\Theta \right)$

where λ_t:Ω→ is the adjoint displacement function of the l^th load case. This is the solution of the linear elastic equations with adjoint force term given in weak form as the linear integral form:

_adj(Ω,δu):=∫_Ω∥σ(u_l)∥^2p−1σ(u_l):σ(δu)  (113)

Since S can be decomposed into a sequence of operations applied to the stress integral ∫_Ω∥σ(u_l)∥^2p this computation breaks down as follows. First, the derivative of this sequence of operations is computed using the chain rule from ordinary calculus. Then, the shape derivative of the stress integral is computed, for which Céa's method is used and described below.

Let Ω_ε denote the variation of Ω generated by the normal speed function Θ. The computation of the desired shape derivative, denoted as DS_Ω·Θ, begins as:

$\begin{array}{cc}{D{S}_{\Omega }·\Theta =\frac{d}{d\epsilon }S\left({\Omega }_{\epsilon }\right){❘}_{\epsilon =0}=\frac{d}{d\epsilon }{\left({\sum }_{l=1}^L\frac{c\left({\tilde{\sigma }}_l\left({\Omega }_{\epsilon }\right)\right)}{c_l}\right)}^{-1}❘}_{\epsilon =0}& \left(114\right)\end{array}$

First, the chain rule from ordinary calculus is used to bring the derivative d/dε onto the stress integral, where then the techniques of shape differentiation take over, as shown below:

${D{S}_{\Omega }·\Theta =-{\left(\sum _{l=1}^L\frac{C\left({\tilde{\sigma }}_l\left({\Omega }_{\epsilon }\right)\right)}{C_l}\right)}^{-2}\sum _{l=1}^L\frac{C^{\prime }\left({\tilde{\sigma }}_l\left(\Omega \right)\right)}{C_l}\frac{d}{d\epsilon }{\tilde{\sigma }}_l\left({\Omega }_{\epsilon }\right){❘}_{l=0}=-{\left[S\left(\Omega \right)\right]}^2{\sum }_{l=1}^L\frac{c^{\prime }\left({\tilde{\sigma }}_l\left(\Omega \right)\right)}{c_l}\frac{{\left[{\tilde{\sigma }}_l\left(\Omega \right)\right]}^{1-2p}}{2p}\frac{d}{d\epsilon }{\int }_{\Omega \epsilon }{\sigma \left(u_{l,\epsilon }\right)}^{2p}❘}_{l=0}$

where u_l,ε is the displacement function for the l^th load case in the domain Ω_ε.

To continue, shape differentiation is used. Céa's method is used to compute the shape derivative of the stress integral remaining in Equation 115. The subscripted l is suppressed in the description below because it is the same for all load cases. The Lagrangian follows:

(Ω,u,λ):=∫Ω∥σ(u)∥^2p+A(Ω,u,λ)+(Ω,λ).  (116)

where

A(Ω,u,λ):=∫_Ωσ(u):e(λ)  (117)

is the integral form appearing in the weak form of the elasticity equations in Ω.

Here, σ(u) is the stress tensor of the displacement u and e(λ) is the “virtual strain” of the “virtual displacement” λ. Thus, (Ω, u, λ) is the “virtual work” associated to the pair of displacements. Also, (Ω, λ) is the integral linear form appearing in the weak form of the elasticity equations in Ω which encodes the “virtual work” done by the applied body and boundary tractions. Now Céa's method proceeds in three steps.

Step 1. Setting the variation of the Lagrangian with respect to λ equal to zero yields the linear elasticity equations satisfied by u. This is by design, since (Ω, u, λ)+(Ω, λ) is precisely the left-hand-side of the weak elasticity equations.

Step 2. Setting the variation of the Lagrangian with respect to u yields a related system of equations for λ called the adjoint state equations. The solution of these equations is called the adjoint state. The weak version of these equations is:

0=∫_Ω2p∥σ(u)∥^2p−1σ(u):σ′(δn)+A(Ω,u,λ)  (118)

for all variations δu of u. Here, σ′(δn) is the derivative of the stress tensor with respect to the displacements. This has a simple form, because stress is linearly related to strain which is in turn linearly related to u—the gradient is a linear operator, or quite simply, σ′=σ. The symmetry of is also used. The above shows that λ satisfied the linear elastic equations but with a new “adjoint force” term encoded in the integral linear form:

_adj(Ω,δu):=∫_Ω2p∥σ(u)∥^2p−1σ(u):σ′(δu)  (119)

Step 3. Céa's method hinges on the following fact, which can be proven: the shape derivative of the Lagrangian, computed by ignoring the shape dependence of u and λ and then plugging in the state for u and the adjoint state for λ, is equal to the shape derivative of ∫_Ω∥σ(u)∥^2p itself. Since ail expressions in the Lagrangian are volumetric integrals, only the formula for the shape derivative of such an integral of a shape-independent function over Ω is needed.

Next, assume that Θ vanishes on the Dirichlet and in the homogeneous Neumann boundary of Ω (otherwise the shape derivative would include other terms). Finally, assume that body forces are absent, for instance by ignoring the effects of gravity in the given load case. Together, these last two assumptions have the effect of dropping the term from the formula for the shape derivative (Equation 116). Therefore, the following results after all assumptions:

$\begin{array}{cc}\frac{d}{d\epsilon }{\int }_{{\Omega }_{\epsilon }}{\sigma \left(u_{\epsilon }\right)}^{2p}{❘}_{\epsilon =0}={\int }_{\partial \Omega }\left({\sigma \left(u\right)}^{2p}+\sigma \left(u\right):e\left(\lambda \right)\right)\Theta & \left(120\right)\end{array}$

where u is the state and λ is the adjoint state.

Strain-Based Safe Life Fatigue

The strain-based safe life calculation also uses Equations 108 and 109, however the C_σ0 parameter becomes C_ε. A strain versus cycles form of the SN-curve is used as input. The strain value is computed from the max principal stress (max_Ωσ^p) according to Neuber's rule, below:

$\begin{array}{cc}\epsilon =\frac{K^2*\sigma }{2E}.& \left(121\right)\end{array}$

where K is the concentration factor specified by a user (defaults to 1) and E is the material's Young's Modulus.

Returning to FIG. 5A, as described above with reference to the Treatment of SN-Curves and FIG. 5B, an expected number of loading cycles for each of the one or more in-use load cases for the physical structure is determined 510 using the maximized stress or strain element and the data relating fatigue strength to loading cycles. A fatigue safety factor inequality constraint for the modeled object is redefined 512 based on a damage fraction calculated from the required number of loading cycles for the modeled object and the expected number of loading cycles for each of the one or more in-use load cases for the physical structure. In some implementations, the inequality constraint can be normalized and adjusted using adaptive controllers, e.g., as described above with reference to FIGS. 4D-4F.

In some implementations, the one or more specifications include two or more specifications of respective different materials from which the physical structure will be manufactured. Accordingly, the data for the specifications include data relating fatigue strength to loading cycles for each of the different materials. As such, the determining 510 includes determining a separate number of expected loading cycles for each of the different materials. Further, the redefining 512 includes calculating a separate fatigue safety factor for each of the different materials based on respective damage fractions calculated from respective ones of the numbers of expected loading cycles for the different materials. From the separate safety factors, a minimum value of the fatigue safety factors for the different materials is used to redefine the fatigue safety factor inequality constraint for the modeled object.

In some implementations, the one or more in-use load cases for the physical structure includes two or more in-use load cases for the physical structure and the one or more design criteria include a required number of loading cycles for the modeled object for each of the two or more in-use load cases for the physical structure. Calculating the separate safety factor for each respective one of the different materials includes summing load-specific damage fractions corresponding to the two or more in-use load cases, wherein each load-specific damage fraction includes a number of expected loading cycles, for one of the different materials and one of the in-use load cases, divided by the required number of loading cycles for the one of the in-use load cases. The separate safety factor is obtained by inverting the sum of the load-specific damage fractions to obtain the separate safety factor. See e.g., Equation 109.

The shape change velocities for an implicit surface in a level-set representation of the three dimensional shape are computed 514 in accordance with at least the fatigue safety factor inequality constraint. The shape change velocities can be computed with respect to other constraints, and proxy derivatives can be computed where a shape derivative is unavailable or not well defined, e.g., as described above with reference to Arbitrary Constraint Handling and FIGS. 4A-4H. Then, the level-set representation is updated 516 using the shape change velocities to produce an updated version of the three dimensional shape of the modeled object. The shape derivative for the Aggregated Fatigue Metric described above with respect to Equations 110-120 can be used to compute the shape change velocity for the fatigue safety factor inequality constraint.

In some implementations, the computing 514 includes computing at least one shape change velocity using an amount determined from a shape derivative formula that approximates a shape derivative of the fatigue safety factor, e.g., the volume shape derivative described above with reference to Equation 105, The formulation is modified to include the stress shape derivative instead of the volume shape derivative.

In some implementations, an importance factor is used for handling the at least one design constraint, including the fatigue safety factor, which is adjusted based on whether or not one or more other constraints were violated in a prior iteration of the iteratively modifying, e.g., as described above with respect to FIG. 4A, and generally with respect to Arbitrary Inequality Constraint Handling and Equations 53-60.

In some implementations, a target value of the volume fraction is adjusted using adaptive controllers, e.g., as described above with reference to FIGS. 4D-4F. Specifically, a target value of the volume fraction or a minimum thickness based inequality constraint is adjusted between an initial target value and a final target value across multiple iterations of the iteratively modifying. The adaptive control is used to stabilize changes made in the amount determined from the shape derivative formula as the target value is adjusted across the multiple iterations.

In some implementations, adaptive control is used to adjust a total contribution of the amount determined from the shape derivative formula to the shape change velocities used in the updating. Adaptive controllers can be used to adjust the contribution of the proxy shape derivative to the total advection velocity, i.e., according to a modified version of the Augmented Lagrangian method for constraint handling, as described above with reference to Equations 102-107.

The performing 506, the finding 508, the determining 510, the redefining 512, the computing 514 and the updating 516 repeat until a check 518 determines that a predefined number of shape modification iterations have been performed or that the generatively designed three dimensional shape of the modeled object in the design space has converged to a stable solution for the one or more design criteria and the one or more in-use load cases.

The check 518 can be a check as described above with reference to FIGS. 3A & 4A. Then, the generatively designed three dimensional shape of the modeled object can be provided, e.g., for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems. The techniques can be combined with hybrid, hollow, and hybrid-hollow methods for topological optimization.

Damage Tolerant Fatigue

Safe Life Fatigue constraints as described above help to ensure that an object does not fail under user-specified load cases by preventing fatigue damage during a specified number of loading cycles. In contrast, damage tolerant fatigue techniques focus on limiting fatigue damage to be within an allowable limit to help ensure the fatigue damage will be detected in the field before the part fails. The key objective is to compute the critical fatigue crack length and ensure that it does not exceed the thickness of the design component during the service inspection interval. Hence, a key requirement of this method is the ability to enforce thickness constraints on topology optimized designs. The thickness constraint can be enforced implicitly, for example by increasing the volume of the design which results in all parts of the design increasing in thickness.

It is noted that, like the Safe Life Fatigue approach, the Damage Tolerant Fatigue techniques described below can be used in any combination with the above-described techniques, e.g., Controlled Convergence, Seeding, and Arbitrary Constraint Handling. Further, in some implementations, the Safe Life Fatigue and Damage Tolerant Fatigue approaches are both available where fatigue design constraints for a given problem are solved with the Safe Life approach in some cases, while in other cases, fatigue design constraints are solved using the following Damage Tolerant Fatigue techniques. An appropriate approach can be selected based on the intended use of a part being designed. For example, in aerospace applications, a part is commonly designed according to a damage tolerant approach. The availability of the two modes together allow for addressing a greater range of design problems.

FIG. 6A shows an example of a process of iteratively modifying a generatively designed three dimensional shape of a modeled object in the design space in accordance with the one or more design criteria, including at least one damage tolerant fatigue constraint. The three dimensional shape includes a level-set representation of an implicit surface and the one or more design criteria include a required number of loading cycles for the modeled object for each of the one or more in-use load cases for the physical structure.

Numerical simulation of the modeled object is performed 608 in accordance with a current version of the three dimensional shape and the one or more in-use load cases to produce a current numerical assessment of a physical response, e.g., a structural response, of the modeled object, as described above with reference to FIGS. 3A, 4A, and 5A.

An expected number of loading cycles for each of the one or more in-use load cases for the physical structure is determined 610 using the current numerical assessment and the measure of thickness to enforce the design criterion that limits the minimum thickness. The point with the maximum strain or stress from the current numerical assessment can be used to determine the expected number of loading cycles. Computing a current thickness and enforcing the minimum thickness is described below.

Thickness Constraint

First, a thickness constraint is computed and enforced. One or more of multiple methods can be used to measure and enforce thickness, e.g., a combination of at least two distinct thickness measures, as described below. Also, it is to be appreciated that the described Damage Tolerant Fatigue techniques can be implemented with density based topology optimization, e.g., using the SIMP method, rather than with a boundary based topology optimization, e.g., the level-set method, as described in connection with FIG. 6A.

FIG. 6B shows graphical representations of geometries 600A, 600B in which thickness of the geometries 600A, 600B are computed according to different measuring techniques. Geometry 600A shows a body 601 measured using a ray casting method. Ray-cast thickness at point x in geometry 600A is defined as h^r(x)=|x−x^r| where x_r is the point where a ray cast in the negative direction of a vector n at point x hits the geometry 600A. Ray casting can be repeated for different directions of the vector n at point x, including when the vector n is normal to the geometry 600A.

By contrast, sphere-fitting method is shown in geometry 600B and involves finding the diameter of the largest sphere that can fit inside the domain while touching the measuring point x in the geometry 600B. This is done by running a bisection algorithm for locating the sphere centroid x_s between points x and x_r. Discrete sampling locations {x₁^p, x₂^p, . . . x_i^p} are defined on the sphere perimeter/surface using the polar angle a. The sphere with origin x_s is deemed to fit inside the domain when the maximum level-set with respect to the domain is below a particular threshold Δh, i.e., the following condition is satisfied:

maxψ(x_i^p≤Δh  (122)

The sphere-fitting thickness is then defined as h^s(x)=2|x−x_s|. In some implementations, thickness is defined as h(x=max{h^s(x), h^r(x)}. If the thickness at a point not on the domain surface is required, it is first projected on to the surface using the normal of the level-set before thickness on the surface is measured. The sphere-fitting method can be repeated multiple times for multiple directions (resulting in multiple x_r points, depending on the selected direction). The largest sphere computed can be used to define the thickness.

A minimum thickness constraint g^t can be defined as

$\begin{array}{cc}{g}^t\left(\Omega \right)-h_{\min }\ge 0.& \left(123\right)\end{array}$ $\begin{array}{c}\min \\ x\in \Gamma \end{array}h\left(x\right)-h_{\min }\ge 0$

where Γ is the domain surface and h_min is the target minimum thickness. An approximate shape derivative is given by:

$\begin{array}{cc}\frac{d{g}^t}{d\Omega }\left(x\right)\approx \{\begin{array}{c}h\left(x\right)-h_{\min }\mathrm{if}h\left(x\right)<h_{\min }\\ 0\mathrm{otherwise}.\end{array}& \left(124\right)\end{array}$

The thickness constraint can be applied using methods described above with respect to Arbitrary Constraint Handling and FIGS. 4A-4H. In some implementations, the importance factor for the thickness constraint is modified by a sigmoid function S^r(ξ) such that the thickness constraint only comes into effect during the latter stages of optimization

p _t ^j ←S ^r(ξ)p_t^j  (125)

The sigmoid function is defined as:

$\begin{array}{cc}{S}^r\left(\xi \right)=\frac{1}{1+e^{-\left(12\xi -6\right)}}& \left(126\right)\end{array}$

with ξ computed as follows

$\begin{array}{cc}\xi \left(t\right)=\frac{t-\hat{t}}{n-\hat{t}}& \left(127\right)\end{array}$

where {circumflex over (t)} identifies the start of thickness constraint application which is taken as the maximum of volume reduction iterations n_v and the first violation of thickness constraint.

FIG. 6C shows an example of a graph tracking curves 601A-E of sigmoid functions at different rates. Curve 601A is curve of a sigmoid function with rate of 1; curve 601B is a curve of a sigmoid function with a rate of 2; curve 601C is a curve of a sigmoid function with a rate of 3; curve 601D is a curve of a sigmoid function with a rate of 4; and curve 601E is a curve of a sigmoid function with a rate of 5.

Although two techniques for measuring thickness are described here, it is understood that any suitable technique for measuring thickness of a body can be used, without loss of generality. Further, the techniques described can be combined, i.e., a combination of two distinct thickness measurements, to improve the results. For example, the two distinct thickness measures can include (i) a first distance measure being a length within the modeled object of a ray cast in a negative normal direction from a surface point of the modeled object, and (ii) a second distance measure being a diameter of a largest sphere that touches the surface point of the modeled object and fits inside the modeled object as determined by checking discrete sampling locations defined on the sphere's surface, such as described above with reference to FIG. 6B.

Damage Tolerant Fatigue Constraint

Next, the damage tolerant fatigue constraint is described. The material's critical fatigue crack length h^d is defined as

$\begin{array}{cc}{h}^d={\left(\frac{1}{\frac{1}{\left(\frac{{\rho }_g}{a^2}-1\right)}-{\mathrm{CI}\left(Y\Delta \sigma \sqrt{\pi }\right)}^{{\rho }_g}\left(\frac{{\rho }_g}{2}-1\right)}\right)}^{\left(\frac{1}{{\rho }_g/2-1}\right)}& \left(128\right)\end{array}$

where C is the number of fatigue cycles, ρ_g is the modulus of the material's fatigue crack growth curve and can be user-specified, e.g., 2.8, a is the initial crack length, which can be set as equal to the minimum defect detection size of the chosen part inspection method, with a default value of a=1 mm, I is a stress intensity factor, e.g., with a default value of I=1.22×10⁻¹²), and Y describes the crack type, e.g., with a default value Y=1. The stress range Δσ is derived from the stress ratio R, e.g., R=0.1 and maximum stress σ_max as the following:

Δσ=(1−R)σ_max.  (129)

The minimum thickness of the material is then given by

$\begin{array}{cc}{h}_m^d=\frac{h^d}{R_T^d}{S}^f& \left(130\right)\end{array}$

where R_T^d is the desired critical thickness to width ratio and S^ƒ is the prescribed design safety factor. The R_T^d multiplier and Y crack type variables are related and can be found in engineering tables describing fatigue crack properties for a given material.

FIG. 6D shows an example of a table 601F describing fatigue crack properties. In the table 601F, the crack type variables Y are shown in relation to corresponding R_T^d multipliers. Y is the geometry correction factor used in the stress intensity factor expression. The correction factor describes the link between the crack length a and the thickness of the feature W. The table 601F shows values of Y and a/W for different materials. The table 601F can be user-specified, with specified values for Y and the ratio a/W for different materials.

The crack length a is set equal to a critical crack length ac, which tends towards infinity, i.e., catastrophic failure of the part. During optimization, this gives a minimum thickness to optimize for, i.e., if the user desires an a/W ratio of 0.1 (i.e., the minimum thickness should be 10 times ac) then the ac value is calculated and a minimum thickness is applied as a multiple of this value. Generally, the value of Y is small and therefore represents a small crack thickness ratio, i.e., a large thickness of the part versus the critical crack length in the feature of the part. This has the effect of increasing the critical crack length in the feature so that it is easily detectable by the chosen inspection techniques.

One way to add an optimization constraint to satisfy damage tolerant fatigue requirements is to use a constant global thickness target. Assume the design is at the fatigue limit where C is equal to the number of fatigue cycles required and σ_max is the corresponding stress from an SN-curve. Equation 128 is used to compute critical crack length h_mi^d for each material. The min thickness constraint (Equation 123) is modified with the target thickness set to the result of Equation 130 using the worst case crack length. The target thickness remains constant throughout the optimization:

$\begin{array}{cc}{h}_{\min }=\underset{i}{\stackrel{n_m}{\min }}{h}_{m_i}^d& \left(131\right)\end{array}$ $\begin{array}{c}\min \\ x\in \Gamma \end{array}h\left(x\right)-\underset{i}{\stackrel{n_m}{\min }}{h}_{m_i}^d\ge 0.$

Alternatively, the damage tolerant fatigue safety factor can be computed for each point x on the surface of the design at each iteration. This uses the stress σ_lc(x) and thickness h(x) at each point to compute the supported number of cycles for each load case C_lc(x):

$\begin{array}{cc}{C}_{\mathrm{lc}}\left(x\right)=\left[\frac{1}{\left(a^{p_g/2-1}\right)}-\frac{1}{{\left(h\left(x\right)R_T^d\right)}^{p_g/2-1}}\right]\left[\text{⁠}\frac{1}{{I\left(Y\left(1-R\right){\sigma }_{\mathrm{lc}}\left(x\right)\sqrt{\pi }\right)}^{p_g}\left(\frac{p_g}{2}-1\right)}\right]& \left(132\right)\end{array}$

The damage tolerant fatigue safety factor of a material S_d can now be worked out similarly to safe life fatigue method, using Miner's rule:

$\begin{array}{cc}{S}^d={\begin{array}{c}\min \\ x\in \Gamma \end{array}\left[{\sum }_{\mathrm{lc}=1}^n\frac{c_{\mathrm{lc}}\left(x\right)}{c_{\mathrm{lc}}^{\mathrm{req}}}\right]}^{-1}& \left(133\right)\end{array}$

where C_lc(x) is the number of cycles obtained from Equation 132 and c_lc^req is the required number of cycles for the given load case lc.

During optimization, the fatigue constraint is defined as S^d−S_T≥0 where S_T^d is the target fatigue safety factor−the lowest of all material safety factors. This can be enforced using methods described above with respect to Arbitrary Constraint Handling and FIGS. 4A-4H. The volume or thickness shape derivative (Equation 124) can be used as a proxy shape derivative, as described above. In order to use the thickness shape derivative, an extra step should be taken of converting the point stress and thickness data corresponding to the critical safety factor 4 into a thickness target h_min using Equations 128 and 130. Note the importance factor should be multiplied by sigmoid type functions, as described above with reference to FIG. 6C.

FIG. 6E shows an example of a process of generating one or more portions of a 3D model of an object to be manufactured, using one or more generative design processes that solve for damage tolerant fatigue constraints on the body.

A computer aided design program obtains 602 a design space for a modeled object for which a corresponding physical structure will be manufactured, one or more design criteria for the modeled object, one or more in-use load cases for the physical structure, and a critical fatigue crack length for a material from which the physical structure will be manufactured. The obtaining 602 can be performed, e.g., as described above with reference to FIGS. 3A, 4A, and 5A. In some implementations, obtaining the critical fatigue crack length for the material includes obtaining one or more specifications of the material from which the physical structure will be manufactured, and calculating the critical fatigue crack length for the material from information in the one or more specifications, e.g., SN-curves and engineering tables as described above with reference to FIGS. 5B and 6D. The information can include a modulus of the material's fatigue crack growth curve.

The program iteratively modifies a generatively designed three dimensional shape of the modeled object in the design space in accordance with the one or more design criteria, the one or more in-use load cases for the physical structure, and the one or more specifications. The iteratively modifying can include modifying both a geometry and a topology of the three dimensional shape of the object, e.g., as described above with reference to FIGS. 3A, 4A, 5A, and 6A. The iteratively modifying includes enforcing 604 a design criterion that limits a minimum thickness of the generatively designed three dimensional shape of the modeled object, the minimum thickness being based on the critical fatigue crack length for the material. The minimum thickness can be a ratio of the critical fatigue crack length, such as a minimum thickness that is 10× the critical fatigue crack length. The thickness can be measured using a variety of techniques, e.g., as described above with reference to FIG. 6B.

The generatively designed three dimensional shape of the modeled object is provided 606 for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

As described above with reference to Handling Arbitrary Inequality Constraints and FIGS. 4A-4H, enforcing a minimum thickness inequality constraint can include, in some implementations, using a volume fraction or a minimum thickness based inequality constraint as a proxy for the design criterion that limits the minimum thickness. Using the volume shape derivative results in the fatigue safety being met correctly, however in some implementations the thickness shape derivative replaces the volume shape derivative. The volume fraction or the minimum thickness based inequality constraint is modified using an importance factor, which is set to zero during an initial stage of the iteratively modifying and is adjusted during a subsequent stage of the iteratively modifying based on whether or not one or more other constraints were violated in a prior iteration of the iteratively modifying.

Also as described above with reference to FIGS. 4D-4F, PID control can be used to adjust the target value of the volume fraction or the minimum thickness based inequality constraint between an initial target value and a final target value across multiple iterations of the iteratively modifying. The PID control can also be used to adjust and stabilize changes made in an amount of modification for the modeled object, determined from evaluation of the volume fraction or the minimum thickness based inequality constraint, as the target value is adjusted across the multiple iterations.

Returning to FIG. 6A, a fatigue safety factor inequality constraint for the modeled object can be redefined 612 based on a damage fraction calculated from the required number of loading cycles for the modeled object and the expected number of loading cycles for each of the one or more in-use load cases for the physical structure.

In some implementations, the one or more in-use load cases for the physical structure includes two or more in-use load cases for the physical structure, and the one or more design criteria includes a required number of loading cycles for the modeled object for each of the two or more in-use load cases for the physical structure. In those implementations, the determining 610 the expected number of loading cycles includes determining a separate number of expected loading cycles for each of multiple points e.g., points identified as critical using one or more techniques. Alternatively, the determining 610 is done for each point x on the surface of the design at each iteration as detailed in Equations 132 and 133 on the implicit surface for each of the two or more in-use load cases.

In some implementations, the redefining 612 the fatigue safety factor inequality constraint includes summing, for each of the multiple points, load-specific damage fractions corresponding to the two or more in-use load cases. Each load-specific damage fraction includes an expected number of loading cycles, for one of the multiple points and one of the in-use load cases, divided by the required number of loading cycles for one of the in-use load cases, to produce a sum of the load-specific damage fractions for each of the multiple points. The sums are inverted, and a minimum value of the inverted sums is used to redefine the fatigue safety factor inequality constraint for the modeled object. See e.g., Equation 109.

Shape change velocities for the implicit surface are computed 614 in accordance with at least the fatigue safety factor inequality constraint. Then, the level-set representation is updated 616 using the shape change velocities to produce an updated version of the three dimensional shape of the modeled object. In some implementations, the computing 614 includes computing at least one shape change velocity using an amount determined from a shape derivative formula that approximates a shape derivative of the fatigue safety factor.

The performing 608, the determining 610, the redefining 612, the computing 614, and the updating 616 is repeated until a check 618 determines that a predefined number of shape modification iterations have been performed or that the generatively designed three dimensional shape of the modeled object in the design space has converged to a stable solution for the one or more design criteria and the one or more in-use load cases, e.g., as described above with reference to FIG. 3A.

Design Digital Twin

The relationship between a manufactured part and a designed part can be made into a digital twin so that the difference between the intended design and physical part is reduced. This is particularly prevalent in additively manufactured parts where effects of thickness and build angle can have significant effect on the performance of the part. It is difficult to account for this in the design process, particularly when applying other optimization techniques, e.g., as described above with reference to FIGS. 3A-6E. Therefore, known relationships between the thickness of a feature and its build angles can be applied on a per-element basis, whereas conventional techniques assume an equal strength across the entire part.

FIG. 7A is an example of a process of iteratively modifying a generatively designed three dimensional shape of a modeled object in the design space in accordance with the one or more design criteria, including a stress constraint. Numerical simulation of the modeled object is performed 708 in accordance with a current version of the three dimensional shape and the one or more in-use load cases to produce a current numerical assessment of a physical response of the modeled object. The physical response can be a structural response of the modeled object under the one or more in-use load cases. The numerical simulation can be performed 708, for example, as described above with reference to FIG. 3A.

A maximized stress or strain element is found 710 for each of the one or more in-use load cases for the physical structure from the current numerical assessment of the physical response of the modeled object. The maximum stress or strain element can be the stress, strain, or both, of a point, location, or region of a modeled object. A description of Fatigue Constraints and calculating stress constraints is described above, in reference to FIGS. 5A and 6A.

Dependency on Feature Size and Build Angle

Additively manufactured components (e.g., metal components) have varying strength for different build angle and component thickness values. The Von Mises stress target at a point x is modified as a function of the thickness h(x) and build angle β(x) which is measured with respect to some constant build direction. Thickness can be measured by any suitable technique or combination of techniques, including the techniques described above with reference to FIG. 6B.

FIG. 7B shows a graphical representation of an example of a build angle measured on an object 700. The build angle at a point x on the surface is computed using

$\beta \left(x\right)=\frac{x}{2}-{\cos }^{-1}\left[n\left(x\right),b\right]$

where n(x) is the normal to the surface at x. An internal location for an object is a location that is not on the surface of the object. For those locations, a build angle can be computed by first projecting the internal location along a predetermined build direction, to a surface location of the three dimensional shape. Then, a build angle is determined for the internal location from an angle of the normal of the surface location with respect to the predetermined build direction, using the above-mentioned formula β(x).

The Von Mises stress constraint is used and takes the following form:

$\begin{array}{cc}\underset{\Omega }{\max }\sigma \left(x\right)-{\sigma }_T\le 0.& \left(134\right)\end{array}$

where σ_T is a fixed stress target for each material. This now varies as a function of thickness and build angle resulting in the stress constraint being modified as

[σ(X)−{tilde over (σ)}_T(h(x),β(x))]≤0 ∀x∈Ω.  (135)

The variability on build angle and thickness is transferred from the stress target to the stress itself. For convenience, and while maintaining the same safety factor:

$\begin{array}{cc}\frac{{\sigma }_T}{{\tilde{\sigma }}_T\left(h\left(x\right),\beta \left(x\right)\right)}\sigma \left(x\right)-{\sigma }_T\le 0\forall x\in \Omega .& \left(136\right)\end{array}$

The constraint can be enforced using methods described above with reference to Arbitrary Constraint Handling and FIGS. 4A-4H, with a shape derivative of a combined distance-normal-curvature penalty used and derived below. In some implementations, the volume shape derivative may be used as the proxy shape derivative, as shown in Equation 105.

Next, a derivative of a shape derivative for a distance-normal-curvature penalty is presented. First, derivations for the shape derivatives of signed distance functions, functions of unit normal vectors, and surface curvature functions are presented, respectively. Then, a derivation of the combined penalty from these three shape derivatives is presented.

Shape Derivative of the Signed Distance Function

Let Ω⊂ be a domain with boundary ∂Ω. In this specification the distance function of ∂Ω, which is the mapping d_∂Ω:→ defined by:

$\begin{array}{cc}{d}_{\partial \Omega }\left(x\right):=\underset{yϵ\partial \Omega }{\min }x-y& \left(137\right)\end{array}$

For each x∈ define the set II_∂Ω(x):={y∈∂Ω:d_∂Ω(x)=∥x−y∥} as the set of y∈Ω for which the minimum above is achieved. In the case where Π_∂Ω(x) consists of exactly one point, that point is known as the nearest-point projection of x onto ∂Ω and is denoted by proj_∂Ω(x). Also assign a sign to the distance function as follows: for x∈Ω, sign_∂Ω(x)=−1 and for x∈\Ω sign_∂Ω(x)=+1. The signed distance function of ∂Ω can now be defined as

sdf _∂Ω(x):=sign_∂Ω(x)*d_∂Ω(x).  (138)

It is easy to see that d_∂Ω is 1-Lipschitz, i.e., its difference quotients are uniformly bounded by one. Thus by Rademacher's theorem, it is differentiable almost everywhere, and this is true regardless of the nature of ∂Ω The medial axis of ∂Ω is defined as the set S⊆ where the square of the distance function is not differentiable. d_∂Ω is never differentiable on ∂Ω whereas d_∂Ω² is differentiable on ∂Ω at least when ∂Ω is a differentiable surface. So the medial axis is defined using d_∂Ω² to exclude all but the exceptional points of ∂Ω. By definition, S has a Lebesgue measure equal to zero. It can be shown that the set of points x∈ for which Π_∂Ω(x) consists of at least two distinct points is a subset of S.

For support in the following description, the following facts in differential geometry are presented. For ease of description, assume that the boundary of Ω is “reasonable” in the sense that its tangent planes exist and vary differentiably so that the curvature of ∂Ω can be defined using the tangential derivatives of the unit normal vector field. Also assume ∂Ω is at least C², meaning that an can be expressed as the union of a finite collection of graphs of C² functions.

Let n_∂Ω:∂Ω→ be the unit outward normal vector field of ∂Ω. Let k₁(y), k₂(y) be the principle curvatures at y∈∂Ω with respect to non and let E₁(y), E₂(y) be an orthogonal basis of the tangent plane of ∂Ω at y aligned with the principle curvature directions at y.

Proposition 1. Let Ω be a domain with C² boundary. Suppose x∉S. Then sdf_∂Ω is twice-differentiable at x. If in addition, x∉∂Ω, then d_∂Ω is also twice-differentiable at x. A point x∈ has a unique projection onto ∂Ω if and only if x∉S. In such a case, d_∂Ω=∥x−proj_∂Ω(x)∥ and the gradient of sdf_∂Ω at x is

$\begin{array}{cc}\nabla {\mathrm{sdf}}_{\partial \Omega }\left(x\right)=n_{\partial \Omega }\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)=\frac{x-{\mathrm{proj}}_{\partial \Omega }\left(x\right)}{{\mathrm{sdf}}_{\partial \Omega }\left(x\right)}\mathrm{when}\notin \partial \Omega & \left(139\right)\end{array}$

For every x∈ and every y∈Π_∂Ω(x), the following holds:

−k₁(proj_∂Ω(x))×sdf_∂Ω(x)≤1 ∀i=1,2  (140)

Moreover, the closure of S consists of all points in S together with all points for which at least one of the inequalities above is strict.

Suppose x is not in the closure of S. Then

$\begin{array}{cc}{\nabla }^2{\mathrm{sdf}}_{\partial \Omega }\left(x\right)={\sum }_{i=1}^2\frac{k_i(\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)}{1+k_i\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)\times {\mathrm{sdf}}_{\partial \Omega }\left(x\right)}{E}_i\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)\otimes E_i\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)& \left(141\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\nabla {\mathrm{proj}}_{\partial \Omega }\left(x\right)={\sum }_{i=1}^2\frac{1}{1+k_i\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)\times {\mathrm{sdf}}_{\partial \Omega }\left(x\right)}{E}_i\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)\otimes E_i\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)& \left(142\right)\end{array}$

Let T_∈:→ be a one-parameter family of transformations of arising as the flow of a vector field Θ:→. This means: T_ε0 satisfies the equations:

$\begin{array}{cc}\frac{\partial T_{\epsilon }}{\partial \epsilon }=\Theta \circ T_{\epsilon }& \left(143\right)\end{array}$ $\begin{array}{cc}{T}_0=\mathrm{Identity}& \left(144\right)\end{array}$

Recall that the Lagrangian derivative of a shape-dependent function f_Ω:→ with respect to this variation is defined as:

$\begin{array}{cc}{f}_{\Omega \left(x\right)}:=\frac{\partial }{\partial \epsilon }{f}_{T_{\epsilon }\left(\Omega \right)}\left(T_{\epsilon }\left(x\right)\right){|}_{\epsilon =0}& \left(145\right)\end{array}$

and the Eulerian derivative of f_Ω with respect to this variation is defined as:

f _Ω′(x):={dot over (f)}Ω(x)−∇f_Ω(x),Θ(x)  (146)

The next result gives the Eulerian derivative of the signed distance function.

Proposition 2. If x does not belong to the closure of S, then the Eulerian derivative of the signed distance function with respect to the variation generated by a vector field Θ at x satisfies:

sdf _∂Ω′(x)=−Θ(proj_∂Ω(x)),n_∂Ω(proj_∂Ω(x))  (147)

A proof follows. First, compute the Lagrangian derivative of the signed distance function with respect to the variation generated by a vector field Θ at x. For this sdf_Tε(∂Ω) (T_ε(x))=∥T_ε(x)−y_ε∥ is used where y_ε minimizes ∥T_ε(x)−y∥ over all y∈T_ε(∂Ω). The result is

[sdf_∂Ω(x)]=n_∂Ω(proj_∂Ω(x)),Θ(x)−Θ(proj_∂Ω(x))  (148)

which is obtained with the help of differentiating the condition sdf_∂Ω(T_ε⁻¹())=0. Similarly, it can be shown that the gradient term in the definition of the Eulerian derivative is n_∂Ω(proj_∂Ω(x)), Θ(x). The desired formula follows.

Consider a shape function of the form:

Φ(Ω):=∫_Ωϕ(x,sdf_∂Ω(x))dx  (149)

where ϕ is a smooth function of its arguments. The shape derivative of ϕ is computed using the standard tools of shape calculus, as well as the formula for the Eulerian derivative of the signed distance function.

Proposition 3. Let T_ε:→ be a one-parameter family of transformations of arising as the flow of a vector field Θ: →. Then

$\begin{array}{cc}\frac{d}{d\epsilon }\varphi \left(T_{\epsilon }\left(\Omega \right)\right){|}_{\epsilon =0}=\text{⁠}-{\int }_{\Omega }{\partial }_2\varphi \left(x,{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right){\Theta }^{\perp }\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)\mathrm{dx}+{\int }_{\partial \Omega }\varphi \left(x,0\right){\Theta }^{\perp }\left(x\right)d\sigma \left(x\right)& \left(150\right)\end{array}$

where Θ^⊥is the normal component of Θ on ∂Ω and ∂_2φ is the partial derivative of c with respect to its second argument. If the Eulerian derivative of a quantity is denoted by a prime, then

$\begin{array}{cc}\begin{array}{c}\frac{d}{d\epsilon }\varphi \left(T_{\epsilon }\left(\Omega \right)\right){|}_{\epsilon =0}={\int }_{\Omega }{\left[\varphi \left(x,{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right)\right]}^{\prime }\mathrm{dx}+\\ {\int }_{\partial \Omega }\varphi \left(x,{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right){\Theta }^{\perp }\left(x\right)d\sigma \left(x\right)\\ ={\int }_{\Omega }{\partial }_2\varphi {\left(x,{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right)\left[{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right]}^{\prime }\mathrm{dx}+\\ {\int }_{\partial \Omega }\varphi \left(x,0\right){\Theta }^{\perp }\left(x\right)d\sigma \left(x\right)\\ =-{\int }_{\Omega }{\partial }_2\varphi \left(x,{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right)\mathrm{dx}+\\ {\int }_{\partial \Omega }\varphi \left(x,0\right){\Theta }^{\perp }\left(x\right)d\sigma \left(x\right).\end{array}& \left(151\right)\end{array}$

Observe that the first integral appearing in the formula above is not in the standard form of the Hadamard-Zolésio structure theorem, namely a surface integral over ∂Ω. The formula does still depend on Θ^⊥ on ∂Ω, so the theorem still holds. Using a change of variables based on the co-area formula, this integral can be expressed in the standard form. The following result is due to Dapogny et al., Geometric Constraints for Shape and Topology Optimization in Architectural Design. Computational Mechanics, Springer Verlag, (2017, 59(6), pp. 933-965).

Proposition 4. Let Ω be a domain with C² boundary. Then

$\begin{array}{cc}{\int }_{\Omega }{\partial }_2\varphi \left(x,{\mathrm{sdf}}_{\partial \Omega }\left(x\right)\right){\Theta }^{\perp }\left({\mathrm{proj}}_{\partial \Omega }\left(x\right)\right)\mathrm{dx}={\int }_{\partial \Omega }\left({\int }_0^{T\left(y\right)}\frac{{\partial }_2\varphi \left(y-\tau n_{\partial \Omega }\left(y\right),\tau \right)}{J\left(y,\tau \right)}d\tau \right){\Theta }^{\perp }\left(\right)d\sigma \left(y\right)& \left(152\right)\end{array}$

where T(y) is the distance to the medial axis of ∂Ω along the ray emanating inwards from y∈∂Ω and J(y, τ) is the Jacobian factor of the change of variables, which is

[J(y,T)]⁻¹=Π_i=1²(1+τk_i(y))  (153)

Note the Jacobian factor above has a particularly nice alternate form. Expanding the product yields:

$\begin{array}{cc}\begin{array}{c}{\left[J\left(y,T\right)\right]}^{-1}=1+\tau \left(k_1\left(y\right)\right)+k_2\left(y\right))+{\tau }^2{k}_1\left(y\right)k_2\left(y\right)\\ =1+\tau H_{\partial \Omega \left(y\right)}+{\tau }^2{k}_{\partial \Omega \left(y\right)}\end{array}& \left(154\right)\end{array}$

where H_∂Ω and K_∂Ω are the mean curvature and the Gauss curvature of ∂Ω respectively. This is because these are the invariants of the second fundamental form of ∂Ω, of which the principal curvatures are the eigenvalues.

Shape Derivative of a Function of the Unit Normal Vector

Let Σ be any oriented, smooth surface embedded in Euclidean space and let Y:→ be a smooth vector field defined on the background Euclidean space and that is completely independent of Σ. Consider shape functions of the form:

Φ(Σ):=∫_Σϕ(x,(Y(x),N_Σ(x))dσ(x))  (155)

where ϕ: ×→ is a smooth function of its arguments. The objects appearing in the integrand above are: the unit normal vector N_Σ of Σ, the surface area element dα of Σ. The shape derivative of Φ at a given shape Σ is calculated with respect to any given shape variation. This means: let T_ε(Σ) be a variation of E generated by a one-parameter family of transformations T_ε of Euclidean space; then the shape derivative with respect to this variation is

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}& \left(156\right)\end{array}$

To calculate the shape derivative of Φ, the transformation T_ε is first described in more detail. There are a number of ways to proceed, and all are equivalent to first order in ε. The velocity method of Hadamard is chosen for the following description because it is well-suited to level set-based shape optimization, but it is understood that other methods can be used, such as those described in M. C. Delfour and J. P. Zolésio, Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, (2^nd Ed. 2011).

Suppose T_ε:→ is the transformation of arising as the flow of a vector field Θ:→. In other words, T_ε satisfies the equations

$\begin{array}{cc}\frac{\partial T_{\epsilon }}{\partial \epsilon }=\Theta \circ T_{\epsilon }& \left(157\right)\end{array}$ $T_0=\mathrm{Identity}.$

∇^∥ is the tangential gradient operator, div^∥ is the tangential divergence operator, and H_Σ is the mean curvature of Σ. Also, ϕ(x,Y(x), N_Σ(x)) as ϕ(x,q(x)) with q(x):=Y(x), N_Σ(x) is written to emphasize the structure of the integrand, and the partial derivative of ϕ is denoted with respect to its second argument as

$\frac{\partial \varphi }{\partial q}.$

Proposition 1. Let Σ be an oriented, smooth surface embedded in Euclidean space and let T_ε be a variation of Σ generated by a vector field Θ:→. The shape derivative of Φ at Σ is given by

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi (T_{\epsilon }(\sum )❘}_{\epsilon =0}={\int }_{\sum }\left(\mathrm{div}(\frac{\partial \varphi }{\partial q}\left(x,q\right)Y\right)+\frac{\partial \varphi \left(x,q\right)}{\partial n}-\varphi \left(x,q\right)H_{\sum }){\Theta }^{\perp }d\sigma & \left(158\right)\end{array}$

where Θ^⊥ is the normal component of Θ on Σ, and

$\frac{\partial }{\partial n}$

means the total directional derivative in the normal direction. A proof follows. Apply the calculus of Eulerian derivatives to obtain:

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}={\int }_{\sum }\frac{\partial \varphi }{\partial q}\left(x,q\right)q^{\prime }d\sigma +{\int }_{\sum }\left(\frac{\partial \varphi \left(x,q\right)}{\partial n}-\varphi \left(s,q\right)H_{\sum }\right){\Theta }^{\perp }d\sigma & \left(159\right)\end{array}$

where the prime denotes the Eulerian derivative.

q′=Y,N_Σ′=Y′,N_Σ+Y′,N_Σ+Y,N_Σ′, =Y^∥,∇^∥Θ^⊥, because Y is shape-independent so its Eulerian derivative vanishes, and N_Σ′=−∇^∥Θ^⊥. Thus

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}=-{\int }_{\sum }\frac{\partial \varphi }{\partial q}\left(x,q\right)\langle Y^{\parallel },{\nabla }^{\parallel }{\Theta }^{\perp }\rangle d\sigma +{\int }_{\sum }\left(\frac{\partial \varphi \left(x,q\right)}{\partial n}-\varphi \left(x,q\right)H_{\sum }\right){\Theta }^{\perp }d\sigma ={\int }_{\sum }\left({\mathrm{div}}^{\parallel }\left(\frac{\partial \varphi }{\partial q}\left(x,q\right)Y^{\parallel }\right)+\frac{\partial \varphi \left(x,q\right)}{\partial n}-\varphi \left(x,q\right)H_{\sum }\right){\Theta }^{\perp }d\sigma & \left(160\right)\end{array}$

after integrating by parts.

Consider a specific application of the formula derived in the Equation 160. That is, the formula to the function ϕ in Dapogny et al., is applied, namely

ϕ(X,N_Σ):=∥Y−N_Σ∥²=2−2X,N_Σ  (161)

where Y is a smooth, shape-independent unit vector field defined on Euclidean space. Since as a function of its argument, ϕ(x,q)=2−2q, then

$\frac{\partial \varphi }{\partial q}$

is the constant function

$\frac{\partial \varphi }{\partial q}\left(q\right)=-2.$

Therefore Proposition 1 yields in this case and:

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}=-2{\int }_{\sum }\left({\mathrm{div}}^{\parallel }\left(Y^{\parallel }\right)+\left(\frac{\partial \langle Y,N_{\sum }\rangle }{\partial n}+\left(1-\langle Y,N_{\sum }\rangle \right)H_{\sum }\right)\right){\Theta }^{\perp }d\sigma .& \left(162\right)\end{array}$

At this point the identity

$\mathrm{div}\left(Y\right)=\frac{\partial Y^{\perp }}{\partial n}-Y^{\perp }{H}_{\sum }+{\mathrm{div}}^{\parallel }\left(Y^{\parallel }\right)$

is used that is valid for any vector field on the background Euclidean space and restricted to Σ. Therefore,

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}=-2{\int }_{\sum }\left(\mathrm{div}\left(Y\right)+H_{\sum }\right){\Theta }^{\perp }d\sigma & \left(163\right)\end{array}$

This formula is almost identical to the formula in the paper to Dapogny, et al. There are two differences: the first is that only the tangential divergence operator appears in paper to Dapogny; the second is that the mean curvature appears with the opposite sign in paper to Dapogny. The first difference is very likely due to an error, whereas the second difference is due to the fact that the paper to Dapogny, uses the opposite sign convention for the mean curvature. In the paper, H_Σ is defined with respect to the inward unit normal, and here it is defined with respect to the outward unit normal. (The former convention has the advantage that the mean curvature of a sphere is positive and equal to 2, while the latter convention has the advantage that the theoretical formulas contain fewer minus signs in total despite the fact that the mean curvature of the unit sphere is −2.)

One can show that the mean curvature of a level set of the form {x∈:F(x)=0} for some level set function F:→ is equal to −div(∇F/∥∇F∥) where ∇F/∥∇F∥ is the unit normal vector field of the level set. Thus if the Y in the formula of Proposition 1 is of this form, the shape derivative can be interpreted as the difference of the mean curvature of Σ and the mean curvature of the level set.

A second example concerns the 2.5D manufacturability penalty. In this case, the penalty function Φ is used as before, but with integrand the function is

ϕ(x,q(x)):=[q(x)]²E(x)  (164)

where q(x):=Y(x),N_Σ(x) as before but with a crucial difference−Y is in fact a constant vector field which represents the milling direction, thus

$\frac{\partial \varphi }{\partial q}\left(q\right)=2qE\left(x\right).$

The function E(x) is a product of offset exponential bump functions. Therefore it follows that

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}={\int }_{\sum }\left(di{v}^{\parallel }\left(2{\mathrm{qEY}}^{\parallel }\right)+\frac{\partial }{\partial n}\left(q^{2}E\right)-q^{2}E{H}_{\sum }\right){\Theta }^{\perp }d\sigma ={\int }_{\sum }\left(2\langle {\nabla }^{\parallel }\left(qE\right),Y^{\parallel }\rangle +\frac{\partial }{\partial n}\left(q^{2}E\right)-q^{2}E{H}_{\sum }\right){\Theta }^{\perp }d\sigma & \left(165\right)\end{array}$

since Y is constant.

Shape Derivative of Surface Curvature Functions

Let Σ be a smooth oriented surface in . Consider shape functions which depend on the surface curvature invariants of Σ (i.e. the mean curvature and the total curvature) of the form:

Φ(Σ):=∫_Σϕ(H_Σ,∥A_Σ∥²)dσ  (166)

where H_Σ is the mean curvature of Σ, ∥A_Σ∥² is the total curvature of Σ, and ϕ is a smooth function of its arguments.

A description follows for computing a shape derivative for shape functions of this form. In other words, it will be shown how to compute the derivative

$\begin{array}{cc}{\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right)❘}_{\epsilon =0}& \left(167\right)\end{array}$

for any variation T_ε(Σ) where T_ε:→ is a one-parameter family of transformations of arising as the flow of a vector field Θ:→. This means: T_ε satisfies the equations

$\begin{array}{cc}\frac{\partial T_{\epsilon }}{\partial \epsilon }=\Theta ·T_{\epsilon }& \left(168\right)\end{array}$

T₀=Identity.

The principal geometric objects associated to the surface Σ are the induced surface metric tensor h_Σ, i.e., the first fundamental form, the unit normal vector field N_Σ, and the second fundamental form A_Σ. A calculation of the derivatives of these components with respect to the variation introduced above is described above (or rather their components with respect to a local basis of the tangent bundle of Σ because these quantities are tensorial).

To carry out the calculations, some preliminary notions are introduced in the context of embedded surfaces in Riemannian geometry. First, a local parametrization of Σ is used near an arbitrary point x∈Σ to parametrize a sufficiently small tubular neighborhood of Σ near x by means of the mapping (ε,x)T_ε(x). Under this parametrization, the ε-direction maps to the vector field Θ. A pair of vector fields {tilde over (E)}₁, {tilde over (E)}₂ is introduced that are linearly independent sufficiently near Σ and tangent to T_ε(Σ) for every ε. This is achieved by pushing the local coordinate basis for the tangent bundle of Σ with respect to the parametrization forward to a sufficiently small tubular neighborhood of Σ. This means: if Σ^∘₁, E^∘₂ denotes this basis, then the following definition is used E_i(y):=DT_ε(E^∘_i(x)) for any y sufficiently near Σ of the form y=T_ε(x) for some x∈Σ. Here DT_ε is the matrix of partial derivatives of T_ε. Second, a vector field N is defined on this tubular neighborhood with the property that if y sufficiently near Σ is of the form y=T_ε(x) for some x∈Σ, then N(y) is the unit normal vector field of T_ε(Σ) at y. The notation X=X^∥+X^⊥N for the orthogonal decomposition of a vector field X defined on this tubular neighborhood with respect to N. Note that if y=T_ε(x) for some x∈Σ, then X^∥(y) is tangent to T_ε(Σ).

Let h_ε and A_ε be the pull-backs of the induced surface metric tensor and second fundamental form of T_ε(Σ), respectively, to Σ. Then their components with respect to the local basis {tilde over (E)}₁, {tilde over (E)}₂ satisfy

[h_ε]_ij={tilde over (E)}_i,{tilde over (E)}_j∘T_ε and [A_ε]_ij=N,∇_{{tilde over (E)}i}{tilde over (E)}_j∘T_ε  (169)

Let [h_ε]^ij denote the components of the inverse of the matrix h_ε. Note that the sub scripted ε is dropped when referring to these quantities on Σ.) Next, the vector field N satisfies

N,N ≡1 and N,{tilde over (E)}_i≡0.  (170)

Finally, by differentiating the defining equations for N, the following is deduced:

∇_{{tilde over (E)}i}N,N≡0 and {tilde over (E)}_i,∇{tilde over (E)}_jN=−[A_ε]_ij  (171)

These two equations uniquely determine the vector field ∇_{{tilde over (E)}i}N. In fact, one can see that these equations establish ∇_{{tilde over (E)}i}N=−Σ_jk[h_ε]^jk[A_ε]_ij{tilde over (E)}_k. on T_ε(Σ) A related result is the following lemma 1.

Lemma 1. The following formula holds:

∇_ΘN=−Σ_klh^kl(A_Σ(E_k,Θ^∥)+E_k(Θ^⊥))E_l.  (172)

Proof. The equation N,N≡1 is differentiated in the Θ-direction to deduce that ∇_ΘN must be tangential. Then the equations N,{tilde over (E)}_i≡1 for i=1, 2 are differentiated in the Θ-direction to deduce

$\begin{array}{cc}\langle {\nabla }_{\Theta }N,E_i\rangle =-\langle N,{\nabla }_{\Theta }{\tilde{E}}_i\rangle =-\langle N,{\nabla }_{{\tilde{E}}_i}\Theta \rangle =-\langle N,{\nabla }_{{\tilde{E}}_i}\left({\Theta }^{\parallel }+{\Theta }^{\perp }N\right)\rangle =-\langle N,{\nabla }_{{\tilde{E}}_i}{\Theta }^{\parallel }+{\tilde{E}}_i\left({\Theta }^{\perp }\right)N+{\Theta }^{\perp }{\nabla }_{{\tilde{E}}_i}N\rangle =-A_{\sum }\left(E_i{\Theta }^{\parallel }\right)-E_i\left({\Theta }^{\perp }\right)& \left(173\right)\end{array}$

for i=1, 2. This implies the result.

Lemma 2. Suppose T_ε is the variation of Σ generated by the vector field Θ and let [h_ε]_ij be the components of the pull-back of the induced surface metric tensor of T_ε(Σ) to Σ with respect to any local basis E₁, E₂ for the tangent bundle of Σ. Then

$\begin{array}{cc}{\frac{\partial {\left[h_{\epsilon }\right]}_{ij}}{\partial \epsilon }❘}_{\epsilon =0}=\langle {\nabla }_{E_i}^{\parallel }{\Theta }^{\parallel },E_j\rangle +\langle {\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel },E_i\rangle -{2\left[A_{\sum }\right]}_{ij}{\Theta }^{\perp }& \left(174\right)\end{array}$

where A_ij are the components of the second fundamental form of Σ with respect to the local basis and ∇^∥ denotes the surface covariant derivative operator.

Proof. The calculations can be performed in any local basis since it can be shown that the results are independent of the basis (this is a foundational principle of differential geometry). Therefore the local basis introduced above is used at an arbitrary x∈Σ. Hence by the definition and properties of covariant derivatives in the following results

$\begin{array}{cc}{\frac{\partial {\left[h_{\epsilon }\right]}_{ij}}{\partial \epsilon }❘}_{\epsilon =0}=\Theta \langle {\tilde{E}}_i,{\tilde{E}}_j\rangle =\langle {\nabla }_{\Theta }{\tilde{E}}_i,E_j\rangle +\mathrm{same}\mathrm{but}\mathrm{with}i\leftrightarrow j=\langle {\nabla }_{{\overline{E}}_i}\Theta ,E_j\rangle +\mathrm{same}\mathrm{but}\mathrm{with}i\leftrightarrow j& \left(175\right)\end{array}$

since the difference between the covariant derivatives above is the Lie bracket [Θ, {tilde over (E)}_i] which vanishes because both θ and {tilde over (E)}_i are push-forwards of coordinate basis vector fields. (In this calculation the notation

$\mathrm{Xf}:={\sum }_k{X}^k\frac{\partial f}{\partial x^k}$

is used, where ƒ:→ is any differentiable function; this notation reflects the standard differential geometric practice of conflating the concept of a vector field X and the directional derivative operator in the X-direction.) Next, using the definition of the second fundamental form of Σ, the following is obtained:

$\begin{array}{cc}\begin{array}{c}\frac{\partial {\left[h_{\epsilon }\right]}_{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}={\tilde{E}}_i\langle \Theta ,{\tilde{E}}_j\rangle -\langle \Theta ,{\nabla }_{{\tilde{E}}_i}{\tilde{E}}_j\rangle +\mathrm{same}\mathrm{but}\mathrm{with}i\leftrightarrow j\\ =E_i\langle {\Theta }^{\parallel },E_j\rangle -\langle {\Theta }^{\parallel },{\nabla }_{E_i}^{\parallel }{E}_j\rangle -{\Theta }^{\perp }\langle N_{\sum },{\nabla }_{{\tilde{E}}_i}{\tilde{E}}_j\rangle +\\ =\mathrm{same}\mathrm{but}\mathrm{with}i\leftrightarrow j\\ \langle {\nabla }_{E_i}^{\parallel }{\Theta }^{\parallel },E_j\rangle +\langle {\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel },E_i\rangle -2{\Theta }^{\perp }{A}_{\mathrm{ij}}\end{array}& \left(176\right)\end{array}$

Lemma 3. Suppose T_ε is the variation of Σ generated by the vector field Θ and let [A_ε]_ij be the components of the pull-back of the second fundamental form of T_ε(Σ) to Σ with respect to the local basis introduced above. Then

$\begin{array}{cc}\frac{\partial {\left[A_{\epsilon }\right]}_{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}=E_i\left(E_j\left({\Theta }^{\perp }\right)\right)-{\nabla }_{E_i}^{\parallel }{E}_j\left({\Theta }^{\perp }\right)-{\left[A_{\sum }^2\right]}_{\mathrm{ij}}{\Theta }^{\perp }+{\nabla }_{E^i}^{\parallel }{A}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)+A_{\sum }\left(E_i,{\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel }\right)+A_{\sum }\left(E_j,{\nabla }_{E_i}^{\parallel }{\Theta }^{\parallel }\right).& \left(177\right)\end{array}$

Proof. Again, the calculations can be performed in any local basis since it can be shown that the results are independent of the basis. Therefore the local basis introduced above is used at an arbitrary x∈Σ. By the definition of the pull-back of the second fundamental form of T_ε(Σ), its components with respect to this local basis satisfy [A_ε]_ij=N,∇_{{tilde over (E)}i}{tilde over (E)}_j∘T_ε. Hence by the definition and properties of covariant derivatives in the following results:

$\begin{array}{cc}\frac{\partial {\left[A_{\epsilon }\right]}_{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}=\Theta \langle N,{\nabla }_{{\tilde{E}}_i}{\tilde{E}}_j\rangle =\langle {\nabla }_{\Theta }N,{\nabla }_{{\tilde{E}}_i}{\tilde{E}}_j\rangle +\langle N,{\nabla }_{\Theta }{\nabla }_{{\tilde{E}}_i}{\tilde{E}}_j\rangle .& \left(178\right)\end{array}$

Using Lemma 1 in the first term of Equation 178 gives

∇_ΘN,∇_{{tilde over (E)}i}{tilde over (E)}_j=A_Σ(∇_Ei^∥E_j,Θ^∥)−∇_Ei^∥E_j(Θ^⊥.  (179)

The following deductions are made for the second term of Equation 178:

$\begin{array}{cc}\begin{array}{c}\langle N,{\nabla }_{\Theta }{\nabla }_{{\tilde{E}}_i}{\tilde{E}}_j\rangle =\langle N,{\nabla }_{{\tilde{E}}_i}{\nabla }_{{\tilde{E}}_j}\Theta \rangle \\ =\langle N,{\nabla }_{{\tilde{E}}_i}\left({\nabla }_{{\tilde{E}}_j}{\Theta }^{\parallel }+E_j\left({\Theta }^{\perp }\right)N+{\Theta }^{\perp }{\nabla }_{{\tilde{E}}_j}N\right)\rangle \\ =\langle N,{\nabla }_{{\tilde{E}}_i}({\nabla }_{{\tilde{E}}_j}^{\parallel }{\Theta }^{\parallel }+A_{\sum }\left(E_j,{\Theta }^{\parallel }\right)N+\\ E_j\left({\Theta }^{\perp }\right)N+{\Theta }^{\perp }{\nabla }_{{\tilde{E}}_j}N)\rangle \\ =A_{\sum }\left(E_i,{\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel }\right)+E_i\left(A_{\sum }\left(E_j,{\Theta }^{\parallel }\right)\right)+\\ E_i\left(E_j\left({\Theta }^{\perp }\right)\right)-{\left[A_{\sum }^2\right]}_{\mathrm{ij}}{\Theta }^{\perp }\end{array}& \left(180\right)\end{array}$

since N,N=1 and the equations uniquely characterizing ∇_{{tilde over (E)}i}N together imply that

N,∇ _{{tilde over (E)} i} ∇_{{tilde over (E)}j}N=−∇_{{tilde over (E)}i}N,∇_{{tilde over (E)}i}N=Σ_klh^klA_ikA_jl=−[A_Σ²]_ij  (181)

Substituting the two derivations into Equation 178 now yields

$\begin{array}{cc}\begin{array}{c}\frac{\partial {\left[A_{\epsilon }\right]}_{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}=E_i\left(E_j\left({\Theta }^{\perp }\right)\right)-{\nabla }_{E_i}^{\parallel }{E}_j\left({\Theta }^{\perp }\right)-{\left[A_{\sum }^2\right]}_{\mathrm{ij}}{\Theta }^{\perp }-\\ A_{\sum }\left({\nabla }_{E_i}^{\parallel }{E}_j,{\Theta }^{\parallel }\right)+A_{\sum }\left(E_i,{\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel }\right)+E_i\left(A_{\sum }\left(E_j,{\Theta }^{\parallel }\right)\right)\\ =E_i\left(E_j\left({\Theta }^{\perp }\right)\right)-{\nabla }_{E_i}^{\parallel }{E}_j\left({\Theta }^{\perp }\right)-{\left[A_{\sum }^2\right]}_{\mathrm{ij}}{\Theta }^{\perp }+\\ {\nabla }_{E_i}^{\parallel }{A}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)+A_{\sum }\left(E_i,{\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel }\right)+A_{\sum }\left(E_j,{\nabla }_{E_i}^{\parallel }{\Theta }^{\parallel }\right)\end{array}& \left(182\right)\end{array}$

using the definition ∇_Ei^∥A_Σ(E_j,Θ^∥)=−−A_Σ(∇_Ei^∥E_j,Θ^∥)−A_Σ(E_j,∇_Ei^∥Θ^∥)+E_i(A_Σ(E_j,Θ^∥)) of the surface covariant derivative of a covariant two-tensor.

The surface curvature invariants of Σ are: the total curvature, i.e., the squared norm of the second fundamental form with respect to the induced surface metric tensor, denoted ∥A_Σ∥²; and the mean curvature, i.e., the trace of the second fundamental form with respect to the induced surface metric tensor, denoted H_Σ. Note that the Gauss curvature or intrinsic curvature of Σ is given in terms of these invariants by the formula G_Σ:=(H_Σ²−∥A_Σ∥²)/2. The derivatives of these quantities are now given with respect to the variation introduced above.

Proposition 4. Suppose T_ε is the variation of Σ generated by the vector field Θ and let H_ε0 be the mean curvature of T_ε(Σ). Then

$\begin{array}{cc}\frac{\partial h_{\epsilon }\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}={\Delta }^{\parallel }{\Theta }^{\perp }+{A_{\sum }}^2{\Theta }^{\perp }+{\sum }_{\mathrm{ij}}{h}^{\mathrm{ij}}{\nabla }_{E_i}^{\parallel }{A}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)& \left(183\right)\end{array}$

where Δ^∥ is the surface Laplace-Beltrami operator.

Proof. Lemma 2 and 3 is applied to the definition H_ε∘T_ε:=[h_ε]^ij[A_ε]_ij. That is,

$\begin{array}{cc}\begin{array}{c}\frac{\partial H_{\epsilon }\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}=\frac{\partial {\left[h_{\epsilon }\right]}^{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}{A}_{\mathrm{ij}}+h_{\mathrm{ij}}\frac{\partial {\left[A_{\epsilon }\right]}_{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}\\ =-h^{\mathrm{ik}}{h}^{\mathrm{jl}}\left(\langle {\nabla }_{E_k}^{\parallel }{\Theta }^{\parallel },E_l\rangle +\langle {\nabla }_{E_l}^{\parallel }{\Theta }^{\parallel },E_k\rangle -2{\Theta }^{\perp }{A}_{\mathrm{kl}}\right)A_{\mathrm{ij}}+\\ \sum _{\mathrm{ij}}{h}^{\mathrm{ij}}\left(E_i\left(E_j\left({\Theta }^{\perp }\right)\right)-{\nabla }_{E_i}^{\parallel }{E}_j\left({\Theta }^{\perp }\right)-{\left[A_{\sum }^2\right]}_{\mathrm{ij}}{\Theta }^{\perp }\right)+\\ \sum _{\mathrm{ij}}{h}^{\mathrm{ij}}\left(\begin{array}{c}{\nabla }_{E_i}^{\parallel }{A}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)+A_{\sum }\left(E_i,{\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel }\right)+\\ A_{\sum }\left(E_j,{\nabla }_{E_i}^{\parallel }{\Theta }^{\parallel }\right)\end{array}\right)\\ ={\nabla }^{\parallel }{\Theta }^{\perp }+{A_{\sum }}^2{\Theta }^{\perp }+\sum _{\mathrm{ij}}{h}^{\mathrm{ij}}{\nabla }_{E_i}^{\parallel }{A}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)\\ ={\Delta }^{\parallel }{\Theta }^{\perp }+{A_{\sum }}^2{\Theta }^{\perp }+\langle {\nabla }^{\parallel }{H}_{\sum },{\Theta }^{\parallel }\rangle \end{array}& \left(184\right)\end{array}$

As required, λ^∥0 is simply identified where it appears and the Codazzi equation is used to replace the ∇^∥A_Σ with the ∇^∥H_Σ term.

Proposition 5. Suppose T_Σ is the variation of Σ generated by the vector field Θ and let ∥A_ε∥² be the total curvature of T_ε(Σ) Then

$\begin{array}{cc}\frac{\partial {A_{\epsilon }}^2\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}=2A_{\sum }:\left({\nabla }^{\parallel }{\nabla }^{\parallel }{\Theta }^{\perp }\right)+2\mathrm{Tr}\left(A_{\sum }^3\right){\Theta }^{\perp }+2{\sum }_{\mathrm{ijkl}}{h}^{\mathrm{ij}}{h^{\mathrm{kl}}\left[A_{\sum }\right]}_{\mathrm{kl}}{A}_{\sum }\left({\nabla }_{E_i}^{\parallel }{A}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)\right)& \left(185\right)\end{array}$

Proof. Apply Lemma 2 and 3 to the definition ∥A_ε∥²∘T_ε:=[h_ε]^ik[h_ε]^jl[A_ε]_ij[A_ε]_kl. The required calculations are similar to those of the proof of the previous lemma.

Variation of Integrated Surface Curvatures

In this section, surface integrals of the surface curvature invariants of Σ are considered. First, the following preliminary result is described.

Lemma 6. Suppose T_ε is the variation of E generated by the vector field Θ and let [h₂₄₉]_ij be the components of the pull-back of the induced surface metric tensor of T_ε(Σ) to Σ with respect to any local basis E₁, E₂ for the tangent bundle of Σ. Let ϕ_Σ:Σ→ be a surface-dependent scalar function that is shape-differentiable. Then

$\begin{array}{cc}\frac{d}{d\epsilon }{\int }_{T_{\epsilon }(\sum )}{\varphi }_{{\sum }_{\epsilon }}d{\sigma }_{\epsilon }{❘}_{\epsilon =0}={\int }_{\sum }\frac{\partial {\varphi }_{\sum \epsilon }\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}d\sigma +{\int }_{\sum }{\varphi }_{\sum }\left({\mathrm{div}}^{\parallel }\left({\Theta }^{\parallel }\right)-H_{\sum }{\Theta }^{\perp }\right)d\sigma & \left(186\right)\end{array}$

where div^∥ is the surface divergence operator and H_Σ is the mean curvature of Σ.

Proof. By the change of variables formula for surface integrals,

$\begin{array}{cc}\begin{array}{c}\frac{d}{d\epsilon }{\int }_{T_{\epsilon }(\sum )}{\varphi }_{{\sum }_{\epsilon }}d{\sigma }_{\epsilon }{❘}_{\epsilon =0}={\int }_{\sum }\frac{1}{\sqrt{\det \left(h\right)}}\frac{\partial }{\partial \epsilon }\\ \left({\varphi }_{{\sum }_{\epsilon }}\circ T_{\epsilon }\sqrt{\det \left(h_{\epsilon }\right)}\right){❘}_{\epsilon =0}d\sigma \\ ={\int }_{\sum }\frac{\partial {\varphi }_{{\sum }_{\epsilon }}\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}d\sigma +\\ {\int }_{\sum }\frac{{\varphi }_{\sum }}{\sqrt{\det \left(h\right)}}\frac{\partial \sqrt{\det \left(h_{\epsilon }\right)}}{\partial \epsilon }{❘}_{\epsilon =0}d\sigma \end{array}& \left(187\right)\end{array}$

To conclude, Lemma 2 and the following calculation is used:

$\begin{array}{cc}\begin{array}{c}\frac{\partial \sqrt{\det }\left(h_{\epsilon }\right)}{\partial \epsilon }{❘}_{\epsilon =0}=\frac{1}{2}\sqrt{\det \left(h\right)}\mathrm{Tr}\left(\sum _{\mathrm{ij}}{h}^{\mathrm{ij}}\frac{\partial {\left[h_{\epsilon }\right]}_{\mathrm{ij}}}{\partial \epsilon }{❘}_{\epsilon =0}\right)\\ \frac{1}{2}\sqrt{\det \left(h\right)}\mathrm{Tr}(\sum _{\mathrm{ij}}{h}^{\mathrm{ij}}(\langle {\nabla }_{E_i}^{\parallel }{\Theta }^{\parallel },E_j\rangle +\\ \langle {\nabla }_{E_j}^{\parallel }{\Theta }^{\parallel },E_i\rangle -2{\Theta }^{\perp }{A}_{\mathrm{ij}}))\\ =\sqrt{\det \left(h\right)}\left({\mathrm{div}}^{\parallel }\left({\Theta }^{\parallel }\right)-H_{\sum }{\Theta }^{\perp }\right)\end{array}.& \left(188\right)\end{array}$

Proposition 7. Let ϕ be a surface integral operator of the form stated in the Introduction. Suppose T_ε is the variation of Σ generated by the vector field Θ. Then

$\begin{array}{cc}\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right){❘}_{\epsilon =0}={\int }_{\sum }{\partial }_1\varphi \left(H_{\sum },{A_{\sum }}^2\right)\frac{\partial H_{{\sum }_{\epsilon }}\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}d\sigma +& \left(189\right)\end{array}$ ${\int }_{\sum }{\partial }_2\varphi \left(H_{\sum },{A_{\sum }}^2\right)\frac{\partial {A_{{\sum }_{\epsilon }}}^2\circ T_{\epsilon }}{\partial \epsilon }{❘}_{\epsilon =0}d\sigma +$ ${\int }_{\sum }\varphi \left(H_{\sum },{A_{\sum }}^2\right)\left({\mathrm{div}}^{\parallel }\left({\Theta }^{\parallel }\right)-H_{\sum }{\Theta }^{\perp }\right)d\sigma .$

Substituting the results of Propositions 4 and 5:

$\begin{array}{cc}\begin{array}{c}\frac{d}{d\epsilon }\Phi \left(T_{\epsilon }(\sum )\right){❘}_{\epsilon =0}={\int }_{\sum }{\partial }_1\varphi \left(H_{\sum },{A_{\sum }}^2\right)\left({\Delta }^{\parallel }{\Theta }^{\perp }+{A_{\sum }}^2{\Theta }^{\perp }\right)d\sigma +\\ {\int }_{\sum }{\partial }_1\varphi \left(H_{\sum },{A_{\sum }}^2\right)\langle {\nabla }^{\parallel }{H}_{\sum },{\Theta }^{\parallel }\rangle d\sigma +\\ 2{\int }_{\sum }{\partial }_2\varphi \left(H_{\sum },{A_{\sum }}^2\right)(A_{\sum }:\left({\nabla }^{\parallel }{\nabla }^{\parallel }{\Theta }^{\perp }\right)+\\ \mathrm{Tr}\left(A_{\sum }^3\right){\Theta }^{\perp })d\sigma +\\ 2{\int }_{\sum }{\partial }_2\varphi \left(H_{\sum },{A_{\sum }}^2\right)\\ {\sum }_{\mathrm{ijkl}}{h}^{\mathrm{ij}}{h^{\mathrm{kl}}\left[A_{\sum }\right]}_{\mathrm{kl}}{A}_{\sum }({\nabla }_E^{\parallel }{\left(,A\right)}_{\sum }\left(E_j,{\Theta }^{\parallel }\right)d\sigma +\\ {\int }_{\sum }\varphi \left(H_{\sum },{A_{\sum }}^2\right)\left({\mathrm{div}}^{\parallel }\left({\Theta }^{\parallel }\right)-H_{\sum }{\Theta }^{\perp }\right)d\sigma \\ ={\int }_{\sum }{\partial }_1\varphi \left(H_{\sum },{A_{\sum }}^2\right)\left({\Delta }^{\parallel }{\Theta }^{\perp }+{A_{\sum }}^2{\Theta }^{\perp }\right)d\sigma +\\ 2{\int }_{\sum }{\partial }_2\varphi \left(H_{\sum },{A_{\sum }}^2\right)(A_{\sum }:\left({\nabla }^{\parallel }{\nabla }^{\parallel }{\Theta }^{\perp }\right)+\\ \mathrm{Tr}\left(A_{\sum }^3\right){\Theta }^{\perp })d\sigma -\\ {\int }_{\sum }\varphi \left(H_{\sum },{A_{\sum }}^2\right)H_{\sum }{\Theta }^{\perp }d\sigma \end{array}& \left(190\right)\end{array}$

since all terms containing Θ^∥ can be shown to equal ∫_Σdiv^∥(ϕ(H_Σ,∥A_Σ∥²)Θ^∥)dσ which vanishes by the generalized Stokes' Theorem.

Shape Derivative of a Combined Distance-Normal Curvature Penalty

Let Ω be a domain in . The general form of a class of shape functions for penalizing geometric features of Ω, based on the signed distance function, as well as the normal vector field and the curvature invariants of the boundary surface ∂Ω. Specifically, consider the shape function:

(Ω):=∫_Ωϕ(χ,d_∂Ω(χ))dχ+∫_∂Ωψ(χ,d_∂Ω(χ)·n_∂Ω(χ)·H_∂Ω(x))dx  (191)

where ϕ and ψ are smooth functions of their arguments, and d_∂Ω, n_∂Ω, H_∂Ω are the signed distance function, the normal vector field, and the mean curvature of ∂Ω, respectively. Most of the geometric performance functions, e.g., surface area, Willmore energy for shape smoothing, thickness penalty, overhang angle penalty can be cast in this form.

The calculation of the shape derivative of Equation 191 breaks up into pieces thanks to the Chain Rule for shape differentiation. Let Ω_ε be a variation of Ω generated by the boundary velocity function Θ:∂Ω→, with tangential and normal parts denoted by Θ^∥ and Θ^⊥ respectively. Then,

$\begin{array}{cc}\frac{d}{d\epsilon }\left({\Omega }_{\epsilon }\right){❘}_{\epsilon =0}={\int }_{\Omega }{\partial }_1\varphi ·\Theta +{\partial }_2\varphi ·{\left[d_{\partial \Omega }\right]}^{·}+\varphi \mathrm{div}\left(\Theta \right)+{\int }_{\partial \Omega }{\partial }_1\psi ·\left({\Theta }^{\parallel }+n_{\partial \Theta }{\Theta }^{\perp }\right)+{\partial }_2\psi ·{\left[d_{\partial \Omega }\right]}^{·}+{\partial }_3\psi ·{\left[n_{\partial \Omega }\right]}^{·}+{\partial }_4\psi ·{\left[H_{\partial \Omega }\right]}^{·}+{\int }_{\partial \Omega }\psi \left({\mathrm{div}}^{\parallel }\left({\Theta }^{\parallel }\right)-H_{\partial \Omega }{\Theta }^{\perp }\right)& \left(192\right)\end{array}$

where ● denotes the so-called Lagrangian derivative and ∂_k denotes taking the partial derivative with respect to the k^th slot of the operand. The Lagrangian derivative of a function ƒ is defined as the derivative of ƒ along the variation induced by Θ. The Lagrangian derivative of a vector field is defined similarly.

Each of the Lagrangian-differentiated expressions above has a formula in terms of Θ^∥ and Θ^⊥ which have been derived above. Moreover, it can be shown that the dependence on Θ^∥ drops out of the final expression entirely (this result is known as the Hadamard-Zolésio Structure Theorem). This is shown by introducing the Eulerian derivative defined by ƒ′:=ƒ^●−∇ƒ·Θ for volumetric integrands, and ƒ′:=ƒ^●−∇^∥ƒ·Θ^∥ for surface integrands. Now thanks to Stokes' Theorem, all the explicit dependence on Θ^∥ disappears:

$\begin{array}{cc}\begin{array}{c}\frac{d}{d\epsilon }\left({\Omega }_{\epsilon }\right){❘}_{\epsilon =0}={\int }_{\Omega }\mathrm{div}\left(\varphi \Theta \right)+{\partial }_2\varphi ·{\left[d_{\partial \Omega }\right]}^{\prime }+{\int }_{\partial \Omega }{\partial }_1\psi ·n_{\partial \Omega }{\Theta }^{\perp }+\\ {\partial }_2\psi ·{\left[d_{\partial \Omega }\right]}^{\prime }+{\partial }_3\psi ·{\left[n_{\partial \Omega }\right]}^{\prime }+{\partial }_4\psi ·{\left[H_{\partial \Omega }\right]}^{\prime }+\\ {\int }_{\partial \Omega }{\mathrm{div}}^{\parallel }\left({\mathrm{\psi \Theta }}^{\parallel }\right)-\psi H_{\partial \Omega }{\Theta }^{\perp }\\ ={\int }_{\partial \Omega }{\mathrm{\varphi \Theta }}^{\perp }+{\int }_{\partial \Omega }\left({\partial }_1\psi ·n_{\partial \Omega }-\psi H_{\partial \Omega }\right){\Theta }^{\perp }+\\ {\int }_{\Omega }{\partial }_2\varphi ·{\left[d_{\partial \Omega }\right]}^{\prime }+{\int }_{\partial \Omega }{\partial }_2\psi ·{\left[d_{\partial \Omega }\right]}^{\prime }+\\ {\partial }_3\psi ·{\left[n_{\partial \Omega }\right]}^{\prime }+{\partial }_4\psi ·{\left[H_{\partial \Omega }\right]}^{\prime }\end{array}.& \left(193\right)\end{array}$

The remaining Eulerian-differentiated terms can be handled using the following equations derived-above. These are:

[d_∂δ]′=−Θ^⊥∘proj_∂δ  (194)

[n_∂δ]′=−∇^∥Θ^⊥  (195)

[H_∂δ]′=Δ^∥Θ^∥+∥A_∂δ∥²Θ^⊥  (196)

where proj_∂δ is the nearest-point projection onto ∂δ, and Δ^∥ is the Laplace-Beltrami operator of ∂δ, and ∥A_∂δ∥² is the squared Frobenius norm of the second fundamental form of ∂δ (a.k.a. the total curvature or sum of the squared principal curvatures). Substituting these into the Equation 193 yields the shape derivative formula.

FIG. 7C shows an example of a process of generating one or more portions of a three dimensional model of an object to be manufactured, using a stress constraint. A computer aided design program obtains 702 a design space for a modeled object, for which a corresponding physical structure will be additively manufactured, design criteria for the modeled object including at least one stress constraint, at least one in-use load case for the physical structure, and a specification of one or more materials from which the physical structure will be additively manufactured.

The specification includes information, e.g., dataset or a function of thickness, build angle, or both, to generate a corresponding strength for the material. The specification indicates multiple strength values for the one or more materials, that are each dependent upon one or both of a thickness of the physical structure to be built using the one or more materials, and a build angle. The object to be manufactured is additively manufactured using the one or more materials.

A generatively designed three dimensional shape of the modeled object is produced 704, including modifying both a geometry of the three dimensional shape and a topology of the three dimensional shape, in accordance with the design criteria, the at least one in-use load case, and the specification of one or more materials. The producing 704 includes varying, during object shape and topology modification of the modeled object, evaluation of the stress constraint at different locations on or in the modeled object in accordance with respective values from the multiple strength values. In some implementations, other constraints are used, e.g., a strain constraint or a displacement constraint. Each strength value corresponds to one or both of the thickness and the build angle at each one of the different locations. Note that the described Design Digital Twin techniques can be implemented with density based topology optimization, e.g., using the SIMP method, rather than a boundary based topology optimization, e.g., the level-set method as described in connection with FIG. 7A. In either case, the generatively designed three dimensional shape of the modeled object is provided 706 for use in additively manufacturing the physical structure.

In some implementations, each strength value corresponds to the thickness or both the thickness and the build angle at each one of the different locations. Further, the producing 704 includes: measuring the thickness at each one of the different locations; and computing a respective strength value based on at least the measured thickness at each one of the different locations. As described above with reference to FIG. 6B, measuring the thickness can include using a measure of thickness for the generatively designed three dimensional shape of the modeled object that is a combination of at least two distinct thickness measures. The two distinct measures can include a combination of at least (i) a first distance measure being a length within the modeled object of a ray cast in a negative normal direction from a surface point of the modeled object, and (ii) a second distance measure being a diameter of a largest sphere that touches the surface point of the modeled object and fits inside the modeled object as determined by checking discrete sampling locations defined on the sphere's surface.

In some implementations, modifying both the geometry and the topology of the three dimensional shape includes enforcing a design criterion that limits a minimum thickness of the generatively designed three dimensional shape of the modeled object, e.g., a minimum thickness being based on the critical fatigue crack length for the material, as described with reference to FIGS. 6A-6E. A minimum thickness can be a ratio of the critical fatigue crack length, such as a minimum thickness that is 10× the critical fatigue crack length.

Returning to FIG. 7A, an expected number of loading cycles is determined 712 for each of the one or more in-use load cases for the physical structure using the maximized stress or strain element and data relating fatigue strength to loading cycles. The data relating fatigue strength to loading cycles can be empirical, i.e., from a table of previously measured strengths for a corresponding material, e.g., as shown in FIG. 6D, or the data can be generated based on a function that generates, as output, a fatigue strength for an input number of loading cycles for a given load case of the object constructed using a corresponding material.

A fatigue safety factor inequality constraint is redefined for the modeled object based on a damage fraction calculated from the required number of loading cycles for the modeled object and the expected number of loading cycles for each of the one or more in-use load cases for the physical structure. Techniques for Handling Arbitrary Constraints, including inequality constraints can be applied, as described above with reference to FIGS. 4A-4H.

Shape change velocities for an implicit surface in a level-set representation of the three dimensional shape are computed 716 in accordance with at least the fatigue safety factor inequality constraint. The level-set representation is updated 718 using the shape change velocities to produce an updated version of the three dimensional shape of the modeled object. In some implementations, computing the shape change velocities includes computing at least one shape change velocity using a gradient determined from a shape derivative of one or both of the thickness and the build angle at each of the different locations, e.g., the shape derivative shown in Equation 193 and derived above.

Alternatively or in combination, i.e., where there are multiple constraints, shape change velocities can be computed using at least one shape change velocity with an amount determined from a shape derivative formula that approximates a shape derivative of one or both of the thickness and the build angle at each of the different locations. One example shape derivative formula is the volume shape derivative formula, described above with reference to Equation 105.

As described above with reference to Arbitrary Inequality Constraints and Equations 91-102, importance factors can be applied to inequality constraints to modify a volume fraction based inequality constraint using a shape derivative formula that includes the volume fraction based inequality constraint. This constraint can also be a minimum thickness inequality constraint described above with reference to Equations 123-127, or a stress based inequality constraints, as described above with reference to FIGS. 5A-6E. The constraints can be modified using an importance factor, which is adjusted based on whether or not one or more other constraints were violated in a prior iteration of the iteratively modifying. An example of the modification of the importance factor is multiplying with a violation multiplier.

The performing 708, the finding 710, the determining 712, the redefining 714, the computing 716 and the updating 718 is repeated until a check 720 determines that a predefined number of shape modification iterations have been performed or that the generatively designed three dimensional shape of the modeled object in the design space has converged to a stable solution for the one or more design criteria and the one or more in-use load cases, as described above with reference to FIG. 3A. Any constraint can be run with design digital twin as long as numerical simulation is run during the process, e.g., any constraint handled as described above with reference to FIGS. 4A-4H. Therefore, digital twin is not exclusive to damage tolerance fatigue. Similarly, any number of constraints could co-exist alongside design digital twin as it only affects the material model used in the numerical simulation.

Singularities and Disconnections

Stress constraints, including Von Mises stress constraints, can be essential for running generative design, but they are hard to implement due to few elements in a Finite Element Analysis model showing very high stresses. A probabilistic method is provided, below, that can mitigate or outright eliminate such high stresses. Although the examples below refer to stress constraints, these techniques can be applied to avoiding singularities and preventing disconnections during generative design for any constraint. Singularities can occur, for example due to sharp re-entrant corners or bad meshing. The process described below can be implemented to be performed automatically as part of a generative design process, e.g., including the techniques described above with reference to FIGS. 3A-7C.

FIG. 8A shows an example of a process of iteratively modifying a generatively designed three dimensional shape of a modeled object in the design space in accordance with the one or more design criteria, while avoiding excessive sudden changes and minimizing the likelihood of disconnections. A computer aided design program obtains 814 a design space for a modeled object, for which a corresponding physical structure will be manufactured, one or more design criteria for the modeled object, and one or more in-use load cases for the physical structure. A generatively designed three dimensional shape of the modeled object is iteratively modified in the design space in accordance with the one or more design criteria and the one or more in-use load cases for the physical structure, e.g., as described above with reference to FIGS. 3A, 4A, 5A, 6A, and 7A.

One solution for avoiding singularities is to replace the maximum stress by a simple percentile value σ^x in Equation 134, modified below.

$\begin{array}{cc}\begin{array}{c}\max \\ \Omega \end{array}\sigma \leftarrow {\sigma }^{x}x\in \left\{0.9,0.95,0.99,0.999\mathrm{etc}\right\}& \left(197\right)\end{array}$

where σ^x denotes the stress of the element in the x^th percentile when elements are sorted in increasing order of σ. However, this results in the maximum stress in the domain max no-being too sensitive to the value of x. Instead, the percentile x is converted into a standard normal deviate z using the inverse of the error function

$\begin{array}{cc}z={\mathrm{erf}}^{-1}\left(x\right)=\frac{1}{2}\sqrt{\pi }\left(x+\frac{\pi }{12}{x}^3+\frac{7{\pi }^2}{480}{x}^5+\frac{127{\pi }^3}{40320}{x}^7+\dots \right)& \left(198\right)\end{array}$

Next, the maximum stress is computed from

$\begin{array}{cc}\begin{array}{c}\max \\ \Omega \end{array}\sigma \leftarrow \mu \left(\sigma \right)+z\chi \left(\sigma \right)& \left(199\right)\end{array}$

where μ(σ) and χ(σ) denote the mean and standard deviation of the stress distribution, both of which can be computed according to any conventional technique.

TABLE 1, below, shows stress singularities for an object that are avoided based on the techniques described in this section:

TABLE 1
		Maximum Stress Based
	Maximum Stress based	on Standard Deviation
Percentile x	on Percentile σx	μ(σ) + erf−1(x) χ(σ)
.95	133.87	123.597
.99	224.08	144.41
.999	500.315	168.56
.9999	815.34	188.84

The percentile based method brings stability to the oscillatory nature of global stress, by smoothing the maximum stress at a point based on the standard deviation of stress values for points on the object in the same percentile.

High Velocity Smoothing

The previous section described regulating the maximum value of constraint values in the domain, i.e., max_Ωgⁱ. The same is true for the constraint derivative dg_t/d_Ω, denoted by dg_t for clarity.

FIG. 8B is a graphical representation of an example of geometry disconnection during optimization. Object 804 is shown at different stage of optimization. Typically, geometric ports can disconnect (shown at 806) from the main body of the design where the load path is taken by the Ersatz material. Having a small region with very high advection velocity causes disconnections and localized shape changes. A combination of velocity clamping and high velocity smoothing can be used to prevent this.

Normalization: all shape derivatives can be normalized such that the magnitude of the maximum value on the surface of the current design F approximates the voxel size Δs

$\begin{array}{cc}{\left(d{g}_t\left(x_K\right)\right)}_{\mathrm{normalised}}=\Delta s\frac{{\mathrm{dg}}_t\left(x_k\right)}{❘{\max }_{x_k\in \Gamma }{\mathrm{dg}}_t\left(x_k\right)❘}& \left(200\right)\end{array}$

This implies an advection time of T=1/Δs is sufficient to prevent the geometry from advecting beyond one voxel. An easy way to achieve this is through velocity clamping.

Velocity clamping: the first step in velocity clamping is computing a reference value for the shape derivative following the percentile approach described above with reference to Equations 197-199.

(dg_t)_ref=μ(dg_t)+erf⁻¹(x)χ(dg_t)  (201)

where μ and χ denote the mean and standard deviation of the shape derivative on the surface and x∈{0.9,0.95,0.99,0.999 etc.} is a user-given percentile. Note that the percentile value x could be much higher for well behaved shape derivatives (without excessive highs) such as strain energy and lower for widely fluctuating shape derivatives, such as stress.

Next, for all grid points x_k in the inner narrow band, i.e., a narrow band width equal to w_nbΔs, velocity values higher than the reference value (dg_t)_ref are clamped while preserving the sign β as follows:

$\begin{array}{cc}{\left({\mathrm{dg}}_t\left(x_k\right)\right)}_{\mathrm{clamped}}=\mathrm{\beta min}\left[❘{\mathrm{dg}}_t\left(x_k\right)❘,❘({\mathrm{dg}}_{t_{\mathrm{ref}}}❘\right]-w_{nb}\Delta s\le \psi \left(x_k\right)\le 0\beta =\frac{{\mathrm{dg}}_t\left(x_k\right)}{\left[{\mathrm{dg}}_t\left(x_k\right)\right]}& \left(202\right)\end{array}$

where ψ(x_l) denotes the level-set of the grid point. However, this results in a sudden change in the velocity profile. Velocity clamping should be done in such a manner that there is a smooth transition in the velocity profile.

High velocity smoothing: The objective in high velocity smoothing is to smooth all high velocities (|dg_t|>|(dg_t)_ref|) such that there is a gradual smoothing of all values above the reference velocity computed from Equation 201. A user given parameter ρ, e.g., ρ=0.85, is used to scale all velocity values less than the reference velocity

$\begin{array}{cc}{\left({\mathrm{dg}}_t\left(x_k\right)\right)}_{\mathrm{smooth}}=\mathrm{\rho \beta }❘{\mathrm{dg}}_t\left(x_k\right)❘\frac{❘{\mathrm{dg}}_t\left(x_n\right)❘}{❘{\left({\mathrm{dg}}_t\right)}_{\mathrm{ref}}❘}k<k_{\mathrm{ref}}& \left(203\right)\end{array}$

where k_ref denotes the position index of the grid point once they are sorted in order of the shape derivative.

The smoothed velocity of all grid points with velocity greater than the reference velocity are found by fitting a cubic polynomial, as follows:

(dg_t(x_k))_smooth=aξ³+bξ²+cξ+d k_ref≤k≤n  (204)

where ξ denotes the shifted index τ=k−k_ref. The unknown coefficients {a, b, c, d} can be found solving the linear system of equations obtained by settings the appropriate boundary conditions at ξ=0 and ξ=ξ_max=n−k_ref

$\begin{array}{cc}\left[\begin{array}{cccc}{\left(0\right)}^3& {\left(0\right)}^2& {\left(0\right)}^1& 1\\ {\left({\xi }_{\max }\right)}^3& {\left({\xi }_{\max }\right)}^2& {\left({\xi }_{\max }\right)}^1& 1\\ 3{\left(0\right)}^2& 2{\left(0\right)}^1& 1& 0\\ 3{\left({\xi }_{\max }\right)}^2& 2{\left({\xi }_{\max }\right)}^1& 1& 0\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\\ d\end{array}\right]=\left[\begin{array}{c}{\mathrm{dg}}_t\left(x_{k_{\mathrm{ref}}}\right)\\ {\mathrm{dg}}_t\left(x_n\right)\\ \nabla {\left({\mathrm{dg}}_t\left(x_{k_{\mathrm{ref}}}\right)\right)}_{\mathrm{smooth}}\\ 0\end{array}\right].& \left(205\right)\end{array}$

The gradient of the smooth velocity at ξ=0 can be obtained as

$\begin{array}{cc}\nabla {\left({\mathrm{dg}}_t\left(x_{k_{\mathrm{ref}}}\right)\right)}_{\mathrm{smooth}}=\rho \nabla \left({\mathrm{dg}}_t\left(x_{k_{\mathrm{ref}}}\right)\right)\frac{❘{\mathrm{dg}}_t\left(x_n\right)❘}{❘{\left({\mathrm{dg}}_t\right)}_{\mathrm{ref}}❘}.& \left(206\right)\end{array}$

Finally, the smoothed normalized velocity is given by

$\begin{array}{cc}{\left({\mathrm{dg}}_t\left({\chi }_k\right)\right)}_{\mathrm{smooth}\mathrm{normalized}}=\mathrm{\Delta s}\frac{{\left({\mathrm{dg}}_t\left(x_k\right)\right)}_{\mathrm{smooth}}}{❘{{\max }_{x_k\in \Gamma }\left({\mathrm{dg}}_t\left(x_k\right)\right)}_{\mathrm{smooth}}❘}.& \left(207\right)\end{array}$

Returning to FIG. 8A, within the iteratively modifying, a numerical simulation of the modeled object in accordance with a current version of the three dimensional shape and the one or more in-use load cases is performed 816 to produce a current numerical assessment of a physical response, e.g., a structural response, of the modeled object. Shape change velocities are computed 818 for an implicit surface in a level-set representation of the three dimensional shape, e.g., as described above with reference to FIG. 3A.

Then, the shape change velocities are changed 820 in accordance with a polynomial function that has been fit to at least a portion of the shape change velocities above a reference velocity. For example, the polynomial function can be a cubic polynomial, as described above with reference to Equation 204. Although the example given is a cubic polynomial, other different-order polynomials can be used.

The level-set representation is updated 822 using the shape change velocities to produce an updated version of the three dimensional shape of the modeled object, e.g., as described above with reference to FIG. 3A. In some implementations, the reference velocity is set based on a mean and a standard deviation of a shape derivative on the implicit surface, e.g., as described above.

The performing 816, the computing 818, the changing 820, and the updating 822 is repeated in the iteratively modifying until check 830 determines that a predefined number of shape modification iterations have been performed or that the generatively designed three dimensional shape of the modeled object in the design space has converged to a stable solution for the one or more design criteria and the one or more in-use load cases, e.g., as described above with reference to FIG. 3A. Finally, the generatively designed three dimensional shape of the modeled object can be provided for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

Backup and Restore

Another technique for preventing geometric disconnections is backing up of all critical data, including the geometric level-set at each iteration, and restoring the data when the result of advection is not desirable. Thus, in some implementations, after the updating 822 but before the repeating, excessive changes can be checked 824 for that were made during the updating. In those cases, the current version of the three dimensional shape can be set 826 as the updated version of the three dimensional shape for a next iteration to undue the excessive change, i.e., the excessive change is undone. Then, shape changes for the next iteration of the iteratively modifying can be slowed 828, as described below. Slowing the changes for the next iteration can include reducing a target volume change for the generatively designed three dimensional shape of the modeled object for the next iteration of the iteratively modifying.

To prevent a repeat of the same undesirable outcome, the volume change Δv_t computed for each iteration t, as described above with reference to FIGS. 4A-4H is modified as follows and using a multiplier β which slows down the convergence by applying a smaller volume change:

Δv_t←βΔv_t.  (208)

This multiplier can be initiated as β=1 at the start and updated in every iteration depending on the outcome of each advection using a fixed increment Δβ(>0):

β←max(0,β−Δβ) if advection undesirable  (209)

β←min(1,β−Δβ) otherwise  (210)

The next task is to classify the outcome of each advection as desirable/undesirable. One way to achieve this is to monitor the change in the Lagrangian, i.e., as described above with reference to Equations 15-25. This can be implemented in a more granular manner by imposing limits on allowable relative changes for objectives. The relative change in an objective/constraint for a given iteration t is defined as

$\begin{array}{cc}{\mathrm{\Delta 𝒥}}_t=❘\frac{{𝒥}_t-{𝒥}_{t-1}}{J_{\mathrm{ref},t}}❘& \left(211\right)\end{array}$ $\begin{array}{cc}{\mathrm{\Delta g}}_t^j=p_t^j❘e_t^j-e_{t-1}^j❘& \left(212\right)\end{array}$

where the constraint error e_t^j is computed using a PID stabilized version of Equation 95 with the importance factor term p_t^j computed using techniques described above with reference to Arbitrary Inequality Constraints and FIGS. 4A-4H. The reference value of the objective, denoted _ref,t, is given by

$\begin{array}{cc}{𝒥}_{\mathrm{ref},t}=\begin{array}{c}t\\ \max ❘i_{}❘\\ i=0\end{array}.& \left(213\right)\end{array}$

The maximum allowable change in objectives is limited to for objective decrease and increase, respectively. The maximum allowable change in constraints is limited to for positive (good), negative (bad) changes in constraints s, respectively. The sign of the constraint inequality needs to be checked to determine the appropriate limit.

Thus, in some implementations, before the performing 816 is done, elements that are generated from the current version of the three dimensional shape for the numerical simulation, but that are partially but not entirely within the implicit surface with Dirichlet boundary conditions are identified and removed 815 before performing the numerical simulation. The element and nodes refer to how the design geometry is represented during the performing 816. In finite element simulation, the domain is replaced by a collection of elements. Each element is defined by a set of nodes along the boundary of the element and, in some cases, inside the element. Generally, the more elements and nodes in the model, the higher the accuracy of the simulation. The set of techniques for creating these nodes and elements is referred to as meshing.

During optimization, the shape of the current domain repeatedly changes. The model can be re-meshed at every iteration for the finite element model to accurately represent the current domain. Since this is computationally expensive, some generative design processes deactivate the elements that lie outside the current domain at each iteration. The representation of the current design can be improved in the numerical simulation model by removing elements not entirely within the implicit surface. The overall accuracy can be improved by considering such cut elements as having partial stiffness.

Further, in some implementations, the checking 824 for the excessive change can include comparing a change in the one or more design criteria resulting from the updating with a predefined limit on an amount of change allowed for the one or more design criteria in a single iteration of the iteratively modifying.

Ersatz Materials

A common practice in the SIMP method of topology optimization is to use variable density materials or ersatz materials. Elements with density ρ=1 are deemed to be inside the domain whereas those with ∈<ρ<1 are on the boundary with ∈ (typically set to a small value such as 0:001) denoting the density of materials outside the design. During finite element analysis, the stiffness of each element is multiplied by the density along with a penalty factor p to penalize intermediate densities

K←ρ ^p K.  (214)

The level-set method often uses only two states (ρ∈{∈,1}) as the boundary and therefore is more clearly defined. The penalty factor is not required as there are no intermediate densities, i.e., K←ρ^pK. However, this sometimes leads to disconnections as elements outside the domain with ρ=∈ can support the load path.

Additionally, it has been observed that having material with ρ=∈ severely affects the buckling safety factor predicted from the finite element model. Ersatz material can be removed completely from the finite element model where there is a significant impact on the predicted buckling factor. This is achieved by grouping all elements with p=E at each iteration and removing all such groups that are not connected to any nodes with Dirichlet boundary conditions.

As noted above, in some implementations of the process shown in FIG. 8A, before the performing 816, elements generated from the current version of the three dimensional shape for the numerical simulation are identified 815 that are partially but not entirely within the implicit surface. Then, densities of the identified elements can be set equal to respective volume fractions of the identified elements, where a volume fraction of an identified element is a fractional amount of the identified element that falls within the implicit surface. In some implementations, the performing 816 can include penalizing stiffness of the identified elements according to the densities, as described below.

As another solution to severe disconnection problems, a tetrahedral-cutting algorithm can be applied, which reinstates the use of ersatz materials for elements on the boundary of the current domain. Essentially, the stiffness of such elements are better approximated by setting the density equal to the volume fraction of the element inside the domain. The stiffness is then penalized according to the density, as in the SIMP method, with a penalty factor of p=1:

$\begin{array}{cc}\rho \left(e_i\right)=\frac{v\left(e_i\cap \Omega \right)}{v\left(e_i\right)}K\left(e_i\right)\leftarrow \rho \left(e_i\right)K\left(e_i\right).& \left(215\right)\end{array}$

Computing the volume of each element inside the domain is not trivial as some elements with valid overlap with the domain may have all nodes outside the elements.

FIG. 8C is a graphical representation of an example of a geometry 808 with simulation elements 810 classified based on intersection of the elements 810 with the geometry 808. The problem is illustrated below with respect to elements 810 e₅⁰, e₆⁰, e₇⁰. Such elements can be subdivided 812 until nodes are on either side of the 0^th iso-contour, i.e.,

$\begin{array}{cc}\begin{array}{c}\min \\ x_i\in e\end{array}\psi \left(x_i\right)\begin{array}{c}\max \\ x_i\in e\end{array}\psi \left(x_i\right)<0& \left(216\right)\end{array}$

where x_i denotes the nodal coordinates of element e. Note that certain elements (e.g., element 810 e₈⁰) may have to be subdivided multiple times until they satisfy the condition of Equation 216. The required depth of subdivision l can be determined by recursively subdividing the element and comparing the total volume inside the domain

$\begin{array}{cc}V\left(e_i^l\bigcap \Omega \right)={\sum }_{j\in C\left(e_i^l\right)}V\left(e_j^{l+1}\bigcap \Omega \right)& \left(217\right)\end{array}$

where C(e_i^l) denotes the indices of the children of element e_i^l. The recursive subdivision can be stopped when the above condition is true. In some implementations, the voxel size Δs is set to half the average edge length of solid elements. Given the smallest feature size≈Δs, only 1 level of subdivision is typically sufficient.

The volume fraction of a tetrahedron element can be cut by the domain. When computing the intersection point of an element edge with the 0^th iso-contour, linear interpolation can be used

$\begin{array}{cc}{\chi }_1^1=\frac{\psi \left(x_1^0\right)x_2^0-\psi \left(x_2^0\right)x_1^0}{\psi \left(x_1^0\right)-\psi \left(x_2^0\right)}& \left(218\right)\end{array}$

FIG. 9 is a schematic diagram of a data processing system including a data processing apparatus 900, which can be programmed as a client or as a server. The data processing apparatus 900 is connected with one or more computers 990 through a network 980. While only one computer is shown in FIG. 9 as the data processing apparatus 900, multiple computers can be used. The data processing apparatus 900 includes various software modules, which can be distributed between an applications layer and an operating system. These can include executable and/or interpretable software programs or libraries, including tools and services of one or more 3D modeling programs 904 that implement the systems and techniques described above. Thus, the 3D modeling program(s) 904 can be CAD program(s) 904 and can implement one or more generative design processes (e.g., using level-set based method(s) for generative design) for topology optimization and numerical simulation operations (finite element analysis (FEA) or other). Further, the program(s) 904 can potentially implement manufacturing control operations (e.g., generating and/or applying toolpath specifications to effect manufacturing of designed objects). The number of software modules used can vary from one implementation to another. Moreover, the software modules can be distributed on one or more data processing apparatus connected by one or more computer networks or other suitable communication networks.

The data processing apparatus 900 also includes hardware or firmware devices including one or more processors 912, one or more additional devices 914, a computer readable medium 916, a communication interface 918, and one or more user interface devices 920. Each processor 912 is capable of processing instructions for execution within the data processing apparatus 900. In some implementations, the processor 912 is a single or multi-threaded processor. Each processor 912 is capable of processing instructions stored on the computer readable medium 916 or on a storage device such as one of the additional devices 914. The data processing apparatus 900 uses the communication interface 919 to communicate with one or more computers 990, for example, over the network 980. Examples of user interface devices 920 include a display, a camera, a speaker, a microphone, a tactile feedback device, a keyboard, a mouse, and VR and/or AR equipment. The data processing apparatus 900 can store instructions that implement operations associated with the program(s) described above, for example, on the computer readable medium 916 or one or more additional devices 914, for example, one or more of a hard disk device, an optical disk device, a tape device, and a solid state memory device.

Embodiments of the subject matter and the functional operations described in this specification can be implemented in digital electronic circuitry, or in computer software, firmware, or hardware, including the structures disclosed in this specification and their structural equivalents, or in combinations of one or more of them. Embodiments of the subject matter described in this specification can be implemented using one or more modules of computer program instructions encoded on a non-transitory computer-readable medium for execution by, or to control the operation of, data processing apparatus. The computer-readable medium can be a manufactured product, such as hard drive in a computer system or an optical disc sold through retail channels, or an embedded system. The computer-readable medium can be acquired separately and later encoded with the one or more modules of computer program instructions, e.g., after delivery of the one or more modules of computer program instructions over a wired or wireless network. The computer-readable medium can be a machine-readable storage device, a machine-readable storage substrate, a memory device, or a combination of one or more of them.

The term “data processing apparatus” encompasses all apparatus, devices, and machines for processing data, including by way of example a programmable processor, a computer, or multiple processors or computers. The apparatus can include, in addition to hardware, code that creates an execution environment for the computer program in question, e.g., code that constitutes processor firmware, a protocol stack, a database management system, an operating system, a runtime environment, or a combination of one or more of them. In addition, the apparatus can employ various different computing model infrastructures, such as web services, distributed computing and grid computing infrastructures.

A computer program (also known as a program, software, software application, script, or code) can be written in any suitable form of programming language, including compiled or interpreted languages, declarative or procedural languages, and it can be deployed in any suitable form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program does not necessarily correspond to a file in a file system. A program can be stored in a portion of a file that holds other programs or data (e.g., one or more scripts stored in a markup language document), in a single file dedicated to the program in question, or in multiple coordinated files (e.g., files that store one or more modules, sub-programs, or portions of code). A computer program can be deployed to be executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a communication network.

The processes and logic flows described in this specification can be performed by one or more programmable processors executing one or more computer programs to perform functions by operating on input data and generating output. The processes and logic flows can also be performed by, and apparatus can also be implemented as, special purpose logic circuitry, e.g., an FPGA (field programmable gate array) or an ASIC (application-specific integrated circuit).

Processors suitable for the execution of a computer program include, by way of example, both general and special purpose microprocessors, and any one or more processors of any kind of digital computer. Generally, a processor will receive instructions and data from a read-only memory or a random access memory or both. The essential elements of a computer are a processor for performing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to receive data from or transfer data to, or both, one or more mass storage devices for storing data, e.g., magnetic, magneto-optical disks, or optical disks. However, a computer need not have such devices. Moreover, a computer can be embedded in another device, e.g., a mobile telephone, a personal digital assistant (PDA), a mobile audio or video player, a game console, a Global Positioning System (GPS) receiver, or a portable storage device (e.g., a universal serial bus (USB) flash drive), to name just a few. Devices suitable for storing computer program instructions and data include all forms of non-volatile memory, media and memory devices, including by way of example semiconductor memory devices, e.g., EPROM (Erasable Programmable Read-Only Memory), EEPROM (Electrically Erasable Programmable Read-Only Memory), and flash memory devices; magnetic disks, e.g., internal hard disks or removable disks; magneto-optical disks; and CD-ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry.

To provide for interaction with a user, embodiments of the subject matter described in this specification can be implemented on a computer having a display device, e.g., an LCD (liquid crystal display) display device, an OLED (organic light emitting diode) display device, or another monitor, for displaying information to the user, and a keyboard and a pointing device, e.g., a mouse or a trackball, by which the user can provide input to the computer. Other kinds of devices can be used to provide for interaction with a user as well; for example, feedback provided to the user can be any suitable form of sensory feedback, e.g., visual feedback, auditory feedback, or tactile feedback; and input from the user can be received in any suitable form, including acoustic, speech, or tactile input.

The computing system can include clients and servers. A client and server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other. Embodiments of the subject matter described in this specification can be implemented in a computing system that includes a back-end component, e.g., as a data server, or that includes a middleware component, e.g., an application server, or that includes a front-end component, e.g., a client computer having a graphical user interface or a browser user interface through which a user can interact with an implementation of the subject matter described is this specification, or any combination of one or more such back-end, middleware, or front-end components. The components of the system can be interconnected by any suitable form or medium of digital data communication, e.g., a communication network. Examples of communication networks include a local area network (“LAN”) and a wide area network (“WAN”), an inter-network (e.g., the Internet), and peer-to-peer networks (e.g., ad hoc peer-to-peer networks).

While this specification contains many implementation details, these should not be construed as limitations on the scope of what is being or may be claimed, but rather as descriptions of features specific to particular embodiments of the disclosed subject matter. Certain features that are described in this specification in the context of separate embodiments can also be implemented in combination in a single embodiment. Conversely, various features that are described in the context of a single embodiment can also be implemented in multiple embodiments separately or in any suitable subcombination. Moreover, although features may be described above as acting in certain combinations and even initially claimed as such, one or more features from a claimed combination can in some cases be excised from the combination, and the claimed combination may be directed to a subcombination or variation of a subcombination.

Similarly, while operations are depicted in the drawings in a particular order, this should not be understood as requiring that such operations be performed in the particular order shown or in sequential order, or that all illustrated operations be performed, to achieve desirable results. In certain circumstances, multitasking and parallel processing may be advantageous. Moreover, the separation of various system components in the embodiments described above should not be understood as requiring such separation in all embodiments, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products.

Thus, particular embodiments of the invention have been described. Other embodiments are within the scope of the following claims. In addition, actions recited in the claims can be performed in a different order and still achieve desirable results.

Claims

1. (canceled)

2. A method comprising: obtaining, by a computer program, one or more load cases and one or more design criteria for a modeled object;iteratively modifying, by the computer program, a three dimensional shape of the modeled object in accordance with the one or more design criteria and the one or more load cases, wherein the iteratively modifying comprises regulating shape change velocities for an implicit surface representation of the three dimensional shape that exceed a reference velocity, wherein the reference velocity is set based on a mean and a standard deviation of a shape derivative on the implicit surface; andproviding, by the computer program, the three dimensional shape of the modeled object for use in manufacturing a physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

3. The method of claim 2, wherein regulating the shape change velocities comprises changing the shape change velocities in accordance with a polynomial function.

4. The method of claim 3, comprising generating the polynomial function by fitting a cubic polynomial to velocity values above the reference velocity.

5. The method of claim 2, wherein the iteratively modifying comprises slowing shape changes for a next iteration of the iteratively modifying in response to an excessive change having been made and undone in a current iteration of the iteratively modifying.

6. The method of claim 5, comprising detecting the excessive change based on a predefined limit on an amount of change allowed for the one or more design criteria in a single iteration of the iteratively modifying.

7. The method of claim 5, wherein slowing the shape changes for the next iteration comprises reducing a target volume change for the three dimensional shape of the modeled object for the next iteration of the iteratively modifying.

8. The method of claim 2, wherein the iteratively modifying comprises: identifying elements that are partially but not entirely within the implicit surface; andremoving any of the identified elements that are not connected to any nodes of a current version of the three dimensional shape with Dirichlet boundary conditions.

9. The method of claim 2, wherein the iteratively modifying comprises: identifying elements that are partially but not entirely within the implicit surface;setting densities of the identified elements equal to respective fractional amounts of each identified element that falls within the implicit surface; andpenalizing stiffness of the identified elements according to the densities during numerical simulation.

10. The method of claim 2, wherein the iteratively modifying comprises computing shape change velocities using an amount determined from a shape derivative formula that approximates a shape derivative for multiple locations on or in the modeled object.

11. The method of claim 10, wherein the shape derivative formula comprises a volume fraction or a stress based inequality constraint that is modified using an importance factor, which is adjusted based on whether or not one or more other constraints were violated in a prior iteration of the iteratively modifying.

12. The method of claim 11, wherein the iteratively modifying comprises: adjusting a target value of the inequality constraint between an initial target value and a final target value across multiple iterations of the iteratively modifying; andusing a proportional-integral-derivative controller to stabilize changes made in the amount determined from the shape derivative formula as the target value is adjusted across the multiple iterations.

13. The method of claim 11, wherein the iteratively modifying comprises using a proportional-integral-derivative controller to adjust a total contribution of the amount determined from the shape derivative formula to shape change velocities used during the iteratively modifying.

14. A system comprising: a first computer communicatively coupled with a network, the first computer comprising a first non-transitory storage medium having instructions of a computer program stored thereon; anda second computer communicatively coupled with the network, the second computer comprising a second non-transitory storage medium having instructions of the computer program stored thereon;wherein the instructions of the computer program are configured to cause the first and second computers to operate cooperatively to obtain one or more load cases and one or more design criteria for a modeled object,iteratively modify a three dimensional shape of the modeled object in accordance with the one or more design criteria and the one or more load cases, wherein the iterative modification comprises regulating shape change velocities for an implicit surface representation of the three dimensional shape that exceed a reference velocity, wherein the reference velocity is set based on a mean and a standard deviation of a shape derivative on the implicit surface, andprovide the three dimensional shape of the modeled object for use in manufacturing a physical structure corresponding to the modeled object using one or more computer-controlled additive manufacturing systems.

15. The system of claim 14, wherein regulating the shape change velocities comprises changing the shape change velocities in accordance with a polynomial function.

16. The system of claim 15, wherein the instructions of the computer program are configured to cause the first and second computers to operate cooperatively to generate the polynomial function by fitting a cubic polynomial to velocity values above the reference velocity.

17. The system of claim 14, wherein the iterative modification comprises slowing shape changes for a next iteration of the iterative modification in response to an excessive change having been made and undone in a current iteration of the iterative modification.

18. The system of claim 17, wherein the instructions of the computer program are configured to cause the first and second computers to operate cooperatively to detect the excessive change based on a predefined limit on an amount of change allowed for the one or more design criteria in a single iteration of the iterative modification.

19. The system of claim 17, wherein slowing the shape changes for the next iteration comprises reducing a target volume change for the three dimensional shape of the modeled object for the next iteration of the iterative modification.

20. The system of claim 14, wherein the iterative modification comprises: identifying elements that are partially but not entirely within the implicit surface; andremoving any of the identified elements that are not connected to any nodes of a current version of the three dimensional shape with Dirichlet boundary conditions.

21. The system of claim 14, wherein the iterative modification comprises: identifying elements that are partially but not entirely within the implicit surface;setting densities of the identified elements equal to respective fractional amounts of each identified element that falls within the implicit surface; andpenalizing stiffness of the identified elements according to the densities during numerical simulation.

22. The system of claim 14, wherein the iterative modification comprises computing shape change velocities using an amount determined from a shape derivative formula that approximates a shape derivative for multiple locations on or in the modeled object.

23. The system of claim 22, wherein the shape derivative formula comprises a volume fraction or a stress based inequality constraint that is modified using an importance factor, which is adjusted based on whether or not one or more other constraints were violated in a prior iteration of the iteratively modifying.

24. The system of claim 23, wherein the iterative modification comprises: adjusting a target value of the inequality constraint between an initial target value and a final target value across multiple iterations of the iteratively modification; andusing a proportional-integral-derivative controller to stabilize changes made in the amount determined from the shape derivative formula as the target value is adjusted across the multiple iterations.

25. The system of claim 24, wherein the iterative modification comprises using a proportional-integral-derivative controller to adjust a total contribution of the amount determined from the shape derivative formula to shape change velocities used during the iteratively modification.