Mechanistic concurrent nested topology design and optimization: bridging advanced manufacturing to optimized product realization

Abstract

A mechanistic nested topological design optimization method for using Convolution-Hierarchical Deep-learning Neural Network Tensor Decomposition (C-HiDeNN-TD) comprises (1) obtaining at least one design variable and at least one design requirement; (2) decomposing the at least one design variable into at least one equation in each of at least two directions; (3) solving the at least one design variable in each of the at least two directions using a C-HiDeNN-TD solver to obtain at least one part design; and (4) outputting an optimal design based on the at least one part design when the optimal design meets the least one design requirement.

Description

Attorney Docket No.: 0116936.295WO2 MECHANISTIC CONCURRENT NESTED TOPOLOGY DESIGN AND OPTIMIZATION: BRIDGING ADVANCED MANUFACTURING TO OPTIMIZED PRODUCT REALIZATION CROSS-REFERENCE TO RELATED PATENT APPLICATION

[0001] This application claims priority to and the benefit of U.S. Provisional Patent Application Serial No.63 / 534,599, filed August 25, 2023, which is incorporated herein in its entirety by reference. FIELD OF INVENTION

[0002] The present disclosure relates to the technical field of Convolution-Hierarchical Deep- learning Neural Network -Tensor Decomposition (C-HiDeNN-TD) which pertains to a method, algorithm, and computer program employed in engineering analysis, and more particularly, to the use of C-HiDeNN-TD to solve a structure, nested lattice, and nested material system. BACKGROUND OF THE INVENTION

[0003] The background description provided herein is for the purpose of generally presenting the context of the present invention. The subject matter discussed in the background of the invention section should not be assumed to be prior art merely as a result of its mention in the background of the invention section.

[0004] Early works in optimal design of multi-scale structures used homogenization-based approaches with small square cells and evaluated equivalent mechanical properties through numerical homogenization. In the late 1990s, focus shifted to "mono-scale" approaches optimizing the distribution of a homogeneous material. Additive manufacturing (AM) has sparked renewed interest in multi-scale structure design, allowing for complex mono-scale and delicate multi-scale structures. Topology optimization (TO) and lattice infill are dominant strategies for lightweight structure design. Combining topology optimization with multi-scale structures holds promise for superior performance. Natural multi-scale structures like bone and bamboo inspire these efforts, despite not fully understanding their specific functions and underlying mechanisms.

[0005] However, current topology optimization tools are limited to a few scales, and the potential of nested topology optimization in engineering design has yet to be realized.

[0006] For real-time control and / or optimization of the design and manufacturing processes, one needs both a fast and differentiable tool. Currently, many works rely on reduced-order method toAttorney Docket No.: 0116936.295WO2 avoid computing a more accurate but computationally intensive solution of the governing equations. However, using of Convolution Hierarchical Deep-learning Neural Network Tensor Decomposition(C-HiDeNN-TD), which enables the higher-order accuracy with arbitrary smoothness of the design features, in topology design and optimization is still absent.

[0007] Therefore, it is desirable to use C-HiDeNN Tensor Decomposition(C-HiDeNN-TD) for mechanistic concurrent nested topology design and optimization for advanced materials system. SUMMARY OF INVENTION In light of the foregoing, this invention discloses a mechanistic nested topological design optimization method for using Convolution-Hierarchical Deep-learning Neural Network Tensor Decomposition (C-HiDeNN-TD). The method comprises (1) obtaining at least one design variable and at least one design requirement; (2) decomposing the at least one design variable into at least one equation in each of at least two directions; (3) solving the at least one design variable in each of the at least two directions using a C-HiDeNN-TD solver to obtain at least one part design; and (4) outputting an optimal design based on the at least one part design when the optimal design meets the least one design requirement.

[0008] In one embodiment, between step (3) and step (4) the method further comprises a step of updating the at least one design variable based on the at least one part design.

[0009] In one embodiment, when the optimal design fails to meet the least one design requirement, repeat the steps (2) to (3) until the optimal design meets the least one design requirement.

[0010] In one embodiment, the at least one equation comprises at least one equilibrium equation.

[0011] In one embodiment, the at least one equation further comprises at least one scale coupling equation.

[0012] In one embodiment, step (3) comprises a lattice level design process and a material level design process.

[0013] In one embodiment, wherein the lattice level design process further comprises (i) setting up an initial value for the at least one design variable; (ii) providing a mode q; (iii) solving the at least one design variable in each of the three directions for the mode q; (iv) obtaining the at least one part design when the at least one design variable meet a convergence requirement.

[0014] In one embodiment, the method further comprises providing a new mode q and repeatingAttorney Docket No.: 0116936.295WO2 the steps (iii) and (iv) when the at least one design variable meet a convergence requirement.

[0015] In one embodiment, the material level design process comprises selecting a material, and obtaining a homogenized material property based on the selected material.

[0016] In one embodiment, the step (3) further comprises a step of analyzing a sensitivity deviation of the at least one part design.

[0017] In another aspect of the invention, a method for graphics processing units (GPU) programming of a mechanistic nested topological design optimization using Convolution- Hierarchical Deep-learning Neural Network Tensor Decomposition (C-HiDeNN-TD) comprises (1) obtaining at least one design variable and at least one design requirement; (2) decomposing the at least one design variable into at least one equation in each of at least two directions; (3) solving the at least one design variable in each of the at least two directions using a C-HiDeNN-TD solver to obtain at least one part design; and (4) outputting an optimal design based on the at least one part design when the optimal design meets the least one design requirement.

[0018] In one embodiment, between step (3) and step (4) the method further comprises a step of updating the at least one design variable based on the at least one part design.

[0019] In one embodiment, when the optimal design fails to meet the least one design requirement, repeat the steps (2) to (3) until the optimal design meets the least one design requirement.

[0020] In one embodiment, the at least one equation comprises at least one equilibrium equation.

[0021] In one embodiment, the at least one equation further comprises at least one scale coupling equation.

[0022] In one embodiment, step (3) comprises a lattice level design process and a material level design process.

[0023] In one embodiment, wherein the lattice level design process further comprises (i) setting up an initial value for the at least one design variable; (ii) providing a mode q; (iii) solving the at least one design variable in each of the three directions for the mode q; (iv) obtaining the at least one part design when the at least one design variable meet a convergence requirement.

[0024] In one embodiment, the method further comprises providing a new mode q and repeating the steps (iii) and (iv) when the at least one design variable meet a convergence requirement.

[0025] In one embodiment, the material level design process comprises selecting a material, and obtaining a homogenized material property based on the selected material.Attorney Docket No.: 0116936.295WO2

[0026] In one embodiment, the step (3) further comprises a step of analyzing a sensitivity deviation of the at least one part design.

[0027] These and other aspects of the present invention will become apparent from the following description of the preferred embodiment taken in conjunction with the following drawings, although variations and modifications therein may be affected without departing from the spirit and scope of the novel concepts of the disclosure. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] The accompanying drawings illustrate one or more embodiments of the invention and together with the written description, serve to explain the principles of the invention. Wherever possible, the same reference numbers are used throughout the drawings to refer to the same or like elements of an embodiment.

[0029] Fig.1 illustrates structure of the HiDeNN for the displacement in a 1D problem with linear elements, in which the whole neural network consists of small neural network blocks shown inside the dashed line box.

[0030] Fig.2 illustrates Illustration of the patch domain and final convolution patch for linear element e. The path size is defined as s = 2; panel (a) shows the patch domain associated with node i – 1; panel (b) shows the patch domain associated with node I; panel (c) shows the final convolution patch associated with element e with the relation.

[0031] Fig.3 illustrates H1norm error vs. computation time plot for a 1D Poisson problem. FEM uses dashed lines, while C-HiDeNN uses solid lines. Curves with the same polynomial order p have the same color. Both FEM and C-HiDeNN are solved via a Jacobi-preconditioned conjugate gradient iterative solver implemented on a small GPU (NVIDIA GeForce GTX 1650, 4GB memory).

[0032] Fig.4 illustrates a diagram of flow chart of C-HiDeNN-TD for topology optimization.

[0033] Fig.5 illustrates differences between the panel (a) SIMP-TO approach, in which an extra filter is needed to control the length-scale and avoid a checkerboard pattern; and panel (b) C- HiDeNN-TD-TO approach, in which no extra filter is needed.

[0034] Fig.6 illustrates a density field with a checkboard pattern; panel (a) shows checkerboard pattern; panel (b) shows a density filter is applied to the density shown in panel (a) as traditional TO approaches. The checkerboard pattern is reduced. Panel (c) shows the C-HiDeNN built-in filter is applied to (a).

[0035] Fig.7 illustrates an example reflecting the different between SIMP approach and C-Attorney Docket No.: 0116936.295WO2 HiDeNN-TD-TO approach; panel (a) shows topology optimization SIMP density plot; panel (b) shows C-HiDeNN-TD-TO density plot with 27 interpolation points in each element.

[0036] Fig.8 illustrates three parameters in the convulation shape functions of one example; panel (a) shows cantilever beam example with a uniform loading on the right side and a clamped boundary condition on the left side. The size is set to M = N = 50 and L = 100; panel (b) shows illustration of the patch size s, and dilation parameter a.

[0037] Fig.9 illustrates one example of convolution polynomial order study. The patch size is set to s = 4. Example (a) uses a small dilation parameter to degenerate into the TD-TO method. From the comparison, polynomials of higher order produce a better objective value with more support structures because of the smoothness of the corresponding solutions.

[0038] Fig.10 illustrates another example of convolution polynomial order study. The patch size is set to s = 4. Example (a) uses a small dilation parameter to degenerate into the TD-TO method. From the comparison, a larger dilation parameter produces a better design and contributes to the design length-scale control.

[0039] Fig.11 illustrates a chart reflecting DoFs vs computational cost between FEM-TO and C-HiDeNN-TD-TO.

[0040] Fig.12 illustrates three cantilever beam cases with different mesh sizes, in which higher resolution TO produces more optimal designs with more supporting structural features.

[0041] Fig.13 illustrates an example of application of C-HiDeNN-TD-TO to a drone design problem.

[0042] Fig.14 illustrates drone design results with volume constraint 0.2 with different design resolutions. Higher smoothness and more detailed structures can be seen for high-resolution design.

[0043] Fig.15 illustrates drone design results with volume constraint 0.12 with different design resolutions.

[0044] Fig.16 illustrates a flowchart of C-HiDeHH-TD for lattice design.

[0045] Fig.17 illustrates lattice design for: panel (a) Maximum bulk modulus; panel (b) Maximum shear modulus; panel (c) Negative Poisson’s ratio design. The resolution of showed examples are 400 × 400 × 400.

[0046] Fig.18 illustrates a flowchart for C-HiDeNN-TD2 two-scale concurrent design

[0047] Fig.19 illustrates diagrams of multi unit cell concurrent drone design problem definition. The initial design domain is divided into different regions as shown in (b) front view and (c) sideAttorney Docket No.: 0116936.295WO2 view.

[0048] Fig.20 illustrates multi unit cell concurrent drone designed structure.

[0049] Fig.21 illustrates offline and online stages of two-stage clustering analysis methods; the left panel shows the offline stage contains three steps as shown in the figure, which will generate a compressed RVE database. The compressed RVE can be then used to predict mechanical responses of the RVE.

[0050] Fig.22 illustrates geometry for an example problem: a two-phase periodic composite represented in 2D with plane strain where blue is a continuous matrix phase and yellow represents inclusions.

[0051] Fig.23 illustrates color contours showing clusters distribution in panel (a) the inclusions and anel (b) the matrix.

[0052] Fig.24 illustrates 3D charts of component-wise magnitude plots for (a) ^^^^ ^^ ^^ ^^^^ ^^ ^^; (b) ^^^^ ^^ ^^; (c) ^^^^ ^^^^ ^^ ^^; spikes along the diagonal direction for all three interaction tensor surface plots suggest self interaction has more contribution than the rest of clusters in cluster-wise stress increment. ^^^^ ^^^^ ^^ ^^and ^^^^ ^^^^ ^^ ^^have the similar magnitude along their diagonal direction due to the homogeneous reference material assumption. ^^^^ ^^^^ ^^ ^^has different magnitudes for matrix and inclusion phase along the diagonal direction, implicitly representing a heterogeneous reference material.

[0053] Fig.25 illustrates 3D charts of component-wise magnitude plots for (a) ^^^^ ^^ ^^ ^^^^ ^^ ^^; (b) ^^^^ ^^ ^^; (c) ^^^^ ^^^^ ^^ ^^.

[0054] Fig.26 illustrates plots of σ^^^^^^vs. ε^^^^^^.

[0055] Fig.27 illustrates plots of ^^^^^^^^vs. ^^^^^^^^.

[0056] Fig.28 illustrates a 3D plot of random samples of strain state; two hundred final states were selected, and four evenly spaced intermediate steps to reach the final states were recorded for a total of NT = 1, 000 samples.

[0057] Fig.29 illustrates histogram of the difference between the FFNN σMpredictions and validation data set for all NV= 150 samples. Most of the predictions made by FFNN has an l2norm less than 2 × 10−5, showing that the FFNN produce an accurate prediction of the stress state.

[0058] Fig.30 illustrates plots reflecting cross-correlation between SCA and FFNN for x-, y-, and shear-directions. The solid black lines are the ground truth: all perfect predictions should lay on those lines. The FFNN-predicted stress compo¬nents lay around the solid black lines, and all three casesAttorney Docket No.: 0116936.295WO2 have correlation coefficients of 1, exhibiting very strong predictive confidence.

[0059] Fig.31 illustrates plots reflecting stress-strain plot with both the FFNN and SCA solutions for the validation data, showing that the FFNN results successfully reproduced SCA results.

[0060] Fig.32 illustrates diagram of two dimensional convolutional neural network. In the training part, the number of threshold η∗= 5. After 4 repeats of convolution, the data will be passed to FFNN. The dimension of input for FFNN becomes Nfl= 74. The hidden layer in FFNN has 74 neurons. Finally, FFNN gives three macro strain components.

[0061] Fig.33 illustrates histogram of the l2norm for the CNN prediction using validation data points. Most of the predictions made by CNN has an l2norm less than 1 × 10−5, showing the CNN produce an accurate prediction of the strain state.

[0062] Fig.34 illustrates plots reflecting loading prediction by CNN vs. loading applied in SCA using validation data points. The solid black lines are the ground truth: all perfect predictions should lay on those lines. All three cases have correlation coefficients higher than 0.99, suggesting the trained CNN can provide a good accuracy in predicting applied strain.

