Systems and methods for resolving numerical instabilities

Abstract

To determine physical behavior of real-world objects.SOLUTION: Using a computer-based model representing a real-world object, a system adds pseudo-constitutive modeling to the computer-based model to alleviate local numerical instabilities caused by ill-conditioned elemental stiffness operators within elements of the model. A computer-based model, representing a real-world object using a plurality of elements, is defined which indicates one or more materials represented by the elements. Equations describing physics-based behaviors of the one or more materials are defined. A stabilization equation that is a function of a non-linear deformation gradient matrix is defined. A simulation is performed of the real-world object, subject to a load, using the defined computer-based model, the defined equations describing physics-based behaviors, and the defined stabilization equation. Performing the simulation includes applying the stabilization equation to each of the plurality of elements.SELECTED DRAWING: Figure 5

Description

[Background technology]

[0001] Many systems and programs are available on the market for designing parts using computer-aided design (CAD) or computer-aided engineering (CAE). These so-called CAD systems allow users to construct and manipulate complex three-dimensional models of objects or assemblies of objects. CAD systems therefore provide a representation of the modeled object using edges or lines, or in certain cases, faces or polygons. The lines, edges, faces, or polygons may be represented in various ways, such as, for example, non-uniform rational B-splines (NURBS).

[0002] These CAD systems manage parts or assemblies of parts modeled by objects, primarily specifications of geometric shapes. In particular, a CAD file contains specifications from which geometric shapes are generated. From the geometric shapes, a representation is generated. The specifications, geometric shapes, and representations may be stored in a single CAD file or multiple CAD files. CAD systems include graphical tools to represent the modeled object to the designer; these tools are dedicated to displaying complex objects. For example, an assembly may contain thousands of parts. CAD systems can be used to manage models of objects stored in electronic files.

[0003] The advent of CAD and CAE systems allows for a wide range of representation possibilities for objects. One such representation is the finite element model (FEM). An FEM, or other CAD, CAE, or computer-based model, may be programmed so that the model has the properties of the underlying object it represents. When an FEM model or other such computer-based model is programmed in such a manner, it can be used to perform simulations of the object it represents. For example, an FEM may be used to represent the interior cavity of a vehicle, an acoustic fluid surrounding a structure, and any number of real-world objects and systems. When a given model represents an object and is programmed accordingly, it can be used to simulate the real-world object itself. For example, an FEM representing a stent may be used to simulate the use of the stent in an actual medical setting.

[0004] Computer-based models can be used to improve the design of the object that the model represents. Design improvements can be identified through the use of computer-based optimization techniques that run a series of simulations to identify modifications to the design of the model and, in turn, the underlying real-world object that the model represents. Summary of the Invention

[0005] Although computer-based methods exist for determining the physics-based, i.e., physical behavior, of real-world objects, these existing methods may provide inaccurate results. Specifically, existing methods may experience divergence and local numerical instability when determining a principal solution for models that undergo large displacements and deformations. The divergence / instability in obtaining a principal solution can result in a failure to complete the analysis, i.e., not reaching a converged solution. Therefore, improvements to existing simulation and modeling methods are needed. Embodiments provide such improvements. Embodiments are directed to designing real-world objects, determining the behavior of real-world objects, and simulations that provide pseudo-constitutive modeling to resolve numerical instabilities in determining a principal solution for a simulation.

[0006] One such exemplary embodiment is a computer-implemented method for determining the physical behavior, i.e., physics-based behavior, of a real-world object. Such an embodiment begins by defining, in a processor's memory, a computer-based model representing the real-world object using a plurality of elements. The defined model describes one or more materials represented by each element of the plurality of elements. The embodiment continues by defining equations describing the physical behavior of the one or more materials and defining a stabilization equation that is a function of a nonlinear deformation gradient matrix. A simulation of the real-world object subjected to a load is then performed using (i) the defined computer-based model, (ii) the defined equations describing the physical behavior, and (iii) the defined stabilization equation. In such an embodiment, performing the simulation includes applying the stabilization equation to each of the plurality of elements. The results of performing the simulation are indicative of the physical behavior, i.e., physics-based behavior, of the real-world object.

[0007] According to one embodiment, the "elements" (i.e., elements of a computer-based model used to represent real-world objects) are tessellated elements of a finite element model (i.e., a computer-based model). Note, however, that embodiments are not limited to utilizing tessellated elements; instead, embodiments may employ any elements, such as triangular or quadrilateral membranes or shells.

[0008] In an embodiment, the computer-implemented method continues further by, for each element of the plurality of elements in the model, associating an artificial internal force with the element based on the defined stabilization equation. In another embodiment, the computer-implemented method continues further by, for each element of the plurality of elements in the model, associating an artificial-based and physics-based behavior (i.e., a combination of artificial-based and physics-based behavior) with the element based on the associated artificial internal force. According to one embodiment, the computer-implemented method performs a simulation of a real-world object subjected to loads using the artificial internal forces associated with each element and the artificial-based and physics-based behavior associated with each element.

[0009] In yet another embodiment, a computer-implemented method performs a simulation of a real-world object using the defined stabilization equation to reduce a condition value associated with a given element of the plurality of elements, where the condition value is a function of a nonlinear deformation gradient matrix.

[0010] In another embodiment, a computer-implemented method updates a computer-based model based on results of performing a simulation, and determines updated physical behavior of the real-world object by performing a simulation of the real-world object using (i) the updated computer-based model, (ii) the defined equations describing the physics-based behavior, and (iii) the defined stabilization equations.

[0011] In another embodiment, the computer-implemented method further includes iterating (i) updating the computer-based model based on the determined updated physical behavior and (ii) determining the updated physical behavior until the determined updated physical behavior satisfies a criterion.

[0012] In some embodiments, the computer-based model is any one or combination of a finite element model, a boundary element method, a finite difference method, a finite volume method, or a discrete element method.

[0013] In embodiments, the real-world object represented by the computer-based model may be any real-world object. For example, in one embodiment, the real-world object is one of an automobile, an industrial facility, an aircraft, a civil structure, a marine device, a medical device, a consumer product, an electronic device, an armored vehicle, or a manufacturing facility. Furthermore, the real-world object represented by the computer-based model may be used in one of an automobile application, an industrial facility application, an aircraft application, a civil structure application, a marine device application, a medical device application, a consumer product application, an electronic device application, an armored vehicle application, or a manufacturing facility application, among other examples. In this manner, embodiments may optimize the design of real-world objects and optimize applications that utilize real-world applications.

[0014] In an embodiment, the corresponding physics-based behavior of the car may be, but is not limited to, deformation of the car's body panels during a collision.

[0015] In an embodiment, the corresponding physics-based behavior of the industrial equipment may be, but is not limited to, deformation of one or more pieces of industrial equipment during operation.

[0016] In an embodiment, the corresponding physics-based behavior of the airplane may be, but is not limited to, deformation of the airplane's wings subjected to aerodynamic forces.

[0017] In an embodiment, the corresponding physics-based behavior of the civil structure may be, but is not limited to, deformation of one or more structural components subjected to a load.

[0018] In an embodiment, the corresponding physics-based behavior of the marine device may be, but is not limited to, deformation of the marine device subjected to varying degrees of hydrostatic pressure.

[0019] In embodiments, the corresponding physics-based behavior of the medical device can be, but is not limited to, the deformation of a syringe needle tip subjected to torque.

[0020] In an embodiment, the corresponding physics-based behavior of the consumer good may be, for example, but not limited to, the deformation of the exemplary consumer good when subjected to deformations associated with normal use of the consumer good.

[0021] In an embodiment, the corresponding physics-based behavior of the electronic device may be, but is not limited to, deformation of the electronic device when subjected to deformation associated with a ground impact.

[0022] In an embodiment, the corresponding physics-based behavior of the armored vehicle may be, but is not limited to, the deformation of the armored vehicle when subjected to deformations associated with absorbing an impact.

[0023] In an embodiment, the corresponding physics-based behavior of the manufacturing equipment may be, but is not limited to, deformation of one or more pieces of the manufacturing equipment during operation.

[0024] Yet another embodiment is directed to a system for determining the physical behavior of a real-world object. In such an embodiment, the system includes a processor and a memory having computer code instructions stored thereon. The processor and memory are configured to use the computer code instructions to cause the system to implement any embodiment or combination of embodiments described herein.

[0025] Another embodiment is directed to a cloud computing implementation for determining the physical behavior of real-world objects. Such an embodiment is directed to a computer program product executed by a server in communication with one or more clients over a network. The computer program product comprises a computer-readable medium having program instructions that, when executed by a processor, cause the processor to implement any embodiment or combination of embodiments described herein. [Brief explanation of the drawings]

[0026] The foregoing will be apparent from the following more particular description of exemplary embodiments, as illustrated in the accompanying drawings, in which like reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale, emphasis instead being placed upon illustrating embodiments.

[0027] [Figure 1] FIG. 1 is a flowchart of an exemplary topology optimization workflow. [Figure 2] Figure 2 shows the representation of bifurcation by large deformation and geometric nonlinear modeling. [Figure 3] Figure 3A shows a two-dimensional (2D) computer-based model of a quadratic four-node element with integration points in an undeformed configuration, according to one embodiment. Figures 3B-3E are plots showing the condition numbers calculated at the element integration points of the 2D quadratic four-node element of Figure 3A as a function of the deformation gradient matrix F for four different types of kinematic deformation. [Figure 4] FIG. 4 is a flowchart of a nonlinear CAE workflow according to one embodiment. [Figure 5] FIG. 5 is a flowchart of a method for determining the physical behavior of a real-world object according to an embodiment. [Figure 6] FIG. 6 illustrates the transformation of material particles P from an undeformed state to a deformed state, according to one embodiment. [Figure 7] FIG. 7 is a flowchart of an optimization method according to one embodiment. [Figure 8] FIG. 8 shows how the embodiment can be incorporated into an existing geometrically nonlinear finite element. [Figure 9] 9A-9D show two-dimensional (2D) and three-dimensional (3D) C-shaped structures without added artificial building material modeling and with artificial building material modeling according to one embodiment. [Figure 10] FIG. 10 shows 2D simulation results determined using existing methods and 2D simulation results determined using an embodiment. [Figure 11]FIG. 11 shows 3D simulation results determined using existing methods and 3D simulation results determined using an embodiment. [Figure 12] FIG. 12 is a flowchart of a sensitivity-based design optimization workflow, according to one embodiment. [Figure 13] Figure 13 shows a benchmark model for geometric nonlinear topology optimization in which both elastic and elastoplastic material properties are defined. [Figure 14] 14A and 14B show the results of a geometrically nonlinear topology stiffness optimization for an elastic physical material solved using existing methods and embodiments, respectively. [Figure 15] FIG. 15 shows two elastic topology optimization designs determined using existing methods and one embodiment. [Figure 16] 16A and 16B show a nonlinear transient 3D topology optimization model. [Figure 17] Figure 17A shows a design density layout of a 3D model that fails optimization performed using existing methods. Figure 17B is an optimization iteration history plot of the 3D model of Figure 17A showing the objective function value of normalized maximum vertical displacement and the normalized volume constraint value as a function of optimization iteration. Figure 17C is a 3D model showing stresses of a failed primal solution of the 3D model of Figure 17A, with shading indicating localized artificial instability. Figure 17D is a 3D model showing stresses of the models of Figures 17A and 17C when numerical stabilization according to one embodiment is applied. [Figure 18]Figure 18A shows the optimized design density layout of the converged topology optimization. Figure 18B is an optimization iteration history plot of the 3D model of Figure 18A showing the objective function value of normalized maximum vertical displacement and the normalized volume constraint value as a function of optimization iteration. Figure 18C depicts a 3D model where shading indicates the stresses of the primary solution of the initial design at the maximum displacement point. Figure 18D depicts a 3D model where shading indicates the stresses of the primary solution of the optimized design at the maximum displacement point to which an embodiment is applied. [Figure 19] FIG. 19 is a simplified block diagram of a computer system for automatically determining the physical behavior of real-world objects, according to one embodiment. [Figure 20] FIG. 20 is a simplified block diagram of a computer network environment in which embodiments may be implemented. DETAILED DESCRIPTION OF THE INVENTION

