Technical Field of the Invention

The present invention concerns, in general, Finite Element Method (FEM)-based Computer- Aided Engineering (CAE) for structural free-form design, and, in particular Adaptive Free-Form Design Optimization.

The present invention finds advantageous, though not exclusive, application in the free-form design of structures for the subsequent Additive Layer Manufacturing (ALM). In fact, the present invention may also find application in the free-form design of structures for their subsequent Layer less Additive (casting techniques), Non-Additive (multi-axis machines, spark-machining, etc.), and mixed Additive-Subtractive Manufacturing.

State of the Art

The emerging additive layer manufacturing provides designers with an enormous, previously unthinkable, variety of shapes for objects. It is based on the addition of material with the "quasi-absence" of tools, thus overcoming the limits of traditional manufacturing based on removal of material, in terms of freedom of the possible objects shapes. Most commercial CAE-FEM softwares have been developed to design objects having a shape complexity limited by subtractive production constraints.

Additive manufacturing considerably broadens the freedom and complexity of the possible objects shapes, but also drastically increases the computational resources necessary for implementing CAE-FEM design. Optimal solutions therefore have to be found in a wider range of permissible configurations.

Free-form topology optimization based on iterative algorithms requires a vast number of FEM analysis iterations that result in a considerable increment of the computational time, especially when large structures are to be designed. For example, in the case of topology optimization performed with the iterative SIMP (Solid Isotropic Material with Penalization) algorithm, in some cases the computational capacities are exceeded with no useful results.

Object and Summary of the Invention

The object of the present invention is to overcome the limits of iterative optimization algorithms, also including those of the aforementioned SIMP algorithm, so as to considerably reduce the computational cost and, at the same time, provide a completely free-form topology optimization.

According to the present invention, a computer-aided system is provided for FEM -based structure design, as claimed in the appended claims.

Brief Description of Drawings

Figure 1 schematically shows a computer-aided system to FEM-based structure design according to the present invention.

Figure 2 shows a block diagram of the operations implemented by the computer-aided system according to the present invention.

Figure 3 shows anisotropic quantities of a generic mesh element.

Figure 4 shows benefits resulting from isotropic and anisotropic mesh adaptations with respect to a uniform mesh. Detailed Description of Preferred Embodiments of the Invention

The present invention will now be described in detail with reference to the accompanying drawings to enable those skilled in the art to embody and use it. Various modifications to the described embodiments will be immediately appreciable to those skilled in the art, and the generic principles described herein can be applied to other embodiments and applications without departing from the scope of the present invention, as defined in the appended claims. Thus, the present invention is not intended to be limited to the disclosed embodiments, but is to be accorded the widest scope consistent with the principles and features described and claimed.

Figure 1 schematically shows a computer-aided system, designated as a whole by reference number 1, for FEM-based structure design for subsequent additive layer manufacturing.

The computer-aided system 1 basically comprises a computer 2 with a user input device, in the example shown comprising a keyboard and a mouse, and a graphical display device, in the example shown in the form of a screen.

Computer 2 basically comprises a processor and an internal and/or external memory device, in which a program for structural design based on the finite element method is stored, and which, when executed by the processor, causes the computer 2 to become programmed to implement an improved SIMP algorithm, hereinafter referrred to as SIMPATY for brevity, comprising the operations hereinafter described with reference to the flow chart shown in Figure 2.

In particular, as shown in Figure 2, the computer 2 is programmed to acquire (block 100) an initial structure design configuration, which can be entered by an operator via the user input device and comprises:

^■ a design domain Ω,

^■ a constrained boundary portio ,

^■ an unloaded boundary portion ,

■ a loaded boundary portion

^■ an applied load f, whether superficial, inertial, concentrated or distributed,

^■ Young's modulus E of the material, and

^■ Poisson ratio v. Initial structure design configuration is then used by the computer 2 in the standard linear elasticity equations to compute a displacement of the structure being designed, when loaded with load f on boundary portion ΓΝ and constrained on boundary portion -

Computer 2 is further programmed to compute (block 200) an initial mesh for the design

domain Ω, by using any mesh generation software available in the literature.

Computer 2 is further programmed to receive (block 300) control parameters, described in more detail hereinafter and indicated as CTOL, MTOL and pmin, intended to control the structure topology optimization, described further on with reference to block 400, and the anisotropic mesh adaptation described further on with reference to blocks 500, 600 and 700.