[0063] Fig.35 illustrates topology optimization setup with FFNN. The FFNN will be used to compute non-linear material responses to drive for a new design. This replaces the constitutive law commonly used for the macroscale with a homogenized response of the microstructure for each point in the macroscale. Mathematically, σM(XM) = ρ(XM)FFFNN(ε(XM),∀XM∈ ΩM

[0064] Fig.36 illustrates amaterial responses in panel (a) vs optimized beam structure with non-linear material responses in panel (b). Comparing a) and b), the joint of all truss members for the non-linear case is located towards the bottom side of the truss structure.

[0065] Fig.37 illustrates topology optimization setup with FFNN and CNN. For each material point within the design zone, the FFNN is used to compute the material response, be it linear or non- linear, considering the effect of mi-crostructure.

[0066] Fig.38 illustrates a coparison of optimized beam structure without damage constraints in panel (a) v. optimized beam structure with damage constraints in panel (b). In (b) there are more trusses to avoid high local stresses that result in damage

[0067] Fig.39 illustrates the general framework for N*M nested lattice design.

[0068] Fig.40 illustrates flow charts of design framework of nested lattice design system for advanced materials.Attorney Docket No.: 0116936.295WO2

[0069] Fig.41 illustrates nested structure associated with bamboo design problem.

[0070] Fig.42 illustrates diagrams of bamboo design for case 1 with a volume constraint of 70%.

[0071] Fig.43 illustrates diagrams of bamboo design for case 2 with a volume constraint of 50%.

[0072] Fig.44 illustrates bamboo design for case 3 with a volume constraint shown in table 4. A different volume fraction constraint at different scale yields to a nested lattice structure with different lattice design.

[0073] Fig.45 illustrates a chart of numerical study of the number of nested concurrent scales. Compared with the 2-scale concurrent design, the axial compliance of 7-scale concurrent design is reduced from 1097.47 to 140.95; the bending compliance is reduced from 2523.898 to 237.8368; the torsion compliance is reduced from 3122.207 to 279.773.

[0074] Fig.46 illustrates design regions of car bumper: (a) initial design view; (b) 3D view; (c) front view.

[0075] Fig.47 illustrates design regions of car bumper in region 1 and region 2 of Fig.24.

[0076] Fig.48 illustrates design regions of car bumper in region 3 of Fig.24. DETAILED DESCRIPTION OF THE INVENTION

[0077] The invention will now be described more fully hereinafter with reference to the accompanying drawings, in which exemplary embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this invention will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art. Like reference numerals refer to like elements throughout.

[0078] The terms used in this specification generally have their ordinary meanings in the art, within the context of the invention, and in the specific context where each term is used. Certain terms that are used to describe the invention are discussed below, or elsewhere in the specification, to provide additional guidance to the practitioner regarding the description of the invention. For convenience, certain terms may be highlighted, for example using italics and / or quotation marks. The use of highlighting has no influence on the scope and meaning of a term; the scope and meaning of a term is the same, in the same context, whether or not it is highlighted. It will be appreciated that same thing can be said in more than one way. Consequently, alternative language and synonyms may be used for any one or more of the terms discussed herein, nor is any special significance to be placed upon whether or not a term is elaborated or discussed herein. Synonyms for certain terms areAttorney Docket No.: 0116936.295WO2 provided. A recital of one or more synonyms does not exclude the use of other synonyms. The use of examples anywhere in this specification including examples of any terms discussed herein is illustrative only, and in no way limits the scope and meaning of the invention or of any exemplified term. Likewise, the invention is not limited to various embodiments given in this specification.

[0079] One of ordinary skill in the art will appreciate that starting materials, biological materials, reagents, synthetic methods, purification methods, analytical methods, assay methods, and biological methods other than those specifically exemplified can be employed in the practice of the invention without resort to undue experimentation. All art-known functional equivalents, of any such materials and methods are intended to be included in this invention. The terms and expressions which have been employed are used as terms of description and not of limitation, and there is no intention that in the use of such terms and expressions of excluding any equivalents of the features shown and described or portions thereof, but it is recognized that various modifications are possible within the scope of the invention claimed. Thus, it should be understood that although the invention has been specifically disclosed by preferred embodiments and optional features, modification and variation of the concepts herein disclosed may be resorted to by those skilled in the art, and that such modifications and variations are considered to be within the scope of this invention as defined by the appended claims.

[0080] Whenever a range is given in the specification, for example, a temperature range, a time range, or a composition or concentration range, all intermediate ranges and subranges, as well as all individual values included in the ranges given are intended to be included in the invention. It will be understood that any subranges or individual values in a range or subrange that are included in the description herein can be excluded from the claims herein.

[0081] It will be understood that, as used in the description herein and throughout the claims that follow, the meaning of “a”, “an”, and “the” includes plural reference unless the context clearly dictates otherwise. Thus, for example, reference to “a cell” includes a plurality of such cells and equivalents thereof known to those skilled in the art. As well, the terms “a” (or “an”), “one or more” and “at least one” can be used interchangeably herein. It is also to be noted that the terms “comprising”, “including”, and “having” can be used interchangeably.

[0082] It will be understood that when an element is referred to as being “on”, “attached” to, “connected” to, “coupled” with, “contacting”, etc., another element, it can be directly on, attached to, connected to, coupled with or contacting the other element or intervening elements may also beAttorney Docket No.: 0116936.295WO2 present. In contrast, when an element is referred to as being, for example, “directly on”, “directly attached” to, “directly connected” to, “directly coupled” with or “directly contacting” another element, there are no intervening elements present. It will also be appreciated by those of skill in the art that references to a structure or feature that is disposed “adjacent” another feature may have portions that overlap or underlie the adjacent feature.

[0083] It will be understood that, although the terms first, second, third etc. may be used herein to describe various elements, components, regions, layers and / or sections, these elements, components, regions, layers and / or sections should not be limited by these terms. These terms are only used to distinguish one element, component, region, layer or section from another element, component, region, layer or section. Thus, a first element, component, region, layer or section discussed below could be termed a second element, component, region, layer or section without departing from the teachings of the invention.

[0084] Furthermore, relative terms, such as “lower” or “bottom” and “upper” or “top,” may be used herein to describe one element’s relationship to another element as illustrated in the Fig. s. It will be understood that relative terms are intended to encompass different orientations of the device in addition to the orientation depicted in the Fig. s. For example, if the device in one of the Fig. s is turned over, elements described as being on the “lower” side of other elements would then be oriented on “upper” sides of the other elements. The exemplary term “lower”, can therefore, encompasses both an orientation of “lower” and “upper,” depending of the particular orientation of the Fig.. Similarly, if the device in one of the Fig. s is turned over, elements described as “below” or “beneath” other elements would then be oriented “above” the other elements. The exemplary terms “below” or “beneath” can, therefore, encompass both an orientation of above and below.

[0085] It will be further understood that the terms “comprises” and / or “comprising”, or “includes” and / or “including”, or “has” and / or “having”, or “carry” and / or “carrying”, or “contain” and / or “containing”, or “involve” and / or “involving”, “characterized by”, and the like are to be open- ended, i.e., to mean including but not limited to. When used in this disclosure, they specify the presence of stated features, regions, integers, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, regions, integers, steps, operations, elements, components, and / or groups thereof.

[0086] Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which thisAttorney Docket No.: 0116936.295WO2 invention belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the relevant art and the invention, and will not be interpreted in an idealized or overly formal sense unless expressly so defined herein.

[0087] As used in the disclosure, “around”, “about”, “approximately” or “substantially” shall generally mean within 20 percent, preferably within 10 percent, and more preferably within 5 percent of a given value or range. Numerical quantities given herein are approximate, meaning that the term “around”, “about”, “approximately” or “substantially” can be inferred if not expressly stated.

[0088] As used in the disclosure, the phrase “at least one of A, B, and C” should be construed to mean a logical (A or B or C), using a non-exclusive logical OR. As used herein, the term “and / or” includes any and all combinations of one or more of the associated listed items.

[0089] The description below is merely illustrative in nature and is in no way intended to limit the invention, its application, or uses. The broad teachings of the invention can be implemented in a variety of forms. Therefore, while this invention includes particular examples, the true scope of the invention should not be so limited since other modifications will become apparent upon a study of the drawings, the specification, and the following claims. For purposes of clarity, the same reference numbers will be used in the drawings to identify similar elements. It should be understood that one or more steps within a method may be executed in different order (or concurrently) without altering the principles of the invention.

[0090] The present invention discloses a new mechanistic reduced-order multiscale method for design and optimization of a structure-nested lattice-material system. C-HiDeNN-Tensor Decomposition (C-HiDeNN-TD) for high-resolution topology optimization

[0091] With the development of machine learning methods, data-driven simulation and design becomes increasing popularity in recent years. Machine learning models can be used as a regression model to learn the input-output relationship such as the material stress-strain response, process- structure-property relations. However, these approaches have two difficulties: (1) it usually needs a very large database for training (2) it cannot predict values outside of the initial database. For the first challenge, the data points can be generated using reduced order. For the second challenge, researchers aim to use neural network to solve the partial derivative equations. However, solving directly the PDEs using neural networks has several issues including slow convergence rates, highly relying on initial conditions. Therefore, some fundamental developments are still needed in neural networkAttorney Docket No.: 0116936.295WO2 solving PDEs. Finite element methods have been developed for decades and are proved effectively in solving PDEs such as mechanical problem, multiphysics, coupled physics. Borrowing the idea of finite element and neural network, Hierarchical Deep Learning Neural Networks (HiDeNN) and Convolution Hierarchical Deep Learning Neural Networks (C-HiDeNN) is developed. 1. HiDeNN theory

[0092] To better understand the C-HiDeNN, a brief review of HiDeNN is presented below. The idea of HiDeNN is to use deep learning neural networks to construct the FE-like shape functions so as to solve partial differential equations for science and engineering problems. For example, for a 1D problem, the HiDeNN shape function on node i can be written as: ^^^^( ^^; ^^^∗^, ^^) ≔ ℱ^^( ^^; ^^, ^^, ^^), (1.1)

[0093] where ℱ^^stands for a deep learning neural network with inputs x, weights ^^, biases ^^, and the ^^. ^^^^denotes the shape function for the node at position ^^^∗^.

[0094] The detailed structure of the neural network blocks can be found in Zhang et al. The whole neural network is assembled by adding all the neural network blocks together, as shown in Fig.1. For a 1D problem, the input for the HiDeNN is the location x on the domain. The output is the displacement at this location ^^HiDeNN(^^)=∑ ^^^^=1^^ℎ, ^^ ^^(^^)=∑ ^^^^=^^1 c ^^^^(^^; ^^^∗^, ^^)^^^^, (1.2)

[0095] where the discretized nodalsolution of the problem, and ^^^^ ^^ ^^ ^^ ^^ ^^is the approximated solution function. After the construction of the HiDeNN shape function, the differentiation and integration of the shape function can be implemented similarly to the finite element method. The problem can be solved for displacement via optimization algorithms. Several advantages of the HiDeNN method include: (1) it is independent of an offline database; it can solve arbitrary problems with different boundary conditions; (2) it can be extended to r-h-p adaptivity to improve the mesh and increase solution accuracy; (3) it can utilize advanced machine learning hardware such as GPU and TPU. 2. Convolution-HiDeNN theory

[0096] To further improve the accuracy of the HiDeNN method, Convolution-HiDeNN (C- HiDeNN) is developed. It can be seen as a generalization of the reproducing kernel particle methods (RKPM). C-HiDeNN is based on a newly developed convolution approximation, which is formed by two functions: (1) a general polynomial function; (2) a local patch function that operates on a set of nodal points adjacent to the node as a convolution operator. Within an element e, the displacementAttorney Docket No.: 0116936.295WO2 can be written as: ^^ ^

[0097] where ^^ℯdenotes the support domain of the piece-wise polynomial associated with element e, and ^^^^^^is the local convolution patch domain for each supporting node i in the polynomialfunction. The and definition of the support domain, patch domain, and final convolution patch can be found in Fig.2. ^^^^(^^)is the polynomial function associated with node i, and ^^^^^^, ^^^^(^^)is the local convolution patch function that considers the contribution from the neighborhood ( ^^^^^^) of node i. The convolution patch function is designed to have a similar role as the kernel function in a convolution neural network. It should be noted that there are three key controlling parameters in the convolution patch functions: 1) the patch size s defining the size of the patch domain; 2) the polynomial order p; 3) the dilation parameter a. The present invention emphasises on their influence on TO design.

[0098] The convolution shape function ^^^^(^^)only changes the formulation of the original shape function ^^^^(^^)without adding extra DoFs. This unique feature can increase the smoothness and accuracy of the approximation while avoiding extra DoFs, resulting in faster computation time for a given accuracy. Fig.3 shows the speedup of C-HiDeNN compared to higher-order FEM for a given level of accuracy. A 1D Poisson's equation is solved, and H1 norm error estimator is computed. Overall, when the polynomial order p is the same, C-HiDeNN is faster than FEM for a given error. Note that higher-order FEM uses higher-order elements while both linear FEM (p=1) and C-HiDeNN use linear elements. That is, quadratic FEM (p =2) and cubic FEM (p =3) use 3-node and 4-node elements, respectively, while linear FEM (p =1) and C-HiDeNN use 2-node elements in a 1D problem. Considering this, C-HiDeNN can achieve accurate solutions orders of magnitude faster than linear FEM (p =1). For example, as illustrated in Fig.3, C-HiDeNNs with (s=1, p =1), (s=2, p=2), and (s=3, p=3) are around 4 times, 100 times, and 300 times faster than linear FEM ( ^^ = 1) to achieve 10−4^^1norm error. This clearly demonstrates the improved efficiency of C-HiDeNN.

[0099] The benefits of C-HiDeNN include:

[0100] C-HiDeNN significantly improves accuracy and convergence rates only by changing the formulation of approximation functions. The DoFs in discretized systems remain unchanged compared to FEM with the same node / mesh used;

[0101] C-HiDeNN can achieve arbitrary orders of smoothness and convergence rates by adjustingAttorney Docket No.: 0116936.295WO2 the different controlling parameters in the convolution patch function;

[0102] C-HiDeNN provides a highly scalable framework that can be accelerated using deep learning libraries and parallel computing hardware such as JAX and TPU.

[0103] The present invention focuses on how to use C-HiDeNN framework for topology optimization, which will be introduced in the following parts.

[0104] First, a brief review of the key equations of topology optimization is presented below.

[0105] Depending on different descriptions of the designed geometry, it can be categorized into but not limited to: Solid Isotropic Material with Penalization (SIMP) method, level-set method, and moving morphable components (MMC) method. This invention focuses on the SIMP method, which defines the structure with a density variable ρ( ^^). The notation used in TO is defined in Table 1. Table 1. Notation used for topology optimization ^^ Objective function ρ Design variable Π Virtual work ∇^^Differential operator ^^ Constitutive relation Ω Design domain Γ Design boundary ^^External loading ^^ External traction ρ^^ ^^ ^^Minimal density ^^(ρ) Design constraints