[0028] A description of an exemplary embodiment follows.

[0029] Many simulations using computer-based models representing real-world objects experience divergence and local numerical instability when determining a primal solution for a model (or portion thereof) that undergoes large displacements and / or deformations. Lack of convergence to a primal solution is particularly common in the field of geometrically nonlinear topology optimization, where stiffness operators between solid and void elements can be significantly ill-conditioned. Typically, the numerical ill-conditioning of finite element stiffness operators for assembled systems is artificially caused by some distorted elements in the model. Embodiments solve this problem by locally stabilizing ill-conditioned portions of the model by adding additional artificial constitutive modeling to the physical constitutive modeling, without changing the overall stiffness of the system.

[0030] Previous attempts to achieve numerical stabilization modify the mathematical description of how materials respond to various loads and / or deformations, but these modifications are not related to the condition number of the element stiffness operator. These previous approaches add an overall stiffness to all low-stiffness (i.e., flexible) elements relative to the stiffness operator of the primal solution. This existing approach often manipulates stiffness scaling using a linear scaling factor, but it does not directly solve the problem of ill-conditioned operators and therefore has several drawbacks. For example, heuristic schemes must be applied, and it is not intuitive to identify how much stiffness should be added to flexible elements to make the stiffness operator well-conditioned. Furthermore, if too little stiffness is added to flexible elements, the stiffness operator will still be ill-conditioned. Conversely, if too much stiffness is added to flexible elements, the overall primal solution is likely to change excessively compared to the original solution of the physical system.

[0031] The problem of local numerical instability is particularly prominent in geometrically nonlinear topology optimization. In typical geometrically nonlinear topology optimization methods, a solid-void material layout is enforced, and typically the stiffness ratio between the void material and the solid material is in the range of 10 -9 ~10 -6 , which is within the range of . As a result, the void elements and intermediate density elements are often ill-conditioned with respect to the stiffness operator and therefore distort very significantly while solving the primal solution. Therefore, the solution for the primal solution may often not be determined. As a result, the entire optimization workflow does not converge, and a converged optimized design cannot be obtained.

[0032] FIG. 1 is a flowchart of an exemplary stabilization approach, method 100, for addressing local artificial bifurcation numerical instability in stiffness operators for geometrically nonlinear topology optimization. Method 100 begins in step 101 by receiving or otherwise defining an initial design for the object's design variables, i.e., material layout. In one embodiment, the initial design is represented in a computer-based model in step 101. Next, a finite element analysis 102 of a principal solution involving geometrically nonlinear modeling is performed using the material layout from step 101. According to one embodiment, the principal solution is the result of the FEM analysis. In an embodiment, specific values ​​may be taken from the principal solution for the optimization formulation, e.g., specific displacements in the model during loading. Subsequently, a sensitivity analysis 103 is applied, and then the design variables are updated 104 based on the sensitivity analysis 103. Thereafter, the constituent materials of the model initially defined from step 101 are updated 105 using the updated material layout from step 104. Next, method 100 checks 106 whether the optimization has converged. If the optimization has not converged, method 100 returns to step 102 and performs another optimization iteration, repeating steps 102, 103, 104, 105, and 106. If it is determined in step 106 that the optimization has resulted in a converged design, method 100 results in a design 107 that is optimized with respect to the design variables, i.e., material layout, and method 100 ends.

[0033] Each step 101-107 of method 100 represents a different possible step in a topology optimization workflow. Note that both existing stabilization approaches and embodiments can be implemented in method 100. For example, the embodiments described herein can be implemented in block 102, which is a step of finite element analysis.

[0034] Element removal has been attempted to improve geometrically nonlinear topology optimization problems with ill-conditioned stiffness operators for void elements and medium-density elements. This existing element removal approach can be implemented in the constituent material update step 105 of the stabilization approach shown in method 100. The element removal stabilization approach depends on the design variables and does not depend on the deformation of the primal solution. Therefore, the existing element removal approach is completely different compared to the embodiment that depends on the deformation of the primal solution.

[0035] Although the element removal approach removes low-stiffness elements from the material layout of the primary solution, the element removal approach has drawbacks. For example, in the element removal method, it is not clear what threshold should be applied to the design variables to determine which elements to remove, thereby removing many elements that do not need to be removed to obtain a well-conditioned stiffness operator. Furthermore, the element removal approach has a slow process of reintroducing previously removed elements during optimization iterations, thereby causing many optimization iterations and resulting in high computational costs.

[0036] The fictitious domain approach has also been used in an attempt to improve geometrically nonlinear topology optimization problems with ill-conditioned stiffness operators for void and medium-density elements. This fictitious domain approach is also implemented in the constituent material update step 102. Existing fictitious domain stabilization approaches depend on the design variables and are independent of the deformation of the primal solution.

[0037] The virtual domain approach converts low-stiffness elements into elements with linear stiffness operators based on the design variables of the material layout. The problem is that when using the virtual domain approach, it is not clear which scheme and parameters within the scheme should be applied to convert the original stiffness operators into linear stiffness operators. Furthermore, for very large deformations, the parameters of the virtual domain scheme must be gradually modified during the optimization iterations, which results in many optimization iterations and high computational costs.

[0038] The relaxed convergence criteria approach is yet another past attempt to improve geometrically nonlinear topology optimization problems that have ill-conditioned stiffness operators for void elements and medium-density elements. The relaxed convergence criteria approach can be implemented in the finite element analysis step 102 of the method 100. The relaxed convergence criteria approach relaxes the convergence criteria for the main solution, but does not solve the problem of ill-conditioned stiffness operators. Furthermore, relaxing the convergence criteria for the main solution has drawbacks. For example, it is unclear which scheme and parameters within the scheme should be applied to relax the convergence criteria for the main solution. Furthermore, when using the relaxed convergence criteria approach, finite element nodes attached to flexible elements can initiate artificial vibrations. The artificial vibrations can cause many optimization iterations, resulting in high computational costs. The artificial vibrations can also merge with nodes of other elements, preventing the main solution from converging.

[0039] Still further, plastic material interpolation has been attempted to improve geometrically nonlinear topology optimization problems that have stiffness operators that are ill-conditioned for void and intermediate density elements. Plastic material interpolation can be implemented in the finite element analysis step 102 of the method 100 of FIG. 1 . The plastic material interpolation approach applies a lower penalization for plastic yielding and / or plastic hardening than for the elastic Young's modulus. The plastic material interpolation approach differs from the embodiments described herein because these penalizations are a function of the design variables and not the deformation of the primal solution. The plastic material interpolation approach has drawbacks. For example, the plastic material interpolation approach can only be applied to structures where parts are made of materials that can be represented by an elastic-plastic material model, as opposed to all feasible constitutive material models.

[0040] The hyperelastic substrate approach has also been attempted to improve geometrically nonlinear topology optimization problems that have ill-conditioned stiffness operators for void elements and intermediate-density elements. This hyperelastic substrate approach can be implemented in step 102 of the finite element analysis. In this approach, hyperelasticity (e.g., a rubber material) is applied to the substrate of the design domain. The hyperelastic substrate approach has drawbacks. For example, it can only be applied to structures whose components are made of materials that can be represented by a hyperelastic model. Therefore, the hyperelastic substrate approach is not suitable for addressing all feasible constitutive material models. Furthermore, the hyperelastic substrate approach adds a general overall stiffness to the void elements, whereas the pseudo-constitutive modeling approach described herein is based directly on the condition number.

[0041] Therefore, existing methods (1) only increase the stiffness of flexible elements and are not linked to the condition number of the stiffness operator; (2) apply to optimization workflows that depend on the design variables and not on the deformation of the primal solution; and (3) do not address the core issue of numerically ill-conditioned stiffness operators but instead attempt to avoid quantifying the impact of numerically ill-conditioned stiffness operators on the primal solution. Thus, while previous methods have attempted to manipulate constitutive modeling to numerically stabilize the global stiffness operator, existing methods do not directly manipulate constitutive modeling in relation to the condition number of the element stiffness operator.

[0042] FIG. 2 illustrates an example of different bifurcations caused by large deformations during geometric nonlinear modeling 201. Specifically, FIG. 2 illustrates the difference between global physical buckling instability 202 and local artificial numerical instability 203 caused by a local ill-conditioned stiffness operator. Global buckling is not addressed by embodiments. Instead, embodiments suppress the problem of local artificial numerical instability 203 caused by a local ill-conditioned stiffness operator. This problem is illustrated in FIG. 2 for a computer-based model 204 having a rigid layout 205 and stress distribution 206. When model 204 is subjected to load 207, this causes significant deformation of element 208 and artificial numerical instability, as shown in a close-up 209 of stress distribution 206 of computer-based model 204.

[0043] Large deformations (deformations large enough to negate linearization of the deformation field) and geometrically nonlinear modeling are important in many applications, such as realistic modeling of real-world objects. When a system / object has small deformations and distortions, equilibrium can be determined in the undeformed geometry while still achieving physically accurate results. In contrast, when a system / object has large deformations and geometrically nonlinear behavior, equilibrium is determined in the deformed geometry to achieve physically accurate results. For such systems, geometric modeling is required to achieve accurate results while simulating large deformations, and large deformations and geometrically nonlinear modeling introduce bifurcation points in the modeling response. For example, these bifurcation points can be important in realistic modeling of pre- and post-buckling 201 and ultimately in accurately predicting the stability of the system. Typically, in a global buckling 201 scenario, these physical bifurcation points result in ill-conditioned stiffness operators, which can be numerically stabilized using quasi-static modeling, enforced displacement loading, and the arc-length method, among other examples, to obtain a physical primal solution.