Computer 2 is further programmed to iteratively repeat the operations of topology optimization and of anisotropic mesh adaptation driven by an anisotropic recovery-based a posteriori error estimator, hereinafter described with reference to blocks 400 to 700, which represent the core of the SIMPATY algorithm, which combines the reliability and computational easiness of an anisotropic recovery-based a posteriori error estimator with the effectiveness of an anisotropic mesh adaptation. As is known, structural topology optimization is a mathematical approach through which the layout of the material of the structure being designed in the design domain Ω is optimized, under the load and constraint conditions specified in the initial design configuration, in such a way that the resulting layout satisfies given design and performance targets. Mesh adaptation is instead a numerical procedure that improves the performance of a finite element solver by modifying the size of the mesh elements, which in 2D are typically of a triangular or square shape, while in 3D they are usually tetrahedral or hexahedral. In particular, the aim of mesh adaptation is to make the mesh elements to be smaller where the phenomenon to be investigated exhibits more complex local characteristics, and to instead use larger mesh elements where the physical solution is regular.

Unlike a standard structured mesh, where the size of a mesh element is uniform over the entire domain, an adapted mesh allows reducing the number of mesh elements, i.e., the size of the finite element algebraic system, for the same solution accuracy, or to increase solution accuracy for the same number of mesh elements.

Independently of the criterion adopted, the computational advantages deriving from the use of adapted meshes are appreciable and now consolidated in the literature. Adapted meshes can be classified into isotropic or anisotropic meshes depending on their geometric characteristics.

Isotropic meshes are formed by very regular and (quasi-) equilateral elements, and the only quantity that changes is the size, i.e. the diameter (see Figure 4, middle picture). Conversely, anisotropic meshes can be characterized by highly elongated elements, thus allowing control of the size, shape and orientation of the elements (see Figure 4, right-hand picture), thereby introducing greater freedom in the computational mesh design. In particular, this flexibility is found to be ideal for describing physical problems characterized by high directionality, as long as the mesh elements are aligned with these directional characteristics, for example, shocks in compressible streams in aerospace applications, multi- material flows with abrupt immiscible interfaces in 3D and ALM printing, and boundary layers in viscous flows around bodies or walls.

Since they only allow adjusting the size of the mesh elements, standard isotropic meshes are not able to capture these directional characteristics, while the adaptation of an anisotropic mesh is a powerful tool for improving the quality and efficiency of finite element procedures. By way of example, Figure 4 shows the effect of adaptation on an isotropic mesh (middle picture) and on an anisotropic mesh (right-hand picture) of an advection-diffusion problem in an L-shaped domain, where the solution shows two inner circular layers and three boundary layers. The accuracy of the solutions based on the two adapted meshes is clearly greater than that based on the uniform mesh. For the same solution accuracy, the number of mesh elements for the isotropic and anisotropic meshes is 24,499 and 5,640, respectively. Mesh adaptation can be implemented through heuristic or theoretical criteria. Heuristic approaches basically employ information related to the derivatives of the numerical solution, such as the solution's gradient or Hessian. Instead, the theoretical approaches are based on sound mathematical tools, known as error estimators, which provide explicit control of the discretization error in terms of the exact solution (a priori error estimators) or of directly computable quantities (a posteriori error estimators).

Popular a posteriori error estimators include those that are recovery-based (based on gradient reconstruction), residual-based (based on the residual associated with the discrete solution), and dual-based (based on the error associated with a physical quantity of interest).

Recovery-based a posteriori error estimators were proposed by O. C. Zienkiewicz and J. Z. Zhu in 1992 and are widely used in the engineering field due to their excellent properties. Recovery-based a posteriori error estimators are not confined to a specific problem, are independent of discrete formulation (except for the selected finite element space), are cheap to compute, easy to implement and work extremely well in practice.

The main idea underlying recovery-based a posteriori error estimators is to compute the discretization error by replacing the gradient of the exact solution with a discrete enriched (or reconstructed) gradient with respect to the gradient of the FEM solution.

However, the recovery-based a posteriori error estimators introduced by O. C. Zienkiewicz and J. Z. Zhu in 1992 were confined to an isotropic context, while their anisotropic counterpart was proposed by S. Micheletti, S. Perotto in Anisotropic adaptation via a Zienkiewicz-Zhu error estimator for 2D elliptic problems, Numerical Mathematics and Advanced Applications 2009, Proceedings of ENUMATH 2009, the 8th European Conference on Numerical Mathematics and Advanced Applications, Uppsala, July 2009, Gunilla Kreiss, Per Lotstedt, Axel Malqvist, Maya Neytcheva Editors, Springer-Verlag Berlin, pp. 645-653, 2010, in the two-dimensional case, and subsequently, in the three-dimensional case by P.E. Farrell, S. Micheletti, S. Perotto in A recovery-based error estimator for anisotropic mesh adaptation in CFD, Bol. Soc. Esp. Mat. Apl., 50, 115-137 (2010), and by P.E. Farrell, S. Micheletti, S. Perotto in An anisotropic Zienkiewicz- Zhu type error estimator for 3D applications, Int. J. Numer. Meth. Engng, 85 (6), 671-692 (2011).

