Computer aided design and manufacturing of generative design shape optimization with limited fatigue damage

Abstract

Methods, systems, and apparatus, including computer programs encoded on media, for computer-aided design of physical structures using generative design processes. One method includes obtaining, by a computer-aided design program, a design space for a modeled object, one or more design criteria for the modeled object, one or more in-use loading conditions, and a critical fatigue crack length for a material used to manufacture the physical structure; iteratively modifying a generatively designed three-dimensional shape of the modeled object in the design space as a function of the critical fatigue crack length for the material, wherein the iteratively modifying includes enforcing a design criterion that limits a minimum thickness of the generatively designed three-dimensional shape, the minimum thickness being based on the critical fatigue crack length for the material.

Description

Background Technology

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

[0002] 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 manufacturing of the physical structure of these objects, for example, using computer numerical control (CNC) manufacturing techniques. Typically, CAD software uses the Boundary Representation (B-Rep) format to store a 3D representation of the geometry of the object being modeled. A B-Rep model is a set of connected surface elements that specify the boundaries between the solid and non-solid parts of the 3D object being modeled. In a B-Rep model (often referred to as B-Rep), the geometry is stored in the computer using smooth and precise mathematical surfaces, which contrasts with the discrete and approximate surfaces of mesh models, which can be difficult to use in CAD programs.

[0003] CAD programs have been integrated with subtractive manufacturing systems and technologies. Subtractive manufacturing refers to any manufacturing process that produces a 3D object from stock material (typically a "blank" or "workpiece" larger than the 3D object) by removing portions of the stock material. Such manufacturing processes typically involve the use of multiple CNC machine tool cutting tools in a series of operations, starting with roughing operations, optionally semi-finishing operations, and finishing operations. Besides CNC machining, other subtractive manufacturing techniques include electrical discharge machining, chemical machining, waterjet machining, and more. In contrast, additive manufacturing (also known as solid freeform manufacturing or 3D printing) refers to any manufacturing process that builds a 3D object from raw materials (typically powders, liquids, suspensions, or molten solids) in a series of layers or cross-sections. Examples of additive manufacturing include fused filament manufacturing (FFF) and selective laser sintering (SLS). Other manufacturing techniques for building 3D objects from raw materials include casting and forging (both hot and cold forging).

[0004] Furthermore, CAD software has been designed to automatically generate 3D geometry (generative design) for a single part or one or more parts in a larger system of parts to be manufactured. This automated generation of 3D geometry is often constrained by a design space specified by the user of the CAD software, and the generation is typically subject to design goals and constraints, which can be defined by the user or another party and imported into the CAD software. Design goals (such as minimizing waste material or the weight of the designed part) can be used to drive the geometry generation process toward better design. Design constraints can include both structural integrity constraints for individual parts (i.e., requirements that the part should not fail under expected structural loads during its service life) and physical constraints imposed by the larger system (i.e., requirements that the part should not interfere with another part in the system during its service life). Examples of design constraints include maximum mass, maximum deflection under load, maximum stress, etc.

[0005] Some CAD software already includes tools for facilitating 3D geometry enhancement using lattices and skins of various sizes, thicknesses, and densities. Latticees consist of beams or struts connected to each other at joints or directly to solid parts, and skins are shell structures that cover or enclose the lattice. Such tools allow 3D parts to be redesigned to reduce weight while still maintaining desired performance characteristics (e.g., stiffness and flexibility). These software tools already utilize various types of lattice topologies that can be used to generate manufacturable lattice structures.

[0006] Furthermore, the inputs to the generative design process may include a set of input solids (B-Rep inputs) specifying the boundary conditions for the generative design process; however, many modern generative design solvers do not operate directly on the exact surface boundary representations of their input solids. Instead, sampling the B-Rep and replacing it with volumetric representations such as level sets or tetrahedral or hexahedral meshes is significantly more convenient and efficient for physical simulations and material synthesis computed by the solver. The set of input solids may include “reserved volumes,” which should always exist in the design and represent interfaces of other parts or locations in the system on which boundary conditions (e.g., mechanical loads and constraints) should be applied. Other regions on which geometry should or should not be generated can also be provided in a similar manner, such as input solids defining “barrier volumes,” which represent regions on which new geometry should not be generated. Summary of the Invention

[0007] This specification describes techniques involving computer-aided design of physical structures using generative design processes, wherein a three-dimensional (3D) model of the physical structure can be generated with fatigue damage of finite size, and wherein the resulting physical structure can then be manufactured using additive manufacturing, subtractive manufacturing, and / or other manufacturing systems and techniques.

[0008] Generally, one or more aspects of the subject matter described in this specification may be embodied in one or more methods, including: obtaining, by a computer-aided design program, a design space of a modeling object upon which the manufacture of a corresponding physical structure will be based, one or more design criteria of the modeling object, one or more in-service load conditions of the physical structure, and the critical fatigue crack length of the material used to manufacture the physical structure; iteratively modifying, by the computer-aided design program, the three-dimensional shape of a generative design of the modeling object in the design space according to the one or more design criteria, the one or more in-service load conditions of the physical structure, and the critical fatigue crack length of the material, wherein the iterative modification includes: enforcing a design criterion that limits the minimum thickness of the three-dimensional shape of the generative design of the modeling object, the minimum thickness being based on the critical fatigue crack length of the material; and providing, by the computer-aided design program, the three-dimensional shape of the generative design of the modeling object for manufacture of the physical structure corresponding to the modeling object using one or more computer-controlled manufacturing systems.

[0009] The one or more design criteria may include the required number of load cycles for the modeled object under each of the one or more in-service load conditions of the physical structure, and the acquisition may include obtaining one or more specifications of the material used to manufacture the physical structure, the one or more specifications including data relating fatigue strength to load cycles; and the iterative modification may include: performing numerical simulations of the modeled object based on the current form of the three-dimensional shape and the one or more in-service load conditions to produce a current numerical evaluation of the physical response of the modeled object; and obtaining a current numerical evaluation of the physical response of the modeled object for at least one of the one or more in-service load conditions of the physical structure. The numerical evaluation identifies the maximum stress or strain element; uses the maximum stress or strain element and the data relating fatigue strength to load cycles to determine the expected number of load cycles for each of at least one of the one or more in-service load conditions of the physical structure; redefines the fatigue safety factor inequality constraint of the modeling object based on a damage fraction calculated according to the required number of load cycles of the modeling object and the expected number of load cycles for each of the one or more in-service load conditions of the physical structure; and calculates the shape change rate of the implicit surface in the level set representation of the three-dimensional shape, at least according to the fatigue safety factor inequality constraint.

[0010] The one or more specifications may include two or more specifications for corresponding different materials used to manufacture the physical structure, the data may include data relating fatigue strength to load cycles for each of the different materials, determining the expected number of load cycles may include determining a separate expected number of load cycles for each of the different materials, and redefining the fatigue safety factor inequality constraint may include: calculating a separate fatigue safety factor for each of the different materials based on a corresponding damage fraction calculated according to the corresponding expected number of load cycles for the different materials; and redefining the fatigue safety factor inequality constraint for the modeled object using the minimum of the fatigue safety factors for the different materials; and wherein calculating the shape change rate may include calculating at least one shape change rate using a gradient determined according to the shape derivative of the fatigue safety factor.

[0011] The search may include calculating the maximum stress value under the load condition in use, based at least on the standard deviation of the stress distribution in the current numerical evaluation of the physical response of the modeled object.

[0012] The enforcement may use a thickness measurement on the three-dimensional shape of the generated design of the modeling object, the thickness measurement being a combination of at least two different thickness measurements.

[0013] The at least two different thickness measurements may include (i) a first distance measurement, the first distance measurement being the length within the modeling object from a surface point of the modeling object projected in the negative normal direction; and (ii) a second distance measurement, the second distance measurement being the diameter of the largest sphere fitted within the modeling object and contacting the surface point of the modeling object, as determined by examining discrete sampling positions defined on the surface of the sphere.

[0014] The enforcement may include using an inequality constraint based on volume fraction or minimum thickness as a proxy for the design criterion limiting the minimum thickness, wherein the inequality constraint based on volume fraction or minimum thickness is modified using an importance coefficient, which is set to zero during the initial phase of the iterative modification and adjusted during subsequent phases of the iterative modification based on whether one or more other constraints were violated in previous iterations of the iterative modification.

[0015] The adjustment may include adjusting the target value of the inequality constraint based on volume fraction or minimum thickness between an initial target value and a final target value during the multiple iterations of the iterative modification; and when adjusting the target value during the multiple iterations, using a proportional-integral-derivative controller to adjust and stabilize the changes made to the amount of modification of the modeling object determined according to the evaluation of the inequality constraint based on volume fraction or minimum thickness.

[0016] Obtaining the critical fatigue crack length of the material may include: obtaining one or more specifications of the material used to manufacture the physical structure; and calculating the critical fatigue crack length of the material based on information in the one or more specifications, the information including the modulus of the fatigue crack growth curve of the material.

[0017] The 3D shape of the generative design of the modeling object may include a level set representation of an implicit surface, the one or more design criteria may include the required number of load cycles for the modeling object under each of the one or more in-use load conditions of the physical structure, and the iterative modification may include: performing a numerical simulation of the modeling object based on the current form of the 3D shape and the one or more in-use load conditions to produce a current numerical evaluation of the physical response of the modeling object; using the current numerical evaluation and the thickness measurement to determine the expected number of load cycles for each of the one or more in-use load conditions of the physical structure to enforce the design criteria that limit the minimum thickness; based on the required number of load cycles for the modeling object and The fatigue safety factor inequality constraint of the modeling object is redefined by calculating the damage score based on the expected number of load cycles for each of the one or more in-use load conditions of the physical structure; the shape change rate of the implicit surface is calculated at least according to the fatigue safety factor inequality constraint; the shape change rate is used to update the level set representation to generate an updated pattern of the three-dimensional shape of the modeling object; and the execution, determination, redefinition, calculation, and update are repeated at least once until a predefined number of shape modification iterations have been performed, or until the three-dimensional shape of the generated design of the modeling 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 conditions.

[0018] The one or more in-use load conditions of the physical structure may include two or more in-use load conditions of the physical structure, the one or more design criteria may include the required number of load cycles for the modeling object under each of the two or more in-use load conditions of the physical structure, determining the expected number of load cycles may include determining the individual expected number of load cycles for each of a plurality of points on the implicit surface under each of the two or more in-use load conditions, and redefining the fatigue safety factor inequality constraint may include: summing load-specific damage fractions corresponding to the two or more in-use load conditions for each of the plurality of points, wherein each load-specific damage fraction includes dividing the expected number of load cycles for one of the plurality of points and one of the in-use load conditions by the required number of load cycles for the one of the in-use load conditions to produce a sum of the load-specific damage fractions for each of the plurality of points; taking the reciprocal of each of the sums of the load-specific damage fractions; and using the minimum of the sum of the reciprocals of the load-specific damage fractions to redefine the fatigue safety factor inequality constraint of the modeling object.

[0019] Calculating the shape change rate may include using a quantity determined according to a shape derivative formula that approximates the shape derivative of the fatigue safety factor to calculate at least one shape change rate.

[0020] The shape derivative formula, which approximates the fatigue safety factor, may include an inequality constraint based on a volume fraction or minimum thickness, which is modified using an importance coefficient adjusted based on whether one or more other constraints were violated in the iteratively modified previous iteration.

[0021] The one or more methods may include: adjusting the target value based on the inequality constraint of volume fraction or minimum thickness between an initial target value and a final target value during the iteratively modified multiple iterations; and when adjusting the target value during the multiple iterations, using a proportional-integral-derivative controller to stabilize the change made to the quantity determined according to the shape derivative formula, and adjusting the overall contribution of the quantity determined according to the shape derivative formula to the rate of shape change used in the update.

[0022] The methods described herein, and others, can be implemented using a non-transitory computer-readable medium encoding a computer-aided design program operable to cause one or more data processing devices to perform one or more methods. In some implementations, a system includes: a non-transitory storage medium storing instructions of a computer-aided design program; and one or more data processing devices configured to execute the instructions of the computer-aided design program to perform the one or more methods. Furthermore, such systems may include additive manufacturing machines or other manufacturing machines, and the one or more data processing devices may be configured to execute the instructions of the computer-aided design program to generate instructions (e.g., toolpath specifications for the additive manufacturing machine) for such a machine based on the three-dimensional model, and to manufacture the physical structure corresponding to the object by the machine using the instructions (e.g., the additive manufacturing machine using the toolpath specifications).

[0023] Specific implementations of the subject matter described in this specification can be implemented to achieve one or more of the following advantages. The design can be optimized based on damage-tolerant fatigue constraints. During topology optimization, fatigue-optimal design requirements, such as fatigue-constrained requirements based on structural loads, can also be met. These requirements can be, for example, the number of cycles for in-service load conditions or a time period converted to the number of cycles. The techniques described in this specification can be applied to objects with a variety of different materials, and individual fatigue constraints can be calculated based on the corresponding material information and specified load conditions. Therefore, by generating a design relative to user-specified fatigue constraints, the resulting design of the object can tolerate a certain amount of damage to avoid severe failure of the corresponding physical object, which could otherwise occur if the object suffers absolute fatigue damage. Thus, a predetermined amount of damage tolerance can be constructed into the physical object from the generative design phase beginning with object creation.

[0024] This process can be performed automatically and in conjunction with other techniques, such as level set methods, solid isotropic material penalty (SIMP) methods, conventional body-fit solvers, backup and recovery for improved stability, load-condition-specific advection for geometry to prevent disconnections, proportional-integral-derivative (PID) controllers, other adaptive controllers and related techniques, including adaptive PID tuning and PID autotuning, the latter of which can be used to satisfy arbitrary design constraints.

[0025] Details of one or more embodiments of the subject matter described herein are set forth in the accompanying drawings and the following description. Other features, aspects, and advantages of the invention will become apparent from the description, drawings, and claims. Attached Figure Description

[0026] Figure 1A An example of a system that can be used to design and manufacture physical structures is shown.

[0027] Figure 1B An example of the process of designing and manufacturing physical structures is shown.

[0028] Figure 2A A graphical representation of an example of a geometric mapping between an initial design configuration and a current design configuration is shown.

[0029] Figure 2B A graphical representation of an example of data mapping from a solid grid to a horizontal set raster is shown.

[0030] Figure 2C An example of the problem setting for a topology optimization process based on level sets is shown; Figure 2D Showing from Figure 2C The problem is set up with a narrowband speed instance; and Figure 2E Showing from Figure 2C The speed after expansion.

[0031] Figure 3A An example of a process is shown that uses one or more generative design processes to generate one or more parts of a 3D model of an object to be manufactured.

[0032] Figure 3B A graphical representation of an instance of an initial inoculation geometry applied to the design space is shown, along with the resulting design of a Michell-type arch problem after performing a generative design process on the design space.

[0033] Figure 3C A graphical representation of an instance of the bubble insertion history for an intermediate design of a Michel type arch problem during the execution of the generative design process is shown.

[0034] Figure 4A An example is shown of a process that uses one or more generative design processes that implement arbitrary constraint processing with controlled convergence to generate one or more parts of a 3D model of an object to be manufactured.

[0035] Figure 4B Examples of graphs showing the gradual decrease and increase of the target constraint value when using B-spline tracking are shown.

[0036] Figure 4C A graphical representation of example variations of the μ value as the target volume is shown, taking into account different voxel sizes.

[0037] Figure 4DAn example is shown of a graph tracking the target volume fraction versus the actual volume fraction during an iterative optimization process with approximate volume control but no adaptive control.

[0038] Figure 4E An example is shown of a graph tracking the target volume fraction versus the actual volume fraction during an iterative optimization process with PID control (but without adaptive PID control).

[0039] Figure 4F An example of the process of adaptively modifying the parameter values ​​of a PID controller is shown.

[0040] Figure 4G Examples of graphs showing different measures used in tracking constraint normalization are shown.

[0041] Figure 4H Examples of graphs showing the convergence history of tracking constraints with and without line search are shown.

[0042] Figure 5A An example is shown of a process that generates one or more parts 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 subject.

[0043] Figure 5B An example of a SN curve tracking fatigue strength (stress) on a body over multiple cycles is shown.

[0044] Figure 6A An example is shown of a process for iteratively modifying the 3D shape of a generative design of a modeling object in the design space according to one or more design criteria, including at least one damage tolerance fatigue constraint.

[0045] Figure 6B A graphical representation of the geometry is shown, where the thickness of the geometry is calculated according to different measurement techniques.

[0046] Figure 6C An example of a graph showing the curves of the sigmoid function at different rates is shown.

[0047] Figure 6D An example of a table describing the properties of fatigue cracks is shown.

[0048] Figure 6E An example is shown of a process that generates one or more parts of a 3D model of an object to be manufactured by using one or more generative design processes that solve for damage tolerance fatigue constraints on the subject.

[0049] Figure 7AIt is an example of a process that iteratively modifies the three-dimensional shape of a generative design by modifying modeling objects in the design space according to one or more design criteria, including stress constraints.

[0050] Figure 7B A graphical representation of an instance of a construction angle measured on an object is shown.

[0051] Figure 7C An example of a process is shown that uses stress constraints to generate one or more parts of a 3D model of an object to be manufactured.

[0052] Figure 8A An example is shown of a process that iteratively modifies the 3D shape of a generative design of modeling objects in the design space according to one or more design criteria, while avoiding excessive abrupt changes and minimizing the possibility of broken connections.

[0053] Figure 8B It is a graphical representation of an instance where the geometric connections are broken during optimization.

[0054] Figure 8C It is a graphical representation of instances of geometric figures with analog elements, which are classified based on the intersection of the elements and the geometric figures.

[0055] Figure 9 A schematic diagram of a data processing system is shown, which includes data processing equipment that can be programmed as a client or a server.

[0056] The same reference numbers and markings in each figure indicate the same elements. Detailed Implementation

[0057] Figure 1A An example of a system 100 for designing and manufacturing physical structures is shown. Computer 110 includes a processor 112 and memory 114, and computer 110 is connectable to a network 140, which may be a private network, a public network, a virtual private network, etc. Processor 112 may be one or more hardware processors, each of which may include multiple processor cores. Memory 114 may include both volatile and non-volatile memory, such as random access memory (RAM) and flash RAM. Computer 110 may include various types of computer storage media and devices, which may include memory 114 to store program instructions that run on processor 112, including one or more computer-aided design (CAD) programs 116 that implement three-dimensional (3D) modeling capabilities and include one or more generative design processes for topology optimization, which are performed using numerical simulation, and include material or microstructure shape optimization techniques, geometric or macrostructure shape optimization techniques, or both (e.g., using one or more level set-based topology optimization processes).

[0058] Numerical simulations performed by one or more CAD programs 116 can model one or more physical properties and can use one or more types of simulations to generate numerical evaluations of the physical response (e.g., structural response) of the modeled object. For example, finite element analysis (FEA) can be used, including linear static FEA, one or more finite difference methods, and one or more material point methods. Furthermore, simulations of physical properties performed by one or more CAD programs 116 can include computational fluid dynamics (CFD), acoustic / noise control, heat conduction, computational injection molding, electrical or electromagnetic flux, and / or material curing (which can be used for phase transitions during the molding process) simulations. Additionally, one or more CAD programs 116 can potentially implement hole and / or fixture generation techniques to support clamping and / or manufacturing control functions during manufacturing.

[0059] As used herein, CAD refers to any suitable program used to design a physical structure that meets design requirements, regardless of whether the CAD program is able to interface with and / or control manufacturing equipment. Therefore, one or more CAD programs 116 may include one or more computer-aided engineering (CAE) programs, one or more computer-aided manufacturing (CAM) programs, etc. One or more CAD programs 116 may run locally on computer 110, remotely on one or more remote computer systems 150 (e.g., one or more server systems of one or more third-party vendors accessible from computer 110 via network 140), or both locally and remotely. Therefore, CAD program 116 may be two or more programs operating collaboratively on two or more separate computer processors, wherein one or more programs 116 operating locally on computer 110 can “offload” processing operations (e.g., generating design and / or numerical simulation operations) by causing one or more programs 116 on one or more computers 150 to perform offloading processing operations.

[0060] One or more CAD programs 116 present a user interface (UI) 122 on a display device 120 of a computer 110, which can be operated using one or more input devices 118 of the computer 110 (e.g., a keyboard and mouse). It should be noted that although in Figure 1AWhile the display device 120 and / or input device 118 are shown as separate devices, they may also be integrated with each other and / or with the computer 110, such as within a tablet computer (e.g., a touchscreen may be an input / output device 118, 120). Furthermore, the computer 110 may include or be part of a virtual reality (VR) and / or augmented reality (AR) system. For example, input / output devices 118, 120 may include VR / AR input gloves 118a and / or VR / AR headset 120a. In any case, the user 160 can interact with one or more CAD programs 116 to create and modify one or more 3D models, which may be stored in one or more 3D model documents 130.

[0061] In the example shown, the initial 3D model 132 is a seed model used as input into the generative design process. In this example, the user 160 has defined the mechanical problem to be solved in the generative design process to generate a new 3D model from the initial 3D model 132. In this case, the defined problem is a Michelson-type arch problem, where the user 160 has specified the domain 134 and the load condition 136. However, this is only one of many possible instances.

[0062] In some implementations, user 160 (or other person or program) may specify the design space of the object to be manufactured, the numerical simulation settings for the object's numerical simulation (e.g., FEA, CFD, acoustic / noise control, thermal conduction, computational injection molding simulation, electrical or electromagnetic flux, material curing, etc.) (e.g., one or more loads and one or more materials), at least one design objective for the object (e.g., minimizing material usage), and at least one design constraint for the object (e.g., volume constraint). In some implementations, inputs for the numerical simulation and generative design process may include: one or more regions of the current 3D model, in which new 3D geometry is generated; one or more load conditions defining one or more loads in one or more different directions borne by the physical structure being designed; one or more materials (e.g., one or more isotropic solid materials of the baseline material model determined as the design space); one or more seed model types used as inputs to the generative design process; one or more generative design processes used; and / or one or more lattice topologies in one or more regions of the design space. The inputs to the design and numerical simulation process can include non-design space, different types of components (e.g., rods, bearings, housings), one or more target manufacturing processes and associated parameters, obstacle geometry to be avoided, retained geometry to be included in the final design, and parameters related to various aspects, such as design resolution, composition type, etc.

[0063] Furthermore, one or more CAD programs 116 provide user interface elements in UI 122 to enable user 160 to specify the various types of inputs described above, and all such inputs (or various subsets thereof) can be used in the generative design and numerical simulation processes described in this document. Additionally, the UI 122 of one or more CAD programs 116 allows user 160 to design parts using conventional 3D modeling functions (to construct an accurate geometric description of the 3D design model), and then use the generative design and simulation processes within a specified design space within one or more portions of the 3D design model. Therefore, as will be understood, many possible types of physical structures can be designed using the systems and techniques described in this document, UI 122 can be used to create a complete mechanical problem definition for the parts to be manufactured, and the generative design and numerical simulation processes can accelerate new product development by achieving improved performance without time-consuming physical testing.

[0064] Furthermore, as described herein, one or more CAD programs 116 implement at least one generative design process, enabling the one or more CAD programs 116 to automatically generate one or more parts (or the entire 3D model) of one or more 3D models based on one or more design objectives and constraints (i.e., design criteria), wherein the geometric design is iteratively optimized based on simulation feedback. It should be noted that, as used herein, “optimization” (or “optimal”) does not mean achieving the best design among all possible designs in all cases, but rather that the best (or near-best) design is selected from a finite set of possible designs that can be generated by available processing resources within an assigned time (e.g., specified by a predefined number of shape modification iterations). Design criteria may be defined by user 160 or another party and imported into one or more CAD programs 116. Design criteria may include structural integrity constraints on individual parts (e.g., requiring that the part should not fail under expected structural loads during its use) and physical constraints imposed by the larger system (e.g., requiring that the part be contained within a specified volume so as not to interfere with one or more other parts in the system during use).