[0044] Large deformation and geometrically nonlinear modeling can excessively distort low-stiffness regions (e.g., elements 208) in the model, which can lead to numerically artificial instabilities 203 due to ill-conditioned local element stiffness operators in finite element simulations. Typically, these low-stiffness regions consist of a few finite elements and have little or no effect on the quality of the predicted global physical buckling modeling and the predicted stability of the system. Often, these numerically artificial instabilities cause large deformation and geometrically nonlinear analyses to fail to achieve convergence on a primal solution. This failure can occur because these numerically artificial instabilities can excessively deform a small number of finite elements, making numerical integration of these finite elements infeasible. Furthermore, some of these finite element stiffness operators are indefinite, or even negatively definite, resulting in artificial local bifurcations and large, unphysical displacements. Furthermore, iterative solver schemes applied to geometrically nonlinear analyses continue to oscillate within a primal solution due to ill-conditioned local stiffness operators, thereby preventing a converged solution.

[0045] The embodiments described herein implement a novel approach that adds artificial constitutive material modeling to make stiffness operators well-conditioned and suppress numerical artificial instabilities.

[0046] Figure 3A shows an undeformed 2D quadratic four-node element 301 (with finite element nodes 301a-d) fully integrated using four integration points 308a-d. Figures 3B-3E are plots 302, 303, 304, and 305, respectively, showing the condition numbers 306a-d of element 301 as a function of the deformation gradient matrix F 307a-d, assuming plane strain. The condition numbers 306a-d are calculated at the four numerical integration points 308a-d for four different kinematic deformation modes of element 301: uniaxial compression (plot 302), simple shear (plot 303), volumetric compression (plot 304), and diagonal compression (plot 305). These plots 302, 303, 304, and 305 each show the condition numbers 306a-d of the stiffness operator for small to large deformation rates 307a-d, respectively. Plots 302, 303, 304, and 305 each have a respective visualization 309a-d of the deformation of element 301 as a function of deformation rate 307a-d.

[0047] Plots 302, 303, and 304 show that the three deformation modes result in numerically reasonable condition numbers 306a, 306b, and 306c for both small and excessive deformations. In contrast, the deformation in plot 305 results in a very high condition number 306d, causing the system to become numerically unstable. Plot 305 also shows that condition number 306d increases significantly when the surface of the deformed element approaches one of the element's integration points (in this example, integration point 308c). This observation can also be confirmed for other element types, numerical integration point approaches (e.g., different locations of the integration points), and deformation modes.

[0048] Many simulation applications, typically from the transportation and traffic sectors (e.g., automobiles), aerospace and defense (e.g., airplanes), life sciences (e.g., medical devices), and high-tech (e.g., mobile phones and laptops) apply large-deformation nonlinear finite element modeling to obtain physically correct results when solving for a primal solution. However, these finite element models sometimes fail during convergence on the primal solution because a very small number of elements out of thousands or even millions of elements have local artificial numerical bifurcations that cause stiffness operators to become ill-conditioned. These local artificial numerical bifurcations can occur both when solving for a primal solution alone and when the primal solution is part of an iterative optimization workflow.

[0049] The embodiments address the problems experienced by the aforementioned industries by introducing an additional artificial pseudo-constitutive model that is added to a geometric CAE model (e.g., a finite element model) to eliminate numerically unstable and ill-conditioned element stiffness operators, where the artificial pseudo-constitutive model is based on the condition number of the element stiffness operator.

[0050] FIG. 4 is a flowchart of a workflow 400 for simulating a real-world object to determine the object's behavior. The user-defined objective function of the optimization workflow is to minimize displacements in elements of a finite element model. Workflow 400 begins in step 401 by defining a computer-based model, e.g., a finite element model representing the object. After defining the model in step 401, method 400 proceeds to steps 402 and 405. In step 402, a physical constitutive model is defined. Meanwhile, in step 405, an artificial constitutive model is added to the finite element model defined in step 401. Next, method 400 adds artificial internal forces in step 406 to the artificial constitutive model from step 405. According to one embodiment, the artificial internal forces added in step 406 are based on stresses calculated using the artificial constitutive model in step 405. Subsequently, method 400 adds an artificial stiffness operator to the artificial constitutive model in step 407 based on the artificial internal forces of step 406. The internal physical forces are added to the physical constitutive model defined in step 402 in step 403. The internal physical forces applied in step 403 are based on stresses calculated using the physical constitutive model from step 402. The physical stiffness operators are added to the internal physical forces from step 403 in step 404. The stiffness operators added in step 404 are based on the internal physical forces from step 403. Both the internal physical forces (applied in step 403) and the artificial internal forces (applied in step 406), and the external forces, are combined into a remainder in step 408. In this manner, step 408 uses the external and internal forces from both the physical system (forces from step 403) and the artificial system (forces from step 406) to combine the remainder. Similarly, both the physical stiffness operators (from step 404) and the artificial stiffness operators (from step 407) are combined into stiffness operator 409. In this manner, step 409 assembles a stiffness operator using both the stiffness operator from the physical system (the stiffness operator from step 404) and the stiffness operator from the artificial system (the stiffness operator from step 407). Method 400 then uses the assembled global stiffness operator 409 to solve for the primal solution and equilibrium of the finite element model in step 410.Workflow 400 then checks whether the main solution has converged in step 411. If the solution has not converged, method 400 returns to steps 402 and 405 and repeats steps 402-411 of workflow 400. If it is determined in step 411 that the main solution has converged, workflow 400 yields a converged main solution in step 412 and method 400 ends.

[0051] In workflow 400, steps 401, 402, 403, 404, 411, and 412 are steps that can be part of a standard nonlinear CAE workflow (a finite element modeling workflow). In addition to steps 401, 402, 403, 404, 411, and 412, method 400 includes a step for eliminating numerical instabilities and ill-conditioned element stiffness operators. Specifically, method 400 adds artificial constitutive modeling in step 405 (adding an artificial pseudo-constitutive model using the condition number of the nonlinear deformation matrix F at the finite element integration points), adds artificial internal forces in step 406, and adds artificial stiffness operators in step 407. Thus, method 400 includes steps 405, 406, and 407 in addition to the steps of an existing CAE workflow. In method 400, in step 408, a residual defining structural equilibrium is assembled using the physical internal forces from step 403 and the artificial internal forces from step 406. The physical stiffness operator 404 and the artificial stiffness operator 407 are assembled at 409 to result in an overall structural stiffness operator. A principal solution is then determined at step 410 using the assembled remainders from step 408 and the assembled stiffness operator from step 409. A check is performed at step 411 to determine if the solution has converged. According to one embodiment, the check at step 411 uses the results / values ​​from the principal solution (determined at step 410) to evaluate the user-defined objective function and constraints. If the solution has not converged, additional solver iterations are performed and method 400 returns to steps 402 and 405; if the solution has converged, method 400 moves to step 412.

[0052] Numerical solvability in finite element analysis often depends on the condition number of the stiffness operator. If the deformation matrix F is ill-conditioned, stress calculations can be numerically unstable and the stiffness operator can be ill-conditioned, often causing divergence of the solution for the principal solution. An embodiment introduces local stability pseudo-constitutive modeling using the condition number of the nonlinear deformation matrix F at the finite element integration points. In one embodiment, a constitutive material law is defined that depends on the condition number F, and this constitutive material law induces stresses that minimize the condition number of the stiffness operator in the finite element setup, as shown in FIG. 4. The applied artificial stresses are thereby introduced only in directions that minimize the condition number of the stiffness operator for a few key elements; the applied artificial stresses do not simply increase the stiffness of a given element in all directions.

[0053] Embodiments, for example, method 400, address the numerical instability of local artificial bifurcations caused by ill-conditioned element stiffness operators, and do not affect or modify the overall physical buckling (physical bifurcations) of the system caused by geometric nonlinear modeling.

[0054] Among other advances, embodiments improve over the prior art by implementing rigorous numerical computation of primal solutions for geometrically nonlinear finite element models with ill-conditioned, low-stiffness elements, enabling embodiments to determine physical primal solutions while requiring fewer solver iterations to obtain the primal solution, thereby reducing computational costs.

[0055] Furthermore, for geometrically nonlinear topology optimization workflows where local numerical instability issues are particularly pronounced, embodiments can enforce a solid-void material layout with a high stiffness ratio between the void material and the solid material. This approach can result in convergence of a dominant solution for every optimization iteration. Furthermore, forcing a solid-void material layout with a high stiffness ratio allows the design to be determined without the dominant solution failing, requiring fewer solver iterations per optimization iteration, thereby resulting in lower computational costs. Furthermore, in embodiments, the heuristic scheme and heuristic parameters are not modified during the optimization iteration depending on the design variables, which results in fewer optimization iterations being performed, and therefore, embodiments can require fewer solver iterations per optimization iteration, resulting in lower computational costs.

[0056] Local artificial numerical instabilities often appear during the deformation process, frequently causing non-convergence in the solution scheme for the finite element residual primal solution. Previous approaches to solving local artificial numerical stability problems result in convergence problems in the primal solution. Therefore, these previous approaches cannot numerically obtain a complete primal solution. Furthermore, previous approaches generate oscillations in the convergence history of the primal solution, and the residual is not satisfied for some elements out of thousands to millions of elements.

[0057] The embodiment directly adds a condition number for the stiffness operator to each finite element in the model, thereby strongly coupling the embodiment to the numerical stability of the finite element model. The embodiment also defines a constitutive material law that depends on the condition number of the local element deformation, which induces stresses that minimize the condition number. Furthermore, the embodiment is enforced in local regions (i.e., a few elements) as needed. This approach is also more general because it is a function of the local deformation and therefore does not require a second field (i.e., a second set of design variables). Therefore, the embodiment ensures stable convergence of the primal solution.

[0058] 5 is a flowchart of a method 500 for determining the physical behavior of a real-world object. Method 500 includes eliminating numerical instabilities and ill-conditioned element stiffness operators, according to one embodiment. Method 500 is computer-implemented and, as such, may be implemented by one or more processors executing computer code instructions.