[0106] A general optimization problem for the SIMP method can be written as ^^ ^^ ^^  ^^(ρ) ^^. ^^. : δΠ=∫Ω ∇^^^^^^^^(ρ)∇^^δ ^^ ^^Ω − ∫Ω^^^^δ ^^ρ ^^Ω − ∫Γ^^^^δ ^^ ^^Γ = 0, (1.4) ρ^^ ^^ ^^≤ ρ ≤ 1, 0 < ρ^^ ^^ ^^≪ 1,Attorney Docket No.: 0116936.295WO2

[0107] where ^^(ρ)is the objective function with design variables ρ. The equilibrium equation is written as δΠ in the principle of virtual work form. The design variable constraint is defined as ρ^^ ^^ ^^≤ ρ ≤ 1, where ρ = ρ^^ ^^ ^^corresponds to no material, and ρ = 1 is solid material. The total volume is ^^(ρ) with ^^∗as the volume constraint. The minimal length-scale, which is the minimum allowable size of structural features in the design, is defined as ^^(ρ) and with constraints ^^(ρ) ≥ ^^^^ ^^ ^^. Other design constraints, such as stress constraints, are defined in ^^ ≤ 0. The optimization equation (1.4) can be solved by optimization algorithms such as the(OC) method or the Method of Moving Asymptotes (MMA). 3. C-HiDeNN-TD for topology optimization

[0108] C-HiDeNN-TD is a reduced-order modeling technique of the C-HiDeNN methods. For structural TO, C-HiDeNN-TD can serve the same role as FEM for solving partial differential equations and enforcing physics constraints. As discussed in the introduction, one of the major difficulties in solving a large number of simultaneous equations using numerical methods is the computational cost of large matrix inversion. To this end, the Proper Generalized Decomposition (PGD) and Tensor Decomposition (TD) based reduced order modeling techniques have been developed. The basic idea is to assume the physical quantities of the problem can be written into a form with a separation of variables. For example, the displacement ^^ = [ ^^, ^^, ^^] for a 3D design problem can be written as: ^^(^^, ^^, ^^)=∑ ^^^^=1 ^^ ^^^^(^^)^^ ^^^^(^^)^^ ^^^^(^^),

[0109] where the index^^, ^^^^, ^^^^are 1Dfunctions to be computed,^^,^^,^^ ^^,^^,^^. the design variable ρcan also be decomposed into ρ(^^, ^^, ^^)=∑ ^^^^=1ρ^^^^(^^)ρ^^^^(^^)ρ^^^^(^^). (1.6)

[0110] To the modes of ^^ and ρ, the index of mode in ^^ is defined as k, and the maximum number ofduring each design iteration, this decomposition can be done repeatedly by the Higher-Order PGD (HOPGD) method.

[0111] Consequently, the original 3D problem is converted to finding the corresponding 1D modes of those ^^^^, ^^^^, ^^^^, ^^^^, ^^^^, ^^^^, ^^^^, ^^^^, ^^^^., and ρ^^, ρ^^, ρ^^. Therefore, the computationalAttorney Docket No.: 0116936.295WO2 complexity is significantly reduced.

[0112] The equilibrium equation in the form of the principle of virtual work can be written as δΠ = ∫Ω ∇^^δ ^^^^(^^ρ)∇^^^^ ^^Ω − ∫Ωδ ^^^^ ^^ρ ^^Ω − ∫Γδ ^^^^^^ ^^Γ = 0, (1.7)

[0113] where ∇^^is the symmetric gradient operator, and ∇^^^^ is the strain in Voigt notation. Considering the separation of variables, the principle of virtual work becomes

[0114] where

[0115] From the equation 1.8, the whole 3D problem has been decomposed into several 1D problems. The C-HiDeNN approximation is used for each separated function. The element-wise approximations of solutions then read

[0116] where values of each separated, , , , , , ^^^^^^, ^^^^^^, ^^^^^^denote the global nodalequation 1.8 can be discretized and rearrangedAttorney Docket No.: 0116936.295WO2 to become

[0117] Finally, thethe above set of equations.

[0118] The directional stiffness matrix of each mode for the separated displacement function is defined as ^^ with the integrated external force ^^. Equation 1.10 can be solved sequentially by looping over each mode of ^^. The number of modes is determined by setting a convergence criterion. After convergence, all modes can be optimized again to increase the accuracy further.

[0119] It is evident that for each iteration in Equation 1.10, the dimension of the directional stiffness matrix is only determined by the number of DoFs in each direction. This results in a drastic reduction in the matrix dimension compared to the traditional FEM. For a problem with a mesh size of 1000 × 1000 × 1000 using linear elements, the DoFs of the problem including design density on each node is (1001 × 1001 × 1001 × 4) → (4,012,012,004).

[0120] Solving such a large number of equations is extremely expensive. However, in the present invention, the number of the DoFs for every one dimensional problem would be only 4004. Ultimately, this can lead to orders of magnitude reduction in the overall computational cost.

[0121] With the new C-HiDeNN-TD solver, the TO flow chart can be redesigned as in Fig.4. The overall procedure is similar to the traditional nested TO. The main difference comes from the solver part. In the design framework of the present invention, the C-HiDeNN-TD is used for solving the equilibrium equation 1.7, whereas traditional TO usually uses FEM. The initial design variable is set to be uniformly distributed and satisfies the volume fraction constraint in the whole design domain. The design variable ρ is decomposed into three directions for a 3D problem or two directions for a 2D problem, as defined in equation 1.6. In one embodiment, the present invention uses HOPGDAttorney Docket No.: 0116936.295WO2 to perform the decomposition of the whole density field. In the flowchart, this step refers to the block "Decompose design variables into three directions". After the decomposition, the decomposed density is sent to the C-HiDeNN-TD solver with the convolution shape function defined in equation 1.9. The design variables are optimized using the Optimal Criterion (OC) method. With each design variable update, a convergence criterion is tested to determine whether or not the design has converged.

[0122] The present invention remarks that the overall design framework is flexible and can be coupled with traditional FEM or other C-HiDeNN methods in a hybrid way. For example, FEM can be used during the first few iterations with a coarse mesh to get some initial design before ultimately switching to C-HiDeNN-TD with a higher resolution to achieve finer structures that perform better. 4. Convolution interpolation for built-in density filtering and higher-order density smoothness without extra DoFs

[0123] This section studies the effect of C-HiDeNN approximation on the optimized density distribution. In TO with the SIMP approach, the design variable, i.e., density field, is considered constant inside each element of the background mesh. Let $\Omega_e$ denote an element region. The SIMP-based density description is defined as ρSIMP(ξ)= ρ^^,   ∀ξ ∈ Ω^^(1.11)

[0124] Consequently,density distributions and checkerboard patterns across multiple elements. Using higher-order elements for a smooth description of the density field is one of the solutions to resolve this issue. However, using higher- order finite elements can result in excessive increases in the computational cost of TO due to the additional DoFs in both the physical equations and the design variable. Another way to mitigate checkerboarding and enforce length-scale control is the filtering technique. Several filter methods have been developed, including sensitivity filter, density filter, Heaviside filter, and averaging filter. However, the extra filtering stage increases the computational burden in the design loop, especially for high-resolution problems where the filter may be applied over millions or billions of design DoFs.

[0125] From the above perspectives, the C-HiDeNN approximation is advantageous in providing a built-in filter to improve the design smoothness without increasing the overall number of design variable DoFs and DoFs of the numerically discretized physical system. Therefore, the proposed C- HiDeNN-TD-TO framework of the present invention leads to a higher-order smooth density description compared to SIMP if the same background mesh (resolution) is used.

[0126] The C-HiDeNN approximated density field can be defined using the nodal values of theAttorney Docket No.: 0116936.295WO2 background mesh. Similar to the displacement field, the nodal density can be denoted by ρ^^, and the approximated density then reads:

[0127] where the support domain ^^ℯand local convolution patch domain ^^ are defined the same as for the displacement approximation. The polynomial function associated with node i is ^^^^(ξ) , and the local convolution patch function is ^^ξ^^,^^^^(ξ). In a general sense, they can be different from the displacement approximation.

[0128] The dilation parameter a, the polynomial order p, and the patch size s serve as built-in filtering parameters to control the minimum design length-scale and remove checkerboard patterns. Fig.5 shows the difference between the SIMP-TO approach and the C-HiDeNN-TD-TO approach. Unlike SIMP-TO, the density is interpolated using C-HiDeNN-TD-TO shape functions and the convolution operation serves as a filter on the density field, which will automatically prevent checkerboarding and incur length-scale control. One example is shown in Fig.6. To illustrate the idea, the present invention generates a density field with a checkerboard pattern, as shown in Fig.6. This is commonly seen if the filter is not applied during the design iteration using SIMP. The present invention applied the density filter and the C-HiDeNN built-in filter on the same density field. The minimum length-scale of the density filter is set to 2. For a fair comparison, the patch size s is set to 2 with a dilation parameter a=2 and polynomial order p=2. Fig.6 shows that the C-HiDeNN built-in filter has the capability to remove a checkerboard pattern similar to a dedicated density filter. It is noted that for some 3D examples with a volume fraction less than 0.3, the C-HiDeNN-TD-TO built- in filter may be insufficient to remove all isolated substructures. More studies are needed, including designing a more stable built-in filter and finding optimal parameter sets.

[0129] The other feature is that the convolution interpolation does not increase the DoFs in the design variables but improves the smoothness of interpolation.

[0130] One example to compare the difference between the SIMP approach and the C-HiDeNN- TD-TO approach is shown in Fig.7, which depicts a cantilever beam design example with a mesh size of 100 × 50 × 50 using SIMP and C-HiDeNN descriptions. As shown in Fig.7 panel (a), the density in the SIMP approach is discontinuous in the design domain. For C-HiDeNN-TD-TO, the density is smooth, with a desired higher-order continuity defined by the C-HiDeNN shape function. To represent the C-HiDeNN-based continuous density function in the 3D domain, the presentAttorney Docket No.: 0116936.295WO2 invention selects a few interpolation points in each element. Fig.7 panel (b) shows the plot of C- HiDeNN-TD-TO with 27 interpolation points in each element. It can be seen that, without adding any additional DoFs in the design, the C-HiDeNN-TD-TO gives a smoother design structure than the SIMP-TO. With more interpolation points added, the plot can be a very smooth structure with details on the surface. It should be noted that although the C-HiDeNN approximation increases the smoothness of the design, it can not lead to the same results as high-resolution design since higher resolutions can result in greater design flexibility creating finer structures with better design objective values. 5. Parametric study of the convolution parameters on the final design

[0131] In this section, we study the effect of the convolution parameters on the final design structures. As mentioned earlier, there are three parameters in the convolution shape functions: the patch size s defining the size of the patch domain; the polynomial order p defining the polynomial order of ^^ξ^^,^^^^(ξ); the dilation parameter a in ^^ξ^^,^^^^(ξ) controlling the size of the window function, as shown inpanel (b).

[0132] As shown in Fig.8, the cantilever beam is used with the left side clamped and a uniform loading on the right side. The background design mesh size is 100 × 50 × 50. The design formulations are defined as: ^^ ^^ ^^ & ^^(ρ) =∫^^^^^^ ^^∫^^Ω Ω + Γ ^^ ^^ ^^Γ ^^. ^^. : δΠ=∫Ω∇^^^^^^^^(^^)∇^^^^ ^^ ^^Ω −∫Ω^^^^δ ^^ρ ^^Ω −∫Γ^^^^δ ^^ ^^Γ = 0, ρ^^ ^^ ^^≤ ρ ≤ 1, 0 < ρ^^ ^^ ^^≪ 1, (1.13) ^^^^ ^^ ^^.

[0133] The objective function is system compliance as defined in equation 1.13. The volume fraction ^^∗is set to 0.5 for the design region. The convolution parameters are defined as follows. The patch size s is defined as the size of the extended region outside the element (bold green region). The dilation parameter a is defined as the tuning parameter of the window function. The polynomial order of the convolution patch function is defined as p. The present invention fixes the patch size s but vary the polynomial order p and the dilation parameter a.

[0134] Fig.9 shows the different results related to the polynomial order p. The patch size is set as s=4. To better understand the synergism of C-HiDeNN with the tensor decomposition, the presentAttorney Docket No.: 0116936.295WO2 invention, in one embodiment, sets Fig.9 panel (a) with a very small and positive a, which produces the Tensor Decomposition-TO (TD-TO) result. From these results, higher order p does give a better design (lower objective function) with more support structures, especially when compared to p=1 (the case of TD-TO). The reason is that the higher-order approximation gives a smoother and more accurate solution of the physics field and significantly improves the sensitivity of the density to displacement. Finally, the optimal design converges with p=3 in this case.

[0135] Fig.10 shows another example to test the dilation parameter a. The polynomial order and patch size are p=3, s=4. From the figure, it is obvious that a large dilation parameter gives a lower objective function. The optimal design (objective function) converges with a =4. Furthermore, the dilation parameter seems to contribute to the length-scale filtering, as a larger a might remove more small support structures in the final design. This can be seen inside the black circles in the figure.

[0136] From the above two comparisons, one can see that the greater the polynomial order, p, and the larger the dilation parameter, a, the more accurate the solved displacements are, leading to better topology-optimized designs. However, it should be noted that there exists an optimal pair of p and a for each problem. Using larger numbers for p and a can result in a higher computational cost. The balance between accuracy and cost should be considered when choosing the optimal convolution parameters. 6. Numerical examples for high-resolution design

[0137] Problem 1: high-resolution cantilever beam

[0138] The first example is still the cantilever example, as illustrated in Fig.8. This time, the efficiency of C-HiDeNN-TD-TO is studied with an increasing number of DoFs (higher resolution). In order to make a comparison to the example in existing studies, the volume fraction is set to be 0.12, and the minimal length-scale parameter ^^^^ ^^ ^^in equation 1.4 is set to be 0.1 of the short edge. The objective function and optimization formulations are same as equation 1.13. Three different designresolutions are considered at this time: 480 × 240 × 240, 960 × 480 × 480, and 1440 × 720 ×720. The dilation parameter a is set to 4 with the polynomial order p=3 and a patch size s=4.

[0139] The computational cost with different DoFs for FEM-TO and the proposed C-HiDeNN- TD-TO are shown in Fig.11. The number of FEM DoFs also considering the density variable on each node. The CPU hours in FEM-TO with the SIMP method are estimated according to work done by N.Aage et al. using supercomputer resources. For the FEM-TO method, the CPU time increases from104 hours to 106 hours as the DoFs increase from 108 to 109. However, for C-HiDeNN-TD-TO, theAttorney Docket No.: 0116936.295WO2 CPU time only increases from 33 hours to 112 hours, with a 6,857 times speedup for the larger mesh when compared to the FEM-TO method. Table 2 lists the detailed numbers of the three cases, including the number of equations, the number of modes in the C-HiDeNN-TD-TO method, the CPU time, and the speed-up relative to FEM-TO.