[0065] Various generative design processes can be used to optimize the shape and topology of at least a portion of a 3D model. Iterative optimization of the geometric design of one or more CAD programs 116 for one or more 3D models involves topology optimization, a lightweight approach where the optimal distribution of materials is determined by minimizing an objective function constrained by design constraints (e.g., volume-constrained structural compliance). Topology optimization can be solved using a variety of numerical methods, which can be broadly categorized into two groups: (1) material or microstructure techniques, and (2) geometric or macrostructure techniques. Microstructure techniques are based on determining the optimal distribution of material density and include the Solid Isotropic Material Penalty (SIMP) method and the homogenization method. In the SIMP method, intermediate material densities are penalized to favor ρ = 0 or ρ = 1, representing voids or solids, respectively. In the homogenization method, intermediate material densities are treated as composite materials.

[0066] In contrast, macrostructural techniques treat materials as homogeneous, and the resulting 3D topology of the modeled object is represented as one or more boundaries between one or more solid regions (with homogeneous material) and one or more void regions (without material) within the design space (also known as the domain or subspace of the domain for topology optimization). During the generative design process, one or more shapes of one or more boundaries are optimized, and the topology changes within the domain as shape optimization is combined with the addition / removal and shrinking / growing / merging of one or more void regions. Therefore, the type of final optimized topology that can be generated from the generative design process using macrostructural techniques can depend significantly on the number and size of voids within the seed geometry and the addition and removal of voids during the optimization process.

[0067] It should be noted that, although in Figure 1A Only one seed model 132 is shown (where this model 132 includes a complex solid region 132A with a plurality of holes 132B surrounding the void region), but it should be understood that for any given generative design process iteration, the generative design process described in this document may employ two or more seed geometries / models to improve the final result of shape and topology optimization. Furthermore, during the shape and topology optimization process, one or more voids may be introduced into the solid domain and / or one or more solids may be introduced into the void domain to improve the final result of shape and topology optimization. Therefore, one or more CAD programs 116 may include various types of available seed geometry and intermediate process geometry introductions, as well as user interface elements that allow user 160 to design their own seed geometry and intermediate process geometry introductions. Similarly, user 160 may run two or more generative design process iterations (saving the results of each iteration) until a preferred generative design is produced.

[0068] In various implementations, as described herein, one or more CAD programs 116 provide a generative design shape optimization process that: (1) utilizes controlled convergence, (2) achieves damage prevention within load cycles, (3) has fatigue damage of finite size, (4) uses a constructed material strength model, and / or (5) has singularity and connectivity breakage prevention. In some implementations, one or more CAD programs 116 implement all of the processes listed above, while in others, one or more CAD programs 116 implement a subset of the processes listed above.

[0069] Once user 160 is satisfied with the generated 3D model of the design, the 3D model can be stored as a 3D model document 130 and / or another representation used to generate the model (e.g., an .STL file for additive manufacturing). This can be done upon request from user 160, or upon request from the user for another action, such as sending the 3D model 132 to an additive manufacturing (AM) machine 170, or other manufacturing machinery that can be directly connected to computer 110 or connected via network 140, as shown in the figure. This may involve post-processing performed on local computer 110 or a cloud service to export the 3D model 132 to an electronic document upon which manufacturing is based. It should be noted that an electronic document (hereinafter referred to as a document for brevity) can be a file, but does not necessarily correspond to a file. A document can be stored as a part of a file that stores other documents, as a single file dedicated to the document in question, or as multiple collaborative files.

[0070] In any case, one or more CAD programs 116 may provide document 135 (with a toolpath specification in an appropriate format) to AM machine 170 to produce a complete structure 138, which includes optimized topology and shape (in this example, an arch design generated for a Michel type arch problem). AM machine 170 may employ one or more additive manufacturing techniques, such as particle techniques (e.g., powder bed melting (PBF), selective laser sintering (SLS), and direct metal laser sintering (DMLS)) and extrusion techniques (e.g., fused deposition modeling (FDM), which may include metal deposition AM). Furthermore, user 160 may save or transfer the 3D model for later use. For example, one or more CAD programs 116 may store document 130 including the generated 3D model.

[0071] In some implementations, one or more subtractive manufacturing (SM) machines 174 (e.g., computer numerical control (CNC) milling machines, such as multi-axis, multi-tool milling machines) may also be used in the manufacturing process. These one or more SM machines 174 may be used to prepare an initial workpiece on which one or more AM machines 170 will operate. In some implementations, a partially complete structure 138 is generated by one or more AM machines 170 and / or using casting methods (e.g., investment casting (IC) using a ceramic shell or sand casting (SC) using a sand core), and then one or more portions of this partially complete structure 138 are removed (e.g., finished) by the CNC machine 174 to form the complete structure. In some implementations, one or more CAD programs 116 may provide the SM machine 174 with corresponding documents 135 (toolpath specifications in an appropriate format, such as CNC programs) for manufacturing the part using various cutting tools, etc. Furthermore, in some implementations, the complete structure 138 is generated integrally using one or more SM machines 174.

[0072] In various implementations, one or more CAD programs 116 of system 100 may implement one or more generative design processes as described in this document. The generative design process seeks optimal geometry, topology, or both. For example, the generative design process seeks optimal geometry in an alternative design by minimizing a performance-related objective function subject to constraints.

[0073] minimize

[0074] Make g i (s, u(s)) = 0, i = 1, ..., n g (2)

[0075] Where s is a vector of design variables related to the geometry of the domain, and u is a vector of state variables (e.g., displacements) dependent on s. Additional constraints (e.g., equilibrium) are represented by the set g. i Represented. For simplicity, equality constraints are assumed here. The mathematical programming method used to minimize (1) can be gradient-based or non-gradient-based. Gradient-based methods (in contrast to non-gradient-based methods) typically use more information related to design sensitivity, such as:

[0076]

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

[0078] Figure 1BAn example of a process for designing and manufacturing a physical structure is illustrated. For example, 180 design variables are obtained from one or more CAD programs 116 to generate a generative 3D model. Different generative design processes can be formulated by using different combinations of design variables, which may include lattice, density field, and level set. In some implementations, design variables may include various types of input, such as input received via UI 122, such as the selection between different generative design synthesis methods available through one or more CAD programs 116 in system 100. In some implementations, available generative design synthesis methods include (1) level set-based topology optimization, which provides a basic level set method for topology optimization, (2) lattice and skin optimization, which provides lattice and skin thickness optimization, (3) hybrid topology optimization, which provides topology optimization with lattice filling, (4) inside-out hybrid topology optimization, wherein the lattice filling exists in the negative space between the topology-optimized design and the original design space, (5) hollow topology optimization, which provides a method for topology optimization using internal hollow regions, and / or (6) hybrid-hollow topology optimization, which provides a method for topology optimization using lattice filling and internal hollow regions. Further details regarding such generative design synthesis methods can be found in PCT / US2019 / 060089, filed November 6, 2019, and published on May 14, 2020 as WO 2020 / 097216, entitled “MACROSTRUCTURE TOPOLOGYGENERATION WITH DISPARATE PHYSICAL SIMULATION FOR COMPUTER AIDED DESIGN AND MANUFACTURING,” which is incorporated herein by reference. Furthermore, available generative design synthesis methods may employ the following combined… Figures 3A to 8C One or more of the systems and technologies described.

[0079] Additional design variables are possible, such as (1) a design space for generating the design geometry, for example, a boundary representation (B-Rep) 3D model designed or loaded into one or more CAD programs 116, which serves as a subspace of the optimization domain of the described generative design process, and / or (2) a set of input solids that specify the boundary conditions for generating the design geometry, for example, using UI 122 selection to specify the B-Rep of one or more subspaces, which are reserved as one or more connection points connecting to one or more other parts in a larger 3D model or one or more separate 3D models. For example, different combinations of design variables may be used by one or more CAD programs 116 in response to input from user 160. For example, user 160 may select different generative design synthesis methods to use within corresponding different design spaces within a single 3D model.

[0080] Other design variables may include settings for numerical simulation, such as the density of elements in the FEA model or a homogenized lattice material representation of a selected lattice topology to be used with the topology-optimized 3D shape of the part that generates the design. Design variables may include various design objectives and constraints, such as those described in this document. Furthermore, functionality to assist users in specifying design variables may be provided, for example, by one or more CAD programs 116. For example, a lattice recommender may provide predictions of suitable lattice settings for a given problem using a single solid-state simulation. In some implementations, a lattice recommender described in PCT Publication No. WO 2017 / 186786 A1, entitled “METHOD AND SYSTEM FOR GENERATING LATTICE RECOMMENDATIONS IN COMPUTER AIDED DESIGN APPLICATIONS”, filed April 26, 2017, is used, which is incorporated herein by reference.

[0081] When generative design variables are specified, for example, one or more CAD programs 116 use one or more generative design processes (e.g., using one or more selected generative design synthesis methods) to generate one or more 3D models. In some implementations, one or more generative design processes may use the described level set method, where s from equations 1, 2, and 3 represents the boundary of a solid region implicitly represented using one or more level sets, which are signed distance values ​​computed on a Cartesian background grid. In level set-based topology optimization methods, the external shape of the structure is represented by a one-dimensional higher-order level set function, and changes in shape and configuration are replaced by changes in the level set function values ​​to obtain an optimal structure. The level set function is a function that indicates whether each part of the design domain for setting the initial structure corresponds to a material domain (material phase) that forms the structure and is occupied by material, a void domain (void phase) that forms voids, or the boundary between these two domains, where a predetermined value between the value representing the material domain and the value representing the void domain indicates the boundary between the material domain and the void domain.

[0082] In some implementations of level-set-based topology optimization methods, one or more octree data structures are used to accurately resolve the geometry. Level-set-based topology optimization typically involves optimizing the shape of the design domain using shape derivatives, which are derivatives of the constraint minimization problem with respect to the shape. Applying shape changes to the level set allows for topology changes during shape modifications. The result 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 accompanied by topology changes. For this type of level-set-based topology optimization, and for variations of this type described in this document, one or more of the following methods can be used.

[0083] Linear elastic topology optimization

[0084] Consider the linear elastic boundary value problem of a solid in the domain Ω:

[0085] In Ω (4)

[0086] u=0 at Γ D (5)

[0087] In Γ N (6)

[0088] Where ∈(u) is the linear strain tensor, D is the fourth-order constitutive tensor, u is the displacement vector, f is the external load vector, and t is the line under the outward normal n at the Neumann boundary Γ. NThe specified traction force. For simplicity, it can be targeted at Γ. D Assume homogeneous Dirichlet boundary conditions. The constrained topology optimization problem can then be categorized as follows:

[0089] Minimize J(Ω, u) (7)

[0090] Subject to In Ω (8)

[0091] u=0 at Γ N Above (9)

[0092] D∈(u)n=t in Γ N (10)

[0093] Minimizing compliance can be used as the objective function.

[0094]

[0095] Figure 2A A graphical representation of an example of a geometric mapping 240 between the initial design configuration and the current design configuration is shown. The solution space in topology optimization can be defined by different perturbations of the geometry within the design space. In this context, a linear mapping 242 can be defined to map a given domain 246Ω to a perturbation domain 248Ωt. Under this mapping, material points with coordinates x∈Ω can be mapped 244 to the following terms.

[0096] x t =x+tδv, t≥0 (12)

[0097] Where δv is a defined constant vector field, and t is a scalar parameter (see [reference]). Figure 2A It should be noted that solving the equations using gradient-based mathematical programming methods involves using the directional derivative of the objective function in the direction of the velocity field δv.

[0098]

[0099] Several methods can be used to obtain the directional derivative of the objective function for use in gradient-based optimization methods. Suitable methods for use with gradient-based optimization methods include direct differentiation, semi-analytical derivatives, adjoint methods, and finite difference methods. Furthermore, the following section combines... Figures 4A to 4F It describes in detail other methods for obtaining the value of the directional derivative of the objective function, including approximation techniques.

[0100] Accompanying method

[0101] Evaluating the shape derivative (Equation 13) may require the directional derivative of the state variable u in the direction of the velocity vector δv. This can be determined using the following chain rule.

[0102]

[0103] However, in some implementations, an adjoint method can be used, which involves the formation of a Lagrangian L(Ω, u, λ) that depends on the domain shape Ω, the displacement field u, and the Lagrangian parameter λ.

[0104]

[0105] The stationarity condition of the Lagrange (i.e., δL(Ω, u, λ) = 0) yields a complete set of shape optimization equations. For example, the accompanying problem of compliance minimization (Equation 11) can be given by considering the variation of the Lagrange with respect to displacement u. After introducing the cost function (Equation 11) and restating the domain term using the divergence theorem as follows...

[0106]

[0107] The corresponding boundary value problem, called the adjoint problem, can be transformed into:

[0108] In Ω (17)

[0109] λ = 0 at Γ D Middle (18)

[0110] In Γ N (19)

[0111] This leads to the determination that λ = -u is a solution to the adjoint problem. This means that the adjoint problem (Equations 17-19) does not need to be explicitly solved for the compliance minimization problem (Equation 11). Such problems are called self-adjoint problems, where the solution to the direct problem also produces an adjoint solution. However, this situation is not common and may require solving different adjoint problems depending on the properties of the direct problem and the objective function. The advantages of using the Lagrange as is include the identities:

[0112]

[0113] This equation allows the shape derivative (Equation 13) to be expressed as a boundary integral of the following form:

[0114]

[0115] Without losing generality, we can assume that some boundary variations are independent of the actual shape optimization. In solid mechanics, boundary variations typically take the following forms:

[0116] δv=0 at Γ D(22)

[0117] δv=0 at Γ N Among them

[0118] δv≠0 at Γ N Where σn=0.

[0119] This means that during shape optimization, the boundary Γ N Only the parts without traction can move freely. In this context, with structural compliance (Equation 11) as a cost function, the Lagrangian (Equation 15) in the direction... The above can be transformed into:

[0120]

[0121] Without restricting δv as described in Equation 22, the variant of the Lagrange can include several additional terms. During the iterative optimization of the shape, the shape derivative (Equation 23) can be used as gradient information. To achieve the maximum reduction of the objective function, the boundary perturbation can be selected as follows.

[0122] δv=-(2u·fD∈(u):∈(u)). (twenty four)

[0123] This boundary perturbation can occur along the normal. The direction is applied, where v is the rate of shape change and is given by the following terms.

[0124]

[0125] Volume control

[0126] Topology optimization using only the compliance minimization objective (Equation 11) can produce an optimal topology covering the entire design space. Therefore, some form of volume constraint is usually required. Furthermore, in some implementations, control over volume changes during topology optimization may be important for several reasons: 1) to enforce volume constraints; 2) to provide the user with control over the topology optimization process, for example, with larger volume changes during the initial iterations and smaller changes during later iterations; and 3) to ensure that arbitrary constraints without shape derivatives are satisfied.

[0127] It should be noted that the presence of shape derivatives with respect to constraints may necessitate modifying the rate of shape change in Equation 25. A modified objective function could be considered, where the volume is penalized by the penalty parameter μ in the following terms:

[0128]

[0129] The corresponding shape derivative (Equation 25) can be given by the following terms:

[0130]

[0131] Where μ is a constant along the boundary. The velocity term in the shape derivative (Equation 25) can now have the following additional term:

[0132] v=-(2u·fD∈(u):∈(u)+μ) (28)

[0133] Augmented Lagrange Method

[0134] In some implementations, the augmented Lagrangian method is used. Some methods (see volume control above) may have limitations, such as difficulties in satisfying the specified volume target. Essentially, the final volume of the design may depend on the value of μ specified in Equation 26. In such cases, the augmented Lagrangian method can be used to achieve the desired design constraint. When the final volume target is V... f In the case of (Ω), the following Lagrange quantities are considered for minimizing compliance:

[0135]

[0136] The shape derivative can be given by the following terms:

[0137]

[0138] The penalty parameters λ and μ can be updated in an increasing sequence, such that they converge to the optimal Lagrange multipliers. In some implementations, one or more heuristics are used to update the penalty parameters.

[0139] Body-fit solver

[0140] In some implementations, one or more body-fitted mesh-based solvers are used. Using such body-fitted mesh-based solvers with level set methods involves mapping data from a solid mesh to a Cartesian grid (it should be noted that the inverse mapping is trivial due to the structured nature of the Cartesian grid). This involves, for example... Figure 2B The two mappings shown in the figure illustrate a graphical representation of an example of a data mapping 260 from a solid grid 262 to a horizontal set raster 264.

[0141] 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 through data averaging. For example, by averaging the data at solid node n... jThe solid mesh element ei data in solid mesh 262 is averaged. Furthermore, a linear shape function can be used to map the data in the solid mesh nodes to voxel grid points. Solid mesh node n j The data at that point can be linearly interpolated to level set grid point gi in level set grid 264. This mapping allows the level set method to be used with complex FEA models solved using a body-fit solver.

[0142] A detailed example of the macroscopic topology optimization process is now described. For simplicity, the compliance minimization problem is used in conjunction with the penalty volume (Equation 26). Furthermore, for all the detailed examples below, it is assumed that FEA is used for numerical simulation for ease of presentation, but other numerical simulation types described above may also be used.

[0143] The initial shape can simply be the design space or the intersection shape of the design space and a suitable seed geometry. An initial level set ψ0 can be created by converting the initial shape into a signed distance field (SDF). This conversion can be done using an implementation in the OpenVDB toolkit. Other methods are also possible. The shape is then iteratively optimized until the objective function has converged. The FEA model used for simulation comprises solid elements in the entirety of one or more solid regions that are being optimized to generate the design space. In each iteration, the constitutive model D of each element e in the FEA model is updated by changing the elastic modulus according to the element's relative position to the current level set. For example, elements outside the level set (called void elements) are given a very low stiffness D. v The constitutive model of the internal elements is set to the stiffness D of the original solid material. s This can be achieved by checking the average set of the element nodes:

[0144]

[0145] Where n j This represents the coordinates of the element nodes. Once the FEA model is synchronized with the current geometry represented by the level set, the boundary value problem (Equations 4-6) can be solved to calculate the advection velocity. Shape changes can be applied by advection of the level set ψ using, for example, the Hamilton-Jacobi equations:

[0146]

[0147] Where v is the shape derivative (Equation 28). It should be noted that one or more heuristic shape update methods can be used, and one or more methods can be used to move the 0th isoline of the level set, taking into account the direction of movement specified by the shape derivative. The linear mapping between FEA nodes and level set grid points (see the body-fit solver above) allows FEA results, such as D∈(u):∈(u), to be transferred to the level set used to calculate the advection velocity v.

[0148] Once the objective function has converged, the contour shaping method for extracting the 0th isomorphic contour line of ψ can be used to extract the surface of the final level set. An example of the algorithm is as follows:

[0149] Level set algorithm for level set-based topology optimization

[0150] Input: (Ω, D) v D s ,μ)

[0151] Output: (Ω) s )

[0152] / / Start the level set ψ from the design space

[0153] 1: ψ0=f SDF (Ω)

[0154] / / Iterate until the convergence tolerance c is met.

[0155] 2: i = 0

[0156] 3: while i = 0 or |J i -J i-1 |>c do

[0157] / / Set the constitutive model of the FEA element to the void D v or solid D s

[0158] 4:

[0159] / / Formulate and solve the FEA problem

[0160] 5.

[0161] 6.uK -1 f

[0162] / / Calculate advection velocity

[0163]

[0164] / / Solve the Hamilton-Jacobi equation and obtain the new level set ψ i+1

[0165] 8:

[0166] 9: ψ i+1 ←ψ i

[0167] / / Calculate the target

[0168] 10:

[0169] / / Incremental iteration

[0170] 11: i←i+1

[0171] 12: end while

[0172] / / Obtain the final level set as a solid region

[0173] 13:

[0174] It should be noted that generalizing this algorithm for arbitrary objectives and constraints may require modifications, such as the following: In Section 6, any accompanying problems required to compute the shape derivatives for each objective and constraint should be addressed (see the accompanying methods above). In Section 7, augmented Lagrangian methods (see the augmented Lagrangian methods above) should be used to combine different shape derivatives to produce a single advection velocity. Further details regarding the implementation of these modifications to the generative design process can be found below in the Augmented Lagrangian Algorithm for Constrained Shape Optimization.

[0175] Geometric shapes

[0176] Suppose ∑ is a smooth, watertight, directed surface in Euclidean space, where the normal field N∑ points out into "free space" outside of ∑. A solid object can be created from ∑ by "thickening" the solid object in the negative normal direction. That is: a small... To define

[0177] Ω h :={x-sN ∑ (x): for all x∈[0, h]} (33)

[0178] To put it bluntly, Ω h It consists of all points sandwiched between ∑ and the offset portion offset from ∑ by a distance h in the negative normal direction, the offset portion being denoted as...

[0179]

[0180] Ω h boundary It consists of two non-intersecting surfaces: ∑ itself and the offset surface ∑ h .

[0181] It should be noted that it is usually called The "outward" unit normal vector field (which originates from Ω) h The defined solid material interior (pointing to the solid material exterior) is equal to ∑N on the surface. ∑ However, compared with ∑ h N on ∑ Pointing in the opposite direction. Its position is in y∈∑ h The exact formula at that point is -N ∑ (proj ∑ (y)), where proj ∑ It is a mapping that takes a point y as its nearest point on ∑ - therefore, for example, if y = x - hN is known ∑ (x) and h is small enough, then proj ∑ (y) = x.

[0182] Deformation of geometric figures

[0183] ∑ can be deformed. The deformation can be generated by a surface normal velocity function, which can be expressed as θ. ⊥ : This transformation can be achieved in several ways: for example, by representing ∑ as the zero level set of the function in the background Euclidean space, and relative to θ. ⊥ (Extensions) utilize advection. The method can be equated to first-order deformation in deformation. Importantly, Σ itself undergoes infinitesimal deformation, which is precisely the surface normal velocity function. Any deformation of Σ will have the concept of a "magnitude" of deformation, which is a positive scalar ε. For example, if the deformation is generated by advection, the magnitude of the deformation corresponds to the advection time. The deformed surface can be obtained from T. ε express.

[0184] Once ∑ is deformed, ∑ h It will also deform. It simply "drags" in a way that keeps the offset distance h within ∑ ε With ∑ h Between the deformation forms. Therefore, the thickened object Ω h The deformation of is entirely determined by the deformation of ∑. It should be noted that it can be shown that ∑ h At point y∈∑ h The infinitesimal variation at that point is exactly -θ ⊥ (proj ∑ (y)), where proj ∑ It is a mapping that takes point y as its nearest point on ∑.

[0185] Optimization based on the steepest descent

[0186] The steepest descent method can be used for a certain objective function: surface → Find the optimal ∑. This is based on the shape Taylor formula, which states that ∑ is relative to a certain θ mentioned above. ⊥ : The generated variant and the change value ε satisfy

