Method and system for numerically simulating foam-like material in finite element analysis

[0001] This application is a continuation-in-part of a U.S. patent application Ser. No. 10/669,047 filed Sep. 23, 2003, the entire contents of which are incorporated herein by reference.

[0002] 1. Field of the Invention

[0003] The present invention generally relates to method, system and software product used in finite element analysis, more particularly to numerical simulation of the structural behaviors of highly compressible material such as foam under loading conditions.

[0004] 2. Description of the Related Art

[0005] Finite element analysis (FEA) is a computerized method widely used in industry to model and solve engineering problems relating to complex systems. FEA derives its name from the manner in which the geometry of the object under consideration is specified. With the advent of the modern digital computer, FEA has been implemented as FEA software. Basically, the FEA software is provided with a model of the geometric description and the associated material properties at each point within the model. In this model, the geometry of the system under analysis is represented by solids, shells and beams of various sizes, which are called elements. The vertices of the elements are referred to as nodes. The model is comprised of a finite number of elements, which are assigned a material name to associate with material properties. The model thus represents the physical space occupied by the object under analysis along with its immediate surroundings. The FEA software then refers to a table in which the properties (e.g., stress-strain constitutive equation, Young's modulus, Poisson's ratio, thermo-conductivity) of each material type are tabulated. Additionally, the conditions at the boundary of the object (i.e., loadings, physical constraints, etc.) are specified. In this fashion a model of the object and its environment is created.

[0006] FEA is becoming increasingly popular with automobile manufacturers for optimizing both the aerodynamic performance and structural integrity of vehicles. Similarly, aircraft manufacturers rely upon FEA to predict airplane performance long before the first prototype is built. Rational design of semiconductor electronic devices is possible with Finite Element Analysis of the electrodynamics, diffusion, and thermodynamics involved in this situation. FEA is utilized to characterize ocean currents and distribution of contaminants. FEA is being applied increasingly to analysis of the production and performance of such consumer goods as ovens, blenders, lighting facilities and many plastic products. In fact, FEA has been employed in as many diverse fields as can be brought to mind, including plastics mold design, modeling of nuclear reactors, analysis of the spot welding process, microwave antenna design, simulating of car crash and biomedical applications such as the design of prosthetic limbs. In short, FEA is utilized to expedite design, maximize productivity and efficiency, and optimize product performance in virtually every stratum of light and heavy industry. This often occurs long before the first prototype is ever developed.

[0007] The finite element analysis method is described in detail by Thomas J. R. Hughes in “The Finite Element Method” (1987), published by Prentice-Hall, Inc., New Jersey, which is incorporated herein by reference in its entirety. Generally, FEA begins by generating a finite element model of a system. In this model, a subject structure is reduced into a number of node points which are connected together to form finite elements. Governing equations of motion are written in a discrete form, where the motions of each node point are the unknown part of the solution. A simulated load or other influence is applied to the system and the resulting effect is analyzed using well known mathematical methods.

[0008] FEA software can be classified into two general types, implicit FEA software and explicit FEA software. Implicit FEA software uses an implicit equation solver to solve a system of coupled linear equations. Such software is generally used to simulate static or quasi-static problems. Explicit FEA software does not solve coupled equations but explicitly solves for each unknown assuming them uncoupled. Explicit FEA software usually uses central difference time integration which requires very small solution cycles or time steps for the method to be stable and accurate. Explicit FEA software is generally used to simulate short duration events where dynamics are important such as impact type events.

[0009] One of the most challenging FEA tasks is to simulate an impact event such as car crash. The highly non-linear behavior of the structural materials must be numerically simulated accurately, realistically and efficiently. A number of materials such as steel, aluminum, foam and rubber that used in an automobile must be included in such FEA software. Many components used in an automobile are made of materials such as foam and rubber, the simulation of the structural responses of these materials becomes very important for the overall accuracy of an analysis.

