SYSTEM AND METHOD POR FINITE ELEMENT BASED

TOPOLOGY OPTIMIZATION

The present application claims priority to United States Provisional Patent Application No. 60/751,500 filed December 19, 2005.

BACKGROUND

The present invention relates generally to design optimization methodologies, and relates in particular to design methodologies for developing designs of structures to support a mass with a minimum amount of material. The design of minimal material (and therefore light-weight) structures is a subject of central importance in the development of a very wide range of products from fighter aircraft to wheel chairs. The design of absolute minimum weight structures requires the identification of an optimal topology of the structural members. For a limited number of loading configurations, analytic methods have been applied to identify optimal topologies. For more general loading configurations however, the optimal topology cannot be determined analytically and numerical procedures become necessary. In recent years, researchers have developed topology optimization schemes that utilize finite element analyses (FEA) to evaluate candidate topologies. Through iterative finite element analyses, these schemes have been shown to converge to known optimal topologies. Various commercially available computer software programs exist for performing FEA.

In particular, the homogenization method has been proposed for generating optimal topologies based on minimization of structural compliance for a given design volume. Homogenization methods include the classic homogenization methods, artificial material models or Solid Isotropic Material Penalization (SIMP) methods. In an alternate approach, referred to as Evolutionary Structural Optimization (ESO), inefficient material in the design domain is removed iteratively. The ESO scheme is a relatively simple algorithm that may be readily implemented in commercial finite element software. Studies have revealed however, that while the ESO method does provide reasonable results in some simple test cases, in more complicated problems it has been shown to fail to provide the optimal topology. The homogenization method is more sophisticated than the ESO but, due to the use of a unique variational principal and the addition of specialized constraint equations, it is not amenible to being incorporated in standard finite element software. Topology optimization methods incorporating features of both these methods are called hybrid methods. In the topology optimization methods described above, the optimal topologies are generated based on design variables computed over an element, which leads to numerical instabilities such as checkerboarding and mesh dependency, and node based approaches have been proposed to address this where material density distribution is defined as a nodal variable. The tedious computations involved in determination of effective Young's modulus and the need for customized finite element codes, pose significant limitations to homogenization methods. Evolutionary optimization methods, where least stressed elements were removed successively from the design domain, have gained popularity as these methods are simpler to implement and have been shown to generate good

topologies utilizing standard finite element codes. Other conventional techniques include a constant weight folly stressed (CFS) method in which the total volume of material is held constant. An advantage of this method is that the density levels of elements are allowed to increase or decrease. Such methods, however, are also not immune to numerical instabilities and other limitations.

There is a need, therefore, for a more flexible, efficient and economical methodology for performing topology optimization analyses of structures.

SUMMARY The invention provides a method for providing an optimal topology for a structure based on a set of design criteria including at least one support point and at least one force to be applied to the structure. The method includes the steps of identifying a plurality of nodes within a structure design domain, and assigning an initial density value to the plurality of nodes. The method also includes the steps of conducting a finite element ^■ analysis on the nodes, determining a stress intensity value for each node, ranking the nodes by relative stress intensity values, and adjusting the density value for each node. The method also includes the step of repeating the steps of conducting a finite element analysis on the nodes, determining the stress intensity value for each node, ranking the nodes by relative stress intensity values, and adjusting the density value for each node until a termination criteria is realized, thereby providing an optimal topology. In accordance with an embodiment, the density values of the nodes are reduced based on the ranking of the nodes, and in accordance with further embodiments, the density values of

the nodes are adjusted based on the ranking of the nodes and based on a set of known density distributions.

BRIEF DESCRIPTION OF THE DRAWINGS The following description may be further understood with reference to the accompanying drawings in which:

Figure 1 shows an illustrative view of the process steps in a topology optimization method in accordance with an embodiment of the invention;

Figure 2 shows an illustrative diagrammatic view of a design domain of a system to be optimized in accordance with an embodiment of the invention;

Figure 3 shows an illustrative diagrammatic view of a minimum weight topology of the domain shown in Figure 2 in accordance with an embodiment of the invention;

Figure 4 shows an illustrative diagrammatic view of near optimal topology of the domain shown in Figure 2 in accordance with an embodiment of the invention; Figure 5 shows an optimization history for a conventional finite element analysis for different numbers of iterations;

Figure 6 shows an optimization history for an evolutionary structural optimization in accordance with an embodiment of the invention for a finite number of iterations;

Figure 7 shows an illustrative diagrammatic view of another design domain of a system to be optimized in accordance with an embodiment of the invention;

Figure 8 shows an illustrative diagrammatic view of a minimum weight topology of the domain shown in Figure 7 in accordance with an embodiment of the invention;

Figure 9 shows an illustrative diagrammatic view of the analytically optimal topology of the domain shown in Figure 7 in accordance with an embodiment of the invention;

Figure 10 shows an illustrative diagrammatic view of another design domain of a system to be optimized in accordance with an embodiment of the invention;

Figure 11 shows an illustrative diagrammatic view of a minimum weight topology of the domain shown in Figure 10 in accordance with an embodiment of the invention;

Figure 12 shows an illustrative diagrammatic view of the analytically optimal topology of the domain shown in Figure 10 in accordance with an embodiment of the invention;

Figure 13 shows an illustrative diagrammatic view of another design domain of a system to be optimized in accordance with an embodiment of the invention;