[0140] The final design structures of the three cases are shown in Fig.12. As the design resolution increases, the overall structure becomes more and more similar. For the same volume fraction and the given length-scale control, the topology design have a convergent result. However, for some regions, such as the region in a dashed circle, the increasing resolution results in more support structures and an improved, lesser objective function. Table 2. Comparison of speed-up in different cases. The CPU time of the FEM-TO method is estimated Resolution FEM DoFs C-HiDeNN- Number of CPU time Speed- TD-TO DoFs modes up ρ ^^ FEA-TO C- (estimated) HiDeNN- TD-TO 480×240×240 111,747,844 964,926 34 334 9,000h 33h 273 960×480×480 889,351,684 3,830,616 36 664 72,000h 64h 1,125 1440×720×720 2,996,363,524 7,273,809 42 841 768,000h 112h 6,857

[0141] Problem 2: Drone design problem

[0142] In the second example, the C-HiDeNN-TD-TO is applied to a drone design problem shown in Fig.13. The design region is a rectangular box with loads on the four top corners and the center of the bottom. For the numerical implementation, the four corners are treated as the fixed boundary conditions in all directions. Furthermore, there is a void region in the design space to leave room for batteries and a control unit, as shown in Fig.13. The design goal is to minimize system compliance under some predefined design constraints on the volume fraction and minimum length- scale control, as defined in equation 1.13. Here the present invention studies two different volume fractions, 0.12 and 0.2, with the minimum length-scale set at 0.03L as shown in Fig.13. The void region is performed by enhancing a passive region in the design loop, which means the void region is forced to have zero density at each design iteration.

[0143] The final TO structures are shown in Figs.14-15 with different design resolutions. In Fig. 14, the drone is designed with a volume fraction constraint of 0.2, whereas in figure Fig.15, the drone is designed with a volume constraint of 0.12. It can be observed that the lower volume fraction willAttorney Docket No.: 0116936.295WO2 produce thinner truss structures. The objective function of design with volume constraint 0.2 is 3.3618, and the objective function with volume constraint 0.12 is 3.4162.

[0144] To summarize, the key features of the C-HiDeNN-TD-TO are:

[0145] Reduces the computational cost, enabling thousands of times speedup for giga-scale topology optimization problems and efficient concurrent multi-scale topology design.

[0146] Improved solution accuracy compared to traditional SIMP methods by using a convolution approximation that leads to higher-order smooth descriptions of structures for a given mesh.

[0147] The convolution approximation provides a built-in length-scale filter for topology optimization and avoids the extra filtering step in traditional topology optimization.

[0148] Local design refinement can be easily achieved due to the location dependence of the convolution shape functions.

[0149] The proposed topology design framework can be further accelerated by GPU and has considerable potential to enable fabrication-specific topology design for additive manufacturing, which usually requires extreme model resolutions to predict the physical shape properties accurately. C-HiDeNN-TD-TO democratizes topology optimization for practical problems by preserving solution accuracy while alleviating the formerly extravagant memory and computation requirements. Examples Example 1 C-HiDeNN-TD2for unit cell design and 2-scale concurrent design

[0150] In one embodiment, the present invention focuses on using C-HiDeNN-TD for lattice material design. Theories are presented and implementation are shown with three examples: maximum bulk modules design, maximum shear modulus design, and negative Poisson's ratio design. In another embodiment, the present invention extends the C-HiDeNN-TD to a 2-scale concurrent design. Concurrent design aims to optimize the two scales together. It can have a better performance, light weight, and multi-functionalities. A drone example with different lattice design regions is used to illustrate the theory. C-HiDeNN-TD for lattice unit cell design

[0151] Lattice structures are manufactured using a repetitive pattern and find numerous applications. When studying these materials, one typically considers homogenized material properties. The challenge that researchers aim to address is how to design a repetitive lattice structure that yields the desired properties in the global homogenized material.

[0152] To start with, the present invention firsts define the design formulations:Attorney Docket No.: 0116936.295WO2

[0153] min c ( ^^^^(ρ)) ^^. ^^. : δΠ=∫Ω∇^^^^^^^^∇^^^^ ^^ ^^Ω − ∫Γ^^^^δ ^^ ^^Γ = 0, ^^ is periodic on Γ ^^ ^^ ^^≤ ≤ <^^ ^^ ^^

[0154] function is related to the homogenized material property. Details are presented in the following examples. The equilibrium equations are similar to the part scale case with an exemption of the body force. The displacement ^^ should be periodic on the boundaries. Other constraints are similar to the part scale. For the homogenization theory, the present invention uses the element mutual energies: ^^^^^^^^^^ ^^=1∫Ω^^^^ ^^ ^^ ^^(ε0^^(^^^^ ^^)− ε^∗^(^^^^ ^^)) ^^Ω , (2.2)

[0155] where the ^^^^^^^^order tensor and can be reduced using assumptions. The domain of lattice unit cell is defined as Ω. The material stiffness in building up the lattice is defined as ^^^^ ^^ ^^ ^^. The eigenstrain is defined as ε0^^(^^^^ ^^)with ε∗^^(^^^^ ^^)as the response.

[0156] To solve the homogenization part, in one embodiment, the present invention extends the C-HiDeNN-TD with enforcing of periodic boundary conditions. Recall the equilibrium equations: ^^. ^^. : δΠ=∫Ω∇^^^^^^^^∇^^^^ ^^ ^^Ω − ∫Γ^^^^δ ^^ ^^Γ = 0, (2.3) ^^ is periodic on Γ

[0157] The solution of using C-HiDeNN-TD is similar to the part scale, where one has ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^, ^^^^^^, ^^having periodic boundary conditions. The solution scheme can be found in Fig.16.

[0158] It has the similar pattern of the part-scale design with the difference in the periodic boundary conditions. In the implementation procedure, it will be a 6 calculations with 6 kinds of eigenstrain. It should be noted that after C-HiDeNN-TD homogenization, it will be sent to a Preconditioned Conjugate Gradients (PCG) procedure to increase the accuracy of the homogenization. After that, the gradient of objective function respect to each voxel can be calculated. Numerical exampleAttorney Docket No.: 0116936.295WO2

[0159] In certain embodiments, three examples are presented to show the C-HiDeNN-TD in designing the lattice material with desired property. The first example is the maximum bulk modulus design, where the objective function is ^^11+ ^^22+ ^^33+ ^^12+ ^^23+ ^^13. The second example is the maximum shear modulus design with the objective function ^^44+ ^^55+ ^^66The third example is the negative Poisson's ratio design, where the objective function is ^^12+ ^^23+ ^^13− ^^+ ^^22+ ^^33). The parameter β is the tuning parameters to ensure the convergence. The 1 to the design loop. The volume fraction of three cases are set to 50% of the unit cell.

[0160] The designed unit cells are shown in the Fig.17. The inside of the maximum bulk modulus design is a hollow sphere shape. The cross section of the maximum shear modulus design is similar to a cross shape. This kind of structure is commonly seen in the engineering applications such as the bridge. The shape of negative Possion's ratio has similar pattern to the existing work. It should be noted that the negative Possion's ratio problem may rely on the initial guess and may require the optimal choice of the tuning parameters.

[0161] The optimized objective function and the operation count can be found in Table 3. For the bulk moduli problem, the optimized bulk moduli can be increased from 0.2370 to 1.5355. For the shear moduli problem, the optimized shear moduli can be increase to 0.1297. For the negative Poisson's ratio problem, it can be reduced to -0.5412. Another improvement is the computational speed. All three cases with the resolution 400 × 400 × 400 can be done with around 50h for a single PC. As a reference, in Aage et al, for a 288 × 288 × 288 unit cell design, the CPU hour can reach 1000 h for 250 iteration steps. The speed up on the unit cell design gives the potentials in concurrent design and nested concurrent design. Table 3. Material constants for matrix and inclusion Objectives Resolution Initial Value Optimized Value Computation time Maximize bulk moduli 400 × 400 × 400 0.2370 1.5355 56.2h Maximize shear moduli 400 × 400 × 400 0.0124 0.1297 54.4h Negative Poisson's ratio 400 × 400 × 400 0.3 -0.5412 53.2h C-HiDeNN-TD2for 2-scale concurrent design.

[0162] In one embodiment, the present invention discloses C-HiDeNN-TD2for concurrent multi- unit cell design. Concurrent multi-scale design is based on the multi-scale theory and aims toAttorney Docket No.: 0116936.295WO2 optimize the two or more scales simultaneously. Considering homogenization theory, the present invention extends the original C-HiDeNN-TD to a two-scale design method called C-HiDeNN-TD2.

[0163] The design flowchart for C-HiDeNN-TD2is shown in Fig.18. The whole design framework is a combination of part scale design and unit cell design.

[0164] The design problem definition of the drone structure is shown in Fig.19. The boundary conditions and the loading type are the same as in the previous example. The difference is in the predefined regions of the design domain. For each region, the present invention assumes it consists of the same unit cells. At the initial design stage, all the unit cells of each region have the same initial set-up with a volume fraction of 0.5.

[0165] The drone structure and the unit cells are simultaneously updated during the concurrent design stage. Results are shown in Fig.20. From the figure, different regions have different unit cell structures that better reduce the objective function. By applying the reinforcing connectivity conditions, the present invention ensures each unit cell is perfectly connected with its neighbors. This example shows one potential application of the C-HiDeNN-TD2method. In Example 3, it can be extended to an n-scale concurrent design framework with high-resolution design routines at each scale. Example 2 Clustering discretization methods for generation of material performance databases in machine learning and design optimization

[0166] In certain embodiments, the present invention discloses how to use clustering discretization methods to generate the material performance database, and how to use the generated database for design and optimization. Mechanical science and engineering can use machine learning. However, data sets have remained relatively scarce; fortunately, known governing equations can supplement these data. Three reduced order methods are summarized and generalized: self-consistent clustering analysis, virtual clustering analysis, and FEM-clustering analysis. These approaches have two-stage structures: unsupervised learning facilitates model complexity reduction and mechanistic equations provide predictions. These predictions define databases appropriate for training neural networks. The feed forward neural network solves forward problems, e.g., replacing constitutive laws or homogenization routines. The convolutional neural network solves inverse problems or is a classifier, e.g., extracting boundary conditions or determining if damage occurs.

[0167] The present invention discloses how these networks are applied, then provide a practical exercise: topology optimization of a structure (a) with non-linear elastic material behavior and (b)Attorney Docket No.: 0116936.295WO2 under a microstructural damage constraint. This results in microstructure-sensitive designs with computational effort only slightly more than for a conventional linear elastic analysis.

[0168] Importantly, the present invention discloses a new way to link the structure scale with the material scale, and provide some ideas how to design such a structure. Overview of two-stage clustering analysis methods

[0169] Self-consistent Clustering Analysis (SCA), its close relatives Virtual Clustering Analysis (VCA) and FEM Clustering Analysis (FCA) are two-stage reduced order modeling approaches consisting of an offline data compression process and an online prediction process. This is concisely illustrated in Fig.21. In the offline stage, the original high-fidelity Representative Volume Element (RVE) represented by voxels or elements is compressed into clusters. In the online stage, macroscopic loading is applied to the reduced order (clustered) RVE. The system of equations describing mechanical response is then solved on only the reduced representation as a boundary value problem. The notation used in this section is summarized in Table 4 and 5. Table 4 Notation used for two-stage clustering analysis methods (1) ^^ Material point ^^′Any other material point ^^Normal to boundary ^^M Macroscale stress tensor ^^M Macroscale strain tensor σ( ^^) Microscale stress tensor ε(^^)Microscale strain tensor ^^Unit eigenstress ^^ Unit eigenstrain Γ( ^^, ^^′)Green's operator Γ0Fourier transform of the periodic Green's operator I Counting index for clusters J Another counting index for clusters χ^^(^^) Characteristic function for cluster I ^^^^ ^^Interaction tensor for discretized Lippmann-Schwinger equation ^^ “Interaction tensor” for FCAAttorney Docket No.: 0116936.295WO2 Ω Material domain ΩIDomain of the Ithcluster (subset of Ω) ^^^^Volume fraction of the I cluster ε̃ Reference material strain ^^Reference material stiffness ^^Compliance matrix for reference material ^^ Strain concentration tensor Table 5 Notation used for two-stage clustering analysis methods (2) ^^0Reference deformation gradient ^^′Deformation concentration tensor λ0, μ0Lame's constants of the homogeneous stiffness tensor ξ Fourier point ℱ, ℱ−1Forward and inverse fast Fourier transform techniques △^ Incremental form of an arbitrary variable ^ ^^^^Residual of the Ithcluster ^^ Jacobiam matrix of ^^ with respect to Δε ^^4Fourth-order identity tensor ^^^^^^^^ ^^Tangent stiffness of the material in the Jthcluster ^^^^Number of clusters ^^^^Number of Fourier points ^^ ^^Number of integration points ^^^^Number of finite elements ^^^^ ^^ ^^ ^^ ^^ ^^Young's Modulus of the matrix ν^^ ^^ ^^ ^^ ^^ ^^Poisson's Ratio of the matrix ^^^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^Young's Modulus of the inclusion ν^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^Poisson's Ratio of the inclusionAttorney Docket No.: 0116936.295WO2 σ^^,  ^^ ^^ ^^ ^^ ^^ ^^Yield strength of the matrix ε̅PEffective plastic strain (matrix)

[0170] To define the boundary value problem, consider a material occupying ^^ ⊂ ^^^^. The goal of homogenization is to find the macroscopic constitutive relation between a macroscopic stress 1 ^^ ^^ ^^

[0171] and aεM=1|^^| ∫ ^^ε( ^^)d ^^ , ^^ ∈ ^^, (3.2)

[0172] where | ^^| is the total volume of the region Ω. The present invention defines mathematically the RVE problem as