[0187] Φ(∑ ε )≈Φ(∑)+εDφ ∑ (θ ⊥ (35)

[0188] This is for sufficiently small ε, where DΦ ∑ (Φ ⊥ Symbolize the shape derivative of Φ at ∑.

[0189] In θ ⊥ The shape derivative formula can be used to select a specific θ. ⊥ This ensures that the shape derivative term in the shape Taylor formula is negative. Therefore, if ∑ is relative to this chosen θ ⊥ If a sufficiently small change in ε occurs, the objective function will decrease. After updating ∑ by performing a variation (e.g., by advection of the level set function representing ∑), the updated shape represents the improvement relative to the objective function. To achieve further improvement, this procedure can be repeated iteratively until convergence occurs.

[0190] Shape derivative of a certain class of objective functions

[0191] Consider an objective function of the following form. Assume Φ0: It is a "volume" shape function (e.g., a shape function that can evaluate a volume domain), such as the average structural compliance as a convex linear combination of a set of load conditions and total mass. This allows the "surface" shape function to be defined by the following terms.

[0192] Φ(Σ) :=Φ0(Ω h (36)

[0193] We can assume that we know how to calculate Φ0 at any field Ω and any variation of Ω with respect to Ω. ∈ The shape derivative, whereby the arbitrary variation is derived from the boundary normal velocity function V ⊥ : The outward unit normal field is generated relative to Ω. That is, Φ0 can be standard, which means the shape derivative satisfies the Hadamard-Zolésio structure theorem and provides the following formula for sufficiently small ε.

[0194] Φ0(Ω ε)≈Φ0(Ω)+εDΦ 0，Ω (θ ⊥ )

[0195] in

[0196] Ω's function G Ω This is called the shape gradient of Φ0 under shape Ω. Ω The exact form depends on Φ0 and can be calculated, for example, with respect to average structural compliance and volume. And it was calculated.

[0197] The shape derivative of Φ can be expressed as the shape derivative of Φ0. This can be achieved by utilizing almost entirely the above formula for the shape derivative of Φ0 (applied to Ω). h ), and taking into account the boundary normal velocity in constituting This is accomplished by the properties of two non-intersecting surfaces. That is: ∑ can have a boundary normal velocity V. ⊥ :=θ ⊥ , and ∑ h It can have boundary normal velocity Then, Φ in ∑ is given by the normal velocity function θ ⊥ The generated variant ∑ ε The next first-order change is given by the following terms.

[0198]

[0199] Therefore, the formula for the desired shape derivative is:

[0200]

[0201] It should be noted that y can be simply used as ∑ h The virtual integral variable is used to emphasize that the two integrals are in different spaces and cannot be combined a priori in any simple way.

[0202] Extracting the descent direction

[0203] Recall that the utility of the exact formula for the shape derivative is that it should allow for a choice of θ in a certain way. ⊥ This causes the shape derivative to become negative, thus reducing the objective function to first order according to the shape Taylor formula. However, how this is achieved for the shape derivative of Φ calculated above is not entirely clear. This is due to the fact that there are two competing terms (e.g., the integral over ∑ and the integral over ∑). h (integral on), and taking θ into account ⊥ It is unclear how these competing factors achieve a balance.

[0204] There are two ways to proceed. The first is to apply the Hilbert space method for extracting the descent direction. The second is to apply variable transformations to ∑ h The integral over ∑ allows it to be expressed as an integral over ∑. This is straightforward: recall that for any point y∈∑ h All of these can be written in the form y = x - hN∑(x); for our purposes now, this should be considered as a transition from ∑ to ∑ h The mapping. This mapping is called n: ∑→∑ h Where n(x) := x-hN ∑ (x). Therefore, through the variable transformation formula used for surface integrals...

[0205]

[0206] Where Jac is the Jacobian determinant of n, and is based on the fact that for all x∈∑, proj ∑ (x-hN ∑ (x))=x.

[0207] The results show that the Jacobian determinant of n can be determined. With some aid from differential geometry, it can be shown that Jac(x) = 1 + hH ∑ (x)+h 2 K ∑ (x), where H ∑ It is the mean curvature of ∑, and K ∑ It is the Gaussian curvature of ∑, thus providing

[0208]

[0209] For DΦ ∑ (θ ⊥ This manipulation provides a good formula for the shape gradient of Φ at ∑, expressed as:

[0210]

[0211] Therefore, θ ⊥ =-G ∑ The choice of DΦ ∑ (θ ⊥ The value is negative, which is required by the steepest descent procedure for updating ∑ toward the optimality of Φ.

[0212] Augmented Lagrangian Algorithm for Constrained Shape Optimization

[0213] introduce

[0214] Algorithms can be developed to solve constrained shape optimization problems where a mixture of shape functions subject to equality and inequality constraints (e.g., volume or compliance) is minimized. These constraints can be represented by shape functions (e.g., target volume, aggregate stress metric, or aggregate displacement metric). Different algorithms can be used to handle non-aggregate (also known as pointwise) constraints, such as the norm of stress at each point of the shape. In general, the optimization problem can be of the form: for u satisfying a linear elastic PDE in Ω... Ω

[0215]

[0216] Subject to G E (Ω,u Ω )=0 (46)

[0217] G I (Ω,u Ω )=0 (47)

[0218] Where Adm is a set of acceptable shapes, F is a shape function with scalar values ​​that is differentiable from the shape, and G... E and G I It is a shape function that is shape-differentiable with vector values, which may depend explicitly on the shape or implicitly on the shape through its elastic response under load.

[0219] Classical augmented Lagrange algorithm

[0220] Equality constraints

[0221] The augmented Lagrangian method, in the classical specification of optimizing scalar functions of vector variables, can be most easily applied to optimization problems of the following form with equality constraints.

[0222]

[0223] Subject to g i (x) = 0, for i = 1, ..., k (48)

[0224] The augmented Lagrangian method is an enhancement of the so-called penalty method. The penalty method considers the penalty objective function. And a sequence c that tends to infinity := c k Then, it can be defined. For any i, it can never be g. i (x k The case where ) = 0. However, since an increasing sequence of c values ​​will penalize the constraints more and more severely, it is expected that x will... k Ultimately, the constraints will be satisfied; in this sense, if x k Converges at x* Then for all i, g satisfies i (x * ) = 0. Additionally, x can be expected to... * It is a solution to equation 48.

[0225] The problem with the penalty method is that, due to the poor condition number, an increasing sequence of c values ​​will cause L to... c Minimizing (x) numerically becomes increasingly difficult. The augmented Lagrange method remedies this problem. This algorithm minimizes the "augmented Lagrange quantity".

[0226]

[0227] It also maintains a series of Lagrange multipliers μ. k and the increasing penalty parameter c k Again, it can be defined. These sequences are updated in a way that makes c k It stabilizes at a large and finite value (thus avoiding the inherent pathologicalness of general penalty methods). This occurs when x... k It converges to the solution x of equation 48 * (Therefore it is optimal and feasible) and μ k Converges on x under the KKT conditions * Lagrange multipliers μ * Since the augmented Lagrangian method iterates over both the original variable x and the "dual" variable μ, this algorithm is an example of the primal-dual optimization method.

[0228] A general framework for the augmented Lagrange algorithm may require methods to increase the penalty parameter and tighten the tolerance.

[0229] 1. Select the final tolerance T 最终 .

[0230] 2. Set k = 0. Start with initial tolerances c0 and μ0. Start under the condition of being subjected to ratio T. 最终 The initial tolerance T0 has fewer constraints.

[0231] 3. When it does not converge to the tolerance T 最终 Internal time:

[0232] a. Applying unconstrained optimization algorithms to the problem When convergence to tolerance T is achieved k Internal stopping algorithm. At this convergence level, x... k+1 The output is the value of x.

[0233] b. Check if constraints are satisfied.

[0234] (i) If in tolerance T k If the internal constraints are not satisfied, then assume c k+1 :=Increase(c k And μ k+1 :=μ k .

[0235] (ii) If in tolerance T k If the internal constraints are satisfied, then according to Update the Lagrange multipliers.

[0236] c. Set tolerance T k+1 :=Tighten(T k ).

[0237] d. Incrementing k.

[0238] Note the characteristic of the above algorithm, which is the update of the Lagrange multipliers. This can be understood as follows: The unconstrained minimum of the augmented Lagrange multiplier at iteration k satisfies... or

[0239]

[0240] Of course, the solution x of optimization equation 48, which is subject to equality constraints * and Lagrange multipliers satisfy

[0241]

[0242] g i (x8) = 0 for i = 1, ..., k (52)

[0243] In the augmented Lagrange algorithm, since it is expected that iteration x k Converges at x * Therefore, it is hoped that it can make and g i (x k The transition from 0 to 0 aligns with this objective. Therefore, the Lagrange multiplier update can be interpreted as... This can be achieved using a fixed-point iterative scheme of the form.

[0244] Inequality constraints

[0245] The ingenious extension of the classic augmented Lagrange algorithm, constrained by equality, to the classic case constrained by inequality depends on two facts. For fact 1, the problem...

[0246]

[0247] Subject to g i(x)≤0 for i=1,...,k. (53)

[0248] Equivalent to the following questions

[0249]

[0250] Subject to For i = 1, ..., k (54)

[0251] The new z variable in equation 54 is called a slack variable, and it satisfies the following at the feasible optimal value: Where x * It is the feasible optimal value of equation 53.

[0252] For fact 2, the augmented Lagrange of equation 54 is

[0253]

[0254] Furthermore, the unconstrained optimizations that occur in the augmented Lagrange algorithm can be decomposed and partially solved as follows:

[0255]

[0256] It introduced This is because elementary calculus techniques can be used to explicitly... Perform minimization on top, and according to The function describes the result.

[0257] As a result, equations 53 subject to inequality constraints can be solved by applying the augmented Lagrangian algorithm (see equality constraints above) to the modified Lagrangian function.

[0258]

[0259] Its gradient (after a certain manipulation) is

[0260]

[0261] Note that, for completeness, the solution for minimizing the z variable appears in equation 56. Substituting this solution yields equation 57. For simplicity, the minimization of the z variable is rewritten in the following form.

[0262]

[0263] Where a, b, and s are fixed parameters, they can be solved as follows.

[0264] The first step is to substitute y: = z 2 And use Instead of the above minimization problem, this new problem is now very simple because the function Φ(y) := a(s+y)+b(s+y) 2 It's a simple parabola, and the problem... Find the minimum value of this parabola within the region y ≥ 0. Therefore, the global constrained minimum is either at the global unconstrained minimum of Φ or at the boundary of the constrained region y = 0, and is small either way. This results in y * = -a / 2b-s, or y * =0; Select Φ from * :=Φ(y * This achieves the minimum result. Therefore, after algebraic rearrangement, it becomes...

[0265]

[0266]

[0267] Using s=g i (x) and a = μ i Sum = c / 2, resulting in equation 57.

[0268] Applications to constrained shape optimization problems

[0269] For simplicity, the augmented Lagrange algorithm can be applied to solve equations 45-47, which have one scalar equality constraint and one scalar inequality constraint. Then, Ω∈Adm can be interpreted as representing that Ω satisfies the following admissibility constraints: each surface of Ω∈Adm contains predefined ports; each Ω∈Adm contains predefined reserved regions; and each Ω∈Adm avoids predefined forbidden regions.

[0270] The classic augmented Lagrangian algorithm (see above for equality constraints) can simulate both equality and inequality constraints, but it is suitable for shapes that satisfy admissibility constraints. The augmented Lagrangian quantity is defined as...

[0271]

[0272] in For simplicity, F and G are not specified. E G I μ for linear elastic PDE Ω The dependence of the solution. To consider only the shape function as a surface or volume integral that is a shape dependency, it is known from general principles that the shape gradient of this type of function L, evaluated under a given shape Ω, is... The expression is dL(Ω): The scalar-valued function. The shape gradient of the augmented Lagrangian in shape Ω is

[0273]

[0274] Where dF(Ω) and dG E (Ω), dG I (Ω) represent the shape gradients of the objective function and the constraint function under shape Ω, respectively.

[0275] It is also known from general principles that updating the shape Ω to reduce L... c The value of L can be maintained while admissibility constraints can be achieved in the following way: relative to L c Shape gradient (e.g., L) c The shape gradient can be appropriately zeroed when an admissibility violation is detected, and the Adalsteinsson-Sethian velocity extension algorithm is used to extend dL. c The value of (Ω) is from Expand to The velocity function formed by the projected extension of the narrow band (Ω) will represent the horizontal set function of the advection at a specific time (determined via a linear search or correlation procedure). This is in c, μ E μ I For any fixed value in the problem min Ω∈Adm L c (Ω,μ E μ I This forms the basis of an unconstrained gradient-based shape optimization algorithm for solving problems.

[0276] The augmented Lagrange algorithm for solving equations 45-47 is now given below. This algorithm requires methods for adding penalty parameters and tightening tolerances.

[0277] 1. Select the final tolerance T 最终

[0278] 2. Set "good" = 0. Start from initial c0. Started at T 最终 The initial tolerance T0 has fewer restrictions.

[0279] 3. Initialize the shape Ω0.

[0280] 4. When it does not converge to the tolerance T 最终 Internal time:

[0281] a. Applying unconstrained gradient-based shape optimization algorithms to the problem Until convergence to tolerance T k Inside. At this convergence level, the output shape is Ω. k+1

[0282] b. Check that the constraints are satisfied.

[0283] (i) If in tolerance T k If the internal constraints are not satisfied, then according to and Assume c k+1 :=Increase(c k (Multiple)

[0284] (ii) If in tolerance T k If the internal constraints are satisfied, then according to and Update the Lagrange multipliers.

[0285] c. Set tolerance T k+1 :=Tighten(T k ).

[0286] d. Incrementing k.

[0287] Pre-adsorption operation

[0288] In some implementations, one or more operations are performed on the velocity field before shape advection, including (1) narrowband velocity limiting, (2) advection prevention using an advection mask, and / or (3) velocity extension. For the first of these, narrowband velocity limiting involves limiting the velocity to a narrow band around the 0th isopleth of the level set. For the second, the advection mask has a value of 0 inside the port (the geometric interface containing the Newman and Dirichlet boundary conditions) and a value of 1 elsewhere in the domain. The shape derivative and advection velocity are multiplied by the advection mask to prevent any advection inside the port. For the third of these, the velocity field should be continuous on both sides of the 0th isopleth of the level set. However, the objective function (typically strain energy) is often only available inside the negative narrowband. Velocity extension projects all velocities inside the domain onto the positive narrowband. This can be done by sampling the velocities inside the domain at points found by moving a distance equal to the level set along the negative normal:

[0289]

[0290] exist Figures 2C to 2E Examples of the first and third operations are shown in the figure.

[0291] Figure 2C This is an example of problem setting 200, which is a topology optimization process based on level sets. Problem setting 200 intends to generate a vehicle mount, where a design space 202 is specified relative to the port 204 (also known as reserved geometry or reserved subspace) of the mount to be generated. Figure 2D and Figure 2E It is also shown in Figure 2C The cross section X is shown in perspective. Figure 2D The narrowband velocity 210 is shown. The number of voxels in the narrowband is determined by w. nb denoted as , and Δs represents the size of the voxel. Figure 2E The effect of velocity expansion relative to the 0th isopleth profile 214 of the level set is shown 212.

[0292] The above methods can be used in conjunction with various types of level set-based topology optimization, which will be discussed in detail below. Figures 3A to 8C To be further described. Return to Figure 1B All generative design processes described in this document can be implemented, for example, in one or more CAD programs 116 to provide both: (1) substantial user control over the generative design process and post-processing of the generative design output to form a final acceptable 3D model of the object; and (2) control functionality for providing a 3D model of the generative design for manufacturing the physical structure corresponding to the object. Therefore, the results of the generative design process can be presented to the user, for example, in a UI 122 on a display device 120, along with an option 190 to accept or reject the design.

[0293] If the design is rejected, then Figure 1B The process can return, for example, to one or more new design variables obtained by one or more CAD programs 116 to generate new generative 3D models. For example, the new design variables may include selecting 180 different void addition and / or removal techniques to change the topology of the 3D shape during shape optimization. For example, in an implementation using a level set-based topology optimization method, one or more voids may be introduced into the material domain based on the topological derivative of the objective function during the structural optimization process to allow changes in topology (configuration), such as introducing holes into the material domain. The following section combines... Figure 3A Further details are described regarding one or more instances of the gap introduction technique. Additionally, the resulting 180 or more design variables may include all the different user inputs described in this document that can influence the generative design process.

[0294] Once the design is not rejected, 190 Figure 1BThe process could involve providing a 3D model (195) by one or more CAD programs 116 for use in manufacturing the 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). Provision 195 may involve sending or saving the 3D model to persistent storage for use in manufacturing the physical structure corresponding to the object using one or more computer-controlled manufacturing systems. In some implementations, provision 195 involves, for example, generating a toolpath specification (195A) using the 3D model by one or more computer-controlled manufacturing systems via one or more CAD programs 116, and manufacturing (195B) at least a portion of the physical structure corresponding to the object via a toolpath specification generated for an additive manufacturing machine, for example, by one or more CAD programs 116 via one or more computer-controlled manufacturing systems.

[0295] It should be noted that the 3D model provided (195) may be a 3D model of (185) generated by a generative design synthesis method or a processed version of the generative design output. For example, in some implementations, the 3D mesh model generated by the generative design synthesis method may be converted into a watertight B-Rep 3D model before providing (195). In some implementations, the generative design output may be post-processed using the system and techniques described in U.S. Patent Application No. 62 / 758,053, filed November 9, 2018, entitled “CONVERSION OF GENERATIVE DESIGN GEOMETRY TOEDITABLE AND WATERTIGHT BOUNDARY REPRESENTATION IN COMPUTER AIDED DESIGN,” which is included in U.S. Patent No. 11,016,470, published May 25, 2021, claiming priority to U.S. Patent Application No. 62 / 758,053 and entitled “CONVERSION OF MESH GEOMETRY TO WATERTIGHT BOUNDARY REPRESENTATION.” Furthermore, in some implementations, the post-processed generative design output can be edited using the systems and techniques described in U.S. Application No. 16 / 186,136, filed November 9, 2018, and published November 5, 2019, as U.S. Patent No. 10,467,807, entitled "FACILITATE DEDITING OF GENERATIVE DESIGN GEOMETRY IN COMPUTER AIDED DESIGN USERINTERFACE". Additionally, the following incorporates... Figures 3A to 8CThe described generative design process can also be implemented using the post-processing, editing, and / or provisioning systems and techniques described above. Ultimately, although described within the context of multiple options provided by the CAD program regarding generative design, each of the generative design processes described in this document can be implemented as a stand-alone generative design process within the CAD program. Therefore, it is not a combination of the following... Figures 3A to 8C All generative design processes described need to be implemented together in any given implementation manner.

[0296] Figure 3A The following text is combined Figures 3A to 8C The process of generating one or more parts of a 3D model of an object to be manufactured, using one or more generative design processes (e.g., as described by...). Figure 1A An instance of one or more CAD programs (116 executed). Figure 3A The process is Figure 1B An instance of process 185 as defined in [the document]. Identify the design space, settings for numerical simulation, and other inputs for the selected generative design process to initiate 300 generative models of macroscopic structural (or geometric) types (e.g., level set representations of 3D models) that can be used with one or more selected generative design processes.

[0297] When the generation process employs the basic level set method for topology optimization, 300 level sets are initiated for the design space. It should be noted that the level set method is an example of macroscopic structural generative modeling techniques, where the generative model represents the object to be designed as one or more boundaries between one or more solid regions (containing homogeneous material) and one or more void regions (containing no material) within the design space. Furthermore, identifying input may involve, for example, receiving user input via UI 122 on display device 120, importing information from another program or third-party source, and / or one or more of these inputs may be predefined in a given implementation.

[0298] The setup for numerical simulations may include one or more physical properties to be simulated and one or more types of simulations to be performed, as discussed above, as well as potential surrogate modeling or other approximation methods. In some implementations, the type of numerical simulation is predefined for all uses of the program or taking into account the specific context on which the generative design process already initiated in the program is based. Furthermore, the setup for numerical simulations includes at least one set of load conditions and / or other physical environment information associated with the type of numerical simulation to be performed.

[0299] The design space can be an initial 3D model or one or more portions of an initial 3D model used as the starting geometry. In some cases, the design space can be defined as the boundary volume of all initial geometries, which is used to specify loads 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 the starting geometry can be automatically set by a genetic algorithm or other process. For example, bubble-shaped holes (e.g., hole 132B) can be placed in a domain and a genetic algorithm can be used to change the bubble size and spacing.

[0300] Inoculation and bubble methods

[0301] initial vaccination

[0302] The design space can be initialized using an inoculation process, where the design space is defined by the design space Ω and the seed geometry Ω as shown below. s Defined by the Boolean intersection between them:

[0303] Ω0=Ω∩Ω s (64)

[0304] Where Ω0 is the inoculation geometry Ω s The initial domain is applied after the initial design space Ω. The seed geometry can have various shapes, such as an array of bubbles or a mesh, which have parametric properties, such as bubble diameter and spacing. The parameters can be user-defined or automatically defined by a seeding process (e.g., bubble seeding) described below.

[0305] Initial seeding allows for more efficient optimization and flexibility in design changes. For example, initial seeding can be defined to avoid local minima and to address the need to restart optimization. Furthermore, initial seeding can facilitate the generation of design changes. In some implementations, the seed geometry is user-defined, for example, based on the final geometry generated by a previously executed design process, a stochastic process that randomly initializes the seed geometry, or based on other factors of interest to the user, or a combination of the foregoing.

[0306] Figure 3B A graphical representation of an instance of an initial inoculation geometry applied to a design space is shown, along with the resulting design for a Michelson-type arch problem after performing a generative design process on the design space. Inoculation geometries 314 and 316 represent different geometries used to initialize the design space and can be user-defined or according to an inoculation process described below, such as bubble inoculation. Inoculation geometry 316 and... Figure 1AThe seed model 132 is the same. Final designs 318 and 320 are generated by applying the same generative design process after the design space has been seeded accordingly based on geometries 314 and 316. Final design 320 is... Figure 1A The complete structure 138 is the same. The resulting final designs 318 and 320 have different geometries and different physical properties, such as different strain energies. Therefore, even if the generative design process is the same, the inoculation geometry can have a significant impact on the final design in different design spaces.

[0307] Furthermore, other inputs may 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, other inputs may include lattice topology, volume fraction, unit size, and thickness. Various types of generative design processes can be used individually or in combination, while Figure 3A Representative processes for all these different types of generative design are illustrated. After initiating a generative model for 300 macroscopic structural types, an iterative process for modifying this generative model begins to meet the design criteria of the physical structure, such as satisfying one or more design constraints and maximizing one or more design objectives.

[0308] Numerical simulations of the physical response of the current model (e.g., a level set representation of an implicit surface of a 3D shape) are performed under one or more defined loading conditions. Generally, the numerical simulation treats each solid region in the current model as a homogeneous solid material and each void region as having no material. However, in some implementations, this treatment can be altered. For example, in hybrid topology optimization, numerical simulations of the 302-modeled object are performed when at least a portion of a solid region is treated as having at least one void (e.g., macroscopic structure generation modeling techniques treat it as a solid, and the numerical simulation treats it as a lattice form partially containing voids) or at least a portion of a void region is treated as having at least one solid (e.g., macroscopic structure generation modeling techniques treat it as voids, and the numerical simulation treats it as a lattice form partially containing solid material). As another example, in the case of hollow topology optimization, numerical simulations of the 302-modeled object are performed when at least a portion of a solid region is treated as having voids (i.e., macroscopic structure generation modeling techniques treat it as a solid, and the numerical simulation treats it as containing hollow regions). Ultimately, in the combination of these two approaches, namely in hybrid hollow topology optimization, numerical simulations of the 302 modeled object are performed when at least a portion of the solid region is considered as both a partially void region and a fully void region (i.e., the macroscopic structure generation modeling technique treats it as a solid, and the numerical simulation treats it as a lattice structure that partially contains the hollow region). Further details regarding hybrid, hollow, and hybrid-hollow topology optimization can be found in PCT / US2019 / 060089, filed November 6, 2019, and published as WO 2020 / 097216 on May 14, 2020.

[0309] The results from the simulation are used to update the current model (e.g., the level set representation) based on the current numerical evaluation of the simulated physical response of the current model. For example, the rate of shape change can be calculated for the implicit surface in the level set representation of the 3D shape of the object being modeled, and the calculated rate of shape change can be used to update the level set representation 304 to produce an updated version of the 3D shape of the object being modeled. It should be noted that in various implementations, as part of the iterative modification of the 3D shape of the object being modeled, such as those described below... Figures 4A to 8C As described, other operations can be performed before, during, or after updating a 304.

[0310] Furthermore, to generate topological changes in the current model of the object, one or more gaps 308 can be inserted into the current model in each iteration of the current model's modification or in a selected modification. For example, one or more bubbles with positions, sizes, and shapes determined according to the current model can be inserted 308. Additionally, in some implementations, the insertion of one or more gaps 308 occurs only during the early portion of the iteration and / or only periodically during the iteration (e.g., at regular gap insertion intervals). Figure 3A In the instance, the insertion of one or more gaps occurs only when the current iteration is less than the predefined gap insertion cutoff value, and only when the current iteration is equal to the gap insertion interval.

[0311] bubble method

[0312] Generally, optional determination 306 and insertion 308 can be performed as part of the bubble method. The bubble method allows changes to the topology of the design space from within, since the default level set method only allows changes from the boundaries. The location of the bubbles is identified using topological derivatives, for example, using a topological shape sensitivity method that correlates the shape derivative with the topological derivative.

[0313] In some implementations, the bubble method is applied under the following characteristics:

[0314] 1. Location: The shape derivative of the Lagrange is used as a proxy for the topological derivative. The shape derivative is described in more detail below with respect to Equation 105.

[0315] 2. Frequency: No bubbles are inserted during each iteration. Instead, bubbles are inserted at user-defined or automatically determined intervals. No bubbles are inserted after the volume reduction iteration has been completed.

[0316] 3. Size: The size of the bubble inserted at a given iteration is based on a user-defined or automatically determined ratio β of the current model volume. b To calculate:

[0317] V(B t )=β b V(Ω t (65)

[0318] Where V(B) t V(Ω) is the volume of the bubble inserted at a certain iteration t, and V(Ω) is the volume of the bubble inserted at a certain iteration t. t ) is the volume of the model at iteration t.

[0319] 4. Shape: As shown below, the shape of the inserted bubble is determined based on the distribution of the topological derivatives of the elements in the current model:

[0320]

[0321] Where e k It 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 bubble is increased over multiple iterations.

[0322] Then, the volume of the element with the lowest shape derivative is accumulated until the necessary bubble volume is reached (i.e., according to equation 64). Bubble B generated at iteration t... t The resulting shape, position, and volume therefore satisfy the following terms:

[0323] B t =e1∪e2∪…∪e k stV(B t )=V(e1)+V(e2)+…+ / (e k (67)

[0324] Return to Figure 3A Optionally, after determining that the current iteration (306) is less than the gap insertion cutoff value and equal to the gap insertion interval, one or more gaps (308) are inserted into the current model. Before performing the numerical simulation (302) or before starting the model (300), the gap insertion cutoff value and the gap insertion interval can be predetermined, for example, by user input or automatically. One or more gaps can be inserted according to a bubble method (e.g., the bubble method described above with reference to equations 65-67).

[0325] Figure 3C A graphical representation of an instance of the bubble insertion history 321 for an intermediate model of a Michell-type arch problem during the execution of the generative design process is shown. Starting with model Ω0 322A at iteration 0, the bubble iteration history 321 includes: model Ω1 322B at iteration 1; model Ω at iteration 10... 10 322C; Model Ω at iteration 11 11 322D (with two bubbles inserted); model Ω at iteration 20. 20 322E; Model Ω at iteration 21 21 322F (a bubble was inserted); model Ω at iteration 30 30 322G; and model Ω at iteration 31. 31 322H (two bubbles were inserted), where the bubble insertion occurred at t = {0, 10, 20, 30}. Also in this example, the bubble volume ratio β... b It is set to 0.05, which means that the volume of each inserted bubble is equal to 5% of the volume of the current domain, i.e., consistent with Equation 65.

[0326] Although Bubble Insertion History 321 illustrates inserting bubbles after the first iteration, bubbles can also be inserted during the first iteration of the design generation process. In some implementations, multiple instances of the bubble method can be performed, where different initial conditions and parameters lead to different design variations. Systems implementing the described techniques can provide the results of design variations for the user to select the preferred variation. In some implementations, the system can automatically select design variations based on predetermined (e.g., user-provided) criteria.

[0327] Return to Figure 3A A convergence check 310 can be performed during each iteration. To determine whether the generated design converges to a stable solution, check 310 identifies situations where all design constraints are satisfied and the design objective has not significantly improved since the last one or more iterations. The numerical simulation 302, update 304, any gap insertion 308, and check 310 process iterate until convergence. Furthermore, in some implementations, the iteration process ends once check 312 shows that a predetermined number of shape modification iterations have been completed. It should be noted that the predetermined number of shape modification iterations can be set high enough to fundamentally guarantee that all design constraints will be satisfied, and the following reference... Figures 4A to 4H The described controlled convergence technique can also fundamentally guarantee convergence to a substantially optimal value for any design objective within a predefined number of iterations (if not before).

[0328] Controlled convergence

[0329] In controlled convergence techniques, a predetermined time period or number of iterations is defined at the start of shape and topology optimization. For each design constraint, a target is specified for each iteration to get closer to the current value of the constraint.

[0330] This offers at least two improvements. First, users can specify a time period or the number of iterations. The convergence rate is then controlled to solve the given problem within the user-specified time period or number of iterations. Therefore, a controlled convergence design process yields a suitable solution, taking into account the resource constraints imposed by the user regarding how much time to allocate for solution generation. Second, the possibility of oscillations and fundamental changes can be reduced by specifying target values ​​for each constraint that are closer to their current values.

[0331] Figure 4AAn example is shown of a process for generating one or more parts of a 3D model of an object to be manufactured using one or more generative design processes that implement arbitrary constraint handling with controlled convergence. Controlled convergence is described first, followed by arbitrary constraint handling. The computer-aided design program obtains the design space of the modeled object upon which the physical structure corresponding to the manufacture will be based, one or more design criteria of the modeled object, and one or more in-use load conditions of the physical structure, wherein one or more design criteria include at least one design constraint. For example, as referenced above... Figure 3A As described, obtaining 416 can be accomplished as part of starting 300. (See above reference.) Figures 3B to 3C As described, the inoculation technique can be applied to design spaces, which can further improve the generative design process.

[0332] Identify 418 iterations. An iteration count can be time or a count of iterations calculated by the user or otherwise. Calculate 420 a series of target values ​​for at least one design constraint, from the initial target value to the final target value, based on the iteration count and the function.

[0333] Assumption Represents a series of constraints, where This is the final objective value at iteration n. The objective value of each constraint to be satisfied in iteration t is defined as...

[0334]

[0335] in And N d (ξ) is a d-order B-spline calculated using recursion.

[0336]

[0337] Where N 0 (ξ)=0.

[0338] Although this description assumes the use of B-splines as the function for calculating 420, the same result can be achieved using any suitable smoothing function. Different smoothing functions include functions of the same class but different orders, such as B-splines of different orders. The choice of smoothing function and its order alters how quickly the final design constraint objective value is achieved within the identified 418 iterations. Other examples of smoothing functions that can be used include polynomials of arbitrary order, Lagrange polynomials, and subdivision curves. In some implementations, the system implementing the described technique may prompt the user with multiple reference points, and in response to receiving a reference point selected by the user, the system may, for example, use interpolation or any suitable technique to generate a curve passing through that point.

[0339] Figure 4BExamples of curves 401A and 401B are shown, respectively, illustrating the gradual decrease and increase of the target constraint value when using B-splines to track the target. (When t ≤ n) v Complete from constraints during iteration From the initial value to the final target The B-spline-based transformation. In iteration n v The final target value is maintained during the period <t≤n. In graphs 401A and 401B, the target value of the constraint is given by the following terms:

[0340]

[0341] Special attention should be paid to dividing the iteration into two parts according to Equation 70. This applies when the iteration is less than or equal to n. v The first part includes the different target values ​​for different iterations within that first part. In the second part (i.e., in n...),... v In the following section, the objective value of the constraint at each iteration is the final objective value.

[0342] Graph 401A shows the constraint target starting from the initial constraint value of 0.9. n v =30, n=50, and d=4 B-spline. Plot 401B shows the constraint objective starting from the initial constraint value of 0.25. n v =80, n=80, and d=3 B-splines. Therefore, various initial and final constraint values, as well as different smoothing functions and iteration amounts, can be used to implement optimization methods.

[0343] Approximate volume control

[0344] In some implementations, calculating the series of objective values ​​for the 420 design constraints includes calculating the objective variation of the volume fraction of the generated 3D shape. Some constraints may have shape gradients or shape derivatives that are not well defined, approximate, or completely undefined. In these cases, surrogate shape derivatives can be calculated to improve accuracy and impose more control.

[0345] Applying a constant value μ to all optimization iterations will lead to convergence to a final volume, which depends on a relative magnitude μ of the velocity v, where the velocity depends on the boundary value problem and is given by the following:

[0346]

[0347] u=0 at Γ D (72)

[0348]

[0349] Where Ω is the domain of the solid, D is the fourth-order constitutive tensor of the solid, u is the displacement vector, and f is the external load vector. Under the outward normal vector n at the Newman boundary Γ N The specified traction force. For simplicity, only Γ D The above assumes homogeneous Dirichlet boundary conditions. Alternatively, variable values ​​μ can be used. t This enables the volume target V to be achieved during the topology optimization iteration t. T，t Assume V t-1 Let represent the volume after the (t-1)th iteration. The expected volume change during iteration t is written as:

[0350] ΔV t =V T，t -V t-1 (74)

[0351] This volume change can be approximated as:

[0352] ΔV t ≈T∫ Γ a(v se +μ)dΓ (75)

[0353] Where a∈{0,1} is the advection mask and T is the time step used in solving the Hamilton-Jacobi equations. It should be noted that when using a smaller time step, the approximation error decreases to zero:

[0354]

[0355] The maximum time step is limited by the following Courant-Friedrichs-Lewy (CFL) conditions.

[0356]

[0357] Where C is a constant, Δs is the voxel size, and |v| 最大 The maximum value of the advection velocity is given by the following:

[0358] v=-(2u·fD∈(u):v(u)+μ (78)

[0359] The variables are defined as they are in equations 71-73.

[0360] Assume v u v l Indicates v se The boundary makes v l ≤v se ≤v uThe maximum speed value is now given by the following items.

[0361] |v| 最大 =v u +μ≥|u l +μ| (79)

[0362] |v| 最大 =-(v l +μ)≥|u u +μ| (80)

[0363] |v| 最大 =-(v u +μ)≥|u l +μ| (81)

[0364] |v| 最大 =v l +μ≥|u u +μ| (82)

[0365] |v| 最大 Substituting the values ​​of T into equation 76, we get:

[0366] If v u +μ≥|u l +μ| (83)

[0367] If -(v l +μ)≥|u u +μ| (84)

[0368] If -(v u +μ)≥|u l +μ| (85)

[0369] If v l +μ≥|u u +μ| (86)

[0370]

[0371] It should be noted that when the body force term is zero in equation 78 (f = 0), the numbers for different cases of μ simplify to equations 83 and 84. This results in the strain energy component of the velocity being positive, i.e., v u v l ∈R + This results in the maximum velocity being limited by either equation 83 or equation 84. Regarding m in equation 77... t The upper limit means that for any volume change ΔV, the upper limit is... tValue not found.

[0372] Figure 4C This shows that the μ value varies with the target volume ΔV, taking into account different voxel sizes Δs. t Graphical representations of example variations. Curve 401C corresponds to voxel size 2, curve 401D to voxel size 3, and curve 401E to voxel size 4. In these example variations, the upper / lower bound of μ is found using a bisection algorithm. Essentially, the upper / lower bound lies within (0 + ΔV). t The iteration begins at 1 / 2 and continues until an effective μ satisfying the corresponding assumptions in equations 79-82 is found from equations 83-86. Alternatively, when using multiple advection steps, an advection time that violates the CFL condition can be used, i.e., T = 1 / |v|. 最大 And calculate μ according to the following terms

[0373]

[0374] 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.

[0375] Return to Figure 4A Computer-aided design programs iteratively modify the 3D shape of modeled objects in the design space based on one or more design criteria and one or more in-use load conditions of the physical structure. Iterative modification includes performing 422 numerical simulations, calculating 424 the shape change rate, and using the shape change rate to update 426 one or more level set representations. Modification of the 3D shape can include both the geometry and the topology of the 3D shape. Performing 422 numerical simulations, for example, as referenced above... Figure 3A As described. Execution 422 includes performing a numerical simulation of the modeled object based on the current morphology of the three-dimensional shape and one or more in-use load conditions to produce a current numerical evaluation of the physical response (e.g., structural response) of the modeled object.

[0376] Calculation 424 involves calculating the rate of shape change of the implicit surface in the level set representation of the 3D shape based on the next (e.g., normalized) objective value from a corresponding series of objective values ​​calculated for each design constraint 420. In some implementations, the described techniques can be applied to density-based methods such as SIMP. The objective values ​​can be normalized, and the techniques for this are described below. The series of objective values ​​begins with an initial objective value and ends with the final objective value of the design constraint. (See above reference...) Figure 4B As described by equation 70, in the first part of the iteration (up to n)v During this period, the target value of 420 is calculated based on a smoothing function (e.g., B-spline). When iteratively modified more than n... v When calculating 424, the next target value used is the final target value of the corresponding design constraint.

[0377] Volume control using an adaptive controller that includes a proportional-integral-derivative (PID) controller.

[0378] Next, we will discuss the use of an adaptive controller for accurate volume control.

[0379] Generally, an adaptive controller is a technique that provides feedback to continuously calculate the error between the desired and measured values, and then applies a correction to a control parameter. A proportional-integral-derivative (PID) controller is a type of adaptive controller that continuously calculates the error between the desired and measured values ​​based on the proportional, integral, and derivative terms of a control parameter. PID controllers can be used to control different parameters that affect shape changes (e.g., the relative contribution of a specific constraint gradient). During iterative modifications, the proportional, integral, and derivative components of the PID controller are adjusted to implicitly mitigate or accelerate shape changes by applying different control values ​​to the controlled parameter in response to oscillations in the generated three-dimensional shape of the design.

[0380] The PID component of the PID controller is also adjusted to achieve a first-stage increase in the measured value in response to repeated successes or failures to satisfy the next target value normalized from the measured value. The increase and decrease of the component can be performed using a multiplier value, which is based on the average deviation of the measured value of the component from the target value. For example, the multiplier could be 1 + abs(deviation). As described below, the PID controller can be used to accurately correct the target volume of the model for each iteration of an optimization process that includes achieving controlled convergence.

[0381] For convenience, the following description will be quoted from [source name]. The volume fractions are represented to indicate the target volume fraction and the actual volume fraction at the end of iteration t, respectively.

[0382]

[0383] Where V0 represents the volume of the design space.

[0384] Considering the maximum number of iterations n and the target final volume fraction The target volume fraction for each iteration can be calculated using the controlled convergence process described in reference equation 68 above:

[0385]

[0386] Next, the value of μ is calculated (Equations 83-86) to achieve this volume target in each iteration. Note the error between the desired volume fraction and the actual volume fraction. This is retained at the end of each process, which is caused by the approximation in equation 75, resulting in a magnitude error V after iteration t. T，t -V t It should be noted that this error has already been considered in Equation 74 when calculating the volume change for iteration t+1; the volume change is the sum of the volume change for the next iteration and the error from the previous iteration.

[0387]

[0388] However, in some cases, this is insufficient to achieve the specified volume target. Therefore, a PID controller is used.

[0389] Figure 4D Examples of graphs 401F and 401G are shown, illustrating the tracking of the target volume fraction versus the actual volume fraction during an iterative optimization process with approximate volume control but no adaptive control. Graph 401F shows the difference between the target and actual volumes of the modeled object as the volume increases over multiple iterations. Graph 401G shows the difference between the target and actual volumes of the modeled object as the volume decreases over multiple iterations. In both cases, the target and actual volumes are inconsistent, and the error increases significantly over longer iterations, potentially leading to disastrous results.

[0390] 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. The error in the volume fraction (between the actual volume and the target volume) at iteration t is defined as follows:

[0391]

[0392] The target volume for iteration t is now calculated using the following terms.

[0393]

[0394] Where K p K i K d ∈R + These are PID parameters. It should be noted that the PID controller is applied to the target volume fraction change, not the target volume, to ensure that the PID parameters are independent of the initial domain volume. Unless otherwise stated, K will be used below. p =1,K i =0.1, K d=0.1. For simplicity, this type of PID controller will be represented below, where the subscript t will be omitted for clarity:

[0395]

[0396] Figure 4E Examples of graphs 401H and 401I tracking the target volume fraction versus the actual volume fraction during an iterative optimization process with PID control (but without adaptive PID control, as described below) are shown. Graph 401H shows the positive volume case (where volume is added to the model through subsequent iterations), while graph 401I shows the negative volume case (where volume is subtracted from the model through subsequent iterations). In both cases, the target volume and the actual volume closely follow each other, with the error being less than that of the same optimization performed without a PID controller, as referenced above. Figure 4D As described.

[0397] The final result of this process is based on the control parameter ΔV. t The change in response parameters Apply control. Furthermore, volume control using a PID controller and general adaptive controller techniques can be combined with all other systems and techniques described in this document, including for handling arbitrary equality and inequality constraints described below, and in conjunction with... Figures 5A to 8C Various other shape and topology optimization techniques are combined and applied as described.

[0398] Adaptive PID Tuning

[0399] In some implementations, to provide good control over a wide range of design problems, K in an adaptive (e.g., PID) controller... p K i K d The parameters are modified in response to the controller's behavior by applying multipliers equal to the average deviation of the iteration results from the target. The three key monitored controller states can be:

[0400] Oscillation: Oscillation is defined as the occurrence of more than two consecutive pairs of successes and failures during the observed time period.

[0401] Repeated failures or successes: A success or failure is considered repeated if it occurs within a predetermined number of iterations equal to a specified maximum number of iterations (e.g., 10%).

[0402] Too many failures or successes: If the relative error between a success or failure and the target exceeds a predetermined threshold, such as 10%, then the success or failure is considered excessive.

[0403] PID control theory states that the integral term should be used to reduce system error, while the derivative term should be used to suppress oscillations. The proportional term guides the convergence speed. In practice, it is often found that the derivative term is difficult to adjust and can actually lead to oscillating behavior (a phenomenon known as "derivative shock").

[0404] Figure 4F An example of an adaptive modification of the parameter values ​​of a PID controller is shown. If oscillation is determined (402), the proportional and integral terms are decreased. If the derivative term is currently zero, it is set to a percentage of the proportional term; otherwise, the derivative term is decreased (404). The goal of these changes is to mitigate volume change. If repeated successes or failures are determined (406), the proportional term is increased, the integral term is decreased, and the derivative term is set to zero (408). The goal of these changes is to slightly increase volume change to eliminate small systematic errors. If too many repeated failures or successes are determined (410), the proportional term is left unaffected, the integral term is increased, and the derivative term is set to zero (412). In this case, the goal is to increase volume change proportionally to significant systematic errors.

[0405] To make the newly adjusted parameters take effect, the normal state can be applied within a set number of iterations after the adjustment. Similarly, the time period considered when determining the controller state can be modified to avoid overcompensation. In some implementations, the design problem is allowed to converge quickly in the initial optimization phase. Therefore, during this interval, only the oscillating state is adapted, while all other states are effectively ignored by the adaptive process.

[0406] Therefore, return to Figure 4A As part of calculation 424, the proportional, integral, and derivative components of the adaptive controller (e.g., a PID controller) can also be adjusted to achieve a second-level increase in shape change in response to a repetition of success or failure with respect to the next (e.g., normalized) target value exceeding a threshold. In this example, the second-level increase in shape change is greater than the first-level increase in shape change.

[0407] Adaptive controllers can also be implemented using any suitable machine learning technique to predict changes in control parameters in response to one or more inputs. For example, a controller can be implemented as a neural network with multiple neural network layers. The neural network may include an input layer that receives one or more inputs, and an output layer that outputs changes in the control parameters in response to the inputs. The neural network may include one or more hidden layers between the input and output layers, each hidden layer applying one or more nonlinear functions to the received inputs at that layer, wherein the functions are weighted according to the learned parameter values.

[0408] An adaptive controller implemented as a neural network can be trained using any suitable supervised learning technique. In the case of supervised learning, the controller can be trained on a training dataset that includes the inputs of a PID controller paired with output changes in response to proportional, integral, and derivative changes in constraint values. The training data may include a subset of all input-output pairs generated by the PID controller, for example, the first N iterations, where N is a predefined value. The training data can be further modified, for example, by adding random oscillations or other variations.

[0409] 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 the constraint changes generated by each PID controller as input and generates a final volume change that takes into account all constraints as output. In implementations using a single machine learning model, the training data for training the model can include a subset of input-output pairs as described above, with or without random perturbations, over multiple iterations.

[0410] In addition to proportional, integral, and derivative terms, the adaptive controller can be trained based on additional inputs to generate constraint values ​​or final volume changes. For example, additional inputs can come from the topology optimization process, such as features extracted from the current shape of the part being designed, or features of the part's current stress, strain, and displacement results. The adaptive controller can receive these additional inputs to learn parameter values ​​to generate more accurate control values ​​or final shape changes and further stabilize the optimization process.

[0411] 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 both. For example, if a particular optimization task fails to converge or becomes unstable, multiple instances of the optimization process can be run for all test cases but for different controller settings, such as different weighted values ​​of model parameters, different hyperparameter values ​​(e.g., learning rate, batch size, etc.), or both. For example, the setting of the instance with the best performance according to a certain performance metric can be used for the adaptive controller and the corresponding optimization process to produce the final design. Multiple instances can be run offline, for example, when the number of test cases is small, or adapted online during the optimization process.

[0412] Although neural networks are given as an example machine learning technique for implementing adaptive controllers, any suitable technique, such as fuzzy logic or any appropriate type of regression model, can be used. Furthermore, as described in this document, adaptive controllers can also have one or more additional inputs (not just measurement errors) that provide more information about the trend and current state of the optimization process.

[0413] Constraint normalization during topology optimization

[0414] As referenced above Figure 4A As described in calculation 424, controlled convergence can be combined with the normalization of constraint values. Constraint values ​​can have different orders of magnitude, leading to ill-conditioned optimization problems and gradients that are difficult to achieve for complex constraints. For example, the initial fatigue safety factor could be approximately 10,000, while the objective value could be 1. In such cases, the sensitivity to gradients decreases as the solution approaches the objective, and a final value of 10 (with an error of 9 / 10,000 relative to the initial value) can be considered acceptable due to relatively small errors. Moving the reference value improves the problem as described below in the normalization algorithm.

[0415] Normalization algorithm

[0416] Input: t i t, n, g T，n T k T i K 阈值 K 当前 K g

[0417] Output: g 参考，t

[0418]

[0419] The following points out some characteristics of the normalization algorithm.

[0420] The reference value g should not be updated in each iteration. 参考，t After setting a new reference, sufficient time should be allowed for the constraints to stabilize. This is achieved using the concept of an inner loop iteration (sections 3 through 22 of the normalization algorithm). The maximum length of the inner loop iteration is determined by... It is given that it can be dynamically adjusted for each design problem, for example, adjusted to ensure that there are at least 6 inner loop iterations when the lower limit of the inner loop length is 5.

[0421] The threshold Δg used to determine whether an existing reference is too far from the current value 阈值 Multiplier of moving average of constraint values (For example, 20) is used for calculation. The length of the moving average calculation is determined by... The threshold is given, for example, 10. If an oscillation is detected during the last inner loop, the threshold is increased; otherwise, the threshold is decreased if the current reference value is higher than the threshold.

[0422] Any change to the reference value is limited by the ratio of the existing value to that used (0 < K). g The maximum change Δg calculated for ≤1 (e.g., 0.6) 阈值 =K g g 参考，t-1 Any reduced reference value should also be greater than the current violation of |g. i -g T，n The multiplier of |. Instances of multipliers are given by K. 当前 =1.2 is given.

[0423] Figure 4G Example graph 401J shows different metrics used in the normalization algorithm as detailed above and in the tracking constraint normalization just described above. Specifically, graph 401J shows how the reference value g is updated over time as the constraint value approaches the next target value. 参考，t The final target value is g. T，n , and g t This indicates that the constraint oscillates with time to approach the final target value g. T，n The current value at that time. Reference value g 参考，t The reference value is greater than both the current and final values ​​of the constraint and does not change with each iteration. The reference value is adjusted over time to approach the final target value without exceeding the current target value.

[0424] In some implementations, the value calculated using the normalization algorithm can be assigned to the target reference value. Furthermore, in cases where a sudden change in the reference value may occur at the beginning of each inner loop iteration of the normalization algorithm, an adaptive controller can be used (e.g., as mentioned above). Figures 4D to 4F The described PID controller is used to implement the applied reference g. 参考，t A smooth transition. An example application for a given iteration t is as follows:

[0425]

[0426]

[0427] Arbitrary constraint processing

[0428] Arbitrary constraint handling is a general approach that can be applied to any type of optimization constraint, even when an exact shape derivative is not available. Whether implementing a surrogate shape derivative or an actual shape derivative, the generative design process enables greater accuracy and control.

[0429] Arbitrary equality constraints

[0430] In the optimization problem, consider equality constraints of the type g(Ω, u(Ω)) = 0, where minimizing compliance is used as the objective function. In the absence of a shape derivative... In such cases, the industry practice is to optimize while monitoring g and terminate optimization when g≈0. This means that when the solution approaches zero as g→0, the final result is a non-convergent solution due to the lack of control.

[0431] However, in the absence of a shape derivative, adaptive controllers can be used to enforce equality constraints. Assume g... i (Ω, u(Ω))j=1,...,n g Let represent a series of equality constraints that have been normalized as described above (e.g., using the normalization algorithm described above), such that Considering the normalized error of each constraint, which can be used to define the error of the constraint for each iteration in a manner similar to Equation 91:

[0432]

[0433] This can be substituted into equation 93 to approximately satisfy constraint g in iteration t. j Required volume change:

[0434]

[0435] Where the constraint is negatively correlated with the volume change, I j =-1, otherwise I gi =1. (This is a partial translation of a Chinese word that doesn't translate directly but can be left as is.) Effectively apply control to achieve response parameters The desired result. When multiple constraints exist, each constraint will recommend a different volume change required to satisfy the constraint. The constraints are grouped into positive and negative values:

[0436]

[0437] Then, the individual volume change applied in this iteration can be calculated using the following terms.

[0438]

[0439] Regarding any nonzero ΔV tThe value of μ in the shape derivative of the volume control can be applied in accordance with the method described above for reference equations 71-87 of the approximate volume control.

[0440] Arbitrary inequality constraints

[0441] The technique described above, referring to arbitrary equality constraints and equations 95-98, causes the volume change to converge to zero when it becomes satisfied with the constraints, i.e., in hour This can cause potential problems when using inequality constraints, because it is possible to... The inequality constraints must be satisfied. Therefore, when using inequality constraints, a slack variable called the importance coefficient can be used to manage the relative contribution of each constraint violation. The importance coefficient mitigates uncontrolled variations in convergence and constraints that interfere with minimizing the objective, when all constraints are assigned equal importance. The importance coefficient modulates the relative importance of different constraints during different iterations of the generative design process.

[0442] Assumption To represent the importance coefficient of each constraint, equation 96 can be modified as follows:

[0443]

[0444] The importance coefficients used are predetermined (e.g., user-provided) importance coefficients. The calculation is as follows:

[0445]

[0446] If the inequality constraint g is violated at iteration t... k ,but otherwise At each iteration, the importance coefficients applied are based on all constraints. The sign is used to update. This allows constraint violations to affect the applied importance coefficient and therefore the volume change ΔV. t .

[0447] To prevent A sudden change occurs when the constraint state changes from violated to not violated (or vice versa). For example, a PID controller can be used, as described in reference equations 96 and 97 above. Figures 4D to 4F The change in the stability importance coefficient is described as follows:

[0448]

[0449] Complex generative design problems can be solved using a combination of techniques described above, such as controlled convergence, combinations of inequality and equality constraints, and the use of PID controllers.

[0450] Modified Augmented Lagrangian Method for Constraint Handling

[0451] The augmented Lagrange algorithm described in Equations 28-30 can be modified as follows to adapt to constraint handling of arbitrary equality and inequality constraints as described above.

[0452] Controlled convergence: The classic augmented Lagrange method will constrain the violation term e. j Calculate the difference between the current value and the final value of the constraint, for example, as follows:

[0453]

[0454] For example, see the above text. Figures 4B to 4C As described, controlled constraints converge. Some constraints may converge earlier than others; therefore, the amount of change made to the design can be gradually reduced, i.e., large changes are made during the initial phase, followed by gradually smaller changes.

[0455] Constraints without shape derivatives: While the shape derivatives for many objectives and constraints can be mathematically calculated using the adjoint method (Equations 14-25), implementing the adjoint method in commercial finite element solvers can be a daunting task. Furthermore, optimization constraints can sometimes be evaluated using a user-provided black-box evaluator. Therefore, arbitrary constraint handling can be referenced above, for example. Figures 4B to 4E The described technique is used to utilize surrogate shape derivatives.

[0456] Precise volume control: Precise volume control is often crucial for obtaining good design outputs from complex engineering examples. The approximate volume control, volume control using adaptive controllers, and the μ calculation method introduced in Equations 71-93 above can be integrated into the augmented Lagrangian method, while using a line search algorithm to further improve accuracy.

[0457] First, the constraints can be divided into two groups. The former includes all constraints affected by volume changes, while the latter includes constraints that are not, for example, minimum / maximum thickness or centroid constraints. The shape derivative from the augmented Lagrange method is modified as follows:

[0458]

[0459] Where the constraint error e j Using the PID stability form of Equation 95, the importance coefficient terms described above in reference equations 99-101, derived from arbitrary inequality constraints, are used for calculation:

[0460]

[0461] When the shape derivative dg j When / dΩ is unavailable, a suitable surrogate shape derivative can be used to approximate the shape derivative. For example, from the set Any such constraint can be approximated using the volume shape derivative as follows:

[0462]

[0463] Where the constraint is negatively correlated with the volume change, I j =-1, otherwise I gi =1. Note that the derivative of the volume in the normal direction is 1.

[0464] Next, μ is calculated using the concepts described above: approximate volume control, volume control using an adaptive controller, and equations 71-93. * .use All constraints are treated as arbitrary constraints, and the target volume change Δv is calculated for each iteration t. t Next, μ is calculated using the method described above with reference to the approximate volume control and equations 71-87. * It should be noted that equation 87 needs to be modified as follows to take into account arbitrary constraints and objectives:

[0465]

[0466] In some cases, pitfalls should be considered: when When this situation is detected, μ should be set to... * Set to 1. This is typically done when the constraint converges toward the objective value, i.e., e. j →Occurs when 0.

[0467] For iteration t, update the augmented Lagrange parameter μ j , λ j An adaptive (e.g., PID) controller can be used for stabilization, as described above relative to... Figure 4F As described. For example, see μ below. j Update rules:

[0468]

[0469] As mentioned earlier, line search algorithms (e.g., gradient descent or Newton's method) can be applied to ensure that the volume change achieved after the advection of the level set is relative to the target volume change Δv. t Within a certain acceptable tolerance, a line search is performed to find the optimal multiplier l. t This causes the advection velocity to change from Provided. Note that disabling line search is equivalent to setting l. t=1.

[0470] Figure 4H Examples of convergence histories for tracking constraints are shown in graphs 414A to 414C with and without wire search. Graph 414A shows the applied volume change versus the target volume change Δv in the case of wire search. t Figure 414B shows the applied volume change compared to the target volume change Δv without line search. t Figure 414C shows the velocity multiplier l during convergence. t The history of.

[0471] Return to Figure 4A In some implementations, one or more design criteria include multiple design constraints, including a first inequality constraint and a second inequality constraint. In these implementations, the design constraints are divided into a first group and a second group. The first group contains all the multiple design constraints affected by volume changes, while the second group contains one or more of the remaining design constraints that are not affected by volume changes. Calculating the 424 shape change rate involves using, for example, an augmented Lagrangian method as described above, which applies adjustment factors to the sum of shape change contributions from the first group but not to the shape change contributions from the second group.

[0472] In some implementations, at least one design constraint that does not have a defined shape gradient includes a first inequality constraint and a second inequality constraint. The first inequality constraint has a first input control parameter for a first proportional-integral-derivative (PID) controller and a first importance coefficient multiplied by the shape change provided by the first PID controller, for example, as described above with reference to arbitrary inequality constraints and equations 99-101. The second inequality constraint has a second input control parameter for a second PID controller and a second importance coefficient multiplied by the shape change provided by the second PID controller.

[0473] The rate of shape change of the 424 implicit surface is calculated based on the surrogate shape gradient, and thus may include adjusting both the first and second importance coefficients based on whether one or more other constraints are violated in the iteratively modified previous iterations. For example, the importance coefficients can be modified by multiplying the coefficients by a multiplier for violations of other constraints in the iteratively modified previous iterations. In some implementations, and as generally referred to above, this is also possible. Figures 4D to 4F Referring to equations 99-101, proportional-integral-derivative control can be used to stabilize the adjustment of the first and second importance coefficients.

[0474] Additionally, as mentioned above... Figure 4A , Figure 4F , Figure 4G As described, the oscillation can be more than two consecutive pairs of successes or failures to satisfy the next target value for normalization. Furthermore, the repetition of successes or failures to satisfy the next target value for normalization can be repetitions occurring within 10% or more of the number of iterations, where the threshold is 10%. In some implementations, adjustments in response to oscillations include decreasing the proportional component, decreasing the integral component, and decreasing or resetting the derivative component. In some implementations, adjustments in response to repetitions of successes or failures include increasing the proportional component, decreasing the integral component, and setting the derivative component to zero. In some implementations, adjustments in response to a relative error exceeding a threshold include increasing the integral component and setting the derivative component to zero.

[0475] As mentioned above, refer to adaptive controllers and Figures 4D to 4F As described, at least one design constraint may not have a defined shape gradient. Therefore, in some implementations, computation 424 includes calculating the rate of shape change of the implicit surface based on the output of a surrogate shape gradient adjusted by adaptive control. Specifically, input control parameters are used, which are measures of the error between a normalized current value of the at least one design constraint without a defined shape gradient and a normalized next target value from the corresponding target value in a series of target values. The measure of error varies with a reference value that changes at least once during iterative modifications, which allows for controlled convergence while smoothly satisfying critical constraints.

[0476] In some implementations, 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 for a first adaptive controller, and the second equality constraint having a second input control parameter for a second adaptive controller, for example, as described above with reference to arbitrary equality constraints and equations 95-98. Then, calculating the shape change rate of the 424 implicit surface based on the surrogate shape gradient includes using the maximum shape change amount provided by the first and second adaptive controllers when neither the first nor the second equality constraint is inversely proportional to the shape change, and using the minimum shape change amount provided by the first and second adaptive controllers when neither the first nor the second equality constraint is proportional to the shape change. When one constraint is inversely proportional to the shape change and at least one constraint is proportional to the shape change, the average shape change amount is used.

[0477] After calculating 424, the shape change rate is used to update 426 one or more level set representations to produce an updated form of the 3D shape of the modeled object. 422, calculating 424, and updating 426 are repeated until check 428 determines that a predefined number of shape modification iterations have been performed, or the generated 3D shape of the modeled object in design space converges to a stable solution for one or more design criteria and one or more in-use load conditions. For example, as referenced above. Figure 3A As described, updates 426 and checks 428 can be performed. The computer-aided design program can then provide a 3D shape of a generative design of the modeled object for use in manufacturing a physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

[0478] Fatigue constraint

[0479] Next, a description of the optimized design based on the provided fatigue constraints on the design body is provided. In addition to the above (e.g., Figures 3A to 4H In addition to the techniques described herein, and the accompanying descriptions of controlled convergence, seeding, and arbitrary constraint handling, the design body can be optimized to meet fatigue constraints, which specify the total expected life of the body or the number of cycles imposed on the body under each load condition. These processes can be performed automatically, in addition to optimization procedures (e.g., generative design procedures as described above).

[0480] Two main methods for assessing this constraint are described: safe life calculation and damage tolerance method. The former concerns preventing fatigue damage by keeping stress below an acceptable threshold. The latter accepts the presence of fatigue damage and aims to ensure that fatigue cracks do not cause serious failure before a specific inspection point.

[0481] Safe life fatigue

[0482] Figure 5A This illustrates an example of a process that generates one or more parts of a 3D model of an object to be manufactured, using one or more generative design processes that solve for one or more safety-life fatigue constraints on a subject. The computer-aided design program obtains the design space of the modeled object upon which the physical structure corresponding to the 504 manufacture will be based, one or more design criteria for the modeled object, one or more in-service load conditions of the physical structure, and one or more specifications of the materials used to manufacture the physical structure. For example, as referenced above... Figure 3A If you complete it as described, you will get a 504 error.

[0483] Safe life optimization is based on the fatigue properties of materials. Test samples are used to obtain SN curves available in material databases. Designers can specify the expected number of load cycles over the service life of a part for each relevant load condition. Their cumulative effect guides the allowable design stress. Typically, stress-based safe life fatigue methods are used. However, in low-cycle loading problems, strain-life fatigue methods may be preferred.

[0484] SN curve processing.

[0485] Figure 5B An example of a SN curve 520 tracking fatigue strength (stress) on a body over multiple cycles is shown. For bodies containing multiple materials, as described in more detail below, multiple SN curves can be used.

[0486] The damage tolerance and safe-life fatigue methods described in this paper can both utilize material SN curves. These curves present computational challenges due to their potentially flat and infinite regions, and should be interpreted carefully during the optimization process.

[0487] SN curves typically place stress on the y-axis and cycles on the x-axis because typical fatigue workflows involve finding the appropriate maximum stress constraint starting from a known limit to the number of cycles. This makes the user input for the SN curve a natural [cycle, stress] point. This point definition is well-suited for damage tolerance optimization methods and therefore only requires sorting the points from lowest to highest cycles. However, in the case of the safe-life fatigue method, the optimization proceeds from stress to cycle data. The simulation of the object being designed provides stress data, for which the expected number of cycles should be found from the SN curve. Therefore, an important aspect of the method is to convert the user-provided data points into [stress, cycle] pairs and sort these pairs from lowest to highest stress. Additionally, when a durability segment is detected—a flat region in the curve at the durability stress (below which the number of cycles is practically infinite)—the set of points corresponding to that segment of the curve is reordered from highest to lowest stress. The curve reading algorithm then easily returns the reported highest cycle value instead of interpolating to infinity.

[0488] When the input data includes data for different materials, a separate expected number of load cycles is provided for each of the different materials. This involves returning the number of load cycles from one or more curves fitted to the set of data points in the plastic and elastic regions that correlate fatigue strength with load cycles; then, returning the highest number of load cycles from the set of data points in the durability region that correlates fatigue strength with load cycles.

[0489] Stress-based safe life fatigue method

[0490] Return to Figure 5A One or more design criteria include the required number of load cycles for the modeled object under each of one or more service-related load conditions of the physical structure. Input can be provided as cycles for each load condition, or as time periods optionally convertible to cycles, such as seconds, minutes, hours, days, weeks, months, years, etc. One or more specifications include data relating fatigue strength to load cycles. Data can be provided as the SN curve of the material. Data can measure stress as an indicator of fatigue stress, but in some implementations, such as when the service-related load conditions require a low number of cycles, strain is used.

[0491] The program iteratively modifies the generated 3D shape of modeling objects in the design space based on one or more design criteria, one or more in-use load conditions of the physical structure, and one or more specifications. Iterative modification may include both modifying the geometry and topology of the object's 3D shape.

[0492] For example, iteratively modifying the current numerical evaluation, which involves performing numerical simulations of the 506 modeled object based on its current 3D shape and one or more in-use load conditions to produce the physical response of the modeled object, as referenced above. Figure 3A and Figure 4A As described.

[0493] For each of one or more in-use load conditions of the physical structure, find the 508-maximized stress or strain element from the current numerical evaluation 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 referenced below... Figures 8A to 8C As described, finding 508 involves calculating the maximum stress value under load conditions in use, based at least on the standard deviation of the stress distribution in the current numerical assessment of the physical response of the modeled object.

[0494] Specifically, during topology optimization, the stress of each element can be calculated at each iteration. Optionally, this can be done for singularities (see below). Figures 8A to 8C After adjustment (as described), the maximum principal stress across all elements under a given load condition can be found, and the corresponding number of cycles can be obtained from one or more SN curves. The ratio of the obtained cycles to the required cycles represents the damage fraction. Using Miner's rule, the overall safety factor is calculated based on the damage fractions for all load conditions. Fatigue constraints must address multi-material problems. The fatigue safety factor for each material is calculated separately, and optimization is guided by the minimum value, as follows:

[0495]

[0496] in The fatigue safety factor of the material is calculated using the following methods:

[0497]

[0498] Among them, the fatigue safety factor S f C is given by the minimum fatigue safety factor of all materials used. σ It is for the maximum principal stress (max) Ω σ p The number of cycles obtained from the SN curve, and It is the required number of cycles for a given load condition lc.

[0499] During optimization, the fatigue constraint is defined as S. f -S T ≥0, where This is the target fatigue safety factor. This can be addressed using the controlled convergence, arbitrary constraint handling, and [other methods mentioned above]. Figures 4A to 4H The method described uses the shape derivative dS derived below. f / dΩ is used to enforce this. In some implementations, such as Equation 105, the volume shape derivative can be used as a surrogate shape derivative.

[0500] Shape derivative of polymer fatigue measure

[0501] Assume Ω is The domain in Ω. In this specification, stress-based fatigue metrics on Ω are independently evaluated at each point x∈Ω relative to multiple load conditions repeatedly applied over time. For topology optimization purposes, these point-by-point measurements are aggregated into a global measurement called the aggregated fatigue metric, which approximates the maximum value of the point-by-point measurements.

[0502] The polymer fatigue metric is presented in the following form. Assumption C: It is the inverse of the SN curve for materials including Ω, and assumes This is the reference value for the l-th load condition. This aggregate fatigue metric allows for the summation of damage caused by cyclic loading. The target value of 1 / (safety factor) is optimized, where the damage the part can withstand is summed with the actual damage from the load (stress or strain and cycle number).

[0503] The polymer fatigue measurement of Ω is defined as:

[0504]

[0505] Where σ l (Ω) is an approximation of the maximum value of the stress tensor in Ω under the l-th loading condition. In other words,

[0506]

[0507] Where u l : It is the displacement function that satisfies the linear elastic equation in Ω with respect to the l-th load condition, and σ(u) l ) is the associated stress tensor.

[0508] The shape derivative of the polymer fatigue metric S will be described in more detail below.

[0509] Proposition 1. The shape function S is relative to the normal velocity Θ: The shape derivative of the variant of the generated shape Ω is given by the following terms:

[0510]

[0511] Where λ t : This is the associated displacement function for the l-th load condition. This is the solution to the linear elastic equation, where the associated force term is given in weak form as a linear integral:

[0512]

[0513] Since S can be decomposed into ∫ applied to the stress integral Ω ||σ(u l )|| 2p The series of operations is broken down as follows. First, the derivatives of this series of operations are calculated using the chain rule from ordinary calculus. Then, the shape derivative of the stress integral is calculated, for which the Céa method is used and is described below.

[0514] Assume Ω ε This represents a variation of Ω generated by the normal velocity function Θ. It is represented as DS. Ω The calculation of the expected shape derivative of Θ begins with:

[0515]

[0516] First, the chain rule from ordinary calculus is used to bring the derivative d / dε to the stress integral, where a shape differentiation technique is then introduced, as follows:

[0517]

[0518] Where u l，ε In the domain Ω ε The displacement function for the l-th load condition.

[0519] Next, shape derivatives are used. The Céa method is used to calculate the shape derivative of the remaining stress integral in Equation 115. The subscript 1 is omitted in the following description because it is the same for all loading conditions. The Lagrange quantities are as follows:

[0520]

[0521] in

[0522] A(Ω, u, λ):=∫ Ω σ(u):e(λ) (117)

[0523] It is the integral form of the elastic equation in the weak form of Ω.

[0524] Here, σ(u) is the stress tensor of displacement u, and e(λ) is the virtual strain of the virtual displacement λ. Therefore, This is the "virtual work" associated with that pair of displacements. Additionally, It is the integral linear form of the elastic equation in the weak form of Ω, which encodes the "virtual work" done by the applied body and boundary traction forces. The Céa method is now performed in three steps.

[0525] Step 1. Setting the variation of the Lagrange relative to λ to zero generates a linear elasticity equation satisfied by u. This is intentionally designed because... It is exactly the left-hand side of the weak elasticity equation.

[0526] Step 2. Assume that the Lagrange variable with respect to u generates a set of related equations called the adjoint state equations with respect to λ. The solutions to these equations are called the adjoint states. The weak form of these equations is: for all variations of u...

[0527] 0=∫ Ω 2p||σ(u)|| 2p-1 σ(u): σ′(δn)+A(Ω, u, λ). (118)

[0528] Here, σ′(δn) is the derivative of the stress tensor with respect to displacement. This has a simple form because stress is linearly related to strain, which in turn is linearly related to u—the gradient is a linear operator, or rather simply σ′ = σ. Also used The symmetry of λ. The above proof demonstrates that λ satisfies the linear elasticity equation, but there exists a new "accompanying force" term encoded in integral linear form:

[0529]

[0530] Step 3. The Céa method depends on the fact that, as can be proven, the shape derivative of the Lagrangian calculated by neglecting the shape dependencies of u and λ and then inserting the states of u and the adjoint states of λ is equal to ∫ Ω ||σ(u)|| 2p The shape derivative of itself. Since all expressions in a Lagrange are volume integrals, we only need the formula for the shape derivative of such integrals for shape-independent functions over Ω.

[0531] Next, it is assumed that Θ vanishes on the Dirichlet and within the homogeneous Newman boundary of Ω (otherwise, the shape derivative would include other terms). Finally, it is assumed that there are no body force terms, for example, by neglecting the effect of gravity under a given loading condition. In summary, these last two assumptions have the effect of removing the shape derivative formula (Equation 116) from the equation. The effect of the term. Therefore, the result after all assumptions is as follows:

[0532]

[0533] Where u is the state and λ is the adjoint state.

[0534] Strain-based safe life fatigue method

[0535] Strain-based safe life calculations also use equations 108 and 109, but C σ The parameter has changed to C ε The strain and cyclic form of the SN curve are used as input. According to Neuberger's law, the strain is calculated from the maximum principal stress (max...). Ω σ p The strain value is calculated as follows:

[0536]

[0537] Where K is the concentration factor specified by the user (default is 1) and E is the Young's modulus of the material.

[0538] Return to Figure 5A As mentioned above, refer to the SN curve processing and Figure 5B As described, the expected number of load cycles for each of one or more service-related load conditions of the physical structure is determined using maximized stress or strain elements and data that correlates fatigue strength with load cycles. The fatigue safety factor inequality constraint for the modeled object is redefined based on the damage fraction calculated according to the required number of load cycles for the modeled object and the expected number of load cycles for each of one or more service-related load conditions of the physical structure. In some implementations, the inequality constraint can be normalized and tuned using an adaptive controller, for example, as referenced above. Figures 4D to 4F As described.

[0539] In some implementations, one or more specifications include two or more specifications for the corresponding different materials used to manufacture the physical structure. Therefore, the data regarding the specifications includes data relating fatigue strength to load cycles for each of the different materials. Thus, determination 510 includes determining a separate expected number of load cycles for each of the different materials. Furthermore, redefining 512 includes calculating a separate fatigue safety factor for each of the different materials based on a corresponding damage fraction calculated from the corresponding expected number of load cycles for each of the different materials. Based on the separate safety factors, the fatigue safety factor inequality constraint of the modeling object is redefined using the minimum of the fatigue safety factors for the different materials.

[0540] In some implementations, one or more service load conditions of the physical structure include two or more service load conditions of the physical structure, and one or more design criteria include the required number of load cycles for the modeled object under each of the two or more service load conditions of the physical structure. Calculating a separate safety factor for each corresponding material in the different materials involves summing load-specific damage fractions corresponding to the two or more service load conditions, where each load-specific damage fraction includes dividing the expected number of load cycles for one of the different materials and one of the service load conditions by the required number of load cycles for that service load condition. The separate safety factor is obtained by taking the reciprocal of the sum of the load-specific damage fractions. See, for example, Equation 109.

[0541] The rate of shape change of the implicit surface in the level set representation of the 514 three-dimensional shape is calculated at least based on the fatigue safety factor inequality constraint. The rate of shape change can be calculated relative to other constraints, and surrogate derivatives can be computed where shape derivatives are unavailable or not well-defined, for example, as described above in the handling of arbitrary constraints and... Figures 4A to 4H As described above. Then, the shape change rate is used to update the 516 level set representation to produce an updated form of the 3D shape of the modeled object. The shape derivative of the aggregate fatigue metric described above with respect to equations 110-120 can be used to calculate the shape change rate of the fatigue safety factor inequality constraint.

[0542] In some implementations, calculation 514 includes calculating at least one rate of shape change using a quantity determined according to a shape derivative formula, which approximates the shape derivative of the fatigue safety factor, for example, the volumetric shape derivative described above with reference to equation 105. The formula is modified to include the stress shape derivative instead of the volumetric shape derivative.

[0543] In some implementations, an importance factor is used to address at least one design constraint, including a fatigue safety factor, for example, as described above relative to... Figure 4A The importance coefficients are adjusted, generally relative to the arbitrary inequality constraint treatment and as described in equations 53-60, based on whether one or more other constraints were violated in the iteratively modified previous iterations.

[0544] In some implementations, an adaptive controller is used to adjust the target value of the volume fraction, for example, as referenced above. Figures 4D to 4F As described. Specifically, in multiple iterations of iterative modification, the target value is adjusted between the initial target value and the final target value based on inequality constraints of volume fraction or minimum thickness. When adjusting the target value in multiple iterations, adaptive control is used to stabilize the changes made to the quantity determined according to the shape derivative formula.

[0545] In some implementations, adaptive control is used to adjust the overall contribution of the quantity determined according to the shape derivative formula to the rate of shape change used in the update. The adaptive controller can be used to adjust the contribution of the surrogate shape derivative to the total advection velocity, i.e., according to a modified form of the augmented Lagrangian method used for constraint handling, as described above with reference to equations 102-107.

[0546] Repeat steps 506, 508, 510, 512, 514, and 516 until check 518 determines that a predefined number of shape modification iterations have been performed, or that the generated 3D shape of the modeling object in the design space has converged to a stable solution for one or more design criteria and one or more in-use load conditions. Check 518 can be as described above. Figure 3A and Figure 4A The described checks. Then, a 3D shape of the generative design of the modeled object can be provided, for example, to manufacture the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems. These techniques can be combined with hybrid, hollow, and hybrid-hollow methods for topology optimization.

[0547] Damage tolerance fatigue method

[0548] As described above, safety-life fatigue constraints help ensure that an object does not fail under user-specified load conditions by preventing fatigue damage during a specified number of load cycles. In contrast, damage-tolerant fatigue techniques focus on limiting fatigue damage within permissible limits to help ensure that fatigue damage is detected in the field before part failure. A key objective is to calculate the critical fatigue crack length and ensure that it does not exceed the thickness of the designed part during the service inspection interval. Therefore, a critical requirement for this approach is the ability to enforce thickness constraints on topology-optimized designs. Thickness constraints can be implicitly enforced, for example, by increasing the volume of the design, resulting in an increase in the thickness of all parts of the design.

[0549] It should be noted that, similar to the safe-life fatigue method, the damage-tolerant fatigue technique described below can be used in any combination with the techniques described above, such as controlled convergence, seeding, and arbitrary constraint handling. Furthermore, in some implementations, both safe-life and damage-tolerant fatigue methods are available, where fatigue design constraints for a given problem are addressed in some cases using the safe-life method, while in others, fatigue design constraints are addressed using the damage-tolerant fatigue technique described below. The appropriate method can be selected based on the intended use of the part being designed. For example, in aerospace applications, parts are typically designed using the damage-tolerant method. The availability of both modes together allows for the solution of a wider range of design problems.

[0550] Figure 6A An example is shown of a process for iteratively modifying the 3D shape of a modeling object in a design space according to one or more design criteria, including at least one damage tolerance fatigue constraint. The 3D shape comprises a level set representation of implicit surfaces, and the one or more design criteria include the required number of load cycles for the modeling object under each of one or more load conditions in use of the physical structure.

[0551] Numerical simulations of the 608 modeling object are performed based on its current 3D shape and one or more in-use load conditions to produce a current numerical assessment of the object's physical response (e.g., structural response), as referenced above. Figure 3A , Figure 4A and Figure 5A As described.

[0552] The expected number of load cycles for each of one or more in-service load conditions in the 610 physical structure is determined using current numerical evaluation and thickness measurements to enforce design criteria limiting the minimum thickness. Points with maximum strain or stress from the current numerical evaluation can be used to determine the expected number of load cycles. The calculation of the current thickness and the enforcement of the minimum thickness are described below.

[0553] Thickness constraint

[0554] First, thickness constraints are calculated and enforced. Thickness can be measured and enforced using one or more of various methods, such as a combination of at least two different thickness measurements, as described below. Furthermore, it should be understood that the described damage-tolerant fatigue technique can be implemented via density-based topology optimization (e.g., using the SIMP method) rather than boundary-based topology optimization (e.g., level set methods), such as in combination with... Figure 6A As described.

[0555] Figure 6B Graphical representations of geometries 600A and 600B are shown, where the thicknesses of geometries 600A and 600B are calculated using different measurement techniques. Geometries 600A show the body 601 measured using a ray casting method. The ray casting thickness at point x in geometry 600A is defined as h. r (x)=|xx r |, where x r This refers to the projection of a ray of light along the negative direction of vector n at point x onto a point in geometry 600A. The ray projection can be repeated for different directions of vector n at point x, including the case where vector n is perpendicular to geometry 600A.

[0556] In contrast, a spherical fitting method is shown in geometry 600B, which involves finding the diameter of the largest sphere that can be fitted within the domain while contacting the measurement point x in geometry 600B. This is accomplished by running a bisection algorithm to find the relationship between point x and x... r The spherical centroid x s Discrete sampling positions are defined on the periphery of the sphere / surface using the polar angle α. With origin x s A sphere is considered to fit within the domain when its maximum level set with respect to the domain is below a specific threshold Δh, i.e., when the following condition is met:

[0557]

[0558] The thickness after spherical fitting is defined as h. s (x)=2|xx s In some implementations, the thickness is defined as h(x = max{h... s (x), h r (x)}. If thickness is required at points not on the domain surface, then before measuring the thickness on the surface, first project those points onto the surface using the normals of the level set. The spherical fitting method can be repeated multiple times for multiple directions (generating multiple x). r (The point depends on the chosen orientation). The calculated maximum sphere can be used to define the thickness.

[0559] Minimum thickness constraint g t It can be defined as

[0560]

[0561] Where Γ is the domain surface and h 最小 This is the minimum thickness of the target. The approximate shape derivative is given by the following terms:

[0562]

[0563] Thickness constraints can be handled using the methods described above for arbitrary constraints and Figures 4A to 4H The described method is applied. In some implementations, the importance coefficient regarding the thickness constraint is determined via the sigmoid function S. r (ξ) is modified so that the thickness constraint only takes effect during the later stages of optimization.

[0564]

[0565] The sigmoid function is defined as:

[0566]

[0567] Where ξ is calculated as follows

[0568]

[0569] in Identify the initial application of thickness constraints, which is considered as volume reduction iteration n. v The maximum iteration and the first violation of the thickness constraint in the process.

[0570] Figure 6C Examples of graphs 601A to 601E tracking the sigmoid function at different rates are shown. Curve 601A is the sigmoid function at rate 1; curve 601B is the sigmoid function at rate 2; curve 601C is the sigmoid function at rate 3; curve 601D is the sigmoid function at rate 4; and curve 601E is the sigmoid function at rate 5.

[0571] Although two techniques for measuring thickness are described herein, it should be understood that any suitable technique for measuring the thickness of a subject may be used without loss of generality. Furthermore, the described techniques can be combined (i.e., a combination of two different thickness measurements) to improve the results. For example, the two different thickness measurements may include (i) a first distance measurement of the length within the modeling object from a surface point of the modeling object projected in the negative normal direction, and (ii) a second distance measurement of the diameter of the largest sphere fitted within the modeling object and contacting a surface point of the modeling object, as referenced above. Figure 6B As described, it is determined by examining discrete sampling positions defined on the curved surface of a sphere.

[0572] Damage tolerance fatigue constraint

[0573] Next, the damage tolerance fatigue constraint is described. The critical fatigue crack length h of the material is discussed. d Defined as

[0574]

[0575] Where C is the fatigue cycle number, ρ g This is the modulus of the fatigue crack growth curve of the material and can be specified by the user, for example, 2.8. 'a' is the initial crack length, which can be set to be equal to the minimum defect detection size of the selected part inspection method; the default value is a = 1 mm. 'I' is the stress intensity coefficient, for example, the default value is I = 1.22 x 10⁻⁶. -12 And Y describes the crack type, for example, the default value is Y = 1. The stress range Δσ is as follows based on the stress ratio R, for example R = 0.1 and the maximum stress σ. 最大 Derivation:

[0576] Δσ=(1-R)σ 最大 (129)

[0577] The minimum thickness of the material is then given by the following item.

[0578]

[0579] in It is the desired critical thickness to width ratio, and S f It is the specified design safety factor. The multiplier and Y-type crack type variables are related and can be found in engineering tables describing the fatigue crack properties of a given material.

[0580] Figure 6D An example of Table 601F, which describes the properties of fatigue cracks, is shown. In Table 601F, information about the corresponding... The multiplier indicates the crack type variable Y. Y is the geometric correction factor used in the stress intensity coefficient expression. The correction factor describes the relationship between the crack length a and the characteristic thickness W. Table 601F shows the values ​​of Y and a / W for different materials. Table 601F can be user-specified, and specified values ​​for Y and the ratio a / W exist for different materials.

[0581] The crack length 'a' is set to be equal to the critical crack length 'ac', which tends towards infinity, i.e., severe part failure. During optimization, this gives the minimum thickness to be optimized; if the user expects an a / W ratio of 0.1 (i.e., the minimum thickness should be 10 times ac), then the value of ac is calculated and the 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 part thickness compared to the critical crack length in the part's features. This has the effect of increasing the critical crack length in the features, making it easier to detect by the selected inspection technique.

[0582] One way to add optimization constraints to meet damage-tolerant fatigue requirements is to use a constant global thickness target. This assumes the design is at its fatigue limit, where C equals the required number of fatigue cycles and σ... 最大 This is the corresponding stress from the SN curve. Equation 128 is used to calculate the critical crack length for each material. The minimum thickness constraint (Equation 123) is modified, where the target thickness is set to the result of Equation 130 using the worst-case crack length. The target thickness remains constant throughout the optimization process.

[0583]

[0584] Alternatively, a damage-tolerant fatigue safety factor can be calculated for each point x on the designed surface at each iteration. This uses the stress σ at each point. lc The number of cycles C supported for each load condition is calculated using (x) and thickness h(x). lc (x):

[0585]

[0586] Damage tolerance and fatigue safety factor S of the material d Now, similar to the safe life fatigue method, the Manner rule can be used to calculate:

[0587]

[0588] Where C lc (x) is the cycle number obtained from equation 132, and It is the required number of cycles for a given load condition lc.

[0589] During optimization, the fatigue constraint is defined as S. d -S T ≥0, where It is the target fatigue safety factor—the lowest safety factor among all material safety factors. This can be used in the above discussion regarding arbitrary constraint treatment and... Figures 4A to 4H The described method enforces this. As described above, the volumetric or thickness shape derivative (Equation 124) can be used as a proxy shape derivative. To use the thickness shape derivative, the following additional steps should be taken: use Equations 128 and 130 to apply the derivative corresponding to the critical safety factor. Point stress and thickness data are converted into thickness target h. 最小 It should be noted that, as referenced above... Figure 6C As described, the importance coefficient should be multiplied by a sigmoid type function.

[0590] Figure 6E An example is shown of a process that generates one or more parts of a 3D model of an object to be manufactured by using one or more generative design processes that solve for damage tolerance fatigue constraints on the subject.

[0591] Computer-aided design programs obtain the physical structure corresponding to the 602 manufacturing process based on the design space of the modeling object, one or more design criteria of the modeling object, one or more in-service load conditions of the physical structure, and the critical fatigue crack length of the material used to manufacture the physical structure. For example, as referenced above. Figure 3A , Figure 4A and Figure 5A As described, an executable 602 can be obtained. In some implementations, obtaining the critical fatigue crack length of the material includes obtaining one or more specifications of the material used to manufacture the physical structure, and based on information in one or more specifications (e.g., as referenced above). Figure 5B and Figure 6D The described SN curve and engineering table are used to calculate the critical fatigue crack length of the material. This information may include the modulus of the material's fatigue crack growth curve.

[0592] The program iteratively modifies the 3D shape of the generated design of the modeling objects in the design space based on one or more design criteria, one or more in-service load conditions of the physical structure, and one or more specifications. For example, as referenced above. Figure 3A , Figure 4A , Figure 5A and Figure 6AAs described, iterative modification may include modifying both the geometry and topology of the object's three-dimensional shape. Iterative modification includes enforcing a design criterion of a minimum thickness for the generative design of the object's three-dimensional shape, constrained by 604 constraints. This minimum thickness is based on the material's critical fatigue crack length. The minimum thickness may be a ratio to the critical fatigue crack length, such as a minimum thickness that is 10 times the critical fatigue crack length. For example, as referenced above... Figure 6B As described, the thickness can be measured using various techniques.

[0593] Provides 3D shapes for generating designs of 606 modeling objects, for use in manufacturing physical structures corresponding to the modeling objects using one or more computer-controlled manufacturing systems.

[0594] As mentioned above, this involves handling arbitrary inequality constraints and... Figures 4A to 4H As described, in some implementations, enforcing minimum thickness inequality constraints may include using inequality constraints based on volume fraction or minimum thickness as a proxy for design criteria limiting the minimum thickness. Using the volume shape derivative results in fatigue safety being correctly satisfied, but in some implementations, the thickness shape derivative is substituted for the volume shape derivative. The volume fraction-based or minimum thickness-based inequality constraints are modified using an importance coefficient, which is set to zero during the initial phase of iterative modification and adjusted during subsequent phases of iterative modification based on whether one or more other constraints were violated in previous iterations of the iterative modification.

[0595] Additionally, as mentioned above... Figures 4D to 4F As described, PID control can be used to adjust the target value based on inequality constraints of volume fraction or minimum thickness between an initial target value and a final target value over multiple iterations of iterative modification. When adjusting the target value over multiple iterations, PID control can also be used to adjust and stabilize changes made to the amount of modification to the modeling object determined by the evaluation of inequality constraints based on volume fraction or minimum thickness.

[0596] Return to Figure 6A The fatigue safety factor inequality constraint of the 612 modeling object can be redefined based on the damage score calculated according to the required number of load cycles of the modeling object and the expected number of load cycles of each of one or more in-use load conditions of the physical structure.

[0597] In some implementations, one or more in-use load conditions of the physical structure include two or more in-use load conditions of the physical structure, and one or more design criteria include the required number of load cycles for the modeled object under each of the two or more in-use load conditions of the physical structure. In these implementations, determining the expected number of load cycles (610) includes determining the individual expected number of load cycles for each of a plurality of points (e.g., points identified as critical points using one or more techniques). Alternatively, as detailed in Equations 132 and 133 for the implicit surfaces under each of the two or more in-use load conditions, 610 is determined for each point x on the surface of the design at each iteration.

[0598] In some implementations, redefining the 612 fatigue safety factor inequality constraint involves summing load-specific damage fractions corresponding to two or more in-service load conditions for each of a plurality of points. Each load-specific damage fraction involves dividing the expected number of load cycles for one of the plurality of points and one of the in-service load conditions by the required number of load cycles for one of the in-service load conditions to produce a sum of load-specific damage fractions for each of the plurality of points. The reciprocal of the sum is taken, and the minimum of the sum of reciprocals is used to redefine the fatigue safety factor inequality constraint for the modeled object. See, for example, Equation 109.

[0599] The shape change rate of the implicit surface 614 is calculated at least based on the fatigue safety factor inequality constraint. Then, the shape change rate is used to update the level set representation 616 to produce an updated form of the 3D shape of the modeled object. In some implementations, calculating 614 includes calculating at least one shape change rate using a quantity determined according to a shape derivative formula that approximates the shape derivative of the fatigue safety factor.

[0600] Repeat steps 608, 610, 612, 614, and 616 until check 618 confirms that a predefined number of shape modification iterations have been performed, or until the generated 3D shape of the modeling object in the design space has converged to a stable solution for one or more design criteria and one or more in-use load conditions, for example, as referenced above. Figure 3A As described.

[0601] Designing digital twins

[0602] The relationship between manufactured and designed parts can be realized as a digital twin, reducing the discrepancy between the intended design and the physical part. This is particularly prevalent in additive manufacturing, where the effects of thickness and build-up angles can significantly impact part performance. This is difficult to consider during the design process, especially when applying other optimization techniques, such as those mentioned above. Figures 3A to 6E As described. Therefore, the known relationship between the thickness of a feature and its construction angle can be applied based on each element, whereas conventional techniques assume equal strength across the entire part.

[0603] Figure 7A This is an example of a process that iteratively modifies the three-dimensional shape of a modeling object in the design space according to one or more design criteria, including stress constraints. A 708 numerical simulation of the modeling object is performed based on the current form of the three-dimensional shape and one or more in-service load conditions to produce a current numerical evaluation of the physical response of the modeling object. The physical response can be the structural response of the modeling object under one or more in-service load conditions. A 708 numerical simulation can be performed, for example, as referenced above. Figure 3A As described.

[0604] For each of one or more in-service load conditions of the physical structure, find the maximum stress or strain element from the current numerical evaluation of the physical response of the modeled object. The maximum stress or strain element can be the stress, strain, or both at a point, location, or region of the modeled object. (See above reference...) Figure 5A and Figure 6A The descriptions of fatigue constraints and computational stress constraints are presented.

[0605] Dependence on feature size and construction angle

[0606] Additively manufactured parts (e.g., metal parts) exhibit different strengths for different build angles and part thickness values. The Von Mises stress target at point x is modified to be a function of thickness h(x) and build angle β(x), which is measured relative to a constant build direction. Thickness can be determined by including the reference above. Figure 6B Any suitable technology or combination of technologies described can be used for measurement.

[0607] Figure 7B A graphical representation of an example of a build angle measured on object 700 is shown. The build angle at point x on the surface is represented using... The construction angle is calculated by considering n(x) as the normal vector of the surface at x. The interior position of the object is a position not on the surface of the object. For these positions, the construction angle can first be calculated by projecting the interior position along a predetermined construction direction onto the surface position of the 3D shape. Then, the construction angle of the interior position is determined using the formula β(x) mentioned above, based on the angle between the normal of the surface position and the predetermined construction direction.

[0608] Von Mies stress constraints are used and take the following form:

[0609]

[0610] Where σ T It is a fixed stress target for each material. Now this will vary with thickness and build angle, resulting in the stress constraint being modified to...

[0611]

[0612] The variability in build angle and thickness is shifted from the stress target to the stress itself. For convenience, while maintaining the same safety factor:

[0613]

[0614] The above reference can be used for arbitrary constraint processing and Figures 4A to 4H The described method enforces constraints, where the shape derivative of the combined distance-normal-curvature penalty is used and derived below. In some implementations, such as Equation 105, the volume shape derivative can be used as a surrogate shape derivative.

[0615] Next, we present the derivatives of the shape derivatives of the distance-normal-curvature penalty. First, we present the derivations of the shape derivatives of the signed distance function, the unit normal function, and the surface curvature function, respectively. Then, we present the derivation of the combined penalty of these three shape derivatives.

[0616] Shape derivative of the signed distance function

[0617] Assumption It is a domain with boundaries. In this specification, as a mapping of The distance function is defined by the following terms:

[0618]

[0619] For each set Define a set of y ∈ Ω in which the above minimum value is achieved. When a point consists of exactly one point, that point is called the point from x to x. The nearest point projection on and by This is indicated. Additionally, the sign is assigned to the distance function as follows: for x∈Ω, And for The signed distance function can now be defined as

[0620]

[0621] It can be easily seen that It is 1-Lipschitz, meaning its difference quotient is uniformly bounded by 1. Therefore, by Rademacher's theorem, it is differentiable almost everywhere, and regardless of... Regardless of its nature, this applies. The central axis is defined as a set The square of the distance function is not differentiable. exist The above is never differentiable, and when When it is a differentiable surface, At least in The axis is differentiable. Therefore, the central axis uses... To define to exclude This includes all points except outliers. By definition, the Lebesgue measure of S is zero. It can be confirmed that the point... The set of is a subset of S, and for the set , It consists of at least two distinct points.

[0622] To support the following description, the following facts about differential geometry are presented. For ease of description, it is assumed that the boundary of Ω is "reasonable" in the sense that it has a tangent and that the tangent is differentiably varied such that... The curvature can be defined using the tangent derivative of the unit normal vector field. It is also assumed that... For at least C 2 This means It can be expressed as C 2 The union of finite sets of the graphs of the functions.

[0623] Assumption yes The unit outward normal vector field. Assume k1(y) and k2(y) are in... relative to The principal curvatures, and assume E1(y) and E2(y) are at y. The cross-section is an orthogonal basis aligned with the principal curvature direction at y.

[0624] Proposition 1. Suppose that Ω is a molecule with C 2 A domain with a boundary. Imagine. Then, It is differentiable twice at x. Furthermore, if but It is also differentiable twice at x. If and only if... At time, point exist It has a unique projection. In this case, and The gradient at x is

[0625]

[0626] For each and each The following are true:

[0627]

[0628] Furthermore, the closure of S consists of all points in S along with all points that strictly obey at least one of the above inequalities.

[0629] Suppose that x is not within the closure of S. Then

[0630]

[0631] and

[0632]

[0633] Assume T ∈ : It is used as a vector field Θ: The flow appeared A family of single-parameter transformations. This means: T ε Satisfy the following equation:

[0634]

[0635] T0 = ​​Identity (144)

[0636] Recall that the shape-dependent function f Ω : The Lagrange derivative with respect to this change is defined as:

[0637]

[0638] And f Ω The Euler derivative with respect to this change is defined as:

[0639]

[0640] The following results give the Euler derivative of the signed distance function.

[0641] Proposition 2. If x is not a closure of S, then the Euler derivative of the signed distance function with respect to the change generated at x by the vector field Θ satisfies:

[0642]

[0643] The proof is as follows. First, calculate the Lagrange derivative of the signed distance function with respect to the change generated at x by the vector field Θ. For this... Use (T) ε (x))=||T ε (x)-y ε ||, where y ε In all Minimize ||T ε (x)-y||. The result is...

[0644]

[0645] This is achieved by using conditions It is obtained by differentiation. Similarly, it can be verified that the gradient term in the definition of the Euler derivative is... The expected formula is as follows.

[0646] Consider the following form of shape function:

[0647]

[0648] Where φ is a smoothing function of its independent variable. The shape derivative of φ is calculated using standard tools of shape calculus and the formula for the Euler derivative with a signed distance function.

[0649] Proposition 3. Assume T ε : It is used as a vector field Θ: The flow appeared A family of single-parameter transformations. Then...

[0650]

[0651] Where Θ ⊥ Is Θ in The normal component on, and It is the partial derivative of φ with respect to its second independent variable. If the Euler derivative of a quantity is expressed in terms of a prime number, then...

[0652]

[0653] It is observed that the first integral appearing in the above formula is not in the standard form of the Hadamard-Zolésio structure theorem, that is, in The surface integral on the surface. The formula is in It does indeed still depend on Θ ⊥Therefore, the theorem still holds. Based on the formula for the complementary area, this integral can be expressed in standard form by using variable changes. The following results should be attributed to Dapogny et al., Geometric Constraints for Shape and Topology Optimization in Architectural Design, Computational Mechanics, Springer Verlag, (2017, 59(6), pp. 933-965).

[0654] Proposition 4. Suppose that Ω is a molecule with C 2 The domain of the boundary. Then

[0655]

[0656] Where T(y) is along the path from... The light radiating inwards The distance from the midline, and J(y, τ) is the Jacobian factor of the variable change, which is...

[0657]

[0658] It should be noted that the Jacobian factors above have particularly good alternative forms. Expanding the product yields:

[0659]

[0660] in and They are The mean curvature and Gaussian curvature. This is because these are... The second fundamental form of invariants, where the principal curvatures are eigenvalues.

[0661] Shape derivative of unit normal vector function

[0662] Suppose that ∑ is any directed smooth surface embedded in Euclidean space, and suppose that Y: It is a smooth vector field defined on the background Euclidean space and completely independent of ∑. Consider a shape function of the following form:

[0663] Φ(∑):=∫ ∑ φ(x, <Y(x),N ∑ (x)>dσ(x)) (155)

[0664] Where φ: It is a smooth function of its independent variable. The object appearing in the above integrand is: the unit normal vector N of ∑. ∑Let ∑ be the surface area element dσ. The shape derivative of Φ with respect to any given shape variation is calculated. This means: assuming T ε (∑) is a transformation T of Euclidean space ε The form of ∑ generated by a family of single parameters; then the shape derivative with respect to this form is

[0665]

[0666] To calculate the shape derivative of Φ, we first describe the transformation T in more detail. ε There are many ways to proceed, and all of them are equivalent to the first order of ε. The Hadamard velocity method has been chosen in the following description because it is well-suited for shape optimization based on level sets; however, it should be understood that other methods, such as those described in MCDelfour and JPZol'esio, Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization (2nd edition, 2011), may be used.

[0667] Imagine T ε : It is used as a vector field Θ: The flow appeared Transformation. In other words, T ε Satisfy the following equations

[0668]

[0669] T0 = ​​identity.

[0670] It is the tangential gradient operator, div || It is the tangential divergence operator, and H ∑ It is the average curvature of ∑. Additionally, φ(x, <Y(x),N ∑ (x)>) is written as φ(x, q(x)) (where q(x) :=<Y(x),N∑(x)> To emphasize the structure of the integrand, and the partial derivative of φ with respect to its second independent variable is expressed as...

[0671] Proposition 1. Suppose that ∑ is a directed, smooth surface embedded in Euclidean space, and assume that T ε It is composed of the vector field Θ: A variant of ∑ is generated. The shape derivative of Φ at ∑ is given by the following terms.

[0672]

[0673] Where Θ ⊥ It is the normal component of Θ on ∑, and Let represent the all-directional derivative along the normal direction. The proof is as follows. Applying the calculus of Euler's derivative, we obtain:

[0674]

[0675] The prime number represents the Euler derivative.

[0676] Because Y is independent of shape, its Euler derivative vanishes, and Therefore, after integrating in parts,

[0677]

[0678] Consider the specific application of the formula derived in Equation 160. That is, apply the formula for the function φ from the paper by Dapogny et al., i.e.

[0679] φ( <X,N ∑ >): =||YN ∑ || 2 =2-2 <X,N ∑ > (161)

[0680] Where Y is a smooth, shape-independent unit vector field defined in Euclidean space. Since its independent variable is the function φ(x, q) = 2 - 2q, then... It is a constant function Therefore, in this case, proposition 1 arises and:

[0681]

[0682] At this point, the identity that is valid for any vector field in the background Euclidean space and is bounded by ∑ is used. therefore,

[0683]

[0684] This formula is almost identical to the one in Dapogny et al.'s paper. There are two differences: first, Dapogny's paper only uses the tangential divergence operator; second, the mean curvature appears with the opposite sign to that in Dapogny's paper. The first difference is likely due to an error, while the second difference stems from the fact that Dapogny's paper uses the opposite sign convention for the mean curvature. In the paper, H... ∑It is defined relative to the inward unit normal vector, while here it is defined relative to the outward unit normal vector. (The former convention has the advantage that the mean curvature of the sphere is positive and equal to 2, while the latter convention has the advantage that the theoretical formula generally contains fewer negative signs, but in fact the mean curvature of a unit sphere is -2.)

[0685] It can be verified that a certain level set function F: of The average curvature of the level set of the form is equal to in It is the unit normal vector field of the level set. Therefore, if Y in the formula of Proposition 1 has this form, the shape derivative can be interpreted as the difference between the mean curvature of ∑ and the mean curvature of the level set.

[0686] The second example involves a 2.5D manufacturability penalty. In this case, the penalty function Φ is used as before, but in the presence of an integrand, the function is...

[0687] φ(x, q(x)) := [q(x)] 2 E(x) (164)

[0688] Where, as before, q(x): = <Y(x),N ∑ (x)>, but there is the following important difference: -Y actually represents a constant vector field in the milling direction, therefore The function E(x) is the product of the offset exponential impulse functions. Therefore, the result is...

[0689]

[0690] Because Y is a constant.

[0691] Shape derivative of surface curvature function

[0692] Assume ∑ is Consider a smooth, directed surface in the region. Consider a shape function of the form that depends on the surface curvature invariants (i.e., mean curvature and total curvature) ∑:

[0693] Φ(∑):=∫ ∑ φ(H ∑ ,||A ∑ || 2 )dσ (166)

[0694] Where H ∑ The mean curvature of ∑, ||A ∑ || 2 It is the total curvature of ∑, and φ is a smooth function of its independent variable.

[0695] The following is a description of how to calculate the shape derivative of this form of shape function. In other words, for any variant T... ε (∑) shows how to calculate the derivative.

[0696]

[0697] Where T ε : It is used as a vector field Θ: The flow appeared A family of single-parameter transformations. This means: T ε Satisfy the following equations

[0698]

[0699] T0 = ​​identity.

[0700] The primary geometric object associated with the surface ∑ is the induced surface metric tensor h. ∑ (That is, the first fundamental form), the unit normal vector field N ∑ Second basic form A ∑ The above describes the calculation of the derivatives of these components with respect to the variants introduced above (rather than the calculation of the derivatives of their components with respect to the local basis of the tangent bundle of ∑, since these quantities are tensors).

[0701] To facilitate computation, some preliminary concepts are introduced within the context of embedded surfaces in Riemannian geometry. First, local parameterization of ∑ is used in the vicinity of any point x∈∑ to facilitate mapping. We parameterize the sufficiently small tubular neighborhood of ∑ around x. Under this parameterization, the ε direction is mapped to the vector field Θ. We introduce a pair of vector fields. The pair of vector fields are sufficiently linearly independent in the vicinity of ∑ and are related to T for every ε. ε (∑) Tangent. This is achieved by pushing the parameterized local coordinate basis of the tangent bundle of ∑ forward to a sufficiently small tubular neighborhood of ∑. This means that if E°1, E°2 represent this basis, then E is defined as follows: i (y): =DT ε (E ° i(x)), for some x∈∑, any y that is sufficiently close to ∑ has y=T ε (x) form. Here, DT ε It is T ε The matrix of partial derivatives. Secondly, the vector field N is defined on this tubular neighborhood and has the following property: if y is sufficiently close to ∑, for some x∈∑, y=T ε If N(y) is of the form (x), then N(y) is T. ε (∑) is the unit normal vector field at point y. Notation: X = X ||+X ⊥ N is the orthogonal decomposition of the vector field X defined on this tubular neighborhood with respect to N. It should be noted that if for some x∈∑, y=T ε (x), then X || (y) and T ε (∑) are tangent.

[0702] Assume h ε and A ε They are T ε The induced surface metric tensor of (∑) and the second fundamental form are pulled back to ∑. Then its components are relative to the local basis. satisfy

[0703]

[0704] Assume [h] ε ] ij Represents matrix h ε The inverse component. Note that when referring to these quantities over ∑, the subscript ε is omitted. Next, the vector field N satisfies

[0705]

[0706] Finally, by differentiating the equation defining N, the following terms are derived:

[0707]

[0708] These two equations uniquely determine the vector field. In fact, it can be seen that these equations in T ε (∑) was established The relevant result is Lemma 1.

[0709] Lemma 1. The following formula holds:

[0710]

[0711] Proof. Applying equation Θ in the direction of Θ...<N,N> ≡1 Differentiate to derive It must be tangential. Then, the equations... Differentiating along the Θ direction for i=1, 2, to derive

[0712]

[0713] Where i = 1, 2. This implies the result.

[0714] Lemma 2. Imagine T ε It is a variant of ∑ generated by the vector field Θ and assumes [h ε ] ijAny local basis E1, E2 relative to the tangent bundle of ∑ will T ε The induced surface metric tensor of (∑) is pulled back to the components of ∑.

[0715]

[0716] Where A ij It is the second fundamental form of ∑, with respect to the components of the local basis, and This represents the surface covariant derivative operator.

[0717] Proof. The computation can be performed on any local basis because it can be proven that the result is independent of the basis (this is a fundamental principle of differential geometry). Therefore, the local basis described above can be used at any x∈∑. Thus, by means of the covariant derivative in… Based on the definitions and properties in [the text], the following results were obtained.

[0718]

[0719] Because the difference between the above covariant derivatives is a vanishing Lie bracket. This is because Θ and Both are forward projections of the coordinate basis vector fields. (In this calculation, the notation used is...) in It is any differentiable function; this notation reflects the standard differential geometry practice of combining the concept of a vector field X with the directional derivative operator in the X direction. Then, using the definition of the second fundamental form of ∑, we obtain the following terms:

[0720]

[0721] Lemma 3. Imagine T ε It is a variant of ∑ generated by the vector field Θ and assumes [A ε ] ij This is relative to the local basis T introduced above. ε The second fundamental form of (∑) is brought back to the components of ∑. Then...

[0722]

[0723] Proof. Again, computation can be performed on any local basis, since it can be shown that the result is independent of the basis. Therefore, the local basis described above can be used at any x∈∑. Through T ε The second fundamental form of (∑) is defined by the pullback, whose components relative to this local basis satisfy... Therefore, through the covariant derivative in Based on the definitions and properties in [the original text], the following results were obtained:

[0724]

[0725] Using Lemma 1 in the first term of equation 178, we obtain...

[0726]

[0727] The following derivation is made for the second term of equation 178:

[0728]

[0729] because<N,N> =1 and uniquely represent together The equation means

[0730]

[0731] Now, substituting the two derivatives into equation 178, we obtain...

[0732]

[0733] This uses the definition of the surface covariant derivative of a covariant bitensor. To proceed.

[0734] The surface curvature invariant of ∑ is: the total curvature, i.e., the square norm of the second fundamental form with respect to the induced surface metric tensor, denoted as ||A ∑ || 2 ; and the mean curvature, i.e., the trace of the second fundamental form relative to the induced surface metric tensor, denoted as H ∑ It should be noted that the Gaussian curvature or intrinsic curvature of ∑ is expressed by the formula for these invariants. Now, we give the derivatives of these quantities with respect to the variations introduced above.

[0735] Proposition 4. Imagine T ε It is a variant of ∑ generated by the vector field Θ and assumes H ε It is T ε The average curvature of (∑). Then

[0736]

[0737] Where Δ || It is the Laplace-Beltrami operator for curved surfaces.

[0738] Proof. Apply Lemmas 2 and 3 to the definition. In other words,

[0739]

[0740] As needed, Δ ||Simply identify it where it appears and use the Codazzi equation to... Item substitution

[0741] Proposition 5. Imagine T ε It is a variant of ∑ generated by the vector field Θ and assumes ||A ε || 2 It is T ε The total curvature of (∑). Then

[0742]

[0743] Proof. Apply Lemmas 2 and 3 to the definition. The required computations are similar to those used to prove the previous lemmas.

[0744] Variation of curvature of integral surface

[0745] In this section, the surface integral of the surface curvature invariant of ∑ is considered. First, the following preliminary results are described.

[0746] Lemma 6. Imagine T ε It is a variant of ∑ generated by the vector field Θ and assumes [h ε ] ij Any local basis E1, E2 relative to the tangent bundle of ∑ will T ε The induced surface metric tensor of (∑) is pulled back to the components of ∑. Assume φ ∑ : It is a surface-dependent scalar function whose shape is differentiable. Then...

[0747]

[0748] div || It is a surface divergence operator and H ∑ It is the average curvature of ∑.

[0749] Proof. Using the formula for the change of variables in a surface integral,

[0750]

[0751] Finally, using Lemma 2 and the following calculations:

[0752]

[0753] Proposition 7. Suppose that φ is a surface integral operator of the form stated in the introduction. Imagine T ε It is a variant of ∑ generated by the vector field Θ. Then...

[0754]

[0755] Substituting the results of propositions 4 and 5:

[0756]

[0757] Because it contains Θ || All terms can be proven by the generalized Stokes' Theorem to be equal to the vanishing ∫ ∑ div || (φ(H ∑ ,||A ∑ || 2 )Θ || )dσ.

[0758] Shape derivative of combined distance-normal vector curvature penalty

[0759] Assume Ω is The domain in Ω. The general form of a certain class of shape functions used to penalize the geometric features of Ω is based on the signed distance function and the boundary surface. The normal vector field and curvature invariants. Specifically, consider the following shape function:

[0760]

[0761] Where φ and ψ are smooth functions of its independent variables, and They are The signed distance function, normal field, and mean curvature. Most geometric performance functions (e.g., surface area, Willmore energy for shape smoothing, thickness penalty, overhang angle penalty) can be transformed in this form.

[0762] Thanks to the chain rule used for shape differentiation, the calculation of the shape derivative of equation 191 is decomposed into several parts. Assume Ω ε It is determined by the boundary velocity function Θ: The generated Ω variant, where the tangential and normal components are respectively derived from Θ || and Θ ⊥ This indicates. Then,

[0763]

[0764] Where · represents the so-called Lagrange derivative and This represents the partial derivative with respect to the k-th slot of the operand. The Lagrange derivative of a function f is defined as the derivative of f along the variation induced by Θ. The Lagrange derivative of a vector field is defined similarly.

[0765] Each of the above Lagrange differential expressions has the Θ derived above. || and Θ ⊥The formula is expressed as follows. Furthermore, it can be verified that for Θ... || The dependency has completely disappeared in the final expression (this result is known as the Hadamard-Zolésio structure theorem). This is confirmed by introducing a function for the volume integrand via... Define and for the integrand function of the surface through The Euler derivative is defined. Now, thanks to Stokes' theorem, for Θ... || All explicit dependencies have disappeared:

[0766]

[0767] The remaining Euler differential terms can be handled using the following equations derived above. These equations are:

[0768]

[0769]

[0770]

[0771] in Is to The projection of the nearest point on, and Δ || yes The Laplace-Beltrami operator, and yes The second fundamental form of the squared Frobenius norm (also known as the total curvature or the sum of the square principal curvatures). Substituting these into equation 193 yields the formula for the shape derivative.

[0772] Figure 7C An example is shown of a process for generating one or more parts of a 3D model of an object to be manufactured using stress constraints. The computer-aided design program obtains the design space of the modeled object upon which the physical structure corresponding to the 702 additive manufacturing will be based, the design criteria of the modeled object including at least one stress constraint, at least one in-use load condition of the physical structure, and the specifications of one or more materials used in the additive manufacturing physical structure.

[0773] The specification includes information, such as a dataset or function of thickness, a build angle, or both, to generate the corresponding strength of the material. The specification indicates multiple strength values ​​for one or more materials, each depending on the thickness of the physical structure to be built using one or more materials, and one or both of the build angle. The object will be manufactured using additive manufacturing with one or more materials.

[0774] Generates the 3D shape of the 704 modeled object through a generative design. This includes modifying both the geometry and topology of the 3D shape according to design criteria, at least one in-service load condition, and specifications for one or more materials. Generating a 704 model involves altering the evaluation of stress constraints at different locations on or within the modeled object based on corresponding values ​​from multiple strength values ​​during object shape and topology modifications. In some implementations, additional constraints, such as strain or displacement constraints, are used. Each strength value corresponds to one or both of the thickness and build angle at each location in the different positions. It should be noted that, when combined... Figure 7A The described design digital twin technique can be achieved through density-based topology optimization (e.g., using the SIMP method) rather than boundary-based topology optimization (e.g., level set method). In either case, the generative design of the 706 modeled object yields a 3D shape for use in additive manufacturing of the physical structure.

[0775] In some implementations, each strength value corresponds to the thickness at each location in different positions, or both the thickness and the construction angle. Furthermore, generating 704 includes: measuring the thickness at each location in different positions; and calculating the corresponding strength value based at least on the measured thickness at each location in different positions. (See above reference...) Figure 6B As described, measuring thickness may include using thickness measurements on the 3D shape of a generative design of a modeling object, said thickness measurement being a combination of at least two different thickness measurements. The two different measurements may include a combination of at least the following: (i) a first distance measurement is the length within the modeling object from a surface point of the modeling object projected in the negative normal direction, and (ii) a second distance measurement is the diameter of the largest sphere fitted within the modeling object and contacting a surface point of the modeling object, as determined by examining discrete sampling positions defined on the surface of the sphere.

[0776] In some implementations, modifying both the geometry and topology of the 3D shape involves enforcing design criteria that restrict the generative design of the 3D shape of the modeled object to a minimum thickness. For example, the minimum thickness might be based on the material's critical fatigue crack length, as referenced in [reference missing]. Figures 6A to 6E As described. The minimum thickness can be a ratio of the critical fatigue crack length, such as the minimum thickness being 10 times the critical fatigue crack length.

[0777] Return to Figure 7A Using maximized stress or strain elements and data relating fatigue strength to load cycles, the expected number of load cycles is determined for each of one or more in-service load conditions of the physical structure. The data relating fatigue strength to load cycles can be empirical, i.e., for example, such as... Figure 6DAs shown, a table of previously measured strengths from the corresponding material, or data that may be generated based on a function that outputs fatigue strength as a result of an input load cycle number for a given load condition of an object constructed using the corresponding material.

[0778] The fatigue safety factor inequality constraint is redefined for the modeled object based on the damage score calculated according to the required number of load cycles for the modeled object and the expected number of load cycles for each of one or more in-service load conditions of the physical structure. (See above reference.) Figures 4A to 4H The techniques described are applicable to handling arbitrary constraints, including inequality constraints.

[0779] The shape change rate of the implicit surface in the level set representation of the 716 three-dimensional shape is calculated at least according to the fatigue safety factor inequality constraint. The shape change rate is used to update the 718 level set representation to produce an updated form of the 3D shape of the modeled object. In some implementations, calculating the shape change rate includes calculating at least one shape change rate using a gradient determined based on one or both of the shape derivatives (e.g., the shape derivatives shown in Equation 193 and derived above) of the thickness and construction angle at each location in different positions.

[0780] Alternatively or in combination, i.e., in the presence of multiple constraints, at least one shape change rate can be used to calculate the shape change rate by a quantity determined according to a shape derivative formula, which approximates the shape derivative of one or both of the thickness and the construction angle at each of the different locations. An example shape derivative formula is the volume shape derivative formula described above with reference to Equation 105.

[0781] As described above with reference to arbitrary inequality constraints and equations 91-102, importance coefficients can be applied to inequality constraints to modify volume fraction-based inequality constraints using the shape derivative formula, which includes volume fraction-based inequality constraints. Such constraints can also be the minimum thickness inequality constraints described above with reference to equations 123-127, or as described above with reference to... Figures 5A to 6E The described constraint is a stress-based inequality constraint. The constraint can be modified using an importance coefficient, which adjusts the importance constraint based on whether one or more other constraints were violated in previous iterations of the iterative modification. An example of modifying the importance coefficient is multiplying it by a violation multiplier.

[0782] Repeat steps 708, 710, 712, 714, 716, and 718 until check 720 confirms that a predefined number of shape modification iterations have been performed, or the generated 3D shape of the modeling object in the design space has converged to a stable solution for one or more design criteria and one or more in-use load conditions, as referenced above. Figure 3AAs described above, any constraints can be run by designing a digital twin, provided that a numerical simulation is run during the process, for example, as referenced above. Figures 4A to 4H The description deals with any constraints. Therefore, digital twins are not exclusive to damage-tolerant fatigue. Similarly, any number of constraints can coexist with a design digital twin, since the design digital twin only affects the material model used in the numerical simulation.

[0783] Singularity and disconnection

[0784] Stress constraints, including von Miss stress constraints, can be essential for generative design, but they are difficult to implement because some elements in the finite element analysis model can exhibit very high stresses. Probabilistic methods for mitigating or eliminating such high stresses are provided below. Although stress constraints are mentioned in the examples below, these techniques can be applied to avoid singularities and prevent disconnections during generative design for any constraints. Singularities can occur, for example, due to sharp concave angles or poor meshing. The procedures described below can be implemented to be performed automatically as part of the generative design process, for example, including the references above. Figures 3A to 7C The technology described.

[0785] Figure 8A This illustrates an example of a process that iteratively modifies the 3D shape of a modeling object in a design space according to one or more design criteria, while avoiding excessive abrupt changes and minimizing the possibility of broken connections. The computer-aided design program obtains the physical structure corresponding to the 814 manufacturing process based on the design space of the modeling object, one or more design criteria of the modeling object, and one or more in-use load conditions of the physical structure. For example, as referenced above... Figure 3A , Figure 4A , Figure 5A , Figure 6A and Figure 7A The description describes the process of iteratively modifying the 3D shape of a generative design of a modeling object in the design space based on one or more design criteria and one or more in-use load conditions of the physical structure.

[0786] One solution to avoid singularities is the following modification using the simple percentile value σ from equation 134. x Replace the maximum stress.

[0787]

[0788] Where σ x This represents the stress of the element at the xth percentile when the elements are sorted in ascending order of σ. However, this leads to the domain max. ΩThe maximum stress in σ is overly sensitive to the value of x. Instead, the percentile x is converted to the standard normal deviation z using the inverse of the error function.

[0789]

[0790] Next, the maximum stress is calculated based on the following:

[0791] Where μ(σ) and χ(σ) represent the mean deviation and standard deviation of the stress distribution, both of which can be calculated using any conventional technique.

[0792] Table 1 below shows the stress singularities avoided by the object based on the techniques described in this section:

[0793]

[0794] Table 1

[0795] Percentile-based methods bring stability to the oscillating nature of global stress by smoothing the maximum stress at a point based on the standard deviation of stress values ​​at points on the same percentile of the object.

[0796] High-speed smoothing

[0797] The previous section described the maximum value of the constraint value in the adjustment domain (i.e., max). Ω g i This also applies to the constrained derivative dg. t / dΩ, for clarity by dg t express.

[0798] Figure 8B This is a graphical representation of an example of a geometric connection breaking during optimization. Object 804 is shown as being at a different stage of optimization. Typically, geometric ports can break off from the main body of the design (shown at 806), in which load paths are occupied by a substitute material (Ersatzmaterial). Small areas with very high advection velocities can cause connection breaks and localized shape changes. A combination of velocity clamping and high-velocity smoothing can be used to prevent this phenomenon.

[0799] Normalization: Normalizes all shape derivatives so that the magnitude of the maximum value on the surface of the current design Γ approximates the voxel size Δs.

[0800]

[0801] This means that the advection time T = 1 / Δs is sufficient to prevent the geometry from advection beyond one voxel. A simple way to achieve this is by using velocity clamping.

[0802] Velocity clamping: The first step in velocity clamping is to calculate a reference value for the shape derivative based on the percentile method described in reference equations 197-199 above.

[0803] (dg t ) 参考 =μ(dg) t )+erf -1 (x)χ(dg t (201)

[0804] Where μ and χ represent the mean deviation and standard deviation of the shape derivative on the surface, and x∈{0.9, 0.95, 0.99, 0.999, etc.} are the percentiles given by the user. It should be noted that the percentile value x may be much higher (but not excessively high) for benign shape derivatives such as strain energy, and may be lower for shape derivatives such as those with large fluctuations in stress.

[0805] Next, for the inner narrow band (i.e., the narrow band width is equal to w) nb All grid points x in Δs) k For values ​​higher than the reference value (dg) t ) 参考 The velocity value is clamped, while retaining the sign β as follows:

[0806]

[0807] Where ψ(x) l () represents the level set of grid points. However, this leads to abrupt changes in the velocity profile. Velocity clamping should be done in a way that allows for a smooth transition in the velocity profile.

[0808] High-speed smoothing: The goal of high-speed smoothing is to smooth all high speeds (|dg) t |>|(dg t ) 参考 This allows for progressive smoothing of all values ​​higher than the reference velocity calculated according to Equation 201. The user-provided parameter ρ (e.g., ρ = 0.85) is used to scale all velocity values ​​lower than the reference velocity.

[0809]

[0810] Where k 参考 This indicates the index of the grid point after it has been sorted by its shape derivative.

[0811] The smoothed velocities of all grid points with velocities greater than the reference velocity are found by fitting a cubic polynomial:

[0812] (dg t (x k ))平滑 =αξ 3 +bξ 2 +cξ+dk 参考 ≤k≤n (204)

[0813] Where ξ represents the offset index ξ = kk 参考 The unknown coefficients {a, b, c, d} can be obtained by: using ξ = 0 and ξ = ξ 最大 =nk 参考 By setting appropriate boundary conditions, the obtained system of linear equations can be solved.

[0814]

[0815] The gradient of the smooth velocity at ξ = 0 can be derived as follows:

[0816]

[0817] Finally, the speed of the smoothed normalization is given by the following terms.

[0818]

[0819] Return to Figure 8A Within iterative modifications, numerical simulations of the 816 modeling object are performed based on the current form of the 3D shape and one or more in-use load conditions to produce a current numerical evaluation of the physical response (e.g., structural response) of the modeling object. For example, as referenced above... Figure 3A The description focuses on calculating the shape change rate of an implicit surface in a level set representation of a three-dimensional shape.

[0820] The shape change rate of 820 is then varied according to a polynomial function that has been fitted to at least a portion of the shape change rate that is higher than a reference rate. For example, as described in reference equation 204 above, the polynomial function can be a cubic polynomial. Although the given example is a cubic polynomial, polynomials of other orders can be used.

[0821] For example, as referenced above Figure 3A The described approach uses the rate of shape change to update the 822-level set representation to produce an update pattern for the 3D shape of the modeled object. In some implementations, such as those described above, the reference rate is set based on the mean and standard deviation of the shape derivative on the implicit surface.

[0822] Repeat steps 816, 818, 820, and 822 in iterative modifications until check 830 determines that a predefined number of shape modification iterations have been performed, or that the generated 3D shape of the modeling object in the design space has converged to a stable solution for one or more design criteria and one or more in-use load conditions, for example, as referenced above. Figure 3A As described. Ultimately, a 3D shape of the modeled object can be provided for generating a design, which can then be used to manufacture the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

[0823] Backup and restore

[0824] Another technique to prevent geometric connections from breaking is to back up all critical data, including the geometric level set, at each iteration and restore the data if the advection results are unsatisfactory. Therefore, in some implementations, excessive changes (824) can be checked in the changes made during the update, after update 822 but before repetition. In these cases, the current form of the 3D shape (826) can be set as the updated form of the 3D shape to undue excessive changes in the next iteration, i.e., to undo the excessive changes. Then, as described below, shape changes in the next iteration (828) can be mitigated. Mitigating changes in the next iteration may include reducing the target volume change of the 3D shape of the generative design of the modeled object in the iteratively modified next iteration.

[0825] To prevent the recurrence of the same undesirable results, the following modification using the multiplier β is referenced above. Figures 4A to 4H The volume change Δv described is calculated for each iteration t. t This slows down convergence by applying a small volume change:

[0826] Δv t ←βΔv t (208)

[0827] This multiplier can be started with β = 1 and updated in each iteration using a fixed increment Δβ (> 0) based on the result of each advection:

[0828] β←max(0, β-Δβ) if advection is not expected (209)

[0829] β←min(1,β+Δβ)Others(210)

[0830] The next task is to classify the outcome of each advection as desired / undesirable. One way to achieve this is to monitor changes in the Lagrangian, i.e., as described above with reference to equations 15-25. This can be achieved in a more refined way by imposing constraints on the permissible relative changes in the objective. The relative change in the objective / constraint given iteration t is defined as...

[0831]

[0832]

[0833] Among them, constraint error Using the PID stability form of Equation 95, and by employing the arbitrary inequality constraints referenced above, and Figures 4A to 4H The importance coefficient term of the described technology calculation To calculate. By The reference value for the target is given by the following items.

[0834]

[0835] The maximum permissible change of the target is limited by the amount of target reduction. and increase The maximum permissible change of a constraint is limited by the positive (benign) change of the constraint. Negative change (malignant) The sign of the constraint inequality needs to be checked to determine the appropriate limit.

[0836] Therefore, in some implementations, before performing 816, elements generated in the numerical simulation based on the current form of the 3D shape but partially, and not entirely, within the implicit surface under Dirichlet boundary conditions are identified and removed, after which the numerical simulation is performed. Elements and nodes refer to how the design geometry is represented during 816. In finite element simulations, the domain is replaced by a set of elements. Each element is defined by a set of nodes along the element's boundary 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 used to create these nodes and elements is called mesh generation.

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

[0838] Furthermore, in some implementations, checking for excessive variation in 824 may include comparing the changes caused by updates under one or more design criteria with predefined limits on the amount of variation allowed for one or more design criteria in a single iteration of iterative modification.

[0839] substitute materials

[0840] A common practice in SIMP methods for topology optimization is to use variable-density materials or substitute materials. Elements with a density ρ = 1 are considered to be inside the domain, while those with ∈ < ρ < 1 are on the boundary, where ∈ (usually set to a small value such as 0:001) represents the density of the material outside the design. During the finite element analysis, the stiffness of each element is multiplied by its density and a penalty factor p is applied to penalize intermediate densities.

[0841] K←ρ p K. (214)

[0842] The level set method often uses only two states (ρ∈{∈, 1}) as boundaries and is therefore more clearly defined. Since there is no intermediate density, i.e., K←ρ p K, therefore no penalty factor is needed. However, this can sometimes lead to disconnections because elements outside the domain with ρ = ∈ can support load paths.

[0843] Furthermore, it has been observed that materials with ρ=∈ significantly affect the stability safety factor predicted by the finite element model. Substitute materials that have a significant impact on the predicted stability factor can be completely removed from the finite element model. This is achieved by grouping all elements with ρ=∈ at each iteration and removing all such groups that are not connected to any nodes under the Dirichlet boundary conditions.

[0844] As mentioned above, in Figure 8A In some implementations of the process shown, prior to execution 816, elements generated in 815 during the numerical simulation based on the current form of the three-dimensional shape that are partially, but not entirely, within the implicit surface are identified. The density of the identified elements can then be set equal to the corresponding volume fraction of the identified elements, where the volume fraction is the number of identified elements falling within the implicit surface. In some implementations, execution 816 may include a density penalty on the stiffness of the identified elements, as described below.

[0845] As another solution to the severe disconnection problem, a tetrahedral cutting algorithm can be applied, which recovers the use of substitute material in elements at the boundaries of the current domain. Essentially, the stiffness of such elements is better approximated by setting the density to be equal to the volume fraction of elements within the domain. As in the SIMP method, the stiffness is then penalized based on the density, where the penalty factor is p = 1.

[0846]

[0847] Calculating the volume of each element within a domain is not straightforward, because all nodes of some elements that effectively overlap with the domain may be outside the elements.

[0848] Figure 8C This is a graphical representation of an instance of a geometric figure 808 having analogous elements 810, which are categorized based on the intersection of element 810 and geometric figure 808. The following discussion is relative to the element... The problem is illustrated. Such elements can be subdivided 812 times until the node is on either side of the 0th iso-contour line, i.e.

[0849]

[0850] Where x i This represents the node coordinates of element e. Note that it may be necessary to specify the coordinates for some elements (e.g., element e). The subdivision is repeated multiple times until the conditions of equation 216 are satisfied. The required depth of subdivision l can be determined by recursively subdividing the elements and comparing the total volume within the domain.

[0851]

[0852] in Represents element The index of the child element. Recursive subdivision can stop when the above condition is met. In some implementations, the voxel size Δs is set to half the average edge length of the solid element. Considering that the minimum feature size ≈ Δs, usually only one level of subdivision is sufficient.

[0853] The volume fraction of a tetrahedral element can be reduced based on the domain. Linear interpolation can be used when calculating the intersection of the element's edge with the 0th isomorphic contour.

[0854]

[0855] Figure 9 This is a schematic diagram of a data processing system including a data processing device 900, which can be programmed as a client or server. The data processing device 900 is connected to one or more computers 990 via a network 980. Although in Figure 9Only one computer is shown as data processing device 900, but multiple computers can be used. Data processing device 900 includes various software modules that can be distributed between the application layer and the operating system. These software modules may include executable and / or interpretable software programs or libraries, including tools and services for implementing one or more 3D modeling programs 904 that implement the systems and techniques described above. Thus, one or more 3D modeling programs 904 may be one or more CAD programs 904 and may implement one or more generative design processes (e.g., generative design using one or more level set-based methods) for topology optimization and numerical simulation operations (finite element analysis (FEA) or other operations). Furthermore, one or more programs 904 may potentially implement manufacturing control operations (e.g., generating and / or applying toolpath specifications to implement the manufacture of the design object). The number of software modules used may vary depending on the implementation. Furthermore, the software modules may be distributed across one or more data processing devices connected via one or more computer networks or other suitable communication networks.

[0856] The data processing apparatus 900 also includes hardware or firmware means, including one or more processors 912, one or more auxiliary 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-threaded 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 auxiliary devices 914. The data processing apparatus 900 uses the communication interface 919, for example, via a network 980, to communicate with one or more computers 990. Examples of user interface devices 920 include displays, cameras, speakers, microphones, haptic feedback devices, keyboards, mice, and VR and / or AR devices. The data processing apparatus 900 may store instructions for implementing operations associated with one or more programs described above on, for example, the computer-readable medium 916 or one or more auxiliary devices 914, such as one or more of the following: hard disk devices, optical disk devices, magnetic tape devices, and solid-state storage devices.

[0857] The embodiments and functional operations of the subject matter described in this specification can be implemented using: digital electronic circuits; or computer software, firmware, or hardware, including the structures disclosed in this specification and their structural equivalents; or combinations thereof. 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 a data processing device or for controlling the operation of the data processing device. The computer-readable medium can be an manufactured product, such as a hard disk drive in a computer system or an optical disc sold through retail channels, or an embedded system. The computer-readable medium may be individually required and later encoded using one or more modules of the computer program instructions, for example, after transmission of one or more modules of the computer program instructions via 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 combinations thereof.

[0858] The term "data processing device" encompasses all devices, apparatuses, and machines used for processing data, including, for example, programmable processors, computers, or multiple processors or computers. In addition to hardware, the device may also include code that creates an execution environment for the computer program under discussion, such as code constituting processor firmware, protocol stacks, database management systems, operating systems, runtime environments, or combinations thereof. Furthermore, the device may employ a variety of different computing model infrastructures, such as web services, distributed computing, and grid computing infrastructures.

[0859] Computer programs (also referred to as programs, software, software applications, scripts, or code) can be written in any suitable programming language, including compiled or interpreted languages, declarative or procedural languages, and can be deployed in any suitable form, including as standalone programs or as modules, components, subroutines, or other units suitable for use in a computing environment. A computer program does not necessarily correspond to a file in a file system. A program may be stored as a part of a file containing other programs or data (e.g., one or more scripts stored in a markup language document), a single file dedicated to the program under discussion, or multiple coordinating files (e.g., files storing one or more modules, subroutines, or code sections). A computer program can be deployed to execute on one or more computers located in one site or distributed across multiple sites and interconnected via a communication network.

[0860] The processes and logic flows described in this specification can be executed by one or more programmable processors, which execute one or more computer programs to perform functions by manipulating input data and generating outputs. The processes and logic flows can also be executed by special-purpose logic circuitry such as FPGAs (Field-Programmable Gate Arrays) or ASICs (Application-Specific Integrated Circuits), and the device can be implemented as any of the above.

[0861] Processors suitable for executing computer programs include, for example, both general-purpose and special-purpose microprocessors, as well as any one or more processors in any kind of digital computer. Generally, a processor receives instructions and data from read-only memory or random access memory, or both. The basic components of a computer are a processor for executing instructions and one or more memory devices for storing instructions and data. Generally, a computer also includes, or is operatively coupled to, one or more mass storage devices for storing data, such as magnetic disks, magneto-optical disks, or optical disks, or both. However, a computer does not need to have such devices. Furthermore, a computer may be embedded in another device, such as a mobile phone, 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. Suitable devices for storing computer program instructions and data include all forms of non-volatile memory, media, and memory devices, including, for example, semiconductor memory devices such as EPROM (Erasable Programmable Read-Only Memory), EEPROM (Electrically Erasable Programmable Read-Only Memory), and flash memory devices; magnetic disks, such as internal hard disks or removable disks; magneto-optical disks; and CD-ROMs and DVD-ROMs. Processors and memory may be supplemented by dedicated logic circuitry or integrated therein.

[0862] To provide interaction with the user, embodiments of the subject matter described in this specification can be implemented on a computer having: a display device, such as an LCD (liquid crystal display), an OLED (organic light-emitting diode) display device, or another display for displaying information to the user; and a keyboard and pointing device, such as a mouse or trackball, that the user can use to provide input to the computer. Other types of devices may also be used to provide interaction with the user; for example, feedback provided to the user may be any suitable form of sensory feedback, such as visual feedback, auditory feedback, or tactile feedback; and input from the user may be received in any suitable form, including sound, speech, or tactile input.

[0863] A computing system may include clients and servers. Clients and servers are typically geographically separated and interact via a communication network. The client-server relationship is established by means of computer programs running on respective computers that have a client-server relationship with each other. Embodiments of the subject matter described in this specification can be implemented in computing systems that include back-end components, such as data servers; or middleware components, such as application servers; or front-end components, such as client computers having a graphical user interface or browser user interface through which a user can interact with an implementation of the subject matter described in this specification; or any combination of one or more such back-end, middleware, or front-end components. Components of the system may be interconnected via any suitable form or medium of digital data communication, such as a communication network. Examples of communication networks include local area networks (“LANs”) and wide area networks (“WANs”), interconnected networks (e.g., the Internet), and peer-to-peer networks (e.g., self-organizing peer-to-peer networks).

[0864] While this specification contains numerous details of implementation methods, these details should not be construed as limiting the scope of the claims or potential claims, but rather as descriptions of features specific to particular embodiments of the disclosed subject matter. Certain features described in this specification within the context of individual embodiments may also be implemented in combination in a single embodiment. Conversely, various features described within the context of a single embodiment may also be implemented individually or in any suitable sub-combination in multiple embodiments. Furthermore, although features may be described above as operating in a particular combination, and even initially claimed in this way, one or more features from a claimed combination may be removed from that combination in some cases, and a claimed combination may refer to a sub-combination or a variation thereof.

[0865] Similarly, although operations are shown in a specific order in the accompanying drawings, this should not be construed as requiring such operations to be performed in the specific order shown or in an ordered order, or to perform all the operations shown, in order to achieve the desired result. In some cases, multitasking and parallel processing may be advantageous. Furthermore, the separation of the various system components in the embodiments described above should not be construed as requiring such separation in all embodiments, and it should be understood that the described program components and systems can generally be integrated together into a single software product or packaged into multiple software products.

[0866] Therefore, specific embodiments of the invention have been described. Other embodiments are also within the scope of the appended claims. Furthermore, the actions described in the claims can be performed in different orders and still achieve the desired result.

Claims

1. A method, the method comprising: The computer-aided design program obtains the design space of the modeling object on which the manufacturing of the corresponding physical structure will be based, one or more design standards of the modeling object, one or more in-use load conditions of the physical structure, and the critical fatigue crack length of the material used to manufacture the physical structure. The computer-aided design program iteratively modifies the three-dimensional shape of the generative design of the modeling object in the design space according to one or more design criteria, one or more in-service load conditions of the physical structure, and the critical fatigue crack length of the material, wherein the iterative modification includes enforcing a design criterion that restricts the minimum thickness of the three-dimensional shape of the generative design of the modeling object, the minimum thickness being based on the critical fatigue crack length of the material; as well as The computer-aided design program provides the three-dimensional shape of the generative design of the modeled object for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

2. The method as described in claim 1, The design criteria mentioned above include the required number of load cycles for the modeling object under each of the one or more load conditions in use of the physical structure; The acquisition further includes obtaining one or more specifications of the material used to manufacture the physical structure, the one or more specifications including data relating fatigue strength to load cycles; and The iterative modification mentioned above also includes: Numerical simulations of the modeled object are performed based on the current form of the three-dimensional shape and one or more in-use load conditions to produce a current numerical assessment of the physical response of the modeled object. For at least one of the one or more in-use load conditions of the physical structure, find the stress or strain element that maximizes the stress or strain from the current numerical evaluation of the physical response of the modeled object. The expected number of load cycles for each of at least one of the one or more in-service load conditions of the physical structure is determined using the maximized stress or strain element and the data relating fatigue strength to load cycles. The fatigue safety factor inequality constraint of the modeling object is redefined based on the damage score calculated according to the required number of load cycles of the modeling object and the expected number of load cycles for each of at least one of the one or more in-use load conditions of the physical structure, wherein the fatigue safety factor inequality constraint is defined as the difference between the fatigue safety factor of the material and the target fatigue safety factor being greater than or equal to zero. The rate of shape change of the implicit surface in the level set representation of the three-dimensional shape is calculated at least according to the fatigue safety factor inequality constraint.

3. The method as described in claim 2, The one or more specifications include two or more specifications for corresponding different materials used to manufacture the physical structure, the data includes data relating fatigue strength to load cycles for each of the different materials, determining the expected number of load cycles includes determining an individual expected number of load cycles for each of the different materials, and redefining the fatigue safety factor inequality constraint includes: A separate fatigue safety factor for each of the different materials is calculated based on the corresponding damage fraction calculated according to the corresponding expected load cycle number in the expected load cycle number of the different materials; as well as The fatigue safety factor inequality constraint of the modeling object is redefined using the minimum value of the fatigue safety factor of the different materials. and The calculation of the shape change rate includes using a gradient determined based on the shape derivative of the fatigue safety factor to calculate at least one shape change rate.

4. The method of claim 2, wherein the lookup includes calculating the maximum stress value under load conditions in use based at least on the standard deviation of the stress distribution in the current numerical assessment of the physical response of the modeled object.

5. The method of claim 1, wherein the enforcement of the three-dimensional shape of the generated design of the modeling object uses a thickness measurement, the thickness measurement being a combination of at least two different thickness measurements.

6. The method of claim 5, wherein the at least two different thickness measurements comprise (i) a first distance measurement, the first distance measurement being the length within the modeling object from a surface point of the modeling object projected in the negative normal direction; and (ii) a second distance measurement, the second distance measurement being the diameter of the largest sphere fitted within the modeling object and contacting the surface point of the modeling object, as determined by examining discrete sampling positions defined on the surface of the sphere.

7. The method of claim 6, wherein the enforcement includes using an inequality constraint based on volume fraction or minimum thickness as a proxy for the design criterion limiting the minimum thickness, wherein, The inequality constraint based on the minimum thickness is defined as the difference between the minimum thickness and the target minimum thickness being greater than or equal to zero, wherein the inequality constraint based on volume fraction or minimum thickness is modified using an importance coefficient, which is set to zero during the initial phase of the iterative modification and adjusted during subsequent phases of the iterative modification based on whether one or more other constraints were violated in previous iterations of the iterative modification.

8. The method of claim 7, wherein the method comprises: In the iteratively modified multiple iterations, the target value based on the inequality constraint of volume fraction or minimum thickness is adjusted between the initial target value and the final target value; as well as When adjusting the target value during the multiple iterations, a proportional-integral-derivative controller is used to adjust and stabilize the changes made to the amount of modification of the modeling object determined based on the evaluation of the inequality constraints based on volume fraction or minimum thickness.

9. The method of claim 1, wherein obtaining the critical fatigue crack length of the material comprises: Obtain one or more specifications for the material used to manufacture the physical structure; as well as The critical fatigue crack length of the material is calculated based on information from one or more specifications, including the modulus of the fatigue crack growth curve of the material.

10. The method of claim 1, wherein the three-dimensional shape of the generative design of the modeling object comprises a level set representation of an implicit surface, the one or more design criteria comprise the required number of load cycles of the modeling object under each of the one or more load conditions in use of the physical structure, and the iterative modification comprises: Numerical simulations of the modeled object are performed based on the current form of the three-dimensional shape and one or more in-use load conditions to produce a current numerical assessment of the physical response of the modeled object; The current numerical evaluation and the thickness measurement are used to determine the expected number of load cycles for each of the one or more in-use load conditions of the physical structure in order to enforce the design criteria that limit the minimum thickness; The fatigue safety factor inequality constraint of the modeling object is redefined based on the damage score calculated according to the required number of load cycles of the modeling object and the expected number of load cycles of each of the one or more in-use load conditions of the physical structure; The rate of shape change of the implicit surface shall be calculated at least according to the fatigue safety factor inequality constraint; The shape change rate is used to update the level set representation to produce an updated pattern of the three-dimensional shape of the modeled object; as well as The execution, determination, redefinition, calculation, and update are repeated at least once until a predefined number of shape modification iterations have been performed, or until the 3D shape of the generated design of the modeling 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 conditions.

11. The method of claim 10, wherein the one or more in-service load conditions of the physical structure include two or more in-service load conditions of the physical structure, the one or more design criteria include the required number of load cycles for the modeling object under each of the two or more in-service load conditions of the physical structure, determining the expected number of load cycles includes determining the individual expected number of load cycles for each of a plurality of points on the implicit surface under each of the two or more in-service load conditions, and redefining the fatigue safety factor inequality constraint includes: For each of the plurality of points, load-specific damage scores corresponding to the two or more in-use load conditions are summed, wherein each load-specific damage score includes dividing the expected number of load cycles for one of the plurality of points and one of the in-use load conditions by the required number of load cycles for the one of the in-use load conditions to produce the sum of the load-specific damage scores for each of the plurality of points; Take the reciprocal of each of the sums of the load-specific damage fractions; and The fatigue safety factor inequality constraint of the modeling object is redefined using the minimum of the sum of the reciprocals of the load-specific damage fractions.

12. The method of claim 11, wherein calculating the shape change rate comprises calculating at least one shape change rate using a quantity determined according to a shape derivative formula, the shape derivative formula being approximately the shape derivative of the fatigue safety factor.

13. The method of claim 12, wherein the shape derivative formula for the shape derivative approximating the fatigue safety factor includes an inequality constraint based on a volume fraction or minimum thickness, the inequality constraint based on a volume fraction or minimum thickness being modified using an importance coefficient, the importance coefficient being adjusted based on whether one or more other constraints were violated in the iteratively modified previous iteration, wherein, The inequality constraint based on minimum thickness is defined as the difference between the minimum thickness and the target minimum thickness being greater than or equal to zero.

14. The method of claim 13, wherein the method comprises: In the iteratively modified multiple iterations, the target value based on the inequality constraint of volume fraction or minimum thickness is adjusted between the initial target value and the final target value; as well as When adjusting the target value during the multiple iterations, a proportional-integral-derivative controller is used to stabilize the changes made to the quantity determined according to the shape derivative formula, and to adjust the overall contribution of the quantity determined according to the shape derivative formula to the rate of shape change used in the update.

15. A system comprising: A non-transitory storage medium, wherein instructions for a computer-aided design program are stored on the non-transitory storage medium; as well as One or more data processing devices, the one or more data processing devices being configured to run the instructions of the computer-aided design program to perform the following operations: The design space of the modeling object on which the manufacturing of the corresponding physical structure will be based, one or more design standards of the modeling object, one or more in-use load conditions of the physical structure, and the critical fatigue crack length of the material used to manufacture the physical structure are obtained. The three-dimensional shape of the generative design of the modeling object in the design space is iteratively modified according to one or more design criteria, one or more in-service load conditions of the physical structure, and the critical fatigue crack length of the material, wherein the one or more data processing devices are configured to run the instructions of the computer-aided design program to iteratively modify the three-dimensional shape of the generative design of the modeling object by: being configured to run the instructions of the computer-aided design program to enforce a design criterion that limits the minimum thickness of the three-dimensional shape of the generative design of the modeling object, the minimum thickness being based on the critical fatigue crack length of the material; as well as Provide the three-dimensional shape of the generative design of the modeled object for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

16. The system of claim 15, the system comprising an additive manufacturing machine, wherein the one or more data processing devices are configured to run the instructions of the computer-aided design program to generate a toolpath specification for the additive manufacturing machine based on a three-dimensional model, and to manufacture the physical structure corresponding to the object by the additive manufacturing machine using the toolpath specification.

17. A non-transitory computer-readable medium encoding a computer-aided design program, the computer-aided design program being operable to cause one or more data processing devices to perform operations including: The computer-aided design program obtains the design space of the modeling object on which the manufacturing of the corresponding physical structure will be based, one or more design standards of the modeling object, one or more in-use load conditions of the physical structure, and the critical fatigue crack length of the material used to manufacture the physical structure. The computer-aided design program iteratively modifies the three-dimensional shape of the generative design of the modeling object in the design space according to one or more design criteria, one or more in-service load conditions of the physical structure, and the critical fatigue crack length of the material, wherein the iterative modification includes enforcing a design criterion that restricts the minimum thickness of the three-dimensional shape of the generative design of the modeling object, the minimum thickness being based on the critical fatigue crack length of the material; as well as The computer-aided design program provides the three-dimensional shape of the generative design of the modeled object for use in manufacturing the physical structure corresponding to the modeled object using one or more computer-controlled manufacturing systems.

18. The non-transitory computer-readable medium of claim 17, wherein the one or more computer-controlled manufacturing systems include an additive manufacturing machine, and the operation comprises: Based on the 3D model, toolpath specifications are generated for the additive manufacturing machine. as well as The physical structure corresponding to the object is manufactured by the additive manufacturing machine using the toolpath specification.

19. The non-transitory computer-readable medium of claim 17, wherein the operation comprises the method of any one of claims 2 to 14.

Images