Structural weak spot analysis
1. A method of determining weak spots in an object, comprising the steps of:
- determining maximal principal stress on each selected point of the object andobtaining a corresponding pressure distribution and selecting the pressure distribution corresponding to the point for which maximal stress is the greatest thereby determining the weak spots on the object;
wherein an objective function is further provided to determine optimization, wherein the objective function comprises;
A system and methods for weak spot analysis. A mesh shape is preprocessed to approximate an input model for an object. The mesh shape is analyzed with modal analysis to identify weak regions. A method and system for determining weak spots in an object. The method and system uses an optimization problem which is solved to determine a pressure distribution on the object maximizing maximal principal stress by solving a set of optimization problems maximizing stress for each point of the object.
|METHODS AND SYSTEM FOR ENABLING PRINTING THREE-DIMENSIONAL OBJECT MODELS|
Patent #US 20110087350A1
Current Assignee3D M.T.P. LTD
Sponsoring Entity3D M.T.P. LTD
|Evaluation method and evaluation apparatus for spot welded portion|
Patent #US 7,543,732 B2
Current AssigneeMitsubishi Jidosha Kogyo Kabushiki Kaisha
Sponsoring EntityMitsubishi Jidosha Kogyo Kabushiki Kaisha
|METHOD FOR THE OPTIMIZATION OF STRESS DISTRIBUTION IN ACOUSTIC THERMOGRAPHY APPLICATIONS|
Patent #US 20090204345A1
Current AssigneeSiemens AG
Sponsoring EntitySiemens AG
|Adjoint-based design variable adaptation|
Patent #US 20060058985A1
Current AssigneeSUPERSONIC AEROSPACE INTERNATIONAL LLC
Sponsoring EntitySUPERSONIC AEROSPACE INTERNATIONAL LLC
|Computer program product for automated fatigue and structural analysis of an element|
Patent #US 6,862,539 B2
Current AssigneeThe Boeing Co.
Sponsoring EntityThe Boeing Co.
|Structural stress analysis|
Patent #US 6,901,809 B2
Current AssigneeBattelle Memorial Institute
Sponsoring EntityBattelle Memorial Institute
|In situ recovery from a hydrocarbon containing formation using conductor-in-conduit heat sources with an electrically conductive material in the overburden|
Patent #US 20040020642A1
Current AssigneeRoyal Dutch Shell Plc
Sponsoring EntityRoyal Dutch Shell Plc
|Method for determining a model for a welding simulation and model thereof|
Patent #US 6,768,974 B1
Current AssigneeCaterpillar Incorporated
Sponsoring EntityCaterpillar Incorporated
|Method and system for determining crack nucleation of a part subject to fretting fatigue|
Patent #US 20030074976A1
Current AssigneeJalees Ahmad
Sponsoring EntityJalees Ahmad
|Integrated system for aircraft crack detection|
Patent #US 5,065,630 A
Current AssigneeVought Aircraft Industries Inc.
Sponsoring EntityGrumman Aerospace Corporation
|INDIVIDUALLY TAILORED SOFT COMPONENTS|
Patent #US 20120232857A1
Current Assignee3shape AS
Sponsoring Entity3shape AS
|STRESSES INDUCED BY RANDOM LOADING|
Patent #US 20120239358A1
Current AssigneeAirbus Operations Limited
Sponsoring EntityAirbus Operations Limited
|METHOD FOR THE PREDICTION OF FATIGUE LIFE FOR WELDED STRUCTURES|
Patent #US 20120259593A1
Current AssigneeDeere Company
Sponsoring EntityDeere Company
|MODAL ASSURANCE CRITERION FOR COMPARING TWO MODE SHAPES OF A MULTI-COMPONENT STRUCTURE|
Patent #US 20130030771A1
Current AssigneeAirbus Group India Private Limited
Sponsoring EntityAirbus Group India Private Limited
|OPTIMIZING THE DESIGN OF PHYSICAL STRUCTURES/OBJECTS|
Patent #US 20140214370A1
Current AssigneeHonda Research Institute Europe GmbH
Sponsoring EntityHonda Research Institute Europe GmbH
|Method of fabricating a mechanical part, including a method of predicting the risks of crack initiation in the part in a “fretting-fatigue” situation|
Patent #US 9,081,922 B2
Current AssigneeAirbus Helicopters
Sponsoring EntityAirbus Helicopters
|Frequency Domain Structural Analysis Of A Product Having Frequency-Dependent Material Property|
Patent #US 20150308936A1
Current AssigneeLivermore Software Technology Corporation
Sponsoring EntityLivermore Software Technology Corporation
|INTELLIGENT PAD FOOT SOIL COMPACTION DEVICES AND METHODS OF USING SAME|
Patent #US 20150316526A1
Current AssigneeColorado School of Mines
Sponsoring EntityColorado School of Mines
|Fabrication of objects with enhanced structural characteristics|
Patent #US 9,308,690 B2
Current AssigneeMakerBot Industries LLC
Sponsoring EntityMakerBot Industries LLC
- 1. A method of determining weak spots in an object, comprising the steps of:
determining maximal principal stress on each selected point of the object and obtaining a corresponding pressure distribution and selecting the pressure distribution corresponding to the point for which maximal stress is the greatest thereby determining the weak spots on the object; wherein an objective function is further provided to determine optimization, wherein the objective function comprises;
- View Dependent Claims (2, 3, 4, 5, 6, 7, 8, 9, 10)
This application is a 371 National Stage of PCT Patent Application No. PCT/US2013/066439, filed Oct. 23, 2013, which claims the priority to U.S. Provisional Application No. 61/855,989,filed May 29, 2013, and U.S. Provisional Patent Application No. 61/718,129, filed Oct. 24, 2012. The entire contents of these applications are hereby incorporated by reference in their entirety.
The United States Government claims certain rights in the invention pursuant to NSF Grant IIS-0,905,502.
This invention is directed to a method and system for weak spot and also structural analysis of a shape or object. More particularly, the invention is directed to a method and system for use in commercial applications wherein mechanical stresses and strains can lead to endurance limit problems and system failure for a selected shape (object) and/or a material, such as for printing processes or any cyclic process that creates a work piece subject to mechanical forces. Further, the invention relates to a method and system for use in commercial applications wherein structural problems in objects are characterized for manufacture or for 3D printing based on geometry and material properties in particular and without making specific assumptions based on loads or manual load setup.
In traditional engineering analyses a typical goal is to determine if stresses in an object or system are within bounds permissible for the material and/or the shape, typically well below the yield stress which is the point beyond which the material undergoes deformation and does not return to its original shape upon release of a load on the material.
Engineering analysis protocols typically assume that an engineer should determine realistic static or dynamic loads on a material part or structure. For example, an engine part typically repeats the same cyclic motion, with a known range of forces applied in a particular way. The loads on a piece of furniture or on a building are even more predictable, and, in the latter case, test loads are often codified. In general this approach and complemented by prototype testing, this methodology works well in mass-production, construction, and small-volume construction of complex machinery (e.g., ships and airplanes). This approach however is known to fail in some instances: for example, a well-known example is of the Citigroup Center building which had to be retrofitted because it was not tested for 45-degree wind loads at the design stage. This illustrates a limitation of testing with a set of predefined test loads. 3D printing applications typically follow a different operational path. Unlike conventional manufacturing methods, 3D printing is primarily used to produce unique or highly customized objects. In a typical scenario, a printing company for such processes will receive a large volume of uploads per day. In order to keep production times and costs down, the company must evaluate whether a design is structurally sound; this needs to be done in a rapid and inexpensive way but often without specific knowledge of design function and likely load distribution. In most cases, it is no longer possible to amortize the expense of engineering analysis by an expert to understand the semantic and plausible loads of each printed object by using standard engineering analysis software.
Computational analysis of structural soundness of mechanical parts and buildings is broadly used, but almost always in the context of known sets of loads. While engineers routinely need to evaluate soundness of structures and mechanisms under worst-case scenarios, in most cases, worst-case loads are designed empirically for specific problems (e.g., construction of buildings to withstand loads from flooding or earthquakes). Automatic methods are less common: an important set of methods in the context of modeling under uncertainty is based on the idea of anti-optimization.
In aerospace engineering, filter-based methods were developed to predict worst-case gusts and turbulence encountered by an airplane describes a model the aircraft'"'"'s response to turbulence as a linear filter'"'"'s response to Gaussian white noise. From this model, a worst-case noise sample and resulting strain are obtained.
Three-dimensional (3D) printing and other types of direct digital manufacturing are rapidly expanding industries that provide easy ways to manufacture highly customized and unique products. The development pipeline for such products is radically different from the conventional manufacturing pipeline: 3D geometric models are designed by users often with little or no manufacturing experience, and sent directly to the printer. Structural analysis on the user side with conventional tools is often unfeasible as it requires specialized training and software. Trial-and-error, the most common approach, is time consuming and expensive.
In the context of analysis tailored for 3D printing applications of the type considered herein, current methods evaluate 3D shapes in two main scenarios to discover structure weakness: applying gravity loads and gripping the shape using two fingers at locations predicted by a heuristic method. This set of fixed usage scenarios is often insufficient to expose the true structure weakness for many printed shapes. Other methods focus on purely geometric ways to evaluate whether a structure is suitable for 3D printing based on empirical rules formulated by the 3D printing industry.
In yet another area, structural stability for simple furniture constructed from rigid planks connected by nails is analyzed at interactive rates. These methods also suggest corrections when shapes with poor stability are detected.
Other recent works address various aspects of computational design for 3D printing. One method provides a pipeline to print objects in a composite material that reproduces desired deformation behavior. To achieve this goal, accurate modeling is done of the nonlinear stress-strain relationship of printing materials and how printed models will respond to imposed loads. The space of deformations is a user-supplied input, and structural soundness of the design with respect to other loads is not considered. While some specialized work on CAD for 3D printing exists, overwhelmingly, standard tools with no or little analysis support are used.
Other work proposes a framework to decompose 3D shapes into smaller parts that can be assembled without compromising the physical functionality of the shape so that larger objects can be printed using printers with a small working volume. A standard finite element simulation to estimate stress of the input shape under gravity in a user specified upright orientation. Other works aim to print articulated models that maintain poses under gravity but do not require manual assembly designs and fits a generic, parametrized printable joint template based on a ball and socket joint. The joint provides enough internal friction and strength to hold poses and survive manipulation, but is tuned to its parameters experimentally instead of using a physically based optimization. Other designs a similar ball and socket joint and a hinge joint. An approximate geometric optimization of stresses is performed by maximizing certain cross-sectional areas of the joint.
3D printing has also been used to reproduce appearance to optimize the layering of base materials in a 3D multi-material printer to print objects whose subsurface scattering best matches an input BSSRDF.
Common single-stage 3D printing processes either deposit the liquid material only in needed places (e.g., FDM) or deposit material in powder form layer by layer, and then fuse or harden it at points inside the object (e.g., stereolithography uses photosensistive polymers, and laser sintering fuses regular polymers by heat).
These processes typically use flexible polymers with large elastic and plastic zones in their stress-strain curves. These polymers rarely break if geometric criteria for printability are satisfied, but they can undergo large plastic deformations.
Printing metal, ceramics, and composite materials often involves multiple stages. For example, the object may be printed layer by layer in metal powder with polymer binder. At the next stage, the binder is cured in a furnace, resulting in a green state part, and at the last stage, the metal is fused in a furnace and extra metal is added. Green state is brittle and has low strength, so parts in this state easily damaged. A simpler multistage process is used for relatively brittle composite materials, e.g, gypsum-based multicolor materials; a second curing stage is used to give the material additional strength. Both the green state and the final material are relatively brittle. Whenever binding polymer is mixed layer-by-layer with a different material, resulting material is likely to be highly anisotropic.
Therefore, both brittle and ductile materials are of importance. The former requires predicting where the material is likely to break, and the latter requires predicting extreme deformations likely to become plastic. Due to the layer-by-layer nature of the printing process, anisotropy is common and needs to be taken into account. Some of the loads even during production stages are hard to predict and quantify.
Consequently, it is desirable to design a system and method capable of automatically identifying “worst case” load scenarios, and evaluating possible stresses under various conditions for a system. A fully automatic system of this type would preferably require as input variables only a knowledge of the shape of the object or material of the system under load, total maximal surface forces applied, and allowed stress ranges. From these input variables, the system would then determine possible ways to distribute the load on the surface of the object or material to achieve maximal stresses, and determine if any of these stresses are outside a safe range. Alternatively, for a given maximal stress, the system would determine a distribution of loads producing this stress with minimal total force applied. It is further desirable to design a system and method capable of identifying “worst case” structural problems in objects, particularly for 3D printing, based only on a knowledge of the shape or geometry of the object or material of the system and material properties only, with no specific assumptions made about loads and manual load setups.
In one embodiment, the present invention combines the basic assumptions set forth above with heuristics allowing us to provide a computationally feasible method and system. In this embodiment of invention the following analyses are preferably utilized:
(a) The high stress in most “bad” cases (i.e. load configurations that require minimal force to exceed a limit) tend to be concentrated in relatively small parts of the object (“weak” regions); (b) weak regions (“wr”) can be identified by modal analysis, i.e. artificial “shaking” of the object at its natural frequencies; and (c) once weak regions are isolated, the problem of determining load distributions which maximize stress at a specific weak region can be approximated by a convex optimization problem evaluation (or even a linear programming problem evaluation), for which a solution is unique and reliable solution method tools are available in the art.
In a preferred methodology, the process can be roughly divided into four steps: preprocessing, modal analysis, weak region extraction and stress optimization. The preprocessing step can be done in a variety of ways, and potentially can be eliminated entirely in selected embodiments.
In its current form the method has certain limitations which can be recognized and engineering applications can be adjusted and safe limits established while still applying the method and system effectively. Weaker region detection appears to work well for a broad range of examples, but it is ultimately a heuristic method. Additional heuristics may be needed to provide more comprehensive coverage. Computing multiple modes can be costly (although it can be reasonably easily performed in parallel). Complex object shapes and sizes may require many modes before all weak regions are found. More generally, the mathematical definition of a weak region can to be refined to agree better with basic engineering evaluation and projections: some stress distributions are known to be more dangerous than others, and more likely to lead to fracture (e.g., compressive vs. extensional stress is not distinguished at this point).
In yet another embodiment of the invention, a method and system has been developed that identifies structural problems in objects designed for 3D printing based on geometry and material properties only, without specific assumptions on loads and manual load setup. The problem is formulated as a constrained optimization problem to determine the “worst” load distribution for a shape that will cause high local stress or large deformations. While in its general form this optimization has a very high computational cost even for relatively small models, an efficient heuristic is demonstrated based on modal analysis and an approximation is provided by a linear programming problem that can solve the problem quickly for the typical size of printed models. The method is validated both computationally and experimentally and demonstrate that it has good predictive power for a number of diverse 3D printed shapes.
In this embodiment, a method and system is therefore provided for approximating the solution of the following problem: From the shape of an object and its material properties, determine the easiest (in terms of minimal applied force) ways to break it.
In order to understand the nature of the invention it is useful to consider the reasons why a suitable 3D model cannot be manufactured by current methods or is likely to fail:
- (1) the dimensions of thin features (walls, cylinder-like features, etc.) are too small for the printing process, resulting in shape fragmentation at the printing stage;
- (2) the strength of the shape is not high enough to withstand gravity, at one of the stages of the printing process;
- (3) the printed shape is likely to be damaged during routine handling during the printing process or shipment;
- (4) the shape breaks during routine use.
In most cases, the first problem is addressed by simple geometric rules in a known way and the second is a straightforward direct simulation problem. The instant invention focus is on the other remaining two problems of the sequence. On the one hand, many 3D printed objects are manufactured with a specific mechanical role in mind, and full evaluation is possible only if sufficient information on expected loads is available. On the other hand, jewelry, toys, art pieces, various types of clothing, and gadget accessories account for a large fraction of products shipped by 3D printing service providers. These objects are often expected to withstand a variety of poorly defined loads (picking up, accidental bending or dropping, forces during shipping, etc.).
To be able to predict structural soundness of a printed object, worst-case loads are sought, within a suitably constrained family, that are most likely to result in damage or undesirable deformations. A direct formulation results in difficult nonlinear and nonconvex optimization problems has been developed with PDF constraints. An approximate method for this search, reducing it to an eigenproblem and a sequence of linear programming problems. The method preferably relies on using eigenmodes of the shape. Modal analysis has proven useful in many contexts. The use of Laplacian eigenmodes of simple shapes for computation predates computers and has a long history in model order reduction for a variety of applications including nonlinear elasticity. In graphics literature, first introduced eigenmodes as a basis suitable for simulation applications, and more recently, a number of deformation-related algorithms based on eigenmode bases have been proposed. At the same time, experimental modal analysis (applying periodic forces with different frequencies and measuring displacements at various points) is broadly used to detect structural damage.
The foregoing summary is illustrative only and is not intended to be in any way limiting. In addition to the illustrative aspects, embodiments, and features described above, further aspects, embodiments, and features will become apparent by reference to the following drawings and the detailed description. These and other objects, advantages and features of the invention, together with the organization and manner of operation therefore, will become apparent from the following detailed description where taken in conjunction with the accompanying drawings wherein like elements have like numerals throughout the several drawings described below.
The foregoing and other features of the present disclosure will become more fully apparent from the following description and appended claims, taken in conjunction with the accompanying drawings described hereinbelow. Understanding that these drawings depict only several embodiments in accordance with the disclosure and are, therefore, not to be considered limiting of its scope, the disclosure will be described with additional specificity and detail through use of the accompanying drawings.
In one embodiment of the invention,
In another embodiment of the invention,
In a preferred system and method of one embodiment of the invention, key elements include separating weak spot identification from determining “worst-case” loads; using a global method to identify multiple weak regions at the same time; and using a constrained surface load model to determine the weakness of each region. The actual choice of method for each step may vary. For the first step (weak region identification), a modal analysis is used, based on one of the elasticity discretizations described in detail hereinafter. This methodology computes deformations corresponding to global smoothly varying force distributions and areas of high stresses are identified as weak regions. One can use a purely geometric method, based on computing local thickness of an object using for example a medial axis decomposition, voxel erosion or any other such conventional technique. Modal analysis is a well-established technique used by engineers and architects for determining the structure response under vibration forces such as strong wind, earthquakes and other known nearly-periodic forces. The modal analysis can be performed using a finite element discretization on a polyhedral mesh conforming to the boundary of the object, a boundary element method, or a meshless discretization approach. With regard to the figures, red (“r”) was the highest stress followed by the other colors and blue (“b”) was lowest stress.
For the second step (determination of worst-case forces), different types of constrained force models can be used: in a preferred implementation positive forces are applied normal to the surface, with maximal per point magnitude and total load constrained. This models typical force behavior in handling and transporting 3D printed objects. One can also consider non-normal forces applied in arbitrary directions, impose further constraints on force density variation, and include volume forces such as gravity. As a target for optimization, different norms of the average, or maximal stress on the weak region can be used, such as von Mises stress, maximal principal stress and other conventional stresses known in the art. As an indicator of failure, various material-dependent yield strength, ultimate compressive and tensile strengths can be used. The application of the method is not restricted to the domain of 3D printing. The method and system can be used to predict structural strength/weakness for any object for which a CAD model (a surface mesh, a spline or subdivision-based boundary representation, a solid model represented by a CSG data structure, voxelized model, level set model or another 3D computer representation for which engineering analysis can be performed is available, and for which a typical load scenario is unknown.
Implementation of the invention is accomplished most readily by a computer system (stand alone or in the cloud) having a non-transitory memory to execute computer software implementing the various steps described throughout the application.
In one most preferred embodiment of this first embodiment, the methodology applies a preprocessing step to produce a tetrahedron mesh (called a “tet” mesh in the rest of this document). The starting point for tet mesh generation, for which a variety of techniques and software are available, is a triangle mesh bounding a region of spacec, which can be easily obtained from most common shape representations. Most tet mesh generation techniques add vertices in the interior of the triangle mesh and construct tetrahedra connecting these vertices. As many additional vertices may be needed in the interior to maintain good tet mesh quality, for large input triangle meshes the process may result in very large tet meshes. To keep the time and memory requirements of modal analysis and force optimization within a reasonable range, the triangle mesh may be simplified before the tet mesh is generated, using one of many available feature-preserving mesh simplification methods. (see
In a next step modal analysis carries out computations of natural frequencies and modal displacement of an input shape. With zero damping and zero external force, a shape deformed according to a modal displacement would vibrate indefinitely at the corresponding natural frequency.
In one example,
Modal analysis gives a sequence of distinct deformations (e.g. modal displacements) that are likely to capture common strong deformations experienced by objects, with no prior assumptions on loads. The modes with lowest eigenvalues require less energy to excite than modes with larger eigenvalues (assuming each modal displacement is normalized). For this reasons the low-eigenvalue modes are more common in our everyday life because the less prominent modes requires too much energy to start or its amplitude is too small to notice.
Modal displacement thus is important for revealing structural strength and weakness of a shape or object. The following sections describe how to exploit this information and determine whether a shape is sufficiently strong for 3D printing (or any other appropriate application involving the same issues for objects under stress and having “weak” regions).
In regard to weak region extraction, structural failures occur in regions where the stress exceeds the maximum stress that the material or shape could tolerate. This can depend on whether the material is brittle, and on its yield strength and/or ultimate tensile strength. The goal of this step is to use modal analysis results to identify and isolate weak regions, which we define as continuous domains where the stress is within a given threshold from its maximal value.
Modal analysis can be used to identify weak regions, for reasons explained hereinbefore. The basic idea of modal analysis is to compute a number of eigenvectors of the stiffness matrix of the system object, which correspond to its “natural” vibration modes. Before converting modal displacement into stress, one normalizes all modal displacements so each of them require the same amount of energy to excite. Once normalized, displacements are converted to stresses using a standard linear elasticity model.
In the step of stress optimization, knowing the weak regions alone is not enough to determine whether a shape can be printed. The stress optimization step determines the maximum stress induced at each weak region with a fixed total force budget. If this stress exceeds the maximal permitted values, the structure may fail if the worst-case loading is applied. In its simplest form, the stress optimization can be formulated as a linear programming (LP) problem where we constrain the total magnitude of forces. Currently, the LP is solved using a well-known cvx package available at (http://cvxr.com/). This cvx computer software package is available under a GPL license, but a variety of other well-known and available packages exist for solving this type of problem.
Further examples of this type of analysis are shown in
In one preferred embodiment of this first embodiment, shown in
As shown in
System 100 may also include a display or output device, an input device such as a key-board, mouse, touch screen or other input device, and may be connected to additional systems via a logical network. Many of the embodiments described herein may be practiced in a networked environment using logical connections to one or more remote computers having processors. Logical connections may include a local area network (LAN) and a wide area network (WAN) that are presented here by way of example and not limitation. Such networking environments are commonplace in office-wide or enterprise-wide computer networks, intranets and the Internet and may use a wide variety of different communication protocols. Those skilled in the art can appreciate that such network computing environments can typically encompass many types of computer system configurations, including personal computers, hand-held devices, multi-processor systems, microprocessor-based or programmable consumer electronics, network PCs, minicomputers, mainframe computers, and the like. Embodiments of the invention may also be practiced in distributed computing environments where tasks are performed by local and remote processing devices that are linked (either by hardwired links, wireless links, or by a combination of hardwired or wireless links) through a communications network. In a distributed computing environment, program modules may be located in both local and remote memory storage devices.
Various embodiments are described in the general context of method steps, which may be implemented in one embodiment by a program product including computer-executable instructions, such as program code, executed by computers in networked environments. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. Computer-executable instructions, associated data structures, and program modules represent examples of program code for executing steps of the methods disclosed herein. The particular sequence of such executable instructions or associated data structures represents examples of corresponding acts for implementing the functions described in such steps.
Software and web implementations of the present invention could be accomplished with standard programming techniques with rule based logic and other logic to accomplish the various database searching steps, correlation steps, comparison steps and decision steps. It should also be noted that the words “component” and “module,” as used herein and in the claims, are intended to encompass implementations using one or more lines of software code, and/or hardware implementations, and/or equipment for receiving manual inputs.
In a second embodiment of the invention, in order to carry out the worst-case structural analysis, a formal description of the problem includes a computationally intractable formulation, but it is needed as a foundation for a practical and advantageous approximate version described hereinafter.
An anisotropic linear material model is used and the linear elasticity equations applied to model object behavior for the purposes of determining weak spots and worst-case force distribution. This model is adequate for some materials used in 3D printing, but nonlinear models may be necessary for others, as discussed in greater detail hereinafter. A distinction is emphasized that should be made between simulation with given loads used to determine precise stress distributions and computation used to determine approximate worst-case loads: lower accuracy is acceptable for the latter. The standard elasticity is included to introduce notation. The stress-strain relationship is linear, and stress is related linearly to displacement:
where ε is the strain tensor, σ is the stress tensor, and u is the displacement. C is the rank-4 elasticity tensor, Cijlm, and the notation C: ε denotes application of this tensor to the strain tensor ε,Σl,mCijlmεlm. The choice and effects of elasticity tensor C are discussed hereinafter in greater detail. An orthotropic material is assumed, for which the tensor Cijlm has up to 9 independent parameters. In a coordinate system aligned with material axes, if we represent C as a 6×6 matrix acting on vectors of components of the symmetric strain tensors [ε11,ε11,ε22,2ε23,2ε31,2ε12], its inverse is given by
Where Yii are directional Young moduli, Gij are shear moduli, and νij/Yi=νji/Yj.
For dynamic linear problems with volume force density F, the equation of motion is
where ρ is the density, and the dot signifies the time derivative. This is primarily directed to static problems, but modal analysis can be used at an intermediate stage, and thus the term ρü is retained.
Equation 2 is subject to boundary conditions: we primarily use a surface force density Fs, which is captured by the condition σn=−Fs on the boundary of the object. If the object is attached to a rigid support, Dirichlet conditions u=0 can be imposed on a part of the boundary.
If the equation of motion (2) is written directly in terms of displacement u, we get
Rigid motion, torque and translation constraints for static problems. If the object is not fixed at least at 3 non-collinear points, an arbitrary force distribution will result in motion of the whole object. The interest is in considering unknown forces, with no assumptions on attachment, and thus we need to be able to eliminate global motion. This is achieved by imposing zero total force and zero total torque constraints, which can be written as,
Displacements enter into this system only in the form Lu, and the operator L has infinitesimal rigid motions in its nullspace. To have a unique solution in u, we impose the zero rigid motion constraint, similar to total torque and force constraints:
In cases of interest for a surface force model, the only volume force is gravity. In all but most extreme cases, gravity does not have a major effect, so we concentrate on surface forces. The forces of interest are restricted in three ways.
Only positive normal forces are allowed: Fs=−pn, where n is the surface normal, and p is pressure, thus ignoring friction. This is an important assumption, as for most situations described hereinbefore, friction forces are likely to be significantly lower than normal forces. At the same time, it is hard to model the bounds on ratios between normal and tangential components accurately in the absence of detailed knowledge of loads and surfaces. Without such bounds, any optimization is likely to produce unrealistic tangential results. Similarly, negative surface forces (e.g., electrostatic attraction), are not likely to play a major role and are excluded.
For the condition of pointwise pressure bounded: p<pmax. If a pressure can be unbounded, an arbitrarily high stress may be produced at a point on the surface. While highly concentrated forces are possible, these are rare, and it is assumed that a realistic bound on surface pressure is available.
For the total force is fixed, again, by increasing the total force, arbitrarily high stresses can be obtained. For example, if our primary target is simulating manual handling situations, one can bound the force by a typical force a human can apply by squeezing, and the maximal pressure is derived from the size of the finger tips.
It remains to specify the objective function. One commonly used measure of interest is maximal principal stress, maxΩ maxi=1,2,3|σi|, where σi are the eigenvalues of the stress tensor. The complete problem of finding the worst-case force distribution satisfying the constraints of our model and optimizing this objective function, has the form,
Maximal principle stress is a suitable measure if one is interested in failure of materials, which occurs when the stress in a direction exceeds a bound. For plastic transition, the norm or some other function of the deviatoric stress,
can be of interest.
One observation can be made when the material is isotropic; and C can be written as YĈ, where Y is the Young modulus, and Ĉ is nondimensional, depending only on the Poisson ratio. Then maximal stress does not depend on Y but only on the Poisson ratio.
Solving this problem yields the worst-case principal stress and, importantly, the pressure distribution on the surface resulting in this stress. The maximal stress makes it possible to evaluate the likelihood of damage during the production process, shipping or use. Examining the pressure distribution makes it possible to evaluate how likely such loads would be and determine how the structure of the object can be strengthened.
All constraints in this problem are observed to be of linear equality and inequality constraints, i.e., the constraints are convex. At the same time, the functional is highly non-linear (in fact, not smooth) and non-convex. Replacing maximal principal stress with another point measure maximized over the surface does not change the nature of the problem.
A brute-force solution can be obtained by solving a sequence of problems in which the objective functional max |σ1|2 is maximized for every point, and then the max is taken over all results. Because we are interested in maximizing the norm, even these simpler per-point problems remain nonconvex and nonlinear. Solving the optimization problem in general form is generally impractical, and due to non-linearities and non-convexity, any optimization is likely to get stuck in local minima.
An extension of the algorithm is directed to optimizing for maximal displacements. The main change is replacing a with u in the functional: maxΩ|u|→max. This formulation is more relevant for flexible materials.
An illustration of an efficient approximate algorithm methodology,
There are two problems to address to make the solution of Eqn. (6) practical: (1) the need to solve an optimization problem for each point of the object to determine which one results in minimal stress; and (2) the non linearity and nonconvexity of each subproblem.
To address the first problem, a modal-analysis based heuristic is used that works remarkably well. The second problem is solved by using tr σ as the linear objective functional. The reasons this substitution is possible for a broad range of cases are discussed in detail below. Modal analysis and weak region is a crucial ingredient of our method, which is used to restrict the part of the object where one needs to maximize the stress or another functional.
Computational modal analysis refers to computing eigenvectors (eigenmodes) ui and eigenvalues λi of L:
Lui=λiui, i=1,2 (7)
In the context of structural analysis, this modal analysis is to predict possible damage or deformations in presence of vibrations. Vibrations, i.e., periodically changing loads are not considered; rather, static or quasi-static loads are considered in the prior art. The following assumptions were used in this second embodiment of the invention:
- Assumption 1: Examining a small number of eigenmodes allows us to find all regions of an object where the stress may be high under arbitrary deformation. While this observation is difficult to prove mathematically, reasonable logic suggests that vibrations of an object at different frequencies will result in high stress in all structurally weak regions of the object. Weak regions are those where high maximal stress can be obtained with low energy density relative to other parts of the object.
To validate this assumption, a brute-force optimization was performed on a number of models (see
The approximate convex problem is concerned with the second important change to the problem by replacing the functional in Eqn. (6) with a functional that can be optimized efficiently and that is minimized by a similar pressure distribution, p, to the original. Focus is thus on the maximal stress, although a similar approach can be used for other functionals. It is observed that almost invariably for any deformation and any compressible material with Poisson ratio v sufficiently different from ½: For points where a principal stress is maximal, other principal stresses are small relative to the principal stress.
Validation was performed of this observation by running simulations with a variety of loads and computing the ratio of the maximal principal stress to |trσ|. Over 36 models tested, the average ratio is 0.96 with standard deviation 0.25.
When this is true, the difference between |σmax|=maxi=1,2,3|σi| and |Σ3i=1σi|, i.e., |trσ| is so small, and one can approximate the maximal principal stress with the absolute value of the trace. Finally, it is noted that problems can be solved separately for maximal compressive stress, approximated by −tr σ, and tensile stress, approximated by tr σ. Thus focus is on the tensile stress.
As weak regions correspond to the highest stress area, and estimated stress tends to have a significantly lower accuracy vs. displacement, a weighted average is used of the stress over each weak region. The choice of weighting, as long as it falls off towards the boundary of the region, has relatively small effect on the result. We choose the L2 norm of the stress computed from the eigenmode as the weight w for averaging the stress trace over each weak region.
The following approximate problem formulation is obtained:
For each eigenmode i, i=1 . . . Mm and each of its weak legions, Dij, j=1 . . . Mr, we solve the following linear programming problem:
Lu=0 on Ω,C:(∇u+∇uT)n=pn on ∂Ω
p≥0,p≤pmax on ∂Ω;∫∂ΩpdA=Ftot (8)
Unlike the original problem, this problem has a unique solution that can be computed efficiently using a convex solver.
A discretization and additional optimizations embodiment solves another part of the problem. The problem can be discretized in the simplest conventional way, using piecewise-linear finite elements. The downside of this approach is that a suitable tetrahedral mesh needs to be generated for each input. For 3D printed models, the task is somewhat simplified: as the cost of printing is dominated by the amount of material used, almost all objects printed in practice are effectively thick shells to the extent this is allowed by the structural requirements. For this reason, the meshing does not increase the number of vertices used to represent the object as much as one would expect.
Let n, be the number of vertices, nb≤n be the number of boundary vertices, and m be the number of elements. The discretized quantities are: p the vector of pressures defined at boundary vertices of dimension nb; and u, the vector of displacements of dimension 3n.
In discrete formulation, we optimize the functional,
wTV D Bu (9)
In this formula, V is a 6m×6m matrix, with the volume of element j repeated 6 times on the diagonal for the 6 components of the stress tensor. D is a 6m×6m block-diagonal matrix. For each element, the corresponding 6×6 block is the rank-4 tensor C in matrix form. B is a 6m×3n applying the FEM discretization (or applying a conventional FV or other like function) of ∇+∇T. Finally, wT is a vector that computes and weights the stress tensor weights, so that wTVx discretizes ∫ΩwtrxdV.
The discretized static elasticity equation combined with boundary conditions takes the form
where K is the standard FEM 3n×3n, stiffness matrix, K=BTV DB. The matrix N is a 3n×nb matrix of components of surface normals, returning per-vertex components of external forces (0 for internal vertices, pn for the boundary), and matrix A is the nb×nb, diagonal vertex area matrix.
The discretized formulation of the total force and torque constraints are:
ΣN Ap=0,ΣT N Ap=0 (11)
where Σ is the 3×3n matrix, summing n 3d vectors concatenated into a 3n vector, and T is 3n×3n block-diagonal matrix computing the torques of the surface force vectors.
Putting all these together, the discretized optimization problem is:
w·(V D Bu)→max w.r.t. u and p
ΣN Ap=0, ΣT N Ap=0
0≤pi≤pmax for all i
where Σs sums scalars on the surface, Σv, sums vectors in the volume Ω, and Tν computes torsion for each point. The total dimension of the problem is nb+3n.
Eliminating displacements embodiment. As most of the degrees of freedom in the system are displacements, but quantities of interest are pressures p, eliminating u results in significant calculational speedups (u can be eliminated even for the displacement maximization problem). The elasticity equation −Ku+Nap=0 is not sufficient for this; it has a nullspace of dimension 6 corresponding rigid motion degrees of freedom, so we need to consider the constraints for zero total rigid motion, Ru=0, where
Rewriting this system in the standard constrained system form,
where λ is the Lagrange multiplier for the constraint Ru=0. It is clear from physical considerations that this system is invertible. Let S be the selection matrix
Then, we can express u as u=STC*−1SNAp. In this form, the objective of Eqn. (12) becomes
The displacement-free optimization problem is,
ΣN Ap=0,ΣT N Ap=0,
While the final system has only sparse constraint matrices, it may appear that computing fT for the objective functional requires inverting C*; we observe however that wTVDBSTC*−1SNA=fT can be rewritten as f=(SNA)Tq, where q is the solution of the equation
In other words, it is sufficient to be able to solve a linear system with matrix C*, and the cost of transforming Eqn. (12) to Eqn. (15) is the cost of a single linear solve.
Finally, for modal analysis, it is observed that the results for isotropic models in particular are well-approximated by simpler eigenanalysis of the Laplacian, which yields a considerable speedup.
Algorithm summary and parameters. The main steps in this second embodiment:
- 1. Compute a tetrahedral mesh Ω for an input triangle mesh.
- 2. Compute Mm, modes using an eigensolver.
- 3. For each mode, find Mr weak regions with highest total energy.
- 4. For each region Di solve the problem Eqn. (14) to obtain worst case pressure candidate pi.
- 5. Solve Lu=0, with boundary pressures specified by pi, to obtain displacements ui, and compute actual maximal principal stress σimax for each weak region.
- 6. Maximal stress is determined as maximum of σimax.
Various conventional available software packages are used including MOSEK to solve the linear programming problem, UMFPACK for linear solves, and ARPACK for computing eigenvectors and eigenvalues. Other well-known conventional packages can be used.
The parameters of the algorithm include Mm, Mr, the choice of threshold 1-ε for weak regions, as well as user-defined maximal pressure ρmax (the latter can be regarded as a part of the definition of the force model).
To determine reasonable choices of Mm and Mr, and the modal analysis has been run for a large number of modes (50); and a large number of weak regions per mode for a collection of objects. For each object and each mode, the weakest region was found and checked in which mode it first appears. Its rank in the list of that mode'"'"'s weak regions was computed and sorted by decreasing energy. The results (
We use ε=0.025 in all cases; the dependence of the size of weak regions for one mode on ε is shown in
Validation of the methodologies were done in several computational and experimental ways. Comparison with direct search for the weakest region. Instead of using the modal analysis stage to identify weak regions and using averaged stress or displacement over weak regions as the target quantity to optimize, the same optimization process can be run directly, treating each tetrahedron as a potential weak region.
The weakness map is defined as a scalar field on the surface mapping each point to the maximal principal stress at this point obtained by approximate optimization. Using our method yields a partial weakness map on the union of all weakness regions we consider.
While speedups are already dramatic for extremely small element counts, the higher asymptotic complexity of brute force causes a rapidly increasing speedup for larger models.
Dependence on tetrahedral mesh resolution. To keep the cost of computation low, especially in the context of interactive applications or processing large number of objects at a printing facility, using coarse tetrahedral meshes is desirable. As
Drop test. To verify the method for brittle materials, a randomized force test was done by dropping printed models onto horizontal pegs. These models were dropped from 1 m high, ensuring a nearly random impact orientation and force application. The test setup is pictured in
The testing results, displayed in
For the objects printed in ductile materials, a different test was performed. The shapes were placed into a cardboard box filled with packaging material and applied pressure to the box'"'"'s exterior. This pressure permanently deformed the models inside. Photographs were obtained of the deformed models in a registered position and compared them to the 3D model from which they were printed. Good agreement was observed with the computed map of maximum displacements, i.e., the map similar to the weakness map, but for the displacement maximization problem.
The approach described was compared to certain prior art, as they aim to predict the loads that a printed model is likely to experience. In this invention a more specific force model was used: pinch grips. The prior art presents an empirical model to predict how the object will be gripped with two fingers. There are many designs for which such a grip does not capture typical use cases or mishaps.
It was observed that the “cup” model herein produces a triangle of forces (perhaps at the expense of higher stress) rather than a pair of opposite forces. One possible reason for this is the pressure bound requiring the force to be distributed over a larger area.
Though our pipeline has not yet reached interactive speeds, it is already fast enough to be included in a 3D printing pipeline, For the sizes of models most commonly sent to 3D printing services (see distribution in
Analyzing the algorithm'"'"'s scaling behavior is complicated by its dependence on structural properties—a separate linear program is run for each weak region that is extracted. To make sense of the timings, they have been separated by stage. Modal analysis and weak region extraction are run only once per model;
Material parameters defining the elasticity tensor C must be measured for each of the 3D printers'"'"' materials. It was observed that the computed maximal stress does not depend on the magnitude of the Young modulus in the isotropic case. However, in the anisotropic case, it does depend on the ratios of directional elasticity moduli, which can be significant (see
Young'"'"'s modulus ratio measurements were taken and are shown for three different 3D printing materials that were used to compute our simulation'"'"'s elasticity parameters. In addition, the extent is described to which various materials match our assumptions on stress strain linearity and what accuracy one can expect from predictions of the maximal stress to tensile strength ratio. In all cases, a Poisson ratio of 0.3 was assumed.
Three materials were used in 3D printing: nylon (PA 2200 by EOS Electro Optical Systems), “sandstone” (zpI50 used in the Zprinter series by 3D Systems), and green state stainless steel (420SS powder bound with proprietary binder used by ExOne). They also represent different classes of materials (brittle vs. ductile, isotropic vs. anisotropic).
To determine their properties a three point bending test was conducted consistent with ASTM standard D5032 ([ASTM 20077]), using the Instron 5960 universal testing machine with a ±1OON load cell and a support span of 40 mm.
Among the three materials tested, green state stainless steel fits the definition of brittle material the best. Stress grows linearly with strain for all samples tested (
Models printed in nylon are known to withstand a large range of deformations.
The most complex material tested was the “sandstone.” Though, like green state metal, it has a relatively low tensile strength, it exhibits a significant plastic region (see
An efficient approximation method and system have been described for determining worst-case loads for a geometric object based on its geometry and material properties only. The method is quite reliable (it relies on a linear solver, an eigensolver, and a convex solver, which all can provide convergence guarantees), efficient, and approximates well worst-case stress and displacement distributions.
Only linear elasticity is considered, and the optimized solution at the second stage may not match reality for certain large plastic deformations. We note, however, that the robustly obtained approximate solution can serve as a starting point for a nonlinear solver that would find a more accurate local-minimum solution nearby. More generally, 3D printed materials exhibit a broad range of complex behaviors, some of which may exhibit considerable variation even for the same printing process. Using computational models reflecting material complexity and uncertainty is an important future direction. From the point of view of robustness of the method, tetrahedral mesh generation is the bottleneck. Meshless techniques may yield a fully robust pipeline using only surface geometry as input.
With respect to the use of substantially any plural and/or singular terms herein, those having skill in the art can translate from the plural to the singular and/or from the singular to the plural as is appropriate to the context and/or application. The various singular/plural permutations may be expressly set forth herein for the sake of clarity.
The foregoing description of illustrative embodiments has been presented for purposes of illustration and of description. It is not intended to be exhaustive or limiting with respect to the precise form disclosed, and modifications and variations are possible in light of the above teachings or may be acquired from practice of the disclosed embodiments. It is intended that the scope of the invention be defined by the claims appended hereto and their equivalents.