[0059] Method 500 begins in step 501 by defining a computer-based model in a memory of a processor implementing the method. The computer-based model defined in step 501 represents a real-world object using a plurality of elements, and indicates one or more materials represented by each element of the plurality of elements of the real-world object. The computer-based model defined in step 501 is any one or combination of a finite element model, a boundary element method, a finite difference method, a finite volume method, or a discrete element method. Furthermore, the real-world object defined in step 501 is any one of an automobile, industrial equipment, an aircraft, a civil structure, a marine device, a medical device, a consumer product, an electronic device, an armored vehicle, or a manufacturing facility, among other examples. Equations describing the physics-based behavior of the one or more materials are then defined in step 502. Following step 502, stabilization equations are defined in step 503. According to one embodiment, the stabilization equations are a function of a nonlinear transformation matrix F. A simulation of the real-world object being subjected to the load is then performed in step 504 using (i) the defined computer-based model from step 501, (ii) the defined equations describing the physics-based behavior from step 502, and (iii) the defined stabilization equations from step 503.

[0060] Embodiments of method 500 further include, for each element of the plurality of elements in the defined model (defined in step 501), associating an artificial internal force with the element based on the defined stabilization equation (defined in step 503). Such embodiments may further include, for each element of the plurality of elements, associating an artificial-based and physics-based behavior (i.e., a combination of artificial-based and physics-based behavior) with the element based on the associated artificial internal force.

[0061] The simulation of the real-world object subjected to the loads performed in step 503 may be performed using artificial internal forces associated with each element, and artificial-based and physics-based behaviors associated with each element. Further, in embodiments of method 500, performing the simulation of the real-world object using the defined stabilization equation (in step 504) reduces a condition value associated with a given element of the plurality of elements. Further, in such embodiments, the condition value is a function of the nonlinear deformation gradient matrix.

[0062] Yet another embodiment of method 500 updates the computer-based model based on the results of the simulation performed in step 504. Such an embodiment may perform a subsequent simulation to determine updated physical behavior of the real-world object. According to one embodiment, this subsequent simulation is performed using the updated computer-based model, the defined equations describing the physics-based behavior, and the defined stabilization equations. Yet another embodiment may iterate, i.e., repeat, (i) updating the computer-based model based on the determined updated physical behavior and (ii) determining the updated physical behavior until the determined updated physical behavior meets user-defined criteria.

[0063] In one embodiment of method 500, a model is defined in step 501 by taking measurements from one or more real-world objects. In this manner, the model defined in step 501 represents the real-world objects for which the measurements were taken. In such an embodiment, method 500 can be used to determine the behavior of the real-world object and to determine design improvements for the real-world object itself. In such an embodiment, the measurements of the object are used in step 501 to define a model, and equations are defined in steps 502 and 503 based on the properties of the real-world object. A simulation is then performed in step 504 as described herein to determine results indicative of the behavior of the real-world object. Based on these results, design changes to the real-world object can be determined, the real-world object itself can be modified according to the determined changes, or the real-world object can be manufactured according to the determined changes. Among other exemplary implementations, such functionality can be used to determine improvements to the crashworthiness of a vehicle. For example, an embodiment can determine during a simulation that a structural member of, say, an automobile will fail when subjected to deformations associated with a 60 mile-per-hour crash. Such embodiments can be iterative, for example, using a model of an automobile having thicker structural members, until a thickness of the structural members that will withstand collision-related deformation is determined. Such embodiments can then continue by manufacturing the automobile using structural members with the determined thickness. This process allows embodiments to determine appropriate design changes and manufacture real-world objects with the determined design without the typical requirement for human trial and error.

[0064] A specific application of geometric nonlinear topology optimization involves modeling intermediate or void elements when the finite element model has stability issues. This embodiment (i) is computationally less expensive than previous methods, and (ii) results in fewer solver iterations and fewer optimization iterations, resulting in increased computational efficiency. Furthermore, the solution of the primal solution does not fail, and a fully optimized design is achieved.

[0065] The present embodiment supports numerical stabilization of condition numbers for (i) all types of constitutive models applied to modeling physical systems, and (ii) all types of shape functions applied to large deformation finite element modeling, thereby enabling the embodiment to be utilized for all finite element types, for example, all physics modeling involving large deformations such as static, quasi-static, and transient modeling, and all types and combinations of numerical solver methods applied to determine the primary finite element solution.

[0066] 6 shows a plot 600 illustrating the path line 605 of a material particle P601 as it moves from its reference configuration 606 (P601) to its deformed configuration 607 (p603). The plot 600 shows that the deformation gradient matrix F608 is defined according to classical continuum theory, and the transformation of material particle P601 from its undeformed state X602 to a deformed p603 has a deformation state x604.

[0067] The structural analysis describes the path taken by material particles within the structure when a load is applied, as seen in plot 600. A material particle P 601 is at position x 602 in a reference configuration 606 and deforms to a new position x 604 during the analysis, with a displacement u(X) 605. Since conservation of mass is assumed, particles cannot disappear, and therefore the particle

[0068]

number

[0069] is mapped as a bijection of

[0070] From this, the infinitesimal vectors in the deformed configuration can be expressed as:

[0071]

number

[0072] The transformation matrix is ​​defined as the deformation gradient matrix F as follows:

[0073]

number

[0074] The embodiment defines a strain energy potential F at a particular material point that correlates with the condition number of the deformation gradient. The motivation is that an ill-conditioned deformation gradient F can be calculated using the right Cauchy-Green deformation tensor C=F T This makes F ill-conditioned. Thus, small changes in the deformation gradient result in large changes in the right-hand Cauchy-Green deformation tensor. As a result, the stresses and nodal forces calculated from the right-hand Cauchy-Green deformation tensor exhibit similar sensitivity to the deformation gradient entries. This ill-conditionedness results in oscillations in the stresses at the material points and, ultimately, in the residual nodal forces of the Newton-Raphson iteration scheme, which hinders convergence. Therefore, the objective of the hyperelastic material model implemented by the embodiment is to prevent the system from reaching a state with ill-conditioned deformation gradients at the integration points.

[0075] Frobenius norm || || F This norm is applied because the condition number is induced by two deviation invariants of the deformation gradient F of the right Cauchy-Green deformation tensor.

[0076]

number

[0077] and

[0078]

number

[0079] and the right Cauchy-Green deformation tensor C, which is chosen to be expressed in terms of:

[0080]

number

[0081]

number

[0082] The condition number F of the deformation gradient is induced by the Frobenius norm and is given by:

[0083]

number

[0084] These Frobenius norms are T F and F -t F -1 can be expressed as the square root of the trace of -T =((...) T ) -1 =((...) -1 ) T and we get:

[0085]

number

[0086] The trace of the matrix product of two matrices A and B is given by Tr(AB)=Tr(BA). Therefore, κ F The last term in (F) can be rewritten as:

[0087]

number

[0088] Two invertible matrices A, B∈R n×n is represented by the following formula: (AB) -1 =B -1 A -1 and κF The last term in (F) can be rewritten as:

[0089]

number

[0090] Definition of the right Cauchy-Green deformation tensor C=F T Using F, this equation gives:

[0091]

number

[0092] Trace of a square matrix A∈R n×n is the sum of its eigenvalues

[0093]

number

[0094] is.

[0095] Furthermore, the eigenvalues ​​A of an invertible matrix and its invertible A -1 The relationship between

[0096]

number

[0097] is.

[0098] This is trivial when considering the eigenvalue decomposition of these two matrices. F Applying this to (F) gives:

[0099]

number

[0100] The eigenvalues ​​C are the principal stretches λ squared in the finite strain theory. i 2 :

[0101]

number

[0102] is.

[0103] The definition of the main stretch is kappa F Applying this to (F) gives:

[0104]

number

[0105] Rearranging the definition of the deviational primary stretch relative to the primary stretch gives:

[0106]

number

[0107] κ F λ in (F) i of

[0108]

number

[0109] Substituting in, we get:

[0110]

number

[0111] invariants

[0112]

number

[0113] and

[0114]

number

[0115] is defined as follows:

[0116]

number

[0117]

number

[0118] κ F (F)

[0119]

number

[0120] and

[0121]

number

[0122] Substituting, we finally get:

[0123]

number

[0124] As a result, the upper bound on the condition number κ_F(C) is estimated as follows:

[0125]

number

[0126] The following points: F (F) is the main stretch λ i The deviation stretch

[0127]

number

[0128] This means that the condition number κ F (F) and kappa F (C) implies that they are independent of volumetric stretch, so they do not affect numerical stability to the constituent material model level. However, excessive volumetric compression or tension can lead to convergence problems due to round-off errors in other parts of the finite element method.

[0129] One embodiment is a condition number κ F (C) shows the possibility of candidate artificial materials that can cope with two deviation invariants.

[0130]

number

[0131] and

[0132]

number

[0133] However, other candidates may also be applied.

[0134] A possible set κ of general artificial building material candidates that addresses condition numbers F (C) is the two deviation invariants

[0135]

number

[0136] and

[0137]

number

[0138] of C, and the third invariant J el of C, which is defined by the following formula:

[0139]

number

[0140] Thereby, one set of artificial constitutive material models (N=1, N=2, etc.) that can be implemented for a condition number based on the strain energy potential is:

[0141]

number

[0142] is.

[0143] During the ceremony,

[0144]

number

[0145] incorporates the condition number into the strain energy potential. The strain energy potential applied in the examples described below is a simplified version of this form:

[0146]

number

[0147] , the following equation is obtained:

[0148]

number

[0149] The disclosed proofs and equations show that embodiments rely on a deformation matrix F that represents the kinematics of a physical system, rather than on a constitutive material model of the physical system. The embodiments thereby support numerical stabilization of physical systems consisting of all types of constitutive models (elastic, plastic, viscoelastic, viscoplastic, hyperelastic, damage, geomechanical, micromechanical, etc.), as well as isotropic, orthotropic, and anisotropic models, and further, the embodiments are not limited to these physical constitutive material models.

[0150] The disclosed proofs and equations also show that the embodiments support all types of shape functions applied in large deformation finite element modeling. Therefore, the embodiments may be used with all finite element types, including 2D or 3D elements, as well as first-, second-, and higher-order elements, triangular, quadrilateral, tetrahedral, hexahedral, pentahedral, pyramidal elements, etc., and continuum, shell, membrane, and beam elements. Furthermore, it should be noted that the embodiments are not limited to these element types. The disclosed proofs and equations also show that the embodiments support all types of physics modeling, including large deformations for static, quasi-static, and transient modeling.

[0151] The disclosed proofs and equations also show that embodiments support any type and combination of solvers, including implicit, explicit, and hybrid methods, operator factorization using direct solvers or iterative solvers with preconditioners, Newton-Raphson or incremental iterative techniques, incremental loading, first order methods (e.g., backward Euler), second order methods (e.g., Newmark beta), or higher order methods (e.g., Runge-Kutta).