[0173] where σ = σ(ε;For the homogenization problem, boundary conditions have to be chosen to satisfy the Hill-Mandel condition. In one embodiment, the periodic boundary conditions are used: ε periodic and σ ⋅ ^^ anti-periodic on ∂ ^^. Continuous and discretized Lippmann-Schwinger equation

[0174] By introducing a reference material with distributed elastic stiffness ^^(^^), it has been shown that the RVE problem with periodic boundary conditions is equivalent to the integral equation: ε( ^^) = ε̃( ^^) −∫′ ′ ′ ′ ′^^ Γ ( ^^, ^^ ): (σ( ^^ ) − ^^ ( ^^ ): ε( ^^ )) d ^^ (3.4)

[0175] andboundary as σ − : ε applied to the reference material; Γ ( ^^, ^^′) is the Green's with the reference ′material. The physical meaning of −Γ^^ ^^ ^^ ^^( ^^, ^^ ) is the ε^^ ^^at material point ^^ if the unit eigenstress ^^^^ ^^is applied at material ^^′, with the components of skldefined by ^^^^^^^^^^= δ^^ ^^δ^^ ^^, where δ^^ ^^and δ^^ ^^arefunctions.

[0176] The integral equation given in Eq.3.4 is known as the Lippmann-Schwinger equation, commonly seen to describe particle scattering in quantum mechanics. There is typically no explicitAttorney Docket No.: 0116936.295WO2 form of Γ ( ^^, ^^′), unless the reference material is homogeneous.

[0177] The domain can be decomposed into several sub-regions, called clusters, distinguished mathematically by the characteristic function χ as shown in Eq.3.5. ^^

[0178] where ^^ = 1,2,3, … ,^^subset of the volume within cluster I. These clusters are defined during the offline stage. This allows one to discretize the Lippmann-Schwinger equation, as: ε^^= ε^^− ∑^^^^ ^^^^ ^^:^^− ^^^^: ε , ∀ ^^ ∈ … , ^^

[0179] ^^ ^^| ^^| ^^variable ^ in the Ithcluster; ^^^^is the volume fraction of the Ithclth|^^|uster where the volume of the I cluster is given by|^^^^|. ^^^^ ^^is the interaction tensor given by Ω ^^^^ ^^=1 ^^^^ ^^| ^^| ∫ ^^ ∫ ^^χ ( ^^)χ^^( ^^′)Γ ( ^^, ^^′)d ^^d ^^′. (3.7)

[0180] Theij in the Ithcluster if the uniform unit eigenstress component kl is applied in the Jthcluster.

[0181]

[0182] Remark 1: If homogeneous reference material is used, as in SCA and VCA, ε^^= ε^^.

[0183] Remark 2: A counterpart of Eq.3.6, similar to the idea of FCA, can be expressed as σ^^= σ^^−∑ ^^^^^^=1 ^^^^ ^^: (ε^^− ^^^^: σ^^), (3.8)

[0184] where ^^ is the is the stress in the referencematerial when applying same as original RVE problem. ε^^− ^^^^: σ^^is the volume average eigenstrain in the Ithcluster. The physical meaning of − ^^^^ ^^^^ ^^ ^^ ^^is the average stress component ij in the Ithcluster if the uniform unit eigenstrain component kl is applied in the Ithcluster. Note that the reference material for FCA is the elastic state of the original RVE, instead of a homogeneous reference material as used by the other two methods here.

[0185] Remark 3: The relationship between ^^^^ ^^and ^^^^ ^^is given by ^^^^: ^^^^ ^^= − ^^^^ ^^: ^^^^. (3.9)

[0186] This can be proven by noting that the effect of applying any unit eigenstrain kl in some cluster J of the reference material is equivalent to that of applying eigenstress − ^^ : ^^^^ ^^in the sameAttorney Docket No.: 0116936.295WO2 cluster.

[0187] In the online stage, SCA solves the incremental, discretized Lippmann-Schwinger equation, Eq.3.6, with arbitrary external loading conditions. These can either be of the fixed strain increment △ ε^^type or be of the fixed stress increment △ σ^^type. Offline: Clustering and Interaction Tensor

[0188] The offline stage consists of three primary steps: (1) data collection, (2) unsupervised learning (e.g. clustering), and (3) pre-computation of the interaction of clusters. Computation of the linear elastic response (data collection) and subsequent clustering based on that response are conducted identically for each of the three two-stage methods presented below: the same voxel mesh and clusters are used in the examples for all three methods. The difference between methods comes in the computation of the interaction tensor.

[0189] Data collection

[0190] Data collection provides information used to construct the reduced representation of the system. It typically involves some computation of the response of a fully resolved system, perhaps with a simplified material model, over a limited set of loading cases. One measure of mechanical response that could be collected is the strain concentration tensor used in micromechanics ^^, defined by ε( ^^) = ^^( ^^): ε^^( ^^), which maps between the far-field or applied strain, ε^^and the strain measured at point in the domain, ε.

[0191] The choice of the strain concentration tensor as given above is often suitable, but it depends on the relevant details of the problem. For example, at finite deformations one might consider using the deformation concentration tensor: ^^′( ^^) =∂ ^^( ^^)∂ ^^0^^ℎ ^^ ^^ ^^ ^^( ^^) is the deformation gradient at any given point ^^ within the domain and ^^0is thedeformationcorresponding to the boundary conditions. Alternatively, if and plastic material responses differ substantially, e.g. if one is isotropic and the other anisotropic, including information about the plastic part of the deformation might be desirable.

[0192] Clustering

[0193] The goal of clustering is to reduce the number of degrees of freedom required to represent the system while minimizing the loss of information about the mechanical response. One way of doing this, as alluded to above, is by grouping material points within the domain of interest. If one can assume that the material response within each group is identical, the evolution of the domain can be determined by solving for the response of each group rather than each material point.Attorney Docket No.: 0116936.295WO2

[0194] Using the data generated during the collection phase, one of the many clustering (or unsupervised learning) techniques might be applied to optimized the domain decomposition. In the following examples, Self-Organizing Maps (SOMs) are employed. This is illustrated in Fig.21, where eight clusters (four in each phase) are constructed. The clustering process assigns each material point with a cluster ID, such that clusters are labelled 1,2, … , ^^^^.

[0195] The k-means clustering method has also been used. More elaborate clustering schemes might also be considered. Ongoing efforts include “adaptive clustering”' schemes that mimic adaptive FE methods in their ability to evolve as deformation progresses, and “enriched” machine learning, whereby a priori information from mechanics about the deformation fields outside the bounds of the data collected in Step 1 are used to guide or bound the unsupervised learning. Another unexplored future direction might use a “feature-based” machine learning that includes microstructure information in addition to mechanics information.

[0196] Interaction tensor

[0197] The interaction tensor describes the impact each cluster has on each of the other clusters. Once the clustering process is completed, the interaction tensor can be explicitly computed. Importantly, the integral part only has to be computed once during the offline stage. Only the results of that calculation are then used for the online stage. Three ways to compute the interaction tensor, one used by each of the methods highlighted here (though these are not exclusive to each), are:

[0198] SCA: Fourier transform for ^^^^ ^^. With periodic boundary conditions and homogeneous reference material, the Green's operator has a simple expression in Fourier space, given by Γδ ξ ξ 0 ξ ξ ξ^0^^ ^^ ^^ (ξ) =^^ ^^ ^^ ^^ λ ^^ ^^ ^^ξ ^^^2μ0|ξ|2−2μ0(λ0+2μ0) |ξ|4(3.10) [0 where Γ^^0^^operator Γ0; λ0and μare the of the homogeneous stiffness tensor; ξ is the Fourier point. Then thebe calculated with ^^^^ ^^=1^^ ^^| ^^| ∫ ^^χ^^( ^^)ℱ−1(ℱ(χ^^)ℱ(Γ0))d ^^ , ∀ ^^, ^^ ∈ {1, … , ^^^^} (3.11)

[0200] is ^^((^^^^)2(^^^^)

[0201] VCA: numerical integration for ^^^^ ^^. With an infinite homogeneous reference material, the Green's operator can be expressed in real space. Numerical integration is the most straightforward method to compute the integral equation in Eq.3.7. The computational complexity is ^^(( ^^^^)2),Attorney Docket No.: 0116936.295WO2 where ^^^^is the number of integration points used. In one embodiment, a fast method is proposed to approximate ^^^^ ^^.

[0202] FCA: finite element method for ^^^^ ^^. Based on the physical interpretation of the interaction tensor, the finite element method can also be used. By applying uniform unit eigenstrain component kl in the Jthcluster, the average stress can be computed for all clusters, resulting in ^^^^ ^^^^ ^^ ^^ ^^for all ^^ =1, … , ^^ ^^. Thus, the computational complexity is ^^ ( ^^ ^^)( ^^ ^^)), where ^^ ^^ is the number of finiteelements used. The tensor ^^^^ ^^is similar to the tensor ^^^^ ^^, although ^^^^ ^^is determined by applying strains rather than stresses. Online: Reduced order response prediction

[0203] Once the interaction tensor database is prepared, the discretized Lippmann-Schwinger equation defined in Eq.3.6 solved in the online stage. An incremental form of Eq.3.6 is given by △ε^^=△ ε^^−∑^^^^^^=1 ^^^^ ^^: (△ σ^^− ^^^^:△ ε^^), ∀ ^^ ∈ {1, … , ^^^^}, (3.12)

[0204] △ σ^^is a function of the incremental strain △ ε^^according to the local material constitutive laws. So the unknowns of Eq.3.12 are the strains in each cluster {△ ε} = {△ ε1, … ,△ ε^^ ^^}. For nonlinear microscale constitutive laws, Eq.3.12 is nonlinear and has to be solved iteratively. Newton's iterative method is used by SCA and VCA, while a different iterative method is used by FCA. The residual form of Eq.3.12 is given by ^^^^=△ ε^^−△ ε^^+∑ ^^^^=^^1 ^^^^ ^^: (△ σ^^− ^^^^:△ ε^^), ∀ ^^ ∈ {1, … , ^^^^} (3.13)

[0205] ^^^^ ^^=∂ ^^^^ ∂△ε^^= δ^^ ^^ ^^4+ ^^^^ ^^:( ^^^^^^ ^^ ^^− ^^ ), ∀ ^^, ^^ ∈ {1, … , ^^^^} (3.14)

[0206] in the Jth^^cluster is ^^^^ ^^ ^^. the complementary energy of the clustering-based system. Comparison of the methods

[0207] In order to compare the accuracy and efficiency of each method, a 2D plane strain model for a two-phase material is constructed and shown in Fig.22. The 2D mesh contains 600 × 600 square pixels. The inclusion area fraction is 51%. Material constants of the matrix and the inclusion are given in Table.6. The yeild suface is von Mises surface as shown in Eq.3.15. The hardening lawAttorney Docket No.: 0116936.295WO2 of the matrix material is given in Eq.3.16. Note that this will be approximated by a non-linear elastic behavior for monotonic loading in future sections. Table 6. Material constants for matrix and inclusion of

[0208] The equivalent von Mises stress is σ̅. The yield stress σ^^, ^^ ^^ ^^ ^^ ^^ ^^is given by the hardening law in Eq.3.16 with equivalent plastic strain̅ε̅ ^ ̅^. Strain from 0 to 0.05 is prescribed in the x-direction and zero strain is enforced in the y- and xy-directions. A direct numerical simulation (DNS) of the microstructure under these loading conditions is performed using the FEM and the effective von Mises stress is recorded for later comparison to the reduced order model results.

[0209] A one-time data-compression of the microstructure is performed following the steps outlined offline section, using the strain concentration tensor. The resulting clustering of the microstructure is shown in Fig.23. The interaction tensor for each method is computed using the aforementioned algorithms in Section regarding offline clustering.

[0210] The different methods, used to compute the interaction tensors, each result in a slightly different form of the tensor. There are strong similarities - after all, the same microstructure and clustering is used - though the details differ. Fig.24 shows a magnitude plot of each of the three methods, where magnitude represents the effect of each stress component in the ^^^^ℎcluster on the corresponding strain component in the ^^^^ℎcluster. The trends of the magnitudes shown in Fig.24 and Fig.25 suggest that inter-cluster interaction is not as strong for VCA as for SCA and FCA, although this difference is relatively minor. The strongest interactions are intra-cluster, as shown by the peaks along the diagonal. FCA has two distinct regions of peaks, corresponding to the set of clusters in the inclusion and in the matrix. The tensor is constructed in an ordered way, which results in these two distinct sets of clusters. This is unlikely to change the overall solution accuracy.Attorney Docket No.: 0116936.295WO2

[0211] Once the microstructural database is created and the interaction tensor has been computed, the online prediction is performed. Fig.26 shows that the ROM results for stress in the x-direction in all three cases are within 5% of the DNS results\footnotemark, with small differences between the three ROMs. Fig.26 panel (b) highlights the slight differences between the three methods: SCA follows the same trend as the DNS, but is slightly softer; the response predicted by VCA is softer still, as a result of the constant reference material assumption; the trend exhibited by FCA is slightly different from the DNS overall, although still quite close in value. Fig.27, ^^^M^ ^^plotted against ^^^M^ ^^, shows that SCA has the best agreement with the DNS solution, and is the only prediction within 5% of the reference solution. Note that VCA has different boundary conditions than the DNS solution, i.e. it uses a fictitious surrounding domain. Some deviation from a DNS result with periodic boundary conditions is therefore expected. A 5% deviation from the reference solution (DNS) is indicated by the blue shaded region in Figs.26-27. Machine learning on databases generated with predictive ROMs Motivation

[0212] In design optimization, one might expect there to be many calls to the material subroutine (e.g. one for each element for each load for each design iteration). Thus running ROMs for RVEs might still be time consuming, compared to a model predicting RVE stress responses given a strain state within milliseconds, or a model providing strain state given microstructure stress contours within milliseconds. To tackle this issue, we propose replacing the ROMs by neural networks trained on the RVE responses computed with SCA for micro-stress and macro-stress. These networks preserve microstructure information and are engineered to:

[0213] Predict macro (homogenized) stress or micro (local) stress given a macro-strain using an FFNN for one RVE. In this case, the FFNN plays the role of traditional material constitutive equations and homogenization (to compute the macro-stress). The present invention denotes these ℱℱ^^ℱ^^^ ^^^ ^^^^ℴand ℱℱℱ ^^ ^^.

[0214] Predict macro-strain given any micro-stress distributions using a convolutional neural network (CNN). This is the inverse of the FFNN. The input is a stress distribution within an RVE domain. The output is the macro strain loading applied to this RVE (the boundary conditions). The present invention denotes this ℱ^^ ^^ ^^.

[0215] Compute RVE damage by comparing local von Mises stress to a stress threshold. A CNN is trained to identify the onset of damage within an RVE. In this case, the CNN acts as a classifier,Attorney Docket No.: 0116936.295WO2 which identifies whether or not the applied macro-strain will cause microstructural damage. The present invention denotes this^^ ^^ ^^ ^^ ^^ ^^.Database generation for machine learning using SCA

[0216] The ROMs presented in this chapter can be used to generate a database for machine learning; in this case, we are using SCA with 8 clusters for demonstration. In general, such databases contain ^^^^training samples and ^^^^validation samples. A database with ^^^^=1,000 strain-stress pairs is computed for a nonlinear elastic material by randomly sampling 200 terminal states with four sub- loading steps each, as shown in Fig.28. Similarly, ^^^^=150 strain-stress pairs (30 final states, plus four intermediate steps) are generated for validation. A monotonically increasing load is assumed and all the final points are confined to a spherical space with a radius of 0.05. Table 7 shows the reduction in time required to build a database achieved by using SCA, versus FEM or FFT for this example case. Without a reduced order model, the generation of microstructure database takes days using FEM or FFT. The database contains ^^^^=1000 pairs of macroscopic strains and stresses εM ^^, σM ^^, and local (micro) stress σ^^(^^), where ^^ = 1,2, … , ^^^^. In this case, the present invention only considers a plane strain problem and do not consider stress in the z direction. SCA reduces the total computational time by two orders of magnitude compared to FFT and by four orders of magnitude compared to FEM. Feed forward neural network with database generated by SCA

[0217] In this case study, an FFNN is trained with data generated using SCA. Using the same microstructure as given in Fig.23, SCA computes the stress responses of the RVE (or any arbitrary RVEs) when a monotonically increasing strain is applied. Since SCA provides efficient evaluation of the stress state, it is convenient to train the FFNN on data made with SCA. The FFNN can then replace SCA by ``learning'' the stress state as a function of the strain state. This would establish a straight-forward relationship between strain and stress for near-instantaneous evaluation of RVE stress responses. Note that although a plastic material is described for the matrix material in the RVE analyzed by SCA, the FFNN described here is limited to non-linear elastic (i.e. path independent) material behavior: the present invention approximates the plastic response with a non-linear elastic one, and focus on monotonic loading. Moreover, in the following design case, the FFNN is trained to predict not only the overall RVE stress responses, but also cluster-wise, local stress responses. This is essentially the same process, but the output layer is size ^^^^× ^^^^( ^^ = 1), with one point for each stress component for each cluster. This second form replicates the non-homogenized results of SCA.Attorney Docket No.: 0116936.295WO2

[0218] Training

[0219] The training procedure for an FFNN can be reformulated as an optimization problem. We define the loss function (or cost function) as Mean Square Error (MSE) for the estimated stress and stress computed by RVE using SCA. Assuming one hidden layer, the optimization formulation is given by: find: ^^^^^^^^=2,   ^^^^^^=2,   ^^^^^^^^=3,   ^^^^^^=3min loss1 ^^ ^^ ^^ ^^( ^^=3) ^ ^^−^^=3, 2(3.17) ^^=3,^^^^=^^ ^ ^^,+^+ ^^^^^

[0220] By, , , be reduced. Note that only training data is used in this

[0221] Usually, the MSE gradually decreases with each training step. To ensure the trained neural network is general enough for all possible input states, some data points called verification data are used to monitor trends in the error. The minimization iterations will terminate before the error of the verification data starts to increase. This will ensure the neural network is able to provide certain extrapolating capability for data points that are not within the training set.

[0222] The FFNN was trained on the database described in the previous section. In this case, an FFNN with one hidden layer and 50 neurons was chosen. In the training procedure, 1,000 samples are used to train the neural network. The Levenberg--Marquardt optimization algorithm is used to reduce the MSE.

[0223] Prediction

[0224] After the training process, a fast evaluation of the stress state during a monotonic loading process was performed. The FFNN used for this predicts the macroscale stress tensor for a given macroscale strain tensor following Eq.3.18(a); similarly, Eq.3.18(b) is used to predict the local (cluster-wise) stress tensors given a macroscale strain tensor. σM ^^(FN)= σ^^= ^^^^^^, ^^= ℱℱℱ ^^ ^^(εM^^, ^^) (3.18a) ^^^^(FN)(^^)= ^^^^= ^^^^^^, ^^(^^)= ℱℱ^ℱ^ ^^^ ^^^ ^^^^ℴ( ^^^M^, ^^) (3.18b)

[0225] To demonstrate validation of the macroscale FFFN, the ^^2norm is used to measure the difference between the overall stress predicted by Eq.3.18(a) and the homogenized SCA results for each sample in the validation data set, as computed by:Attorney Docket No.: 0116936.295WO2

[0226] To validate the trained FFNN, another 30 final strain states (unknown during the training) with five load steps each (including the final state) were selected. Fig.29 shows a histogram of difference measured with the ^^2-norm, Eq.3.19, for these new ^^^^= 150 strain-stress pairs of in the validation data set. Most of the test samples have a very low ^^2-norm, which shows that the FFNN is well trained. Fig.30 shows the SCA stress data for each stress component plotted against the FFNN predictions; a perfect match has a slope of one. The associated cross-correlation statistic is one: the FFNN solution match the SCA solution perfectly. At each load step the stress predictions of the FFNN and SCA match, as shown in Fig.31.

[0227] This case study illustrates a convenient workflow that used the reduced order modeling approach to generate a rich microstructure response database for training an FFNN, which is then used for generating fast predictions of the RVE responses. Note that although the validation for the FFNN for the homogenized stress-strain relationship is given in detail here, a similar process has been used for the relationship between macroscale loading and microscale (cluster-wise, or local) stresses, as given in Eq.3.18(b).

[0228] The present invention proposes that the FFNNs shown here can be used in a design optimization process, such as topology optimization or microstructural design, where a fast and accurate material responses prediction is desired. However, note that the material is non-linear elastic and / or under monotonic loading. If plasticity and loading / unloading are considered, a different FFNN setup or a different neural network may be required. A speed comparison of running 150 samples with SCA and FFNN is given in Table 7, where the speedup for online prediction of σMis 10,000 for the FFNN over SCA. This idea will be explored further using both the FFNN predictions (Eq. 3.18(a)(b)) and convolutional neural networks in later sections. Table 7. Comparison of time required to run 150 samples using SCA and FFNN. Method Total time (s) Online speedup over SCA SCA(4+4 clusters) offline: 13.0 + online: 3 × 102---- FFNN training: 60 + prediction: 3 × 10−210,000 Convolutional neural network for boundary condition identification with database generated by SCA

[0229] Using the CNN illustrated in Fig.32, a mapping has been established between local stress distribution and applied external strain on the microstructure.Attorney Docket No.: 0116936.295WO2

[0230] Training

[0231] This CNN was trained on the same database of 1,000 SCA results as was the FFNN. For the CNN, the input is the micro-scale stress at each point (voxel) of the RVE, given by the cluster- wise results of SCA, σ(α, β). Each point in the RVE contains the three 2D stress components, like the RGB channels used for images. The output of the CNN is the macroscale strain ε^^that was applied as the loading conditions and caused the observed stresses. The training can be written as an optimization problem, as given in Eq.3.20. The equations for training a CNN with padding, convolution, and pooling layers repeated ^^^^ ^^ ^^ ^^times is given by: find: ^^^^^^^^, ^^^^^^(l=2,3... ^^^^), in FFNN (3.20) = 1,2 … ^^^^ ^^ ^^ ^^), in CNN^^− ε∗ ^^, 2^^ ))