[0010] Traditionally the structural responses of highly compressible material such as foam have been numerically simulated using the Ogden strain energy function W. A brief summary of the Ogden strain energy function and corresponding engineering/nominal stress σ₀ and Cauthy/true stress σ functions are listed in FIG. 1. More details of the Ogden energy function is described by R. W. Ogden in Chapter 7 of the book titled: “Non-linear Elastic Deformations” (1984), published by Ellis Horwood Limited, United Kingdom, which is incorporated herein by reference in its entirety. A number of commercially available FEA software includes these approaches to simulate foam-like material. For example, LS-DYNA, a general purpose three-dimensional non-linear large-deformation FEA software from Livermore Software Technology Corporation, is capable of simulating foam-like material using the Ogden energy function.

[0011] Today, there are a number of practical problems associated with the simulation of foam-like material in FEA. To implement the Ogden energy function properly in the FEA software requires engineers to spend a tremendous amount of time to prepare experimental data and then convert them into a set of coefficients to fit a polynomial Ogden function for FEA software. Due to highly non-linear characteristics of this polynomial function, the inexactly fitted function has often resulted. This leads to a lengthy iterative trial-and-error process of modifying the input coefficients to match the behavior of foam-like material in the real world. It is therefore desirable to have a new method to numerically simulate foam-like material more efficiently and effectively.

[0012] This section is for the purpose of summarizing some aspects of the present invention and to briefly introduce some preferred embodiments. Simplifications or omissions may be made to avoid obscuring the purpose of the section. Such simplifications or omissions are not intended to limit the scope of the present invention.

[0013] The present invention discloses a system, method and software product to numerically simulate highly compressible material such as foam with a pragmatic approach. Under the assumption of uniaxial loading and isotropic material, a method of calculating stress function f(λ) is developed as shown in FIGS. 2 and 3. The method eliminates the requirement of fitting a polynomial function; instead f(λ) is calculated using only a set of the stress-strain curve from a uniaxial loading test. Based upon the data of this stress-strain curve, the stress function f(λ) is precomputed for a range of equally spaced stretch ratios λ. These results are stored in a computer's RAM as a lookup table and retrieved later during the solution phase of FEA. The present invention does not require time-consuming trial-and-error process of fitting polynomial coefficients. This is a huge advantage over the existing method. The stresses corresponding to a particular strain or stretch ratio can be interpolated rapidly from the pre-computed table directly during the lengthy solution phase.

[0014] According to one aspect of the present invention, during the solution phase of FEA, the stresses at each integration point of the finite elements is computed for each time step. In order to compute the stresses, the eigensolution of a 3×3 stretch tensor is computed first. The resulting eigenvalues are the three main stretch ratios for the element. The corresponding engineering/nominal stress is then looked up or interpolated. These nominal stresses are converted into the true stresses and transformed to the global coordinate system. If the rate dependent viscoelastic effect is needed, the nominal stresses are looked up from a set of curves instead of a single curve.

[0015] Other objects, features, and advantages of the present invention will become apparent upon examining the following detailed description of an embodiment thereof, taken in conjunction with the attached drawings.

[0023] In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention. However, it will become obvious to those skilled in the art that the present invention may be practiced without these specific details. The descriptions and representations herein are the common means used by those experienced or skilled in the art to most effectively convey the substance of their work to others skilled in the art. In other instances, well-known methods, procedures, components, and circuitry have not been described in detail to avoid unnecessarily obscuring aspects of the present invention.

[0024] Reference herein to “one embodiment” or “an embodiment” means that a particular feature, structure, or characteristic described in connection with the embodiment can be included in at least one embodiment of the invention. The appearances of the phrase “in one embodiment” in various places in the specification are not necessarily all referring to the same embodiment, nor are separate or alternative embodiments mutually exclusive of other embodiments. Further, the order of blocks in process flowcharts or diagrams representing one or more embodiments of the invention do not inherently indicate any particular order nor imply any limitations in the invention.

[0025] To facilitate the description of the present invention, it deems necessary to provide definitions for some terms that will be used throughout the disclosure herein. It should be noted that the definitions following are to facilitate the understanding and describe the present invention according to an embodiment. The definitions may appear to include some limitations with respect to the embodiment, the actual meaning of the terms has applicability well beyond such embodiment, which can be appreciated by those skilled in the art:

[0026] FEA stands for Finite Element Analysis.

[0027] Implicit FEA refers to Kx'F, where K is the global stiffness matrix, x is the unknown displacement array and F is the global force array.

[0028] Explicit FEA refers to Ma=F, where M is the mass array, a is the acceleration array and F is the global force array.

[0029] Time domain analysis refers to a FEA that simulates a physical phenomenon in time domain.

[0030] Dynamic analysis refers to a FEA that accounts for the mass and inertia effects of a structure. In general, there are two classes of dynamic analysis: time domain and frequency domain.

[0031] Solution cycle and cycle are used interchangeably to refer to the numbered solution states starting with cycle 0 at time 0.

[0032] The time step refers to an interval of time between subsequent cycles.

[0033] Loading Condition is defined as the external load acting on a structure.

[0034] Uniaxial Load is either tension or compression in one direction.

[0035] Longitudinal Direction is the direction, in which the uniaxial force applies to a material.

[0036] Lateral Direction is the orthogonal direction to the longitudinal direction.

[0037] Tension is a force to pull or stretch a material.

[0038] Compression is a force to squeeze or press a material.

[0039] Stress is defined as the intensity of force or force per unit area of the material.

[0040] Stretch Ratio λ is defined as the ratio between the deformed length and the original length of a material.

[0041] Strain ε is a non-dimensional quantity representing the elongation or shrinkage per unit length of a material under loading condition. The Strain and Stretch Ratio is related with the following equation: ε=λ−1.

[0042] Strain Rate is defined as the amount of strain that develops in one unit of time.

[0043] Poisson's Ratio ν is defined as a function of the strain in lateral direction to the strain in the longitudinal direction. Strain Energy is potential energy generated from the force acting on a material. In the elastic case, the strain energy is defined as W=0.5*ε*σ.

[0044] Deformation Gradient Matrix/Tensor is defined as the local nature of the deformation of a solid body.

[0045] Principal Stretch Ratio is defined as the stretch ratio in the principal loading direction.

[0046] Eigensolution is a mathematical method to determine the natural vibration mode shapes and amplitudes of a structure.

[0047] Eigenvalues are the results obtained from the eigensolution.

[0048] Elastic material is a material that recovers to its original shape after unloading (i.e., the load is removed).

[0049] Hyperelastic material is a material has one-to-one relationship between stress and strain for a large deformation. The stress value can be determined by an energy function such as Stress function f(λ).

[0050] Quasistatic loading is defined as the static force dominates comparing to inertia and damping forces.

[0051] Rate dependent under dynamic loading means that under dynamic loading the material is no longer quasistatic; for the same values of strain, different values of stress may be obtained depending upon the loading velocity or strain rate. Rate dependent under dynamic loading and viscous material are interchangeably used.

[0052] Nominal stress or Engineering stress, interchangeably used herein, is defined as the force divided by undeformed cross section area.

[0053] True stress or Cauchy stress, interchangeably used herein, is defined as the force divided by current or deformed cross section area.

[0054] Local coordinate system or element coordinate system, interchangeably used herein, is a coordinate system used to define element structural responses such as true stress and nominal stress.

[0055] Global coordinate system is a coordinate system used to define the structure in FEA.

[0056] Embodiments of the present invention are discussed herein with reference to FIGS. 1-7. However, those skilled in the art will readily appreciate that the detailed description given herein with respect to these figures is for explanatory purposes as the invention extends beyond these limited embodiments.