[0152] FIG. 7 is a flowchart of a workflow 700 for simulating a real-world object to determine the object's behavior. Workflow 700 is a theoretical implementation of workflow 400 described herein above in connection with FIG. 4. Workflow 700 begins in step 701 by defining a computer-based model, e.g., a finite element model. After defining the model in step 701, method 700 proceeds to steps 702 and 705. In step 705, an artificial constitutive model is added to the finite element model defined in step 701. Next, in step 706, artificial internal forces are applied to the artificial constitutive model defined in step 705. The artificial internal forces added in step 706 include forces based on stresses calculated using the artificial constitutive model of step 705. Embodiments then add artificial stiffness operators to the artificial constitutive model in step 707 based on the artificial internal forces of step 706. In another path of workflow 770, in step 702, a physical constitutive model is added to the finite element model defined in 701. The internal physical forces are then added to the physical constitutive model at step 703. The internal physical forces are based on the stresses calculated using the physical constitutive model from step 702. To continue, a physical stiffness operator is added to the internal physical forces (added at step 703) at step 704. The physical stiffness operator added at step 704 is based on the internal physical forces added at step 703. Both the internal physical forces applied at step 703 and the artificial internal forces applied at step 706 are assembled into a remainder at step 708. Similarly, both the physical stiffness operator from step 704 and the artificial stiffness operator from step 707 are assembled into a stiffness operator at step 709. Then, at step 710, the method 700 solves the primal solution and equilibrium of the finite element model using the global stiffness operator assembled at step 709. The workflow 700 then checks whether the primal solution has converged at step 711. If the main solution has not converged, method 700 returns to steps 702 and 705 and repeats steps 702-711. If step 711 determines that the solution has converged, workflow 700 yields a converged main solution in step 712 and method 700 ends.

[0153] Workflow 700 of FIG. 7 is a theoretical implementation of workflow 400 described above in connection with FIG. 4. In workflow 700, steps 701, 702, 703, 704, 711, and 712 are steps that may be part of a standard nonlinear CAE workflow (finite element modeling workflow). In addition to steps 701, 702, 703, 704, 711, and 712, method 700 includes steps 705-710 for eliminating numerical instabilities and ill-conditioned element stiffness operators. Specifically, method 700 adds artificial pseudo-constitutive modeling in step 705 (using the condition number of the nonlinear deformation matrix F at the element integration points), adds artificial internal forces in step 706, and adds artificial stiffness operators in step 707. In method 700, a residual defining structural equilibrium is assembled in step 708 using the physical internal forces from step 703 and the artificial internal forces from step 706. The physical stiffness operator 704 and the artificial stiffness operator 707 are assembled in step 709 to result in an overall structural stiffness operator. A primal solution is then determined in step 710 using the assembled remainders from step 708 and the assembled stiffness operator from step 709. A check is performed in step 711 to determine if the solution has converged. If the solution has not converged, additional solver iterations are performed and the method 700 returns to steps 702 and 705; if the solution has converged, the method 700 moves to step 712.

[0154] The assembled remainder {R} of the system equilibrium defined in step 708 is the sum of the external forces {P} (e.g., applied loads), the internal forces {I} of the physical system, 人工的} (from step 703), and the internal forces {I 物理的} (from step 706). The remainder {R} assembled in 708 may also include other forces such as contact forces, inertial forces, and thermal loads, among other examples.

[0155] internal force {I 物理的} and {I 人工的} is a function of the field {U} of the primal solution and the element forces {i e} is determined by summing

[0156]

number

[0157] In the formula, [D 物理的 ] and [D 人工的 ] are the constitutive modeling of the physical system and the added artificial constitutive modeling, respectively. The physical constitutive modeling can also be a function of the state variables {α}, e.g., for plasticity.

[0158] Generally, the embodiments make the overall physical internal force significantly greater than the artificial internal force (I 物理的} >> {I 人工的}), the overall principal solution {U} is not changed due to the added artificial constitutive modeling. However, for a small number of artificially ill-conditioned low-stiffness elements, embodiments add artificial internal forces due to the artificial stresses of the added artificial constitutive modeling. Thereby, for these few elements at the element level, the element's physical internal forces are lower than the artificial element forces ({i 物理的,e} < {i 人工的,e}), but the overall primal solution {U} remains unchanged.

[0159] Therefore, the global stiffness operator [K] constructed in step 709 and used in step 710 to determine the geometrically nonlinear primal solution {U} that satisfies the remainder {R}={0} is derived using partial derivatives with respect to the primal solution {U} as follows:

[0160]

number

[0161] During the ceremony,

[0162]

number

[0163] The embodiment uses the physical element stiffness operator [k 物理的,e ] for some important elements, we use an artificial element stiffness operator [k 人工的,e ], the artificial element forces {i 人工的,e} and do not simply increase the stiffness in all directions of a given element.

[0164] The artificial element stiffness operator [k 人工的,e ] is the physical stiffness operator [K 物理的 ] does not introduce so much stiffness that the overall physical buckling (bifurcation) limit given by

[0165] The exemplary workflow 400 shown in FIG. 4 and the exemplary workflow 700 shown in FIG. 7 can be implemented in the source code of a finite element (or CAE) system. The embodiments can also be implemented in a non-intrusive manner if the existing finite element system already supports one or more artificial constitutive material models based on invariants for large deformation modeling. The non-intrusive manner was chosen for the implementation results described herein below, where the embodiments were implemented using the finite element solver Abaqus® provided by the applicant. An example of this non-intrusive implementation is shown in FIG. 8, in which each finite element of a physical model is duplicated into a new artificial finite element using an artificial constitutive material model. The duplicated artificial finite element has a new element ID but maintains the node IDs of the duplicated physical finite element. The artificial constitutive material modeling is thereby consistently added to the physical constitutive material modeling. The assembled finite element system superimposes the behavior of the original physical system and the added artificial stabilization to enforce well-conditioned global operators. The embodiments implemented in Abaqus® support all arbitrary element types.

[0166] Figure 8 shows an exemplary embodiment implemented on an existing geometrically nonlinear finite element solver. In Figure 8, finite element 801 represents the physical modeling of an element, finite element 802 represents the artificial modeling of the element, and finite element 803 represents the assembly of the physical finite element modeling 801 and the artificial finite element modeling 803. Each finite element of the physical model 801 is duplicated into a new artificial finite element 802 using an artificial constituent material model. The duplicated artificial finite element has a new element ID 805 (the edge boundary of the finite element) but maintains the node IDs (e.g., node 804) of the duplicated physical finite element, as shown in 803.

[0167] 9A-9D show exemplary C-shaped structures in both 2D models 901-902 and 3D models 903-904 used to verify the operation of embodiments, e.g., method 400. To verify the operation of embodiments, artificial constituent material modeling 905a-b was added within C-shaped structures 902 and 904 to stabilize the numerical instability of the finite elements in void regions 905a-b, which have significantly lower stiffness. The strictly solid (without artificial constituent material modeling) models 901 and 903 are reference models for comparison with the solid-void (with artificial constituent material modeling) models 902 and 904. Comparing the simulation results of models 901 and 903 (with open voids) with those of models 902 and 904 verifies the embodiment, as the void regions 905a-b do not affect the overall system behavior and response compared to the strictly solid models 901 and 903.

[0168] This validation was performed using 2D finite element modeling (models 901 and 902) and 3D finite element modeling (models 903 and 904). To determine these results (discussed below in connection with FIGS. 10 and 11), models 901-904 were fully constrained on the left side and loaded at two corners on the right side with forces F1=0.02 906 and F2=0.03 907 for 2D models 901 and 902, and F1=0.2 908 and F2=0.3 909 for 3D models 903 and 904. Models 901-904 were left undimensioned, using dimensions 10x10x1 and thickness 1 for 2D models 901 and 902, and dimensions 10x10x10 for 3D models 903 and 904.

[0169] The solid materials in Figures 9A-9D (fabricated with a C-shaped structure) have a Young's modulus E s = 1 and Poisson's ratio ν = 0.3. There are two different versions of the model: strictly solid models 901 and 903, and combined solid-void models 902 and 904. Solid models 901 and 903 contain only C-shaped solid regions. Solid-void models 902 and 904 consist of solid C-shaped regions, which in turn enclose significantly softer void regions 905a-b. The void materials 905a-b have a significantly lower Young's modulus, E v =10 -9 void regions should not affect the behavior and response of the overall system, but finite elements in the void regions are prone to numerical instability due to localized artificial numerical instabilities and bifurcations. Accordingly, Figures 9A-D show that embodiments can numerically stabilize the void elements of solid-void models 902 and 904 without substantially changing the determined behavior and response of the overall system. Therefore, the solutions and responses determined using solid models 901 and 903 are reference solutions for system behavior to verify that the solutions and responses determined using embodiments of models 902 and 904 match.

[0170] FIG. 10 shows 2D structures 1001a-d, 1002a-d, and 1003a-d used in simulations to verify four different element type embodiments (structures 1001a-d correspond to structure 901, structures 1002a-d and 1003a-d correspond to structure 902, and the embodiment was applied to structures 1003a-d). The four element types used were triangular 3-node first-order plane stress element (CPS3) 1004a, triangular 6-node second-order plane stress element (CPS6) 1004b, quadrilateral 4-node first-order plane stress element (CPS4) 1004c, and quadrilateral 4-node first-order plane stress element with reduced-order integration by hourglass control (CPS4) 1004d. The first column 1005 shows the strictly solid models 1001a-d as reference solutions for the principal solutions of the four element types 1004a-d. Second column 1006 shows the structural response of solid-void models 1002a-d without numerical stabilization. For models 1002a-d where numerical stabilization is not applied, the finite element solving principal solution for four element types 1004a-d fails numerically before the final deformation is determined (unlike the results shown in column 1007). Third column 1007 shows the structural response of solid-void models 1003a-d of four element types 1004a-d where numerical stabilization is applied in accordance with embodiments. The structural response in column 1007 numerically and substantially matches the final deformation shown in column 1005. Thus, application of embodiments to 2D models does not modify the results when there is no instability, but when there is instability, e.g., in column 1006, embodiments successfully determine the simulation results.

[0171] Figure 11 shows 3D structures 1101a-c, 1102a-c, and 1103a-c used in simulations to verify the implementation of three different element types (structures 1101a-c correspond to structure 903, structures 1102a-c and 1103a-c correspond to structure 904, and the implementation applies to structures 1103a-c). The three element types are tetrahedral 4-node linear element (C3D4) 1104a, hexahedral 8-node linear element (C3D8) 1104b, and hexahedral 8-node linear element with reduced-order integration with hourglass control (C3D8R) 1104c. The first column 1105 shows the strictly solid models 1101a-c as reference solutions for the primal solution. Second column 1106 shows the structural response of solid-void models 1102a-c where numerical stabilization was not implemented, and therefore the finite element solution of the primary solution failed numerically (unlike the results shown in column 1107) before the final deformation was determined. Third column 1107 shows the structural response of solid-void models 1103a-c where numerical stabilization was implemented in accordance with embodiments described herein. The results in column 1107 for models 1103a-c numerically and substantially match the results in column 1107 for strictly solid models 1101a-c. Thus, application of embodiments to the 3D models does not modify the results when there is no instability, but when there is instability, e.g., in column 1106, embodiments successfully determine the simulation results.