[0232] nested operations shown in Fig.32. It is defined as ℱ η ^^ ^^ ^^(σ^^(α, β))= ^^η^^ℴℴℓ ^^ ^^ℊ^^^η^ℴ ^^ ^^( ^^η^^ ^^ ^^ ^^ ^^ ^^ℊ((σ^^(α, β)))),

[0233] which terms in thetraining problem are define as follows. The weight and bias in the FFNN are ^^^^^^^^and ^^^^^^, and ^^ζ,κ,ηξ,ψ, ^^κ,ηare the weights and biases in the convolution operations. The truth is ε∗ ^^, ^^and theε^^, ^^. The number of training samples is defined as ^^^^, and ^^^^(^^ = ^^^^)is the number ofneurons output layer. In this case, ^^^^(^^ = ^^^^)is three. The sampling is indexed by s.

[0234] The inputs to the CNN are the three stress arrays corresponding to the components of stress given by σ^^(α, β). The outputs are three strain values ε^^^^, ^^(^^ = 1,2,3). Just as for the FFNN, the mean squared error (MSE) is reduced gradually step by one of several optimization algorithms.

[0235] Prediction

[0236] Just as with the FFNN, the present invention can define the function ℱ^^ ^^ ^^(σ^^(α, β)) that describes the operation performed by the trained CNN, in this case ℱ^^ ^^ ^^(σ^^(α, β)) = εM ^^. (3.21)Attorney Docket No.: 0116936.295WO2

[0237] In Fig.3.13, a histogram for the error between predicted ε^^and the validation data set (again, the ^^^^= 150 samples generated previously) is computed using Eq.3.22. In Fig.34, the correlations between the CNN prediction of applied strain and the reference solution of the three strain components is provided, using the same validation data sets as mentioned before. These show that the CNN can effectively map the stress contour to the applied external strain. Such a map may play an important role in linking microstructure information with macroscale information, e.g. connecting microstructure failure strength to a macroscale strain state. This will assist the inverse design problem where the optimum loading state is inferred using local information. ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^^^^^^^ ^^=^^ ^^ ^^(σ^^(α, β)) − ε^^, ^^( ^^ ^^ ^^), ^^ = 1,2, … , ^^^^(3.22)generated by SCA

[0238] The application of a CNN to material microstructure predictions is not limited to the sample problem shown above. Another example use of a CNN is as a microstructure classifier; this is similar to its common application in image classification. By using a microstructure mesh and an applied strain on the microstructure as the input, a CNN can be trained to predict whether the microstructure will become damaged.

[0239] Training

[0240] The training procedure for the classification CNN is given by find: ^^^^^^^^, ^^^^^^(l=2,3... ^^^^), in FFNN (3.23) ^^^^, ^^, ^^^^, ^^, ^^^^, ^^( ^^ = 1,2 … ^^^^ ^^ ^^ ^^ ^^ ^^),( ^^ = 1,2 … ^^^^ ^^ ^^ ^^), in CNN cross entropy=1∑^^^^^(−( ^^∗ ^^log( ^^^^) + (1 − ^^∗ ^^) log(1 − ^^^^)))^^^^ ^=1where: ^^^^= ℱℱℱ ^^ ^^(ℱ^^ℓ ^^ ^^ ^^ℯ ^^(ℱ^^^^^^^^^^ℴ ^^ ^^(… ℱ^2^ ^^ ^^(ℱ^1^ ^^ ^^( ^^^^( ^^, ^^))) … ))) nested )))

[0241] The output of damaged, 1 for damaged. Since thein terms of the cross entropy between the truth value ^^∗ ^^and the predicted value ^^^^. Cross entropy, or log loss, is widely used to measure the performance of a classification model. The CNN in this section is trained on the data extracted from the database used for the previous NNs. The database for training the classifier consists of pairs of micro stress distributions and damage indicators. In this case, a critical- stress-based damage criterion is used to decide whether an RVE is damaged or not: if any von MisesAttorney Docket No.: 0116936.295WO2 stress σ(^^)exceeds a critical von Mises stress σ∗the RVE is considered damaged.

[0242] Prediction

[0243] Once trained, the CNN can predict whether the applied loading will result in microstructure damage without carrying out the full RVE simulation: (3.24)

[0244] In next section, such a CNN will be used in a microstructure-based topology optimization example to illustrate the effect of a microstructural damage constraint on the optimized structure. Microstructure-based, multiscale topology optimization using neural networks

[0245] In certain embodiments, the present invention illustrates how FFNNs and CNNs might be used in topology optimization to achieve microstructure-based design. This differentiates the current approach from classical topology optimization which typically uses simple constitutive relationships. The present invention proposes to compress the RVE response database into an FFNN for forward prediction of RVE stress responses, where it will act similarly to a traditional homogenized constitutive model. However, because a database of RVE responses is used, no functional form of the homogenization is required. The RVE microstructure damage responses is represented with a trained FFNN+CNN; this introduces microstructure damage, linked with the applied strain state of the RVE. By using well-trained FFNNs and CNNs, two different optimization problems are defined as below:

[0246] Topology optimization with a material constitutive law extracted from an FFNN trained on the stress-strain relationship of a given microstructure. In this case, a non-linear material behavior during topology optimization is used to achieve a design that is durable under extreme loading conditions where the material response enters the non-linear region.

[0247] Topology optimization with constraints defined by the FFNN+CNN framework to identify microstructure damage and thereby design durable (damage aware) structures. In this case, the microstructure damage acts as an extra constraint to the topology optimization formulation to achieve a design that alleviates or avoids possible local microstructure damage.

[0248] In these two example problems, the design zone is described with a 60 × 30 mesh of rectangular, linear elements. Each element is a 1 ^^ ^^ × 1 ^^ ^^ square. The elastic material properties are from the homogenized SCA results: E=200 MPa and ν = 0.27. An approximated non-linear elastic material response is extracted from the RVE simulation. As described in the following sections. Note that while plasticity is used, our FFNN is only valid for a non-linear elastic approximation of theAttorney Docket No.: 0116936.295WO2 material response. A point load of 75 N is applied at right bottom corner. The desired volume fraction of the optimized part is set as 0.35 of the original design zone. For the second case with damage the critical stress is defined as 0.7 MPa. Topology optimization with FFNN

[0249] The formulation of microstructure sensitive topology optimization with an FFNN is shown in Eq.3.25. The objective function is defined as the overall strain energy of the structure, Φ, which is to be minimized. A subtle but important difference in the present example in this work is that the FFNN is used to replace the usual linear elastic material response with a nonlinear one. By using a nonlinear material responses database, a new avenue for data-driven material and structure design is illustrated. This is depicted graphically in Fig.35, which shows the use of an FFNN to generate microstructure-based stress-strain response within a topology design framework.(3.25)

[0250] ^^ is the body force, ^^ is the applied traction on the boundary of the design zone, and ^^ is the local displacement in the design zone. ^^∗is the admissible displacement field, and ^^^^defines the boundaries of the macrostructure (the design region). σ(^^^^)and σ(^^^^)are macro stress and strain. The density for each macro mesh is ρ( ^^^^). The desired volume of material remaining in the design zone is ^^∗and defined as 0.35, and ^^(ρ( ^^^^)) isin the design zone. The density of location ^^^^is ρ(^^^^)in the design zone, and the homogenized stress is σ(^^^^), defined as the product of ρ(^^^^)and ℱℱℱ ^^ ^^(ε(^^^^)). Here, ℱℱℱ ^^ ^^represents a trained FFNN that will generate stress predictions based on given strain input ε( ^^^^).

[0251] In this case, thethe RVE responses. Hence, σ^^is directly approximated by the FFNN. In truth, constraints of theused for optimization require aAttorney Docket No.: 0116936.295WO2 functional description of the behavior (this limitation will be addressed in future work), thus the effective von Mises RVE stress versus effective strain curve is approximated using an exponential ̅^^̅= 0.7784 ∗ ^^22.18∗̅ε̅ ^ ̅^̅− 0.8071 ∗ ^^−378∗̅ε̅ ^ ̅^̅. This hyperelastic material definition was fit to the FFNN results. This ad-hoc approach ensures stability of the optimization process by minimizing the overall strain energy using condition-based optimization.