Figure 14 shows an illustrative diagrammatic view of a minimum weight topology of the domain shown in Figure 13 in accordance with an embodiment of the invention; Figure 15 shows an illustrative diagrammatic view of the analytically optimal topology of the domain shown in Figure 13 in accordance with an embodiment of the invention;

Figure 16 shows an illustrative diagrammatic view of another design domain of a system to be optimized in accordance with an embodiment of the invention; Figure 17 shows an illustrative diagrammatic view of a minimum weight topology of the domain shown in Figure 16 in accordance with an embodiment of the invention;

Figure 18 shows an illustrative diagrammatic view of the analytically optimal topology of the domain shown in Figure 16 in accordance with an embodiment of the invention;

Figure 19 shows an illustrative view of the process steps in a topology optimization method in accordance with a further embodiment of the invention;

Figure 20 shows an illustrative graphical representation of a transition process from initial to final probability distribution in a system in accordance with an embodiment of the invention;

Figure 21 shows an illustrative graphical representation of a transition process from initial to final cumulative probability distribution in a system in accordance with an embodiment of the invention;

Figure 22 shows an alternative illustrative graphical representation of a transition process from initial to final probability distribution in a system in accordance with an embodiment of the invention; Figure 23 shows an alternative illustrative graphical representation of a transition process from initial to final cumulative probability distribution in a system in accordance with an embodiment of the invention;

Figure 24 shows an illustrative. diagrammatic view of a minimum weight topology of a domain as shown in Figure 2 for a dual strength material in accordance with an embodiment of the invention;

Figure 25 shows an illustrative diagrammatic view of the analytically optimal topology of the domain shown in Figure 2 for a dual strength material in accordance with an embodiment of the invention;

Figure 26 shows an illustrative diagrammatic view of a minimum weight topology of a domain similar to that shown in Figure 13 for a dual strength material in accordance with an embodiment of the invention;

Figure 27 shows an illustrative diagrammatic view of the analytically optimal topology of a domain similar to that shown in Figure 13 for a dual strength material in accordance with an embodiment of the invention;

Figure 28 shows an illustrative diagrammatic view of a test case expected topology for a methodology of the invention for use on a three-dimensional structure;

Figure 29 shows an illustrative diagrammatic view of a predicted topology using an adaptive topology optimization program in accordance with an embodiment of the invention.

Figure 30 shows illustrative diagrammatic views of classifications used in the generation of top and bottom triangulated surfaces in accordance with an embodiment of the invention; Figure 31 shows illustrative diagrammatic views of classifications used in the generation of lateral triangulated surfaces in accordance with an embodiment of the invention;

Figure 32 shows an illustrative diagrammatic view of a minimum weight topology of the domain shown in Figure 16 for use as input to the geometry extraction algorithm in accordance with an embodiment of the invention; and

Figure 33 shows an illustrative diagrammatic view of a stereolithography (STL) model generated from the image in Figure 32 in accordance with an embodiment of the invention.

The drawings are for illustrative purposes only.

DETAILED DESCRIPTION OF THE ILLUSTRATED EMBODIMENTS

The invention provides topology optimization algorithms that may be implemented as supplemental or front end routines to standard commercial programs, allowing the designer to use finite element codes that are familiar, efficient, and currently available from many different sources. In an example, a methodolgy of the invention may be implemented with the finite element analysis program, Abaqus, sold by Abaqus, Inc. of Providence, Rhode Island, through the use of a set of FORTRAN programming language user subroutines. With this implentation, a designer familiar with using Abaqus may easily perform topology optimization without the need for acquiring and learning new software. It is expected that this scheme may also be easily implemented in other popular finite element software.

The adaptive topology optimization schemes of the invention provide robust topology optimization tools for use in identifying optimal, minimum weight designs of frame-like structures. Ih methodologies of the invention, a user that is already proficient in using a particular commerical program will have little difficulty in performing optimization studies. The optimization method may be applied to address optimization problems that have not been previously studied with other finite element based optimization procedures. The optimization method may be applied to consider the case of dual-material-property structures in which tensile members may be constructed using a different material than the compression members. Alternatively, the method may be applied for cases where a material whose tensile properties are different than the

compressive properties. Preliminary studies have shown that the optimization method may effectively identify optimal topologies for dual-material-property structures. Also, the method may be applied for structures that experience multiple alternative load configurations. Methodologies of the invention may further be employed for topology optimization even when different structural member types have different material properties. This expands the application of the methods to the latest generation of continuous carbon fiber reinforced plastic and ceramic fiber reinforced metal alloy structures. Other generalizations of the procedures include extension to three dimensions, optimization of multiply loaded structures, and extraction of precise geometry for component manufacture.

In various embodiments, the invention provides a node based evolutionary procedure that is immune to numerical instabilities, is easy to implement, and incorporates some attributes of homogenization procedures while retaining computational simplicity. Optimal topologies may be generated by a material enhancement and reduction algorithm. In a procedure of the invention, nodal densities are gradually adjusted based on nodal stress, strain or strain energy levels determined through iterative finite element analyses. The procedure is continued until the discussed volume is achieved. Many finite element based topology optimization methods treat design variables such as material density, stress, etc. as an element variable. Such methods are referred to as element based methods. Element based methods estimate design variables such as density, stress, strain, strain energy density, etc. by averaging at the element centroid or at

Gauss points. Such treatments do not guarantee inter-element continuity and may result in a structure that has elements with alternating high and low density areas referred to as checkerboarding.