[0172] The numerical experiments and findings in Figures 10 and 11 validate the implementation of different element types (1004a-d and 1104a-c) and 2D and 3D elements. Furthermore, the implementation has no practical impact on the modeling quality of the predicted overall physical response and behavior. This is illustrated by the similarity of the results in columns 1005 and 1007, and columns 1105 and 1107. Second, the implementation allows successful simulations involving large deformations involving elements with different stiffnesses, in contrast to methods that do not apply stabilization. This is illustrated in Figure 10, where the simulation in column 1006 fails and the simulation in column 1007 succeeds, and in Figure 11, where the simulation in column 1106 fails and the simulation in column 1107 succeeds.

[0173] An application of the embodiment is in the context of geometric nonlinear optimization, where the condition number of the stiffness operator can change significantly from one optimization iteration to another. The problem of local numerical instability is particularly pronounced for geometric nonlinear topology optimization workflows because they enforce solid-void material layouts with high stiffness ratios between the void and solid materials. Figure 12 shows an example in which an embodiment of artificial constitutive modeling for numerical stabilization of void elements (flexible elements) can be implemented in a geometric nonlinear topology optimization workflow 1200.

[0174] 12 illustrates a workflow for a sensitivity-based design optimization process 1200, according to one embodiment. The design implementation workflow 1200 performs artificial configuration numerical stabilization 1202, as described herein above with respect to workflows 400, 500, and 700 of FIGS. 4, 5, and 7, respectively, along with additional steps 1203-1208.

[0175] The optimization workflow 1200 begins with an initial design described by an initial set of design variables in step 1201. Next, in step 1202, a primal solution including artificial constitutive modeling for numerical stabilization is determined for each optimization design iteration. Each optimization design iteration begins by taking the initial design variables from 1201 and adding the artificial constitutive model to the physical constitutive model in step 1202a. In step 1202b, the physical and artificial models from step 1202a are used to assemble residual and stiffness operators. Next, step 1202c uses the assembled global stiffness operators from step 1202b to solve the primal solution and equilibrium of the CAE model (finite element model). Step 1202d checks whether the primal solution has converged. If step 1202d determines that the solution has converged, the iteration results in the converged primal solution moving to step 1202e. If, at step 1202d, it is determined that the main solution has not converged, the workflow returns to step 1202a.

[0176] After substeps 1202a-e, a converged primal solution is determined and method 1200 moves to step 1203, where values ​​of the objective function and constraints are determined. Non-limiting examples of values ​​of the objective function and constraints include direct measurements of the model (e.g., mass, center of gravity, etc.) or measurements determined by the results of the primal solution (e.g., stiffness, displacement, force, stress, strain, etc.).

[0177] Next, objective function sensitivities and constraints are calculated in step 1204. In one embodiment, objective function sensitivities and constraints based on the primal solution are calculated in step 1204 using direct or adjoint methods and implemented in the Abaqus® kernel.

[0178] The optimization problem is then solved in step 1205 using mathematical programming, i.e., optimization calculations. According to one embodiment, the mathematical programming is strictly based on the values ​​of user-defined design goals relative to an objective function and constraints. It is further noted that in mathematics, computer science, and operations research, mathematical programming, alternatively named mathematical optimization or simply optimization, is the process of selecting the best solution (with respect to some criteria) from several available alternatives. Embodiments of method 1200 may use any such mathematical programming known in the art.

[0179] A new physical model for the next optimization iteration is generated in step 1206 based on the design variable values ​​found in step 1205. Note that in an embodiment, the design variables (determined in step 1205) and the physical model variables (updated in step 1206) may be the same as, for example, the thickness design variable for size optimization. However, the design variables and physical model variables may also be different, for example, for density topology optimization, where the design variable is relative density that maps to physical density and stiffness.

[0180] Step 1207 checks whether the optimization has converged. If the optimization has not converged, method 1200 returns to step 1202 and a new optimization iteration is initiated. If the optimization has converged, the optimization workflow 1200 moves to step 1208 and outputs the final design. According to one embodiment, the optimization workflow is converged when the convergence criteria for the change in the object is less than, for example, 0.1%, and the change in the design variables is less than, for example, 0.5%.

[0181] Figure 13 shows a benchmark model 1300 for geometric nonlinear topology optimization, of which two versions have been defined: one that applies strictly elastic material properties 1302, and one that applies elastoplastic material properties, i.e., elastic material properties 1302 and plastic material properties 1303. Model 1300 is clamped on the right side and has a load P 1301 applied in the y-direction at the midpoint of the right side. Model 1300 consists of 3,600 quadrilateral four-node first-order plane strain elements (CPE4).

[0182] An exemplary implementation applies optimization workflow 1200 to model 1300 for topology optimization. The optimization objective function minimizes the displacement in the load direction at force P 1301 subject to a 50% relative mass constraint for the design domain, thereby maximizing the stiffness of model 1300 for a given mass. The design domain is defined as model 1300 as a whole, where each finite element is a design variable in the topology optimization to determine a new optimized conceptual material layout in the design domain.

[0183] 14A and 14B show geometric nonlinear topology optimization results for the stiffness optimization problem for two versions of model 1300 described herein above in connection with FIG. 13. Optimization workflow 1200 was applied to determine the various optimization results in FIG. 14. Specifically, the results shown in columns 1403a-b were generated without applying numerical stabilization in step 1202. The results in columns 1404a-b were generated by applying numerical stabilization through the use of artificial constitutive modeling in step 1202. The results in columns 1403a-b and 1404a-b show different designs obtained for different load magnitudes. Numerical experiments show that for the geometrically nonlinear topology optimization results of Figures 14A and 14B, the strictly elastic optimization design 1401 (i.e., results generated using a version of model 1300 that applies strictly elastic material properties 1302) fails to converge at load magnitudes higher than P = 100 kN without applying numerical stabilization (column 1403a), and the elasto-plastic optimization design 1402 (i.e., results generated using a version of model 1300 that applies elasto-plastic material properties, i.e., elastic material properties 1302 and plastic material properties 1303) fails to converge at load magnitudes higher than P = 7.5 kN without applying numerical stabilization (column 1403b). In contrast, all optimization cases that applied numerical stabilization by adding artificial constitutive material modeling according to an embodiment of step 1202 converge in finding an optimized design (columns 1404a-b).

[0184] Figure 15 shows two elastic topology optimization designs 1501 and 1502, where the models were subjected to a load magnitude of P = 500 kN. Figure 15 further illustrates the 500 kN load results from columns 1403a and 1404b of Figure 14, where numerical stabilization was not applied (columns 1403a and 1501) and where numerical stabilization was applied by adding artificial constitutive modeling (columns 1404a and 1502). Figure 15 shows that the optimization workflow without numerical stabilization 1501 fails during the optimization iterations to reach a primary solution because the finite elements are numerically unstable and ill-conditioned (as shown in zoom 1503). Therefore, the overall optimization workflow without stabilization fails to converge. The optimization workflow with numerical stabilization by adding artificial constitutive modeling 1502 converges to an optimized design. Furthermore, the optimization iteration history is shown in plots 1504a-b of the added artificial constitutive modeling.

[0185] Additionally, FIG. 15 includes plots 1504a and 1504b of the optimization iteration history. Specifically, plot 1504a shows displacement 1508 compared to optimization iteration 1507a, and plot 1504b shows relative mass 1509 compared to displacement 1507 for optimization iteration 1507b. Plots 1504a-b demonstrate smooth convergence in objectives and also demonstrate the need for fewer optimization iterations, and therefore lower computational costs. FIG. 15 also shows the stress distributions for the physical constitutive material model 1505 and the added artificial constitutive modeling 1506. The stress scaling is different for the two stress distributions 1505 and 1506 shown. The solid portion of the optimized design primarily transmits physical stresses, while the stresses from the added artificial constitutive modeling are present only in locations and directions necessary to ensure well-conditioned global operators and numerical stability.

[0186] The optimization results illustrated in FIGS. 14 and 15 lead to the following conclusions regarding embodiments. First, embodiments support numerical stabilization of physical systems consisting of both elastic and elastoplastic constitutive material models. More generally, embodiments are suitable for all types of physical constitutive models. Second, embodiments add only stresses, and thereby stiffness, at locations and directions necessary for numerical stabilization. Embodiments thereby locally stabilize ill-conditioned portions of a model by adding additional artificial constitutive modeling to the physical constitutive model without changing the overall stiffness of the system. As noted above, many previous approaches add global stiffness to the model at all locations and directions, thereby adding stiffness to the model even at locations and directions where additional stiffness is not necessary. Third, when embodiments are applied to geometrically nonlinear topology optimization workflows, they are independent of design variables but strictly depend on the deformation of the principal solution of the structural layout during optimization iterations. Therefore, heuristic schemes and heuristic parameters depending on the design variables are not modified during optimization iterations. This requires significantly fewer optimization iterations. Because the embodiments require fewer iterations, they have lower computational costs and significantly faster and smoother optimization convergence compared to previous and benchmark references. Furthermore, the embodiments have shorter run times to determine an optimized design. The embodiments outperform all previous numerical stabilization methods in terms of the number of optimization iterations to achieve an optimized geometric nonlinear design. Furthermore, the embodiments prevent the optimization from failing to converge when some elements of the physical system are numerically unstable and ill-conditioned relative to the principal solution of structural equilibrium in a given optimization iteration.

[0187] 16A and 16B illustrate a geometrically nonlinear transient 3D topology optimization in which elastic-plastic material properties 1605 are defined for a finite element model 1600.

[0188] In Figures 16A and 16B, model 1600 includes a flexible beam 1601 positioned on two rigid supports 1602a-b, with the center of beam 1601 exposed to a moving rigid punch 1603. Supports 1602a-b are fully constrained, while punch 1603 is free to move vertically with an initial velocity v = 10 m / s and mass m = 200 kg. Beam 1601 is unconstrained. Contact between the four bodies (beam 1601, supports 1602a-b, and punch 1603) has a friction coefficient of μ = 0.1. Elastic-plastic material properties 1605 are defined for the material of beam 1601, while the two lower supports 1602a-b and punch 1603 are all assumed to be rigid. In one implementation, model 1600 is reduced by a factor of four by applying two planes of symmetry, as shown in Figure 16B. The quarter beam is discretized using 13,482 C3D8 rectilinear hexahedral elements. Figure 16A represents the full model 1600a, and Figure 16B represents the reduced model 1600b by applying symmetry.