[0252] The optimized beam structures are illustrated in Fig.36 panels (a) and (b) for linear elastic material and non-linear elastic material, respectively. Necessary results are provided in Table 8. Two structures show substantial difference in the final shape of the structure. This means the material non- linearity plays an important role in topology optimization, where a new truss structure is realized in order to ensure minimal strain energy. The result suggests the importance of considerable of microstructure-based material non-linearity into structure optimization. It may also possible to include microstructure variation, such as different particle volume fractions, into the structure. In short, FFNN provides an alternative for a data-driven microstructure-based topology optimization, where the microstructural effect can be incorporated into the process and to achieve different designs that meet the design criteria. Table 8. Results of the FFNN-based non-linear-elastic optimization problem. Note that substantial speedup is achieved while retaining the accuracy and microstructural basis of the concurrent approaches. This would provide further speed advantages in 3D linear FFNN FE-SCA FE-FE material (nonlinear) concurrent* concurrent* Initial compliance (strain energy) (N ⋅ 12.6 (6.3) 20.0 (10.0) --- --- cm) Optimized compliance (strain energy) (N 28.0 (14.0) 20.0 (10.0) --- --- ⋅ cm) Database generation + training (s) 0 313 --- --- Optimization calculation time (s) 338 472 23328 220 x 106Factor of speed-up over FE-FE --- 280255 9431 --- No. of iterations 14 18 --- --- Topology optimization with constraints defined by FFNN+CNN

[0253] Topology optimization results may have high stress concentration zones. If not treated properly, the concentration zones may cause unexpected damage and affect the function of the structure. A traditional way to address this is to add stress or strain constraints to optimization. Previous authors have focused on optimization with stress concentration and singularities. In theseAttorney Docket No.: 0116936.295WO2 previous studies, most damage criteria are only related to the macroscopic material model and do not consider micro structure mechanical behavior. In this work, an FFNN+CNN framework trained by the database generated with SCA.

[0254] In order to do this, the FFNN defined in Eq.3.18(b) that performs the operation σ(^^)= ^^ ^ ^ ^is used. This predicts local, rather than homogenized, stresses in the RVE. These localserve as input to the CNN outlined previously, which indicates whether or not damage has occurred. The optimization formulations are thus: minimize: Φ=∫^^m ^^ ^^ dΩ^^+ ∫^^m ^^ ^^d ^^^^, ∀ ^^ ∈ ^^∗^^^^: macrodomain ^^^^: coordinate in macro domain ^^: coordinate in micro domain 0255] where Φ=^^^^=1( ^^ ^^^^ ^^[0 ∑ ( )) ^^^^^^^^^^^^^^is the compliance of the overall structure, ^^ is the number (typically p=3), ^^^^is the displacement for eachelement, and ^^0is the element stiffness. The averaged strain response of the RVE is ε^^. The local stresses within the RVE are defined as σ, and are calculated with a trained FFNN ℱℱ^^ℱ^^^ ^^^ ^^^^ℴ. Other variables are the same as defined above for the FFNN-only optimization. For ^^ RVE ^^^^(^^),the output value of ℱ^^ℓ ^^ ^^ ^^ ^^ ^^ ^^^^ ^^ ^^should be 0, which represents a non-damaged state.

[0256] Since database only gives a criterion, sensitivity analysis is not preferred in this case. Thehere follows a refined optimally criteria (OC) method. The optimization program structure is based on the 99-line topology optimization code written by O. Sigmund (and thus the problem is similar, though the implementation is not identical to the FFNN-based optimization above). During the density update, for each element, three strain components will be passed to the FFNN+CNN. The FFNN+CNN will determine whether the current strain state is acceptable by assessing the local microstructural response. This is shown schematically in Fig.37. If an element has been damaged, the density of this element will be increased by applying a penaltyAttorney Docket No.: 0116936.295WO2 factor, while the densities of the rest of the elements will be decreased to satisfy the volume constraint.

[0257] Similarly to the FFNN example, Fig.38 panel (a) shows the reference case: a linear elastic optimization without a damage constraint. The final compliance is 30 N ⋅ cm. Fig.38 panel (b) shows the optimized design with a microstructure-based damage constraint. The final compliance is 31 N ⋅ cm. Notice that while the compliance is quite similar, the design under a microstructural damage constraint resembles a more conventional truss structure; the optimization has avoided sharp angles and has fewer beams that give rise to stress concentrations likely to result in microstructure-driven damage.

[0258] A summary of the design variables and important parameters related to the optimization is provided in Table 9. The simple examples above show the potential application of an FFNN+CNN database generated by clustering reduce order methods. However, the optimization is just based on an artificially defined optimality criterion. The algorithm may not be stable for all kinds of problems. In the future, we should study the sensitivity and singularity of constraints based on an FFNN+CNN database. Table 9. Results of the FFNN+CNN constraint optimization problem. Note that substantial speedup is achieved while retaining the accuracy and microstructural basis of the concurrent approaches. This would provide further speed advantages in 3D. linear FFNN + FE-SCA FE-FE material CNN concurrent* concurrent* Initial compliance (strain energy) (N ⋅ 295 (148) 295 (148) --- --- cm) Optimized compliance (strain energy) (N 30.0 (15.0) 31.0 (15.5) --- --- ⋅ cm) Database generation + training (s) 0 512 --- --- Optimization calculation time (s) 12.6 14.5 69674 660 x 106Factor of speed-up over FE-FE --- 45 x 1069473 --- No. of iterations 53 61 --- ---

[0259] Two challenges with current approaches to machine learning methods in the mechanical science of materials are: (1) the database generation time and effort are extensive, and (2) the application of machine learning is not well developed or understood by the community. This chapter covers several different topics related to these challenges:

[0260] The present invention has outlined, related, and compared three different clustering-Attorney Docket No.: 0116936.295WO2 discretization methods (SCA, VCA, and FCA) that rely on unsupervised learning for order reduction and the solution of mechanistic governing equations for prediction.

[0261] One of these methods, SCA, was used to develop an example material behavior database suitable for training neural networks. This approach to database development substantially reduces the effort required to acquire the information upon which neural networks may be trained.

[0262] The basic operations, and how these combinations to make predictions of mechanical responses, in an FFNN were outlined. This includes the role of weights, biases and activation functions as well as the description of the training stage of the neural network as a minimization problem using notation common within the mechanical sciences.

[0263] A similar description of convolutional neural networks was developed for two different possible applications: (1) to solve inverse problems where the boundary conditions need to be identified from a known stress distribution and (2) as a classifier to identify if damage will occur within a microstructure given a known stress distribution.

[0264] Two microstructure-sensitive topology optimizations are demonstrated. In the first case, the material response at the microscale derived from the FFNN results, and used to perform design against a load that causes the material to behave in a non-linear elastic way. In the second case, a material damage constrain is added to the optimization, where the CNN is used to identify if damage has occur on the microscale and penalize the design accordingly. Example 3 Mechanistic nested topological optimization design theory for structure, nested lattice, and material selection

[0265] In certain embodiments, the present invention extends the developed C-HiDeNN-TD to solve a structure, nested lattice, and nested material system. Nested lattice designs are ubiquitous in nature that not only serve the structural purpose but also provide multifunctionality. However, realizing such design in engineering applications is challenging due to the high computational effort required for such nested lattice design.

[0266] In one embodiment, the present invention proposes a method for N-scale concurrent nested topology structural design of which for each lattice structure n=1,..., N, the structure can be subdivided into d=1,.., D regions with functional design for applications. The proposed concurrent nested lattice design theory and its implementation through the C-HiDeNN-TD approaches have enormous potential in engineering product development. This will also advance the field of manufacturing science as realizing those high-resolution nested lattice designs requires significantAttorney Docket No.: 0116936.295WO2 advancement in manufacturing. Mechanistic nested topological optimization design theory

[0267] The scales are defined below. The notation table used is shown in Table 10.

[0268] To have a clear understanding of the problem, the present invention first discloses the general framework. As shown in Fig.39, a general design consists of nested scales with multimaterial. All the lattices and sub-lattices design domain is expressed by a general notation[^^][ ^^] ^^−1, ^^ where n is the number of lattice design scales, d is the design index of each n scale, andmchoice or design. To keep track of the nested design, we introduce a index j which tracks the design index of scale n -1. n=0 is considered the starting part-scale design domain without any lattice design and n =1 is the same part-scale design domain with nested lattice design obtained in each optimization iteration. And m=0 denotes the lattice design, and the materials are selected in the lower scales for that particular design. All the notations are summarized in the Table 10. Based on the application of the loading and boundary conditions, the initial design domain can be divided into d design domain in scale n =1. In the lower scales of n, the design pattern can repeat itself (similar to self-similar fractal design) or show non-repeating lattice design. Our objective is to for a given number of scales n, optimize the design domain using nested lattice design d with the materials selection m. The present invention starts with defining each scales, design variables, equality and inequality constraints, objective functions and optimization formulation to achieve such nested lattice design to solve them concurrently. To make the equations concise and easy to be understood, we only show [d] instead of [ ^^]^^−1, ^^. In the numerical implementation procedure, one should use a variable to storage the relation of different [d]s.Table 10. Notation used for C-HiDeNN-TDN x SCAM theory[^][ ^^][ ^^] ^^−1, ^^[ General notation^^]n Index of scales ^^ = 0,1,2 … , ^^ + 1 m Index of material m=0,1,2..., M. m=0 is the lattice design d Index of designed unit cell j Design index of n -1 scale ρ[[^^],[ ^^]^^] ( ^^, ^^, ^^)Design variables for d-th structure / lattice at scale n with material m ^^[ ^^],[ ^^]PK1 stress for d-th design region at scale n ^^[ ^^],[ ^^]Deformation gradient for d-th design region at scale n ^^[ ^^],[ ^^]Displacement for d-th design region at scale n ^^[ ^^],[ ^^]Coordinates for d-th design region at scale n ∂Ω[^^^^],[ ^^]Essential boundary condition on scale nAttorney Docket No.: 0116936.295WO2 ^^ ^^[^^^^],[ ^^]Natural boundary condition on scale n Δ^^[[^^]^^]Increment of deformation gradient of material m at scale n Γ[[^^] [ ^^] ′[ ^^]^^] ( ^^[ ^^] , ^^[ ^^] )Green function represents the strain at ^^[[^^]^^]contributed by an external stress at ^^[′[ ^^]^^]^^[0,[( ^^−1)]^^] ( ^^′)The homogeneous isotropic linear elastic reference material property Equilibrium and scale coupling equations

[0269] In one embodiment, the present invention first defines the equilibrium equations and coupling equations for each scale. To give a general theory, one use finite deformation formulations with PK1 stress ^^ and deformation gradient ^^.

[0270] For scale n=1, the equilibrium equations are defined as:

[0271] At scale n=1: ∇ ⋅ ^^[ ^^=1]( ^^[ ^^=1], ρ[[^^=1,2… ^^],[ ^^]^^=0,1,2… ^^]) + ^^[ ^^=1]= 0,  ∀ ^^[ ^^=1]∈ Ω[ ^^=1](4.1) and firstPK stress are defined as ^^[ ^^=1]and ^^[ ^^=1]. The coordinate system is defined as ^^[ ^^=1]. The essential boundary condition is defined as ∂Ω[^^^^=1]and the natural boundary is defined[ ^^=1]^^. It is obvious that for scale one, all the design^^],[ ^^]^^]in the current scale and lower contributed tothe equilibrium equations. In the present invention assumes each voxel in scalen=1 has a homogenized its corresponding lower scales. For nonlinear material property problem, it is a stress response correspondence to the strain. For linear material property problem, it is a set of homogenized material parameters. The present invention will present this after defining the rest equilibrium equations in the lower scales.

[0273] As for each design region in lower lattice scales, the equilibrium equations are similar, we define the equations based on a general design region[^^],[^^].

[0274] The equilibrium equations are defined as:

[0275] At scale n, design d:Attorney Docket No.: 0116936.295WO2^^[ ^^],[ ^^]= ^^,, ^^,is anti-periodic.

[0276] The first PK stress and the deformation gradient is similar to the scale n =1. The design variables that influence the current equilibrium equations at scale n for design d is defined as [^^= ^^, ^^+1,… ^^],[ ^^]. There is no body force term as the lattice only takes the boundary conditions ashomogenization theory. The essential boundary condition is defined as periodic and the natural boundary is defined as anti-periodic.

[0277] For the material set m=1,2...M, there are various homogenization theory the describing the material response to the external loading. The present invention uses the Lippmann–schwinger homogenization in finite deformation to describe the equilibrium equations in the the material level. Assume the material m is on scale n and the upper scale is n -1. The equitions are:

[0278] At scale n, material m: Δ ^^[[^^]^^]( ^^[[^^]^^]) = Δ ^^[^^^^−]1−∫Ω[^^] [^^]Γ[[^^]^^]( ^^[[^^]^^], ^^[′[ ^^]^^]) : (Δ ^^[[^^]^^]( ^^′) − ^^[0,[( ^^−1)]^^]( ^^′): Δ ^^[[^^]^^]( ^^′)) ^^ ^^′green function is defined as Γ[ ^^]( ^^[ ^^], ^^′[ ^^]), which s[ ^^][^^] [ ^^] [ ^^]train at ^^[ ^^]contributed by an external stress at ^^[′^^^]^. ^^[0,[(^^]the homogeneous isotropic linear elastic referencematerial

[0280] equations for each scale including the material micro structure scale, the present invention defines the scale coupling equations. Please note the scale coupling equations only applies to a scenario where for a specific voxel, its lower scale is a lattice. For the voxel with a material lower scale, one can directly use the Lippmann–schwinger equation. For the voxel with a lattice at lower scale, the present invention use Hill’s lemma to define the relations. Assume a voxel at design [ ^^], [ ^^] with location ^^∗[ ^^],[ ^^], and the lower scale is a lattice with index[^^ + 1], [ ^^]. The coupling equations are defined as:^^[ ^^],[ ^^]( ^^∗[ ^^],[ ^^]) =1 [ ^^+1],[ ^^] [ ^^+1],[ ^^] [ ] [ ]|Ω[ ^^],[ ^^]|∫Ω[ ^^+1],[ ^^]^^ ( ^^ ) ^^ ^^^^+1 , ^^. (4.4)Attorney Docket No.: 0116936.295WO2

[0281] The deformation gradient at lower scale is defined as ^^[ ^^+1],[ ^^], whereas the deformation gradient at upper scale on location ^^∗[ ^^],[ ^^]is defined as ^^[ ^^],[ ^^]( ^^∗[ ^^],[ ^^]). It can be calculated by the integral and average over the whole domain in the lower scale. Optimization formulation and solution scheme

[0282] As defined previously, the design variable is defined[ ^^][ ^^]. For a general purpose, the objective function can be defined as a weighted sum of variousfunctions: ^^^^^^^^^^^[ ^^],[ ^^]

[0283] Where ^^^^,^^,^^formulation ℱ[ ^^],[ ^^]refers to the design objective for the design region d at scale n, on material m. For thisonly have material selection. There is no difference for objective functions of different material. Therefore sub-objective function can be expressed as ℱ[ ^^],[ ^^]. The weight factor ^^[ ^^]also disappears. It should be noted that there is no requirement tofunction for every scale with every design unit cell. For some cases, for example, concurrent design to find the optimal compliance at part scale, the objective function is only ℱ[ ^^=1],[ ^^=1]. The design of lower scale lattice and the selection of material are guided by this single objective function at part scale.

[0284] Apart from the objective function, the design usually has equality constraint and inequality constraints.

[0285] The equality constraint can be written as: ^^[ ][ ] [ ][[^^ , ^^ ^^ , ^^] [ ^^],[ ^^] [ ^^],[ ]^^]= ^^[ ^^]+ ^^[ ^^]+ ℎ[ ^^]^^, (4.6)

[0286] where ^^[ ^^],[ ^^]are the equilibrium conditions. The volume constraints are ^^[[^^],[ ^^]^^]. All other equality constraints are defined as ℎ[ ][^^ ,[ ^^]^^].

[0287] The inequality defined as: ^^[[^^],[ ^^] [, ∗[ ^^],[ ^^] [ ^^],[ ] [ ][ ]^^]= ^^[ ^^]ρ[ ^^]+ ^^[^^^^]+ ^^[^^ , ^^^^], (4.7)

[0288] where the constraints are [∗[ ][ ]^^ℎ ^^^^ , ^^^^]. The manufacturing constraints are ^^[,^^]. Other inequality constraints are [][ ]^^^^ , ^^Finally, the design formulations can beform:Attorney Docket No.: 0116936.295WO2 Find ρ[^^^^],[ ^^], λ[[^^],[ ^^]^^], κ[[^^],[ ^^]^^], ^^[[^^],[ ^^]^^]such that (4.8) Satisfy^^]Minimize: ℒ (ρ[ ,[ ^^] [ ^^],[ ^^][] , ^^[ ^^] )^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ )).

[0290] The Lagrange parameters are,. are used e inequality constraints ^^[ ]for th^^ ,[ ^^]. The above can using existing optimization methods such as themethod (OC) or the Method of Moving Asymptotes (MMA).

[0291] After the introduction of the general theory, the present invention discloses the solution scheme. In the implementation part, the present invention makes several simplifications: 1) only focus on small strain problems; 2) only consider material selection instead of material design.

[0292] First, the solution of part scale using C-HiDeNN-TD is the same as the theory introduced in Example 1. The solution of the lattice scale is the same as the theory introduced in Example 2.

[0293] the present invention discloses the sensitivity derivation for the system compliance. To start with, for scale n=1, the present invention gives the sensitivity respect to the part scale design ρ[ ^^=1],[ ^^=1]. As the present invention only considers material selection, the notation m is not considered in the sensitivity. After the concurrent design, the present invention runs several steps for material selection to choose the optimal material candidate that can decrease the objective function. For scale 1, the sensitivity can be written as: At scale n=1: ∂ ^^ ∂ρ[ ^^=1],[ ^^=1](4.9) ]. as ρ[ ^^=1],. The homogenized material property for each point at scale n=1 is ^^ (ρ[[^^=2,3.. ^^],[ ^^]^^]) and is related to all the lower scales and material. The strain ε[ ^^=1],[ ^^=1]in scale n=1.

[0295] For the lower scale of each lattice, onea general form: At scale d, design d:Attorney Docket No.: 0116936.295WO2 ∂ ^^ ∂ρ[ ^^],[ ^^](4.10) ]

[0296] The sensitivity between ^^^^=1and design variable ρ[ ^^],[ ^^]satisfies the chain rule and one has: ∂ ^^[ ^^−1],[ ^^] 1 ^^ ^^ ^^−1^^],defined as ^^[ ^^+1],[ ^^]. It can be calculated by the homogenization of the lower scale. The test strain is defined as ε( ^^0[ ^^],[ ^^]) and the response strain refers to ε( ^^[ ^^],[ ^^]). In some embodiments, the presenthas few design regions in each scale n, and most of the n scale nested lattices are under similar loading, the sitivity can be simplify by the assumption∂ ^ [^^−1],[^^] [ ^^=1]sen^ ∂ ^^∂ρ[ ^^],[ ^^]=∂ρ[ ^^],[ ^^]. It should be noted that even under this assumption, the homogenization partmaintained nested to ensure a concurrent nested design. After the concurrent optimization of the nested lattice, the material selection can be added to further optimize the objective function.

[0298] The design framework is shown above. It has three parts corresponding to the three level of structure, lattice, and material. The part level and lattice level are similar to the 2-scale concurrent design with a difference in looping over d and n. The material level is used to characterize the material properties of the entire design. For the material selection part, it can also be used to further increase the structural performance. Problem 1: bamboo design

[0299] In one embodiment, the present invention uses the example of bamboo design to illustrate the concept of nested lattice design. Bamboo is a widely observed natural material that exhibits a nested structure. Due to its distinctive composition, bamboo possesses a lightweight design while maintaining significant strength. Another example of nested structures can be found in the Willis Tower, which held the title of the world's tallest building from 1973 to 1998. The tower's tube-like design, reminiscent of bamboo, contributed to its impressive height.

[0300] The design is shown in Fig.41. The design region has several scales: (1) for the part scale n=1, the design is a 250 × 25 × 25 voxel region, and it has a uniform axial force on the top and aAttorney Docket No.: 0116936.295WO2 fixed region on the bottom; (2) for the nested lattice regions, each lattice is a 25 × 25 × 25 box. The design objective is to design a bamboo structure that can take most of the axial force with nested lattice structures. The design variables are part scale structure with 6 scales of nested lattice. The volume constraint of the part scale is set to 70%.

[0301] To investigative the effectiveness of the volume constraints on the design, the present invention discloses three cases: (1) case 1: all unit cells have a volume fraction constraint of 0.7; (2) case 2: all unit cells have a volume fraction constraint of 0.5; (3) case 3: all unit cells have a different volume fraction constraint shown in Table 11. Table 11. Volume fraction of each scale for case 3 Scale n = 1 n = 2 n = 3 n = 4 n = 5 n = 6 n = 7 Volume Fraction 50% 73.76% 72.48% 87.20% 78.56% 72.48% 83.36%

[0302] The nested design of case is shown in Fig.42. To enhance virtualization purposes, the present invention, in one embodiment, chooses to display only the cross-section of the design. In doing so, it is observed that at the part scale, the design automatically generates a rounded shape that closely resembles the actual structure of bamboo. This observation suggests that naturally occurring optimal designs may exhibit similar patterns to those derived through numerical optimization.

[0303] The design of case 2 is shown in Fig.43. After lowering the volume fraction, the design will becomes a rectangle frame with a rectangle hole. This kind of pattern has been widely seen in a truss structures. Similar to case 1, every lattice level has the similar design pattern.

[0304] The results of case 1 and case 2 shows if the constraint of each scale are same, it has the trend to become a fractal design. In case 3, the present invention defines different volume fractions at each scale. The result is shown in Fig.44. C-HiDeNN-TDNx SCAMhas the capability to go beyond fractal design by enabling different design constraints at different scales.

[0305] To understand the influence of the number of nested scales, the present invention adds one additional numerical study.6 cases are defined: 1) 2-scale concurrent design (2) 3-scle concurrent design (6) 7-scale concurrent design. All the cases are constraint with a same volume fraction. The result is shown in Fig.45.

[0306] From the plot, with the increasing number of the concurrent scales, the axial compliance, bending compliance, and the torsion compliance can be decreased. The first few increments of the scales highly reduce the objective function. This is because the first few increments mostly optimize on a large size compared with the lower scales. From the bamboo example, it can be seen that the theory has the capability of enabling a nested n-scale design. By applying different constraints, it canAttorney Docket No.: 0116936.295WO2 go beyond the fractal design. The idea of different design region d and different material m will be discussed in the next example. Problem 2: car bumper design

[0307] In certain embodiments, the present invention discloses using a car bumper problem to demonstrate how to use the nested concurrent topology theory with advanced materials for a car bumper design. To start with, the car bumper is the front part of a vehicle. It has two main functionalities: 1) shock absorption, 2) maintain a good stiffness with a lower compliance. As lower compliance and shock absorption are usually contradictory design. In this problem, the present invention divides the car bumper into different regions.

[0308] For region 1, the design target is the maximum energy absorption. On scale n=1, one only assumes it is a solid material. On scale n=2, one designs the structure to have a negative Poisson's ratio. On the material level, the present invention chooses the material with maximum toughness. For region 2 and 3, on scale n=1, the present invention has the part scale design. On scale n=2 and scale n=3, the present invention defines both regions have a nested lattice. The volume constraint for each unit cell is defined as 50%.

[0309] The lattice design of region 1 is shown in Fig 47. For scale n=2, d=1, it is a negative Poisson's ratio design. There is no scale n=3 for region 1. For scale n=4, the present invention chooses the material mirco-structure with the maximum toughness. It should be noted the material micro-structures are just for illustration purpose. For region 2, the design is a three-scale nested lattice design. The present invention chooses the material micro-structure with the maximum stiffness.

[0310] The design of region 3 is shown in Fig.48. It features nested structures that share similarities with the bamboo problem. However, there is a key difference: unlike the symmetric structure of bamboo, the designed pillar region in this case is non-symmetric. Therefore it can be seen for n=3,d=3, the pillar region has non-symmetric holes. Most of the holes are closer to the impact region. The region near the impact would have a lower bending moment. As a result, more material tends to fill the distant region where the highest stress values are anticipated. This observation aligns with real bamboo structures, where the outer regions are composed of stronger materials and lattices.

[0311] In these embodiments, the present invention expands upon the previously developed C- HiDeNN-TD method, introducing a nested topology design theory with advanced material systems.

[0312] The foregoing description of the exemplary embodiments of the invention has been presented only for the purposes of illustration and description and is not intended to be exhaustive orAttorney Docket No.: 0116936.295WO2 to limit the invention to the precise forms disclosed. Many modifications and variations are possible in light of the above teaching.

[0313] The embodiments were chosen and described in order to explain the principles of the invention and their practical application so as to enable others skilled in the art to utilize the invention and various embodiments and with various modifications as are suited to the particular use contemplated. Alternative embodiments will become apparent to those skilled in the art to which the invention pertains without departing from its spirit and scope. Accordingly, the scope of the invention is defined by the appended claims rather than the foregoing description and the exemplary embodiments described therein.

[0314] Some references, which may include patents, patent applications and various publications, are cited and discussed in the description of this invention. The citation and / or discussion of such references is provided merely to clarify the description of the invention and is not an admission that any such reference is “prior art” to the invention described herein. All references cited and discussed in this specification are incorporated herein by reference in their entireties and to the same extent as if each reference was individually incorporated by reference.

Claims

Attorney Docket No.: 0116936.295WO2 CLAIMS What is claimed is:

1. A mechanistic nested topological design optimization method for using Convolution- Hierarchical Deep-learning Neural Network Tensor Decomposition (C-HiDeNN-TD), comprising: (1) obtaining at least one design variable and at least one design requirement; (2) decomposing the at least one design variable into at least one equation in each of at least two directions; (3) solving the at least one design variable in each of the at least two directions using a C- HiDeNN-TD solver to obtain at least one part design; and (4) outputting an optimal design based on the at least one part design when the optimal design meets the least one design requirement.

2. The method according to claim 1, wherein between step (3) and step (4) the method further comprises a step of updating the at least one design variable based on the at least one part design.

3. The method according to claim 1, wherein when the optimal design fails to meet the least one design requirement, repeat the steps (2) to (3) until the optimal design meets the least one design requirement.

4. The method according to claim 2, wherein the at least one equation comprises at least one equilibrium equation.

5. The method according to claim 4, wherein the at least one equation further comprises at least one scale coupling equation.

6. The method according to claim 3, wherein step (3) comprises a lattice level design process and a material level design process.

7. The method according to claim 6, wherein the lattice level design process further comprises (i) setting up an initial value for the at least one design variable;Attorney Docket No.: 0116936.295WO2 (ii) providing a mode q; (iii) solving the at least one design variable in each of the at least two directions for the mode q; (iv) obtaining the at least one part design when the at least one design variable meet a convergence requirement.

8. The method according to claim 7 further comprises providing a new mode q and repeating the steps (iii) and (iv) when the at least one design variable meet a convergence requirement.

9. The method according to claim 8, wherein the material level design process comprises selecting a material, and obtaining a homogenized material property based on the selected material.

10. The method according to claim 9, wherein the step (3) further comprises a step of analyzing a sensitivity deviation of the at least one part design.

11. A method for graphics processing units (GPU) programming of a mechanistic nested topological design optimization using Convolution-Hierarchical Deep-learning Neural Network Tensor Decomposition (C-HiDeNN-TD), comprising: (1) obtaining at least one design variable and at least one design requirement; (2) decomposing the at least one design variable into at least one equation in each of at least two directions; (3) solving the at least one design variable in each of the at least two directions using a C- HiDeNN-TD solver to obtain at least one part design; and (4) outputting an optimal design based on the at least one part design when the optimal design meets the least one design requirement.

12. The method according to claim 11, wherein between step (3) and step (4) the method further comprises a step of updating the at least one design variable based on the at least one part design.Attorney Docket No.: 0116936.295WO2 13. The method according to claim 11, wherein when the optimal design fails to meet the least one design requirement, repeat the steps (2) to (3) until the optimal design meets the least one design requirement.

14. The method according to claim 12, wherein the at least one equation comprises at least one equilibrium equation.

15. The method according to claim 14, wherein the at least one equation further comprises at least one scale coupling equation.

16. The method according to claim 13, wherein step (3) comprises a lattice level design process and a material level design process.

17. The method according to claim 16, wherein the lattice level design process further comprises (i) setting up an initial value for the at least one design variable; (ii) providing a mode q; (iii) solving the at least one design variable in each of the at least two directions for the mode q; (iv) obtaining the at least one part design when the at least one design variable meet a convergence requirement.

18. The method according to claim 17 further comprises providing a new mode q and repeating the steps (iii) and (iv) when the at least one design variable meet a convergence requirement.

19. The method according to claim 18, wherein the material level design process comprises selecting a material, and obtaining a homogenized material property based on the selected material.

20. The method according to claim 19, wherein the step (4) further comprises a step of analyzing a sensitivity deviation of the at least one part design.

Images