The present invention provides a node based approach in accordance with an embodiment of the invention in which material density is specified at nodes. In a node based approach, design variables are defined as nodal variables and are interpolated within each element. This ensures inter-element continuity and as a result, prevents numerical instabilities. In order to update nodal density values, stress intensities (such as nodal stresses, strains or strain energies) are computed by extrapolation of integration point stresses to the nodes and averaging stress components from adjacent elements. This stress averaging approach is commonly used hi finite element post processing of stresses. Material properties required at element Gauss points are computed using appropriate interpolation functions.

With reference to Figure I ₃ an optimization algorithm used in an embodiment of me invention involves the following steps. The process begins (step 100) and a finite element model of a design domain is generated with prescribed loading and boundary conditions (step 102). The material density at all nodes is then defined as a field variable, and each node is assigned an initial density value of unity (step 104). Displacement finite element analysis is then conducted (step 106), wherein at material integration points the effective modulus is computed according to the relationship that E = E _o p where E is the effective modulus, E ₀ , initial modulus and p is the material density. The principle stress values O ₁ and σ ₂ are then computed and averaged at each node (step 108), and the strain energy density values are then averaged at each node as

A penalization percentage β is then defined (step 110). This parameter represents the fraction of nodes whose density is to he incrementally reduced. The number of nodes to be penalized in the current iteration is defined as Np ¹ = βN _node for the first iteration

Np ¹ = Np ^1"1 + (N _ω ode - Np ^1"1 ) β for subsequent iterations. where Np ¹ is the number of nodes to be penalized in the current iteration, N _nOde is the total number of nodes in the finite element model, and Np ^1"1 , the number of nodes penalized in the previous iteration. The nodal strain energy values are then sorted in the order of increasing magnitude as Us ¹ . The limiting value of strain energy U _LOW is then chosen (step 112) as the strain energy of the N _P ¹ th node in the sorted list, and all nodes with U' < UL _OW are identified as nodes to be penalized. The nodal density of the nodes identified

for penalization is then updated (step 114) according to p _mw ={p _j —&) and the density