[0189] In one implementation, the optimization workflow in 1200 is applied to model 1600 for topology optimization. The optimization objective function is to minimize the maximum vertical displacement over time of reference point RP 1604 relative to punch 1603 at 25 discrete time points t1∈{0.001, 0.002,...,0.025}. Additionally, a 30% relative mass constraint is enforced on the design domain, whereby the dynamic stiffness is maximized for a given mass. The design domain for topology optimization is defined as the entire flexible beam 1601, where each finite element is a design variable in the topology optimization to determine a new optimized conceptual material layout.

[0190] The results of the geometric nonlinear topology optimization for the dynamic stiffness optimization problem described herein above in connection with FIGS. 16 and 16B are shown in FIGS. 17A-17D and 18A-18D.

[0191] FIG. 17A shows a 3D model 1701 representing the density design layout of a failed optimization iteration (specifically, iteration 22) in an implementation where numerical stabilization was not applied. FIG. 17B is a plot 1702 representing the optimization iteration history, showing the normalized maximum vertical displacement objective function value 1706 and the normalized volume constraint value 1707 as a function of optimization iteration 1708. FIG. 17C shows 3D model 1703, the same 3D model from FIGS. 16A-16B, with the same flexible beam 1601, rigid punch 1603, and support 1602a. The shading in model 1703 indicates stresses at the final convergence time increment of the failed principal solution at optimization iteration 22, where a localized artificial instability phenomenon 1704 was exhibited. The phenomenon 1704 ultimately caused the failure of the principal solution, and thus the complete failure of optimization iteration 22, when numerical stabilization was not applied. FIG. 17D shows 3D model 1705, where the shading indicates stresses. Model 1705 is the same 3D model 1701 from Figures 16A-16B, with the same model 1701 of Figure 17A, the same flexible beam 1601, rigid punch 1603, and support 1602a. Furthermore, model 1705 is determined at the same time as model 1703 of Figure 17C, but model 1705 was determined by applying numerical stabilization using the addition of artificial constitutive modeling at 1202. By adding artificial constitutive modeling, no local artificial instabilities and local buckling are observed (as shown in Figure 1709), and the principal solution fully converges.

[0192] The results from FIG. 17C show that when numerical stabilization is not applied, the main solution does not converge at optimization iteration 22. A local instability and buckling phenomenon 1704 in one of the low-density regions is observed, causing the main analysis to fail. Consequently, the full optimization fails at optimization iteration 22. In contrast, as shown in FIG. 17D, when numerical stabilization according to an embodiment is applied, complete convergence of the main solution is achieved for the same model. Thus, using an embodiment, the analysis of the main solution fully converges, and no signs of local instability or buckling are observed during the analysis 1709.

[0193] FIG. 18A shows a 3D model 1801 with an optimized design density layout determined using method 1200. In this implementation, the layout of model 1801 was determined at optimization iteration number 91. FIG. 18B is a plot 1802 representing an optimization iteration history showing objective function values ​​1804 of normalized maximum vertical displacement and normalized volume constraint values ​​1805 as a function of optimization iteration 1806. FIG. 18C shows 3D model 1802, where the shading represents the stresses of the principal solution of the initial design at optimization iteration number 0, the time point with maximum displacement. Model 1802 is the same 3D model from FIGS. 16A-16B with the same flexible beam 1601, rigid punch 1603, and support 1602a. FIG. 18D shows 3D model 1803, where the shading represents the stresses of the principal solution of the optimized design at optimization iteration number 91, the time point with maximum displacement.

[0194] The results in Figures 18A-18D show that application of the embodiment for nonlinear topology optimization leads to a fully converged method for the solid-void topology optimization design, and the dominant solution does not fail during the optimization iteration process. Furthermore, when compared to plots 1702 and 1802 in Figures 17B and 18B, respectively, a smoother convergence behavior of the optimization iterations due to numerical stabilization is observed.

[0195] Embodiments can be used in modeling areas including geometric nonlinear modeling and large deformation continuum finite element modeling, where numerically ill-conditioned stiffness operators lead to numerically artificial instabilities.

[0196] Embodiments require fewer solver iterations to obtain a converged primal solution. Additionally, embodiments can ensure that a converged primal solution is obtained by suppressing numerical artifacts of global operators in the finite element model.

[0197] It should be noted that although embodiments are described herein with respect to structural finite element models, the embodiments are not limited to structural finite element models, but may also be applied to other domains, such as multiphysics finite element models solved using finite element methods, such as finite element models representing thermal, fluid, and / or electromagnetic designs, among other examples.

[0198] Embodiments can be used for geometrically nonlinear topology optimization during conceptual design. In such applications, void and intermediate density elements are often ill-conditioned with respect to stiffness operators, causing the primal solution in the topology optimization workflow to converge slowly or not at all, causing the entire optimization workflow to fail.

[0199] Embodiments can be utilized to force faster convergence of the main solution at each optimization iteration, independent of the design layout.Embodiments of the present technology can force convergence of the main solution to a solution at each optimization iteration, independent of the design layout.

[0200] Advantageously, embodiments can directly link the condition numbers of the element stiffness operators to the constitutive model.

[0201] FIG. 19 is a simplified block diagram of a computer-based system 1920 that can be used to determine the physical behavior of real-world objects according to any of the various embodiments described herein. The system 1920 includes a bus 1923. The bus 1923 serves as an interconnection between the various components of the system 1920. Connected to the bus 1923 is an input / output device interface 1926 for connecting various input and output devices, such as a keyboard, mouse, display, and speakers, to the system 1920. A central processing unit (CPU) 1922 is connected to the bus 1923 and provides for the execution of computer instructions implementing the embodiments. The memory 1925 provides volatile storage for data used to execute computer instructions implementing the embodiments described herein, such as methods 400, 500, 700, and 1200 described above in connection with FIGS. 4, 5, 7, and 12, respectively. The storage device 1924 provides non-volatile storage for software instructions, such as an operating system (not shown) and embodiment configurations. The system 1920 also includes a network interface 1921 for connecting to any of a variety of networks known in the art, including wide area networks (WANs) and local area networks (LANs).

[0202] It should be understood that the exemplary embodiments described herein may be implemented in many different ways. In some instances, the various methods and systems described herein may each be implemented by a physical, virtual, or hybrid general-purpose computer, such as computer system 1920 or a computer network environment, such as computer environment 2020 described herein below in connection with FIG. 20. Computer system 1920 may be converted into a machine that performs the methods described herein, for example, by loading software instructions into either memory 1925 or non-volatile storage 1924 for execution by CPU 1922. Those skilled in the art should further appreciate that system 1920 and its various components may be configured to implement any embodiment or combination of embodiments described herein. Furthermore, system 1920 may implement the various embodiments described herein utilizing any combination of hardware, software, and firmware modules operably coupled internally or externally to system 1920. Furthermore, system 1920 may be communicatively coupled to or embedded within manufacturing equipment to control the equipment to create physical objects, as described herein.

[0203] 20 illustrates a computer network environment 2020 in which embodiments of the present invention may be implemented. In the computer network environment 2020, a server 2021 is linked to clients 2023a-n via a communications network 2022. The environment 2020 may be used to enable the clients 2023a-n, alone or in combination with the server 2021, to perform any of the methods described herein. As non-limiting examples, the computer network environment 2020 may provide cloud computing embodiments, software as a service (SAAS) embodiments, etc.

[0204] The embodiments or aspects thereof may be implemented in the form of hardware, firmware, or software. If implemented in software, the software may be stored on any non-transitory computer-readable medium configured to enable a processor to load the software, or a subset of its instructions. The processor is then configured to execute the instructions to operate a device or cause it to operate in a method described herein.

[0205] Furthermore, firmware, software, routines, or instructions may be described herein as performing certain operations and / or functions of a data processor, although it will be understood that such descriptions contained herein are merely for convenience and that such operations actually result from a computing device, processor, controller, or other device executing the firmware, software, routines, instructions, etc.

[0206] It will be understood that the flow diagrams, block diagrams, and network diagrams may include more or fewer elements, may be arranged differently, or may be represented differently, but it will also be understood that a particular implementation may implement the block diagrams and network diagrams, and the number of block diagrams and network diagrams illustrating the implementation of an embodiment, in a particular way.

[0207] Accordingly, further embodiments may also be implemented in various computer architectures, physical computers, virtual computers, cloud computers, and / or some combination thereof, and therefore the data processors described herein are intended to be illustrative only and not limiting of the embodiments.

[0208] While exemplary embodiments have been particularly shown and described, those skilled in the art will understand that various changes in form and details can be made therein without departing from the scope of the embodiments encompassed by the appended claims.

[0209] The teachings of all patents, published applications, and references cited herein are incorporated by reference in their entirety.

[0210] References [1] Schillinger, D., Duster, A., Rank, E. (2012). The hp-d-adaptive finite cell method for geometrically nonlinear problems of solid mechanics. International Journal for Numerical Methods in Engineering. 89:1171-1202.

[0211] [2] Sigmund, O., Maute, K. (2013). Topology optimization approaches. Structural and Multidisciplinary Optimization, 48(6):1031-1055.

[0212] [3] Sigmund, O. (2022). On benchmarking and good scientific practice in topology optimization. Structural and Multidisciplinary Optimization, 65, 315

[0213] [4] Bruns, TE, Tortorelli, DA (2003). An Element Removal and Reintroduction Strategy for the Topology Optimization of Structures and Compliant Mechanisms. International Journal for Numerical Methods in Engineering.57:1413-1430.

[0214] [5] Bruns, T.E., Sigmund, O. Tortorelli, D.A. (2002). Numerical Methods for the Topology Optimization of Structures that Exhibit Snap Through. International Journal for Numerical Methods in Engineering. 55:1215-1237.

[0215] [6] Bruns, T.E. (2006). Zero density lower bounds in topology optimization. Computer Methods in Applied Mechanics & Engineering.196: 566-578.

[0216] [7] Behrou, R., Lotfi, R., Carstensen, J., Ferrari, F., Guest J. (2021). Revisiting element removal for density-based structural topology optimization with reintroduction by Heaviside projection. Computer Methods in Applied Mechanics & Engineering. 380:113799.

[0217] [8] Wang, F., Lazarov, B. S., Sigmund, O., and Jensen, J. S. (2014). Interpolation scheme for fictitious domain techniques and topology optimization of finite strain elastic problems. Computer Methods in Applied Mechanics & Engineering. 276(7):453-472.

[0218] [9] Dalklint, A., Wallin, M., Tortorelli. D.A. (2021). Structural stability and artificial buckling modes in topology optimization. Structural and Multidisciplinary Optimization. 64(4):1751-1763.