Instead, the combination of a topology optimization with an anisotropic mesh adaptation driven by heuristic criteria based on data provided by the Hessian of the design variable (density of the material) combined with the Hessian of the cost functional with respect to the design variable, and where all the quantities involved in the Hessians are preliminarily filtered, is proposed by K. E. Jensen in Anisotropic mesh adaptation and topology optimization in three dimensions, J. Mech. Design, doi: 10.1115/1.4032266, by K. E. Jensen in Solving stress and compliance constrained volume minimization using anisotropic mesh adaptation, the method of moving asymptotes and a global p-norm, Struct. Multidisc. Optim., doi: 10.1007/s00158-016-1439-9, by K. E. Jensen in Optimising Stress Constrained Structural Optimisation, Research Note, in 23rd International Meshing Roundtable (IMR23), 2014, and by K. E. Jensen and G. Gorman in Anisotropic Mesh Adaptation, the Method of Moving Asymptotes and the global p-norm Stress Constraint, arXiv preprint a rXiv:1410.8104.

Unlike the solutions proposed in the aforementioned literature, the present invention proposes a synergetic combination of a topology optimization with an anisotropic mesh adaptation driven by an anisotropic recovery-based a posteriori error estimator, rather than by a heuristic approach, and without necessarily performing any kind of filtering, in order to develop a CAE-FEM design tool aimed at free-form design, with lower computational costs with respect to those associated with other adaptation methods and with a rigorous procedure from the theoretical standpoint.

The present invention differs from what proposed in the aforementioned articles of S. Micheletti and S. Perotto, Anisotropic adaptation via a Zienkiewicz-Zhu error estimator for 2D elliptic problems, A recovery-based error estimator for anisotropic mesh adaptation in CFD, and An anisotropic Zienkiewicz-Zhu type error estimator for 3D applications, since in the estimator proposed in these articles for the two-dimensional and three-dimensional cases, respectively, is used for the first time in the present invention to drive a topology optimization procedure. No mention at all of a topology optimization is made in these articles.

The present invention further differs from what is proposed by K.E. Jensen in the aforementioned articles Anisotropic mesh adaptation and topology optimization in three dimensions, Solving stress and compliance constrained volume minimization using anisotropic mesh adaptation, the method of moving asymptotes and a global p-norm, Optimising Stress Constrained Structural Optimisation, and Anisotropic Mesh Adaptation, the Method of Moving Asymptotes and the global p-norm Stress Constraint, since the anisotropic mesh adaptation is driven by a theoretical tool, namely a recovery-based a posteriori error estimator, rather than by heuristic criteria.