[0057] Referring now to the drawings, in which like numerals refer to like parts throughout several views. FIG. 1 shows the structural response functions in the original Ogden highly compressible material model and derivation of nominal stress formula as function of longitudinal stretch ratio λ₁ based on the assumption of uniaxial loading and isotropic material. Ogden strain energy function W for compressible material is listed in 110. In these functions, λ₁, λ₂ and λ₃ represent the stretch ratios in three principal directions of a compressible material. The variable J is the volumetric ratio or the deformed volume divided by the undeformed volume. The material coefficient n is related to Poisson's ratio of the material. Coefficients μ and α are the material constants. The classic strain energy function W in 110 is the foundation to numerically simulate highly compressible material using a digital computer. The formula for true stress σ is listed in 120 and nominal stress σ₀ is listed in 130. Under the uniaxial tension or compressible loading, the coefficient n is determined in the formulas listed in 140. Under the isotropic material assumption, the two lateral stretch ratios λ₂ and λ₃ are equal to each other. The relationship between stretch ratio λ₁ in the uniaxial loading direction and lateral stretch ratio λ₃ is also listed in 140. Substitute λ₂ and λ₃ with λ₁ in the nominal stress formula in 130, the nominal stress formula can be written as a function of the longitudinal stretch ratio λ₁ as shown in 150. It is evident that element stresses can easily be calculated, once the coefficients μ and α have been determined; however, the determination of these coefficients for the polynomial function is very difficult in practice.

[0058] FIGS. 2 and 3 in combination show the derivation of stress function f(λ). The original stress function f(λ) in polynomial form is defined in 210. Substitute f(λ) into the nominal stress σ₀ formula in 150, the resulting formula is shown in 220. In 230, Poisson's ratio ν is defined as function of material constant n. After substituting ν in formula 220, the new formula is listed as formula 240. Formula 250 is a result of replacing λ with λ^−ν in formula 240. The next formula 260 is similarly obtained by substituting λ^(−υ)2 into λ in formula 240. The process repeats for the rest of higher order formulas using the same substitution method.

[0059] By rearranging formulas 240, 250, 260 and the rest of higher order formulas, the resulting equation for stress function f(λ) is rewritten as listed in 310. The equation 310 is rewritten as a different form equation as shown in 320. Wherein the value of stress function f(λ) for the stretch ratio of interest λ equals to the summation of a sequence of λ^(−υ)j σ₀(λ^1−υjj), where j is an integer related to the j-th term of the sequence, ν is Poisson's ratio of the compressible material, and σ₀(λ^[−υ]j) is the stress value at stretch ratio λ^(−υ)j. Using the relationship between strain ε and stretch ratio λ in 330, the equation 320 transforms to equation 340, which can be evaluated solely with the stress values corresponding to the strains from the strain-stress curve obtained via a uniaxial tension/compression test of the material of interest. Finally, the nominal stress σ₀ and true stresses σ are rewritten and dependent only on the stress function f(λ) listed in formulas 350 and 360, respectively

[0060] Referring now to FIG. 4, an engineering or nominal stress-strain curve 420 is illustrated. The curve represents a set of experimental data for a compressible material under uniaxial loading condition; both tension and compression tests results are included. According to one embodiment of the present invention, the engineer needs only to prepare the experimental data such as this exemplary curve 420 as input to simulate compressible material in FEA software. In order to capture the stretch ratio of interest, an exemplary input requires a range of strain values shown in 430. If the viscous effect of compressible material is under consideration, the exemplary curve 420 may be replaced by a family of stress-strain curves resulting from dynamic experiments.

[0061] With reference now to FIG. 5, it shows a flow chart 500 for computing the infinite series in equation 340. Because the order of magnitude decreases drastically from one term to the next, the sequence converges very rapidly. According to one embodiment of the present invention, the flow chart 500 of evaluation of function f(λ) is summarized in FIG. 5. At 510, f(λ) is assigned a value equals to λ multiplied by σ₀(λ−1) for a given λ. The stress σ₀ value at strain ε₀ or λ−1 is from engineering test strain-stress data (e.g., the curve defined in the user input phase of FEA software). At 520, λ is stored into a variable λ_old. A new variable λ_new is stored as λ_old^−ν at 530, where ν is Poisson's ratio of the compressible material. At 540, a comparing test is performed for the absolute value of (λ_new−1) being less than or equal to a threshold value ξ. In one embodiment, ξ is set to 0.01. If the test succeeds, the computation of f(λ) has finished, the rest of the terms in the infinite series is too small to affect the final result of the computation. If the test fails, another stress value at strain λ_new−1 is multiplied by λ_new and accumulated into function f(λ) at 550. At 560, the value of λ_new is stored into λ_old. The process goes back to 530 until the computation finishes.