[0219]

[10] Bluhm, G.L., Sigmund, O., Poulios K. (2021)。Internal contact modeling for finite strain topology optimization. Computational Mechanics. 67:1099-1114.

[0220]

[11] Wallin, M., Ivarsson, N., and Tortorelli, D.A. (2018). Stiffness optimization of non-linear elastic structures. Computer Methods in Applied Mechanics & Engineering. 330:292-307.

[0221]

[12] Buhl, T., Pedersen, C.B.W., Sigmund, O. (2000). Stiffness design of geometrically nonlinear structures using topology optimization. Structural and Multidisciplinary Optimization. 19:93-104.

[0222]

[13] Pedersen, C.B.W., Buhl, T., Sigmund, O. (2001). Topology Synthesis of Large displacement compliant mechanisms. International Journal for Numerical Methods in Engineering. 50:2683-2705.

[0223]

[14] Pedersen, C.B.W., Fleck, N.A., Ananthasuresh, G.K. (2006). Design of a compliant mechanism to modify an actuator characteristic to deliver a constant output force. ASME Journal of Mechanical Design, 128:1101-1112.

[0224]

[15] Maute, K., Schwarz, S., Ramm, E. (1998). Adaptive topology optimization of elastoplastic structures. Structural optimization. 15:81-91.

[0225]

[16] Wallin, M., Jonsson, V., Wingren, E. (2016). Topology optimization based on finite strain plasticity. Structural and Multidisciplinary Optimization. 54:783-793.

[0226]

[17] Pedersen, C.B.W. (2002). Revisiting Topology Optimization of Continuum Structures with Elastoplastic Response. In Proceeding 15th Nordic Seminar on Computational Mechanics, Aalborg, Denmark.

[0227]

[18] Lahuerta, R.D., Simoes, E.T., Campello, E.M.B., Pimenta, P.M., Silva E.C.N. (2013). Towards the stabilization of the low density elements in topology optimization with large deformation. Computational Mechanics. 52(4): 779-797.

[0228]

[19] Ortigosa, R., Ruiz, D., Gil A.J., Donoso, A., Bellido, J.C. (2020) A stabilization approach for topology optimization of hyperelastic structures with the SIMP method. Computer Methods in Applied Mechanics and Engineering. 364:112924.

[0229]

[20] Ortigosa, R., Martinez-Frutos, J., Gil A.J., Herrero-Perez, D. (2019). A new stabilization approach for level-set based topology optimization of hyperelastic materials. Structural and Multidisciplinary Optimization. 60:2343-2371.

[0230]

[21] Abaqus. (2022). SIMULIA User Assistance. Dassault Systemes.

[0231]

[22] Riks, E. (2008). On the Purpose and Limitations of Buckling Analysis. 2nd International Conference on Buckling and Post-buckling Behavior of Composite Laminated Shell Structures, Braunschweig, Germany.

[0232]

[23] Zienkiewicz, O.C., Taylor, R.L. and Fox, D. D. The Finite Element Method for Solid and Structural Mechanics, Seventh Edition, Elsevier, (2014)

[0233]

[24] Bathe, K.J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, N. J., (1982).

[0234]

[25] Bergstrom, J. (2015). In Mechanics of Solid Polymers: Theory and Computational Modeling. Elsevier.

[0235]

[26] Tosca. (2022). SIMULIA User Assistance. Dassault Systemes.

[0236]

[27] Kleiber, M., Antunez, H., Hien, T. D. and Kowalczyk, P. (1997). Parameter Sensitivity in Nonlinear Mechanics, Theory and Finite Element Computations, John Wiley and Sons.

[0237]

[28] Choi, K. K. and Kim, N. H (2005), Structural Sensitivity Analysis and Optimization 2: Nonlinear Systems and Applications、Springer-Verlag New York Inc.

[0238]

[29] Michaleris, P., Tortorelli, D. A. and Vidal, C. A. (1994). Tangent Operators and Design Sensitivity Formulations for Transient Non-linear Coupled Problems with Applications to Elastoplasticity, International Journal for Numerical Methods in Engineering 37: 2471-2499.

Claims

1. 1. A computer-implemented method for determining physical behavior of a real-world object, comprising: defining, in a memory of the processor, a computer-based model representing the real-world object using a plurality of elements, the defined model indicating one or more materials represented by each element of the plurality of elements; defining equations that describe the physics-based behavior of the one or more materials; defining a stabilization equation that is a function of a nonlinear deformation gradient matrix; for each element of the plurality of elements in the model, associating artificial internal forces with the element based on the defined stabilization equations, and associating artificial-based and physics-based behaviors with the element based on the associated artificial internal forces; performing a simulation of the real-world object subjected to a load using (i) the defined computer-based model, (ii) the defined equations describing physics-based behavior, and (iii) the defined stabilization equations, wherein performing the simulation includes applying the stabilization equations to each of the plurality of elements, and wherein results of performing the simulation are indicative of the physical behavior of the real-world object.

2. For each element of the plurality of elements in the model: The computer-implemented method of claim 1 , further comprising associating internal physical forces with the elements based on the defined equations that describe physics-based behavior.

3. For each element of the plurality of elements, The computer-implemented method of claim 2 , further comprising assembling a remainder using the associated artificial internal forces and the associated physical internal forces, the assembled remainder defining a structural equilibrium.

4. 2. The computer-implemented method of claim 1, wherein performing the simulation of the real-world object subjected to the load further comprises using the artificial internal forces associated with each element and the artificial-based and physics-based behaviors associated with each element.

5. 2. The computer-implemented method of claim 1, wherein performing the simulation of the real-world object using the defined stabilization equation reduces a condition value associated with a given element of the plurality of elements, the condition value being a function of the nonlinear deformation gradient matrix.

6. updating the computer-based model based on the results of performing the simulation; and 2. The computer-implemented method of claim 1, further comprising: determining updated physical behavior of the real-world object by performing a simulation of the real-world object using (i) the updated computer-based model, (ii) the defined equations describing physics-based behavior, and (iii) the defined stabilization equations.

7. 7. The computer-implemented method of claim 6, further comprising iterating (i) updating the computer-based model based on the determined updated physical behavior and (ii) determining the updated physical behavior until the determined updated physical behavior satisfies a criterion.

8. The computer-implemented method of claim 1 , wherein the computer-based model is any one or combination of a finite element model, a boundary element method, a finite difference method, a finite volume method, or a discrete element method.

9. 10. The computer-implemented method of claim 1, wherein the real-world object represented by the computer-based model is one of an automobile, industrial equipment, an aircraft, a civil structure, a marine device, a medical device, a consumer product, an electronic device, an armored vehicle, or a manufacturing facility.

10. 1. A system for determining physical behavior of a real-world object, comprising: a processor; a memory having computer code instructions stored thereon, said processor and said memory using said computer code instructions to cause said system to: defining in the memory a computer-based model representing the real-world object using a plurality of elements, the defined model indicating one or more materials represented by each element of the plurality of elements; defining equations that describe the physics-based behavior of the one or more materials; defining a stabilization equation that is a function of a nonlinear deformation gradient matrix; for each element of the plurality of elements in the model, associating artificial internal forces with the element based on the defined stabilization equations, and associating artificial-based and physics-based behaviors with the element based on the associated artificial internal forces; 1. A system configured to: perform a simulation of the real-world object subjected to a load using (i) the defined computer-based model, (ii) the defined equations describing physics-based behavior, and (iii) the defined stabilization equations, wherein performing the simulation includes applying the stabilization equations to each of the plurality of elements, and wherein results of performing the simulation are indicative of the physical behavior of the real-world object.

11. The processor and the memory use the computer code instructions to, for each element of the plurality of elements in the model, The system of claim 10 , further configured to associate internal physical forces with the elements based on the defined equations that define physics-based behavior.

12. The processor and the memory use the computer code instructions to instruct the system, for each element of the plurality of elements, The system of claim 11 , further configured to assemble a residue using the associated artificial internal forces and the associated physical internal forces, the assembled residue defining a structural equilibrium.

13. 11. The system of claim 10, wherein the processor and the memory are configured to use the computer code instructions to cause the system to use the artificial internal forces associated with each element and the artificial-based and physics-based behaviors associated with each element when performing the simulation.

14. 11. The system of claim 10, wherein performing the simulation of the real-world object using the defined stabilization equation reduces a condition value associated with a given element of the plurality of elements, the condition value being a function of the nonlinear deformation gradient matrix.

15. The processor and the memory use the computer code instructions to cause the system to: updating the computer-based model based on the results of performing the simulation; and 11. The system of claim 10, further configured to: determine updated physical behavior of the real-world object by performing a simulation of the real-world object using (i) the updated computer-based model, (ii) the defined equations describing physics-based behavior, and (iii) the defined stabilization equations.

16. The processor and the memory use the computer code instructions to cause the system to:

16. The system of claim 15, further configured to iterate (i) updating the computer-based model based on the determined updated physical behavior, and (ii) determining the updated physical behavior, until the determined updated physical behavior satisfies a criterion.

17. The system of claim 10 , wherein the computer-based model is any one or combination of a finite element model, a boundary element method, a finite difference method, a finite volume method, or a discrete element method.

18. 11. The system of claim 10, wherein the real-world object represented by the computer-based model is one of an automobile, industrial equipment, an aircraft, a civil structure, a marine device, a medical device, a consumer product, an electronic device, an armored vehicle, or a manufacturing facility.

19. 1. A computer program for determining physical behavior of real world objects, said computer program executed by a server in communication with one or more clients over a network, said computer program comprising: and a computer-readable medium having program instructions that, when executed by a processor, cause the processor to: defining, in a memory of the processor, a computer-based model representing the real-world object using a plurality of elements, the defined model indicating one or more materials represented by each element of the plurality of elements; defining equations that describe the physics-based behavior of the one or more materials; defining a stabilization equation that is a function of a nonlinear deformation gradient matrix; for each element of the plurality of elements in the model, associating artificial internal forces with the element based on the defined stabilization equations, and associating artificial-based and physics-based behaviors with the element based on the associated artificial internal forces; 1. A computer program comprising computer instructions for performing a simulation of the real-world object subjected to a load using (i) the defined computer-based model, (ii) the defined equations describing physics-based behavior, and (iii) the defined stabilization equations, wherein performing the simulation includes applying the stabilization equations to each of the plurality of elements, and wherein results of performing the simulation are indicative of the physical behavior of the real-world object.

20. 20. The computer program product of claim 19, wherein performing the simulation of the real-world object using the defined stabilization equation reduces a condition value associated with a given element of the plurality of elements, the condition value being a function of the nonlinear deformation gradient matrix.

Images