The advantages deriving from the combination proposed in the present invention between a recovery-based a posteriori error estimator and an anisotropic mesh adaptation procedure driven by this estimator, compared to the procedure proposed by K.E. Jensen based on anisotropic mesh adaptation driven by a heuristic estimator, have been verified in the test case shown in Figure 9 of the article of K.E. Jensen Anisotropic mesh adaptation and topology optimization in three dimensions. In particular, the comparison has been carried out for the same number of nodes in the final adapted grid (3,833 nodes for the SIMPATY algorithm versus 3,121 nodes for K.E. Jensen's algorithm), and minimising a target function defined, in this case, by the compliance, namely the yieldingness of the final structure. The output structure of SIMPATY algorithm was characterized by a compliance equal to 1.31144, versus a compliance computed by K.E. Jensen's algorithm equal to 1.5529, so showing that the SIMPATY algorithm results in a more performant structure, namely more rigid.

Furthermore, the comparison shows the extensive difference between the two algorithms in terms of execution times, on computers with similar characteristics. In particular, the execution time of the SIMPATY algorithm was equal to 0.9 hours, versus the 3.3 hours of K.E. Jensen's algorithm.

To implement the SIMPATY algorithm, the computer 2 is first programmed to compute an optimized topology (block 400) of the structure being designed by properly implementing the aforementioned iterative SIMP algorithm. As is known, this is based on a density function p that may range in and which represents the distribution of material in the structure

(Pmin = material absent (empty), 1 = material present (full)), and on the basis of which Young's modulus E is replaced by an effective modulus where p = 3 is a penalty exponent. The result of the topology optimization is represented by an optimized density function p ₀ ptm. In an alternative embodiment, instead of the iterative SIMP algorithm, the topology optimization could be performed by using any other topology optimization algorithm based on a density function p or on its generalization (known in the literature as phase field). The computer 2 is further programmed to first compute (block 500) an anisotropic recovery- based a posteriori error estimator η that quantifies the error between the gradient of the exact density p and the gradient of its FEM approximation, and then compute a metric M based on this anisotropic estimator, as better described hereinafter. Both η and M are then used to compute a new adapted anisotropic mesh.

In particular, the mathematical tool that drives the adaptation of the anisotropic mesh is represented by a global anisotropic recovery-based a posteriori error estimator η that collects the local anisotropic recovery-based a posteriori error estimators ηk of all the elements K of the mesh Th, with:

The local anisotropic recovery-based a posteriori error estimator ηκ associated with the element K of mesh Th tracks the relevant anisotropic geometrical characteristics of the element K, together with the typical structure of an anisotropic recovery-based a posteriori error estimator.

As proposed in the aforementioned article Anisotropic adaptation via a Zienkiewicz-Zhu error estimator for 2D elliptic problems, the local anisotropic recovery-based a posteriori error estimator ηκ may be computed as follows:

where is the positive semi-definite symmetric matrix:

where:

is a generic vector-valued function, and

is the patch of mesh elements associated with K.

The anisotropic characteristics of element K of the mesh T _h are represented by the quantities

, which define the length and direction of the semi-axes of the ellipse

circumscribing element K, with (see Figure 3).

The "recovery-based" nature of the local anisotropic estimator ηk is instead represented by the quantit i.e by the (recovered error) difference between the

gradient of the discrete density _h and a suitable reconstruction of this gradient, confined to the patch of elements The technique proposed in the aforementioned Anisotropic adaptation via a Zienkiewicz-Zhu error estimator for 2D elliptic problems is adopted for the reconstruction P of the gradient

i.e. the weighted average of the area of the discrete gradients associated with the triangles of the patch Δ _κ is considered. At each iteration k of the iterative procedure, starting from mesh a metric defined

by new values of is then computed, based on the local anisotropic estimators ηκ,

which minimizes the number of mesh elements while at the same time ensuring a certain accuracy. The new values of the quantities are computed on the basis of an error equidistribution criterion, according to which ηκ is approximately constant on every element of the new adapted mesh combined with the minimization of the number of elements in the mes which guarantees a certain accuracy MTOL on the discretization error (in

practice, on η).

Computer 2 is further programmed to compute an adapted anisotropic mesh based on

the previously computed metric and using a suitable mesh generator (block 600).

Computer 2 is further programmed to monitor (block 700) the convergence of the SIMPATY algorithm based on the outcome of two checks, the first of which is for checking the maximum number k _max of iterations, while the second is for checking stagnation in the iterative procedure, interrupting the iterative procedure when the relative change in mesh cardinality (number of mesh elements) in successive iterations is less than a threshold CTOL, thus returning a topologically optimized structure model ready for use in additive layer manufacturing and defined by an adapted anisotropic mesh T _h and an optimized density function

Pseudo-code of the SIMPATY algorithm is shown below. The input data consists of the aforementioned CTOL, MTOL, and parameters, where CTOL is the mesh

cardinality stagnation control parameter, MTOL defines the discretization error accuracy, k _wax represents the maximum number of iterations, p _m m represents the minimum density value of the structure that identifies the absence of material (empty), and T¾ represents the initial mesh of the design domain Ω computed in block 200.

Topology optimization described with reference to block 400 is implemented by the Interior Point Optimizer IPOPT, which solves a generic constrained optimization problem.

Input data for the IPOPT function has been chosen as indicated in http://www.coin- or.org/lpopt/documentation/nodelO.html. even if, in principle, the operator can choose other values or resort to other optimization tools.

Instead, the anisotropic mesh adaptation described with reference to block 600 is implemented via the adapt function.

The two kernels IPOPT and adapt are linked by the key quantity which defines the FEM

discrete density, piecewise linear, at generic iteration it. SIMPATY algorithm thus iteratively alternates topology optimization and anisotropic mesh adaptation driven by an anisotropic recovery-based a posteriori error estimator until either mesh cardinality stagnation or the maximum number of iterations is reached.

The SIMPATY algorithm enables overcoming the limits of the SIMP algorithm extremely well, considerably reducing the computational cost by a factor of roughly three in terms of degrees of freedom, at the same time providing a completely free-form topology optimization.