[0062] An exemplary FEA software implementation 600 of the present invention is shown in FIG. 6. At the input phase 610 of FEA software, at least one set of engineering or nominal stress-strain data for rubber material under uniaxial loading test is read in as a group of corresponding stress-strain pairs. If the loading condition is quasistatic, only one set of stress-strain curve is required. If the rate dependent dynamic loading condition is under consideration, then a family of stress-strain curves representing different strain rate are required. Based on the input stress-strain constitutive curve, a plurality of f(λ) is then generated for a range of λ at regular intervals at the initialization phase 620 using the process listed in flow chart 500 as shown in FIG. 5. These pre-computed f(λ) are stored in computer's RAM as a stress function table for later use in the solution phase. Next in 630, the FEA software calculates structural responses at each integration point of every element for each time step. According to one embodiment, for the element with compressible material such as foam, a subroutine handling foam material is called. Within the subroutine, the deformation gradient matrix is gathered first at 640. For a solid element, either the left or the right 3×3 stretch tensors may represent the deformation gradient matrix. The deformation gradient matrix may be corresponded to the stretch ratio λ or the stretch ratio squared λ² depending upon the solution scheme employed. At 650, an eigensolution is solved for the deformation gradient tensor. The principal stretch ratio can be calculated from the resulting eigenvalues at 660. With these stretch ratios, at 670, the nominal stresses are then calculated using the corresponding f(λ) values stored in the initialization phase. Because all f(λ) values of interest are pre-computed and stored as tabulated form, it is very computational efficient to perform a table lookup and interpolation. In the case of rate dependent dynamic loading, strain rate at the element integration point is used as basis to select two closest curves to determine the corresponding rate dependent stress. The true stresses are calculated and transformed into the global coordinate system at 690.

[0063] With reference now to FIG. 7, a block diagram illustrates a computing device 700 in which the present invention may be implemented, and in which code or instructions implementing the processes of the present invention may be located. The exemplary computer system in FIG. 7 is discussed only for descriptive purposes. It should not be considered a limitation of the invention. Although the following descriptions related to a particular computer system, the concepts apply equally to other computer systems that are dissimilar to that shown in FIG. 7.

[0064] Computer system 700 includes a processor 710 and main random access memory (RAM) 720 connecting to a local bus 705 through a bridge 715. Additional connections to local bus 705 may be made through direct component interconnection or through add-in boards. In the depicted example, network adapter 725, small computer system interface (SCSI) adapter 730, and expansion bus interface 735 are directly connected to local bus 705. In contrast, audio adapter 740, graphics adapter 745, and video adapter 750 are connected to local bus 705 by add-in boards inserted into expansion slots. Expansion bus interface 735 provides a connection for a keyboard and mouse adapter 755, modem 760, and additional memory 765. SCSI adapter 730 provides a connection for hard disk drive 770, tape drive 775, and CD-ROM drive 780.

[0065] In order to communicate with other computer systems via a network, the computer system 700 connects to the network via network adapter 725. The network, Internet or intranet, connects multiple network devices utilizing general purpose communication lines.

[0066] Those of ordinary skill in the art will appreciate that the hardware shown in FIG. 7 may vary depending on the implementation. Other internal hardware or peripheral devices, such as flash ROM (or equivalent nonvolatile memory) or optical disk drives and the like, may be used in addition to or in lieu of the hardware depicted in FIG. 7. Also, the processes of the present invention may be applied to a multiprocessor computer system. In general, Computer system 700 is controlled and coordinated by operating system (OS) software, which performs tasks such as process scheduling, memory management, networking and I/O services. Exemplary OS includes Linux™, Microsoft Windows™.

[0067] Although an exemplary embodiment of invention has been disclosed, it will be apparent to those skilled in the art that various changes and modifications may be made to achieve the advantage of the invention. It will be obvious to those skilled in the art that some components may be substituted with another component providing same function. The appended claims cover the present invention.