value is increased for each of the remaining nodes according to p _mw =(P _f +(x) . To

maintain density values within the range 10 ^"4 < p < 1, nodal densities are adjusted according to the relationship:

if Pmw ≤ ¹ O ^"4 > setp _mw =10 ^" *

tf PNEW ^{≥ 1} ^ ^Set PNEiV = ¹ where p _Nm is the new density value for the next iteration.

The finite element analysis is then repeated (step 116) using the new nodal density distribution, and if a termination criteria is satisfied (step 118) then the program ends (step 120), otherwise the program returns to step 108 above.

To demonstrate the performance of this topology optimization algorithm four structural topologies are generated and are compared against known solutions. For the first test case a topology is also generated using a conventional scheme to demonstrate the improved computational effort and effectiveness of a methodology of the invention in identifying an optimal solution.

Figure 2 shows the design domain and boundary conditions for loading of a cantilever beam. The points of support of the initial volume 10 (of 20 mm by 14.14 mm) are shown at 12 and 14, and the applied force is shown at 16. Figure 3 shows at 20 the analytic optimal topology for this problem that provides the minimum weight. The topology shown at 22 in Figure 4 is slightly heavier (less than 1%) than the topology shown in Figure 2 but is near optimal for strength. For the finite element based scheme, the domain was discretized using 7100 bilinear quadrilateral elements. The optimization procedure seeks uniform strain energy in the domain, with volume reduction as the termination criteria. The topology corresponding to 54% volume reduction is shown in Figure 5. In particular, the topology following ten iterations is shown at 30, twenty iterations is shown at 32, thirty five iterations is shown at 34, forty five iterations is shown at 36, and fifty iterations is shown at 38. It is clear that in less than fifty finite element iterations, the optimal topology is identified.

These results compared with a typical element based material removal scheme for the same design using ESO procedures as shown in Figure 6 where the topology after one hundred iterations is shown at 40, two hundred iterations is shown at 42, four hundred iterations is shown at 44, six hundred iterations is shown at 46 and eight hundred eighty five iterations is shown at 48. It may be seen in Figure 6 that despite requiring many

more finite element iterations, the ESO method fails to identify the optimal topology and convergence to the near optimal solution shown in Figure 5.

As shown in Figure 7, a volume 50 (of 160 m hy 100 m) is supported by supports 52 and 54, and is subjected to a load force as shown at 56. The design domain of Figure 7 results in a centerfan design in accordance with an embodiment of the invention. Using a volume reduction of 71% and a design domain modeled using 12000 bilinear quadrilateral elements, the desired volume reduction is achieved in 150 iterations corresponding to volume reduction of 71%. The resulting structural topology is shown at 60 in Figure 8 and agrees with the analytic result shown at 62 in Figure 9. As shown in Figure 10, a volume 70 (of 100 m by 100 m) may be supported as shown at 72 and 74, and subjected to a load force as shown at 76. Using a target volume reduction of 83% and a design domain modeled using 19200 bilinear quadrilateral elements the desired volume reduction was achieved in 85 iterations. As in the previous example, the numerically generated topology as shown at 78 in Figure 11 agrees with the analytic result as shown at 80 in Figure 80.

As shown in Figure 13, a volume 84 (of 300 m by 300 m) may be supported as shown at 86 and 88, and subjected to a load as shown at 90. The height to span ratio may be, for example, 2.5. Using a target volume of 32% of the initial domain volume, the topology shown at 92 is shown in Figure 14 was achieved in 130 iterations using 14400 elements. Figure 15 shows at 94 the analytic minimum-weight topology for this design domain.

As shown in Figure 16, the domain with boundary conditions and loading provides a volume 130 (of 25 mm by 100 mm) that may be supported as shown at 132

and 134, and subjected to a load as shown at 136. The beam has a width to depth ratio of 4.8. The domain was modeled using 21000 quadratic bilinear elements. For this test case optimization procedure was terminated when a volume reduction of 36% was achieved. The resulting topology is shown at 140 in Figure 17. Analytic minimum-weight topology corresponding to this test case is shown at 142 in Figure 18.

In accordance with an embodiment therefore, the invention provides a node based structural topology optimization scheme. The illustrative examples demonstrate the ability of the proposed scheme to generate optimal topologies and are shown to be immune to numerical instabilities such as checkerboarding. The topologies generated using the current scheme were in very good agreement with minimum-weight topologies generated using analytic optimality criteria methods.

The performance of the new scheme was explored through several numerical examples. It is noted that the desired volume reduction was achieved in all cases with 50- 150 finite element runs. From the test case of Figures 2 - 4, a significant reduction in computational effort as compared to the ESO method is evident. It was observed that with recommended parameters for the ESO scheme, the number of finite element analysis required to achieve a prescribed volume reduction is significantly higher than that of the new scheme. From Figures 5 and 6, it is clear that the topologies generated by ESO and current scheme may be different. For this case, the ESO method generated the near optimal solution (Figure 4) while the current scheme successfully identified the optimal layout (Figure 3).

The evolutionary algorithm of certain embodiments of the invention is guided by two user defined parameters, α and β. Through parametric studies it was noted that with

values 0.1 and 0.01 for α and β, respectively successfully produced optimal topologies. The parameter α controls the degree of penalization applied to each node. It has been observed that a value of α above 0.40 causes rapid volume reduction. Rapid reduction in volume may result it large stress, strain and strain energy density variations, which can lead to non-optimal topologies and can cause slender members to disappear from the domain leading to near optimal solutions. Obtaining an approximate form of the optimal structural topology using higher values of α and β however, may yield some insights to the designer with very few number of finite element analysis.

From the above examples it may be concluded that the proposed scheme is immune to phenomena such as checker boarding and islanding in certain embodiments. It was observed however, that the proposed scheme is sensitive to relative sizes of structural members in optimal solutions. Mesh sensitivity analysis reveals that a very coarse mesh can cause slender structural members in optimal solutions to disappear, leading to near optimal solutions. It should be noted however, that the in such cases the difference in the material volume between optimal and near optimal solutions is negligible. The novel scheme with high mesh resolutions is capable of identifying minimum- weight structures as can be seen from the test cases. The topologies of structural members are identified from material density contours varying from 1 to 10 ^"4 . From the density contours it is evident that proposed scheme drives the design domain to fully dense material and voids, but it may be noted from the density contours the presence of intermediate density contours at the edges of fully dense structural members. This is an outcome of the progressive material penalization scheme used here and presents no difficulty in the interpretation of the topologies generated.

For the above examples, the termination criteria was the prescribed volume reduction, but the current scheme may be easily extended to handle other constraints such as prescribed deflections, compliance, stress levels, multiple loading cases, etc. The presence of optimality criteria along with perimeter controls can make the proposed scheme robust, capable of identifying global optimal solutions.

In accordance with another embodiment, the invention provides a further topology optimization methodology that involves continuous material redistribution and is based on sorted nodal stress, strain or strain energies. In such a topology optimization procedure and with reference to Figure 19, the process begins (step 200) and the first step (step 202) is to define the problem. The problem definition includes specification of the two or three dimensional space that is available for the structure. In particular, the design domain, support locations and types, applied loading, volume of optimal structure, and desired number of iterations i _max are all defined.

For cases with infinite or semi-infinite design domains, a finite size region that extends well beyond the expected topology is specified. The location and type of support conditions must also be specified, as well as the location, direction and magnitudes of any applied mechanical loads that the structure must support. Finally, the problem definition must specify the desired final volume, or mass, of the optimal structure.

A finite element model of the design domain is developed (step 204) consistent with the problem definition parameters (design domain) by discretizing the region into a finite number of elements and nodes and imposing load and support conditions. The domain is typically, though not necessarily, discretized into an array of regular shaped

elements. Load and support conditions consistent with the problem definition are applied to the model.

The desired final mass of the structure is specified at the beginning of the analysis and is held fixed through the course of the numerical iterations. This material mass is initially distributed uniformly throughout the design domain resulting in a uniform, partially-dense material. For first finite element iteration (/= 1), distribute mass uniformly throughout the design domain by assigning (step 206) to each node in the finite element model a relative nodal density of:

" ^• ' ^γ ' ^{l v} ' where Vf is the final structural volume and V _D is the volume of the partially-dense design domain. Since all of the nodes are initially assigned this relative density, the distribution of material densities may be described by the probability distribution function, f ₀ , given by

/XP) = S(P-P ₀ )

where 5 is the Dirac delta function and p is the relative material density (0 < p < 1). The corresponding cumulative distribution function, F ₀₅ is given by

F ₀ (P) = H(P^p ₀ )

where H is the Heaviside step function.

A finite element analysis is then performed (step 208) using an elastic modulus that depends on relative density through the relation

E = E ₄ P

The results of the finite element analysis include stress, strain and strain energy fields.

The nodes are then sorted (step 210) according to stress, strain or strain energy. The nodes are then ranked (step 212) based on stress, strain or strain energy, and nodal densities are assigned according to a prescribed family of statistical distribution functions that gradually transition from a unimodal distribution, where all nodes have the same relative density, to bimodal, where nodes are either fully dense (/7=1) or void (p=p _mm,

where p^ _n «1). The process then repeats steps 208 — 212 until i = z _max (step 214), whereupon the process ends (step 218). The process may also optionally include a step of adjusting (step 216) nodal density for stability as discussed in more detail below.

The desired final material distribution is characterized by distinct regions of fully dense material (p= 1) and regions that have zero relative density. The geometry of the fully dense regions represent the optimized structural topology. For the finite element calculations, a material density of zero would result in zero material stiffness and a singular global stiffness matrix. For these calculations, these regions are assigned a very- small relative density, p _m m « 1. The resulting final material distribution can be described by the probability distribution function, ff, given by

ff (p) = Q-Po) S(P -PnJ + P ₀ S(p -i)

and the corresponding final cumulative distribution is given by

F _f (p) = (l-p _o ) H(p - _Paia ) + p ₀ H(P -I)

- IS -

Note that for both the initial and final material distributions, the specified material mass is maintained.

For intermediate iterations, families of probability and cumulative probability distribution functions can be established which provide a smooth transition from the initial distribution for f _o (p) and F _o (p) above to the final distributions for f _f (p) and

F _j . (/?) . Figures 20 and 21 illustrate the gradual transition from the initial to the final

probability and cumulative probability distributions respectively. Conditions for the probability distribution function are that all material points have densities between zero and one, leading to the requirement

The requirement that the total mass of material be held constant is satisfied by the condition that the average density remain constant according to

Finally, a scheme is required to incrementally transition from the initial to the final material distribution.

One example of a family of distribution functions is given by the uniform distribution. For this family of distribution functions, Figure 20 shows at 146 the transition from initial to final probability distributions, and Figure 21 shows at 148 the transition from initial to final cumulative probability distribution.

Another example of a family of distribution functions is given by the beta distribution

where r and s are adjustable parameters and

F[r + s)

where F is the gamma function. The corresponding cumulative distribution function, also known as the incomplete beta function, is given by

F(p) = βjp,r.s) = - * føy ¹ (L-PT ¹ dp' ^ V = ^s ) o

This choice of material distribution functions satisfies this for any combination of parameters r and s. It can be shown that the requirement for constant average density p _a is satisfied if

To provide a smooth transition from the initial to the final material distribution, the following scheme can be imposed. A non-dimensional time parameter, t, is defined such that £ = 0 corresponds to the initial material distribution and t — 1 corresponds to the final material distribution. To provide a smooth transition from the initial unimodal

distribution where all material points are assigned the initial density to the final bimodal distribution where all material points are fully dense or have essentially zero density, an

appropriate function r(t) must be specified. This is achieved by determining r(t), O ≤ t≤ 1, such that

where Ff is the final bimodal cumulative distribution. The function Hf) can be determined numerically during the finite element iterations numerically using standard root finding and numerical integration algorithms. Alternatively, the function r(t, p _o ) can be prescribed analytically or established a priori with values stored in a lookup table; for example, using a Fortran DATA statement. During the finite element iterations, a bicubic interpolation algorithm is typically used to determine the desired values of r{t, p ₀ ). With this scheme, for the family of distribution functions described by the beta distribution, Figure 22 shows at 150 the transition from initial to final probability distributions, and Figure 23 shows at 152 the transition from initial to final cumulative probability distribution.

At each, finite element iteration, the desired r{t) is used to assign nodal densities according to the current density distribution function. The nodal density assignment is performed based on the sorted nodal stresses, strains or strain energies computed from the previous iteration. Nodes with relatively low stress, strain or strain energy are assigned reduced nodal densities and nodes with relatively high stress, strain or strain energy are assigned increased nodal densities. Through direct assignment of nodal densities, the

desired progression of density distributions is enforced. In computing the element stiffness matrices, the nodal density field is interpolated to give the Young's modulus, E, at each Gauss point according to the relation

E = E _d p

where E _d is the fully dense Young's modulus. Since the density of each node is reassigned at each finite element iteration, convergence to the final topology can be achieved in relatively few finite element iterations.

In the material redistribution scheme described above, nodal densities are assigned based on global nodal strain energy distributions with no consideration of the previous density history of a given node. As a result, if too few iterations are performed, numerical instabilities can occur. To avoid such instabilities, a scheme similar to the Euler method for the numerical solution of first order ordinary differential equations is used (step 216). In this scheme, the density at node n at iteration i+1 is computed using

where δp" is computed from the gradient of p" at iteration i. This gradient is estimated

by fitting a quadratic interpolation through the points p'V/, p",- andp"*ι ₊ j where /?" ^*• ₊ / is the new nodal density obtained from the strain energy sorting algorithm discussed above. The resulting nodal densities are updated according to

Through this scheme, stable convergence to the final topology can be achieved with relatively few iterations. The topology optimization scheme based on nodal strain energy sorting has been implemented using the commercial finite element code, Abaqus. It is expected that the above process should require fewer required iterations. Several different families of intermediate density distributions have been evaluated and have been shown to perform effectively. Families of intermediate density distributions that reliably predict optimal topologies may be identified for certain applications. Once an appropriate family of functions is identified, numerous two and three dimensional test cases may be evaluated to validate the robustness of the methodology.

In accordance with a further embodiment, the invention provides a methodology for dual strength or dual stiffness (dual property) materials. In this scheme, a fictitious modulus field is introduced such that the optimality criteria for dual property materials are satisfied. The design domain is initially taken to be fully dense. Through iterative finite element analyses, nodal densities are progressively adjusted based on current nodal strain energy levels. After each iteration, a new modulus field is established and the process is repeated, until the domain satisfies optimality criteria for minimum- weight or a prescribed volume reduction is achieved. Illustrative examples are presented in which the material domain converges to discrete structural members that are either in the state of pure tension or compression and are stressed to the tensile and compressive strengths, respectively. The resulting structural topologies satisfy optimality criteria for dual strength minimum- weight topologies.

For example, optimality criteria for generating minimum- volume topologies for dual strength materials require that tensile and compressive members follow the

directions of maximum and minimum strain (ε-r, εc) respectively in a virtual strain field which is compatible with the support constraints and which satisfies

where, σj, £τ are the tensile strength, tensile strain and , σr, sc are the compressive stress and strain respectively.

Introducing a fictitious modulus for members in compression and tension, this requirement reduces to

Zl = ^l

E _j , E _Q

with Eγ, E _C being modulus in tension and compression respectively. If it is assumed that the tensile strength is related to compressive strength as σ _τ = nσ _c

combining the preceding equations provides the relation between the fictitious moduli as

^ = n ²

E _c

It may be inferred that once the design domain is reduced to regions of uniaxial tension and compression with stress levels reaching the respective strength limits, the resulting structural configuration will correspond to the minimum- volume topology. Initially the domain is assumed to be fully dense. Convergence to the optimal solution requires that regions in the state of uniaxial tension should take the modulus value of E _T and regions in the state of uniaxial compression take on a modulus value of Ec ■ During the course of the finite element iterations, certain regions experience a biaxial stress state.

In these regions, appropriate intermediate modulus values must be assigned. A procedure for determining an appropriate modulus value for regions of biaxial stressing is described below.

Let σj and σ ₂ denote the principal stresses averaged at a given node with σj > σ ₂ . If σj > σ _? > 0, the node is predominantly tensile and is assigned a modulus, ET. Similarly 0 > σj > σo, the node is taken to be in a compressive state and is assigned a modulus, Ec- For nodes where σj > 0 > σ? , an effective strain energy is defined as

This strain energy density is used to define the following effective modulus for regions in a state of biaxial tension/compression

With this formulation, the fictitious modulus varies continuously as the state of stress in a given node varies in the course of finite element iterations. Converged topologies consist of tensile and compressive members where the tensile members have a modulus of E _T and compressive members have a modulus Ec. Regions of biaxial stress state essentially disappear during the course of the iteration. The continuous fictitious modulus and density fields associated with the node based approach are shown to provide optimal topologies with no numerical instabilities. The above scheme can readily be extended to three dimensions. The optimization algorithm used for the present embodiment is the process discussed above with reference to Figure 1 except that in step 106, p is the density at

any given iteration with E _τ and E _c being the modulus in tension and compression

respectively, and that step 108 further includes establishing a fictitious modulus based on optimality criteria from the equation above for Eφ The process repeats as in Figure 1 until convergence criteria are met.

The optimization procedure was implemented using a MATLAB code which generates input file data for a finite element analysis code ABAQUS system based on the above mentioned optimization algorithm. After each finite element iteration, stress values were extracted from output files by the MATLAB code. Based on nodal stresses strain energy values were computed and material density values are updated. Optimal structural topologies are generated for two cantilever designs using the current scheme. The material considered for both test cases are taken to be dual strength with, the ratio of tensile strength to compressive strength being three (n=3). Topologies generated using current scheme are compared against analytical layouts generated using matrix operator methods. The optimization procedure was terminated when prescribed volume reduction was achieved. The solutions obtained were optimal as the design domain was driven by optimality criteria and the design domain reduces to regions of pure tension and compression with respective fictitious modulus.

The first example for this embodiment considered involved a cantilever loading with dimensions, loading and boundary conditions as shown in Figure 2. The structural topology corresponding to a 60% volume reduction is shown in Figure 5 (after 50 iterations) for the case of material having the same strength in tension and compression. For the case of dual strength material with n=3, the optimal topology is shown at 160 in Figure 24. The theoretical layout is shown at 162 in Figure 25. Comparing these numerical results with theoretical layouts for these cases reveals excellent agreement.

In another dual strength cantilever example, a volume was considered with a span to height ratio of 2.5. The design parameters were similar to those shown in Figure 13 with the distance between supports (86 and 88) being 100 m, and the distance between the applied force (90) and a line connecting the supports is 250 m. Using a volume reduction of 68% yields topologies similar to that shown in Figure 14 for the case of single strength material, but yields a topology as shown at 164 in Figure 26 for the case of dual strength materials. Comparing these results to the theoretical layouts to this case as shown at 166 in Figure 27 reveals excellent agreement.

The methodology of the present embodiment, therefore, may be developed for structures composed of members having different properties in tension and compression. The optimization procedure proposed here is readily implemented. By specifying n=l. optimal structural topologies may be generated, which permits the proposed scheme to be generally applied. It has been observed that optimization of material having different strength in tension and compression increased computational time compared to the material having same strength in tension and compression. This is explained by the fact that establishing the fictitious modulus field requires more equilibrium iterations.

One significant advantage of the proposed scheme over other finite element based topology optimization algorithms is that the parameters that control the optimization procedure are limited to two. The default values of the two parameters α and β used for the illustrative examples were α = 0.10 and β=0.01. Through multiple optimization runs it was observed that higher values of these parameters cause large density variation in a given step, which in turn may cause significant variation in strain energy levels. Large variations in strain energy densities can result in non-optimal solutions. A lower value of

α and β however, increases the computational effort due to a reduced material removal rate.

The illustrative examples considered here converged to the prescribed volume reduction in less than 150 finite element iterations. For the test cases presented here it may be noted that for the same design loads the volume of material for optimal structural topologies of material having same strength in tension and compression and dual strength material was less than 3%. The robustness of the proposed scheme is established by the excellent agreement observed between the analytic solutions and optimal topologies generated using the proposed scheme. To investigate the procedure used to establish fictitious modulus, contour plots where generated for the modulus field corresponding to a converged solution. Tensile and compressive structural members converge to moduli values of Ex and Ec respectively. The node based approach ensures a continuous distribution of moduli in the design domain.

Topology optimization methodologies therefore, are provided for design domains composed of dual strength materials. The optimization procedure was based on optimality criteria for minimum-weight structures, where the material considered have different strength in tension and compression. A node-based approach immune to various numerical instabilities was used to establish both the fictitious modulus distribution and the material density distribution. From the illustrative examples it is concluded that the proposed optimization procedure is robust in identifying minimum- weight topologies for dual strength material domain. The use of a commercial finite element analysis code makes the proposed scheme highly versatile.

These dual strength material demonstration cases are performed using the algorithm described in Figure 1. Alternatively, the algorithm described in Figure 19 may also be employed, in which case step 208 further includes establishing a fictitious moduli based on optimality criteria from the equation above for E _e ff. In accordance with further embodiments, the method may be extended to the design of three dimensional structures as well. In this case, the design domain is taken to be three-dimensional, typically but not necessarily rectangular in shape and discretized using hexahedral brick finite elements. Alternatively, tetrahedral finite elements can be utilized. As an example, Figure 28 shows at 170 the analytic solution for a three dimensional generalization of the two dimensional centerfan shown in 80 in Figure 12. In the three dimensional generalization, a concentrated vertical load is applied at the base of a semi-infinite design domain with eight equally spaced roller supports along a circular path that prevents motion in the vertical and circumferential directions. The corresponding result from the finite element based topology optimization is shown in 172 of Figure 29.

In yet further embodiments, the invention may be applied to the case of multiple loading. Extension of the above methods to the case of multiple loadings is achieved by considering the stress, strain or strain energy distributions associated with each loading. Again, the design domain is initially taken to be partially dense such that the desired final mass of material is uniformly distributed throughout the design domain. During each finite element iteration, material is gradually redistributed using a family of statistical distribution functions which evolve from the initial unimodal distribution to the final binαodal distribution. The material stiffness is taken to be proportional to the material

density such that MIy dense regions have the actual material stiffness and material with zero density have essentially zero stiffness. In each iteration, the stress, strain or strain energy at each node is computed. The nodal stress, strain or strain energies are then sorted and the material density is adjusted according to the current desired density distribution. As shown above for cases of single loading, several cases have been examined in which the optimal topology is correctly identified.

For the case of multiple loadings, the material redistribution scheme is applied after each of the separate load configurations have been applied independently to the design domain with the current density distribution. For each node, the maximum stress, strain or strain energy experienced for each load is saved and the nodal stresses, strains or strain energies are sorted. Then, the material redistribution scheme is applied to the sorted maximum nodal stress, strain or strain energies. In this manner, if a region experiences high stress, strain or strain energy during any of the loadings, its density will be increased. As a result, over the course of the iterations, structural members evolve such that all loadings are supported.

In the single loading procedure, a single load step is applied at time, t=0 ₃ and the material redistribution is performed using time dependent material properties. At 1=0, the unimodel density distribution is imposed and at time t=l, the final bimodal distribution is imposed. For the case of multiple loadings, the time stepping scheme is modified such that numerous load steps are applied. For example, if three load configurations are to be considered, Step 1 is taken to be the first loading, Step 2 is taken to be the second loading, Step 3 is taken to be the third loading, Step 4 is taken to be the first loading for the second iteration, etc. In this example, at the end of every third load step, the

maximum stress, strain or strain energy experienced by each of the nodes during the previous three load steps is sorted and the algorithm is applied to determine the density distribution for the next three load steps.

Validation of the procedure described above was achieved by applying the procedure for test cases where the optimal topology for multiple loadings is known. In certain cases, the optimal topology may consist of either two or three member truss structures, depending on the load orientations, θ _P1 and θp ₂ . These predictions are based on a failure criteria in which the material is taken to be elastic - perfectly plastic and all members reach the yield stress, σy. For three member trusses, since yielding of one member does not correspond to collapse of the structure, the collapse load is higher than the elastic limit. Examination of cases where three members are predicted to be optimal, reveals that lighter two member structures may be found if the failure criterion is taken to be initial yielding. For example, consider the cantilever beam case where P ₁ — P ₂ = P and

θpi = -θp2 = 60°. Using plastic collapse as the failure criterion predicts an optimal three member truss. If however, initial member yielding is taken to be the failure criteria, the volume of material required to support the load for a two member truss is V = 2.585 PL / σy with members oriented at ± 36.36°. This volume is smaller than that for the three

member truss where members are oriented at 45°, 0°, -45° where V = 2.688 PL /σy. Note that for the two member solution, each member has half the total volume of the structure. For the three member solution, the member volume fractions that minimize the required volume are Vi = V ₃ = 0.462, v ₂ = 0.076. It is of interest to note that if plastic collapse is taken to be the failure criterion, the three bar truss topology gives a value of V = 2.564 PL / σy,. Since the finite element based optimization procedure does not consider plastic

deformation and is designed to keep stresses at or below the elastic limit, it is expected that topologies predicted by the finite element based scheme will correspond to analytic topologies based on the elastic limit. Based on the observations discussed above, it is expected that the optimal solution is a two member truss where the members are oriented

at ±α. Due to symmetry of the loadings, the required cross-sectional area of the two members is identical and, since the structure is statically determinate, can be computed from force equilibrium. The required volume to support the load, V, is non-

dimensionalized as V - V l {PL I σ _γ ) . It can be shown that for a given member

orientation, α, the non-dimensional volume required to support both loads is given by

— _ cos θ _p sin θ _p cos a sin a cos a The multiple loading procedure discussed above was applied to the above test case for various load angles 0° < θp < 90°. A 71 x 142 element rectangular design domain was selected. Two hundred load steps (100 topology optimization iterations) were used for these cases. The desired final volume was taken to be 5% of the design domain. Examination of the final density contours reveals that in all cases, two member truss topologies were identified. To estimate the member angle, α, the vertical coordinate, y, of the node closest to the center of each member was extracted from the finite element results. The member angle was then estimated using the formula

a RJ tan ^"1 {y/ L) . These results provide excellent correlation with the theoretical result.

Also, the optimal truss volumes computed with the member angles determined from the finite element topology optimization results are compared to the theoretical values. Again, excellent correlation is observed. The results of this study demonstrate that

W

extension of the optimization procedure to cases involving multiple loading effectively predicts optimal topologies for the test case considered. .

In accordance with further embodiments, methodologies of the invention may be employed in a system for extracting a precise geometry for component manufacture. Coupling of finite element based topology optimization with precise geometry extraction algorithms and rapid manufacturing technologies provides a comprehensive tool set in the development of lightweight component designs. For illustrative purposes, an algorithm is developed in which the results of 2-D topology optimization are used as input. The extraction algorithm defines smooth surfaces that are then saved in a standardized CAD format.

The 2-D finite element topology optimization results consist of 2-D density field data as characterized by nodal densities. For model data extraction, the finite element optimization code has been adapted to create three output files. The node definition file contains the node numbers and x- and y-coordinate data stored, in an N _no d _e x 3 array where N _nOd e is the number of nodes in the finite element model. The element file definition contains the element connectivity in an N _e i _em x 4 array where N _e i _em is the number of elements in the finite element model. The nodal density data file contains the nodal densities at each finite element iteration stored in an N _s tep x N _no d _e array where N _s t_p is the number of finite element iterations. Density fields for all iterations are saved to provide for visualization of the incremental evolution of the structure.

Construction of a 3-D surface model from the 2-D field results is achieved in a two step process. First, the top and bottom surfaces of the model are constructed by defining triangulated surfaces on an element by element basis. Following the marching

cube classification procedure for 3-D surface modeling, each 4 node element is classified by a binary designation in which each node is assigned a value of 1 if its density is above a specified threshold, and 0 if it is below the density threshold. With this method of classification, 16 different surface types are possible as shown at 300 - 330 in Figure 30. These classifications may be represented by 4 digit binary numbers that consist of the four nodal assignments. The dense surfaces can be defined as triangles. For the lateral edges, each element is imagined to be extruded in the thickness direction to a specified thickness. The lateral edges are then classified in a manner similar to that imposed for the top and bottom surfaces. In this case, however, only 4 surface types as shown at 332- 338 in Figure 31 need to be defined. With these surface classifications, the 2-D finite element based optimization results for each 4 node element are used to create triangulated surfaces on 3-D hexahedral elements.

After identifying the triangulated surfaces for the hexahedral elements. The exact location of the triangle vertices are determined through linear interpolation of the nodal densities using the threshold density for the desired surface. Since the density field results often consist of discontinuous variations in nodal density within elements along the edge of the converged domain, generation of smooth surfaces is enhanced by imposing a nodal density smoothing. The following smoothing algorithm can be imposed in which the nodal density of each interior node is replaced with a weighted average of densities of adjacent nodes.

- 1 1 / \ 1 / λ

Pf ⁼ ^P>i ⁺ g [P'-hJ ⁺ PMJ íPU-X +Aj ₊ I J +— (λ-U-I + /Wi +A ₊ Ij ₊ I +A-U ₊ I J

To investigate the procedure for geometry extraction, the case of a simply supported beam with a central load applied to the top of the mid-span is considered. The

optimal topology is known from analytic optimality criteria. A 2 -D finite element based optimization procedure is applied and a contour image providing regions of fully dense material and regions of essentially voids. These results are shown at 340 in Figure 32 and are then used as input for the geometry extraction algorithm to produce a 3-D triangulated surfaces that encloses the fully dense regions. These triangulated surfaces are then saved in a standard CAD format (.STL) that was originally developed for the rapid prototyping technology, stereolithography. STL files are now used routinely as input for a wide variety of rapid manufacturing processes. The resulting STL file for the example case is shown at 342 in Figure 33. Those skilled in the art will appreciate that numerous modifications and variations may be made to the above disclosed embodiments without departing from the spirit and scope of the present invention.