Method and system for simulating fracture propagation in brittle materials

A system, method and software product for simulating fracture propagation in brittle materials is disclosed. Using the fracture energy release rates and the critical separations of the three fundamental fracture modes obtained independently from experiments, and a user defined normalized curve representing the traction-separation law, the present invention simulates the fracture propagation in a non-linear dynamic analysis. By using either dimensional or dimensionless effective separation approach, the actual tractions are calculated at each solution cycle based on user's choice. In addition, loading, unloading and reloading conditions are also continuously traced to ensure the proper constitutive equations are used. Three fracture modes are seamlessly mixed and integrated.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention generally relates to a method, system and software product used in three-dimensional non-linear finite element analysis of a structure, more particularly to simulating fracture propagation in brittle materials.

2. Description of the Related Art

Finite element analysis (FEA) is a computerized method widely used in industry to model and solve engineering problems relating to complex systems such as three-dimensional non-linear structural design and analysis. 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 certain points 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 are created.

FEA is becoming increasingly popular with automobile manufacturers for optimizing both the aerodynamic performance and structural integrity of vehicles. Similarly, aircraft manufacturers rely on 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.

On the most challenging FEA tasks is to simulate fracture (e.g., crack, micro-crack) propagation in brittle materials such as glasses, ceramics, and hard composites. Fracture usually begins when stress applied to a material is concentrated at the tip of a micro-crack. When the stress exceeds a critical value, atomic bonds begin to break, elastic energy is released, and new surface is created as the crack propagates in the material. Brittle fracture is not only an annoying everyday experience or a safety hazard, but also an important technological process for the shaping of hard materials. Controlling the brittle fracture of flint-stone was the crucial step into the stone-age and polishing silicon wafers of 300 mm diameter with tolerable height variations of only a few atom spacing is a technological challenge today. Engineers at the beginning of the last century started to investigate brittle fracture processes and soon realized that the mechanical stress in the solid is concentrated at the crack tip. This stress concentration increases with increasing sharpness of the crack. In a brittle material, the crack tip is atomically sharp and, therefore, the material must sustain very high stresses exceeding the nominal fracture strength of engineering materials.

Fracture can be categorized into three fundamental modes:FIG. 1Ashows an opening or normal mode which is designated as Mode-I100A;FIG. 1Bshows a forward shear or sliding mode designated as Mode-II100B; andFIG. 1Cshows a transverse shear or tearing mode designated as Mode-III100C. These three fracture modes can occur separately or in any combination. Fractures in which two or more modes were operative are termed mixed-mode fractures. Today, all of the existing fracture propagation simulation schemes have problems. Some schemes require same fracture energy release rate for all three modes, some assume a simplified linear, bilinear or trilinear traction-separation law, some do not differentiate Mode-II and Mode-III, while others do not enforce the irreversible condition.

Therefore, there is a need to have a set of general cohesive laws that include all three fracture modes interacting with each other for simulating fracture in brittle material.

SUMMARY OF THE INVENTION

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.

The present invention discloses a system, method and software product for simulating fracture propagation in brittle materials. According to one aspect, fracture energy release rate and critical separation of each of the three fundamental modes is obtained and defined. A curve representing the normalized traction-separation law (TSL) is also defined. The area under the normalized TSL curve is calculated. The maximum traction for each of the modes is calculated using the fracture energy release rate and the critical separation of each individual mode. Depending upon which traction-separation law is selected, a set of constitutive equations for traction and separation is derived accordingly for loading and unloading conditions. The actual traction can be calculated at each solution cycle. The loading and the unloading condition are to be monitored and tracked in order to simulate the fracture propagation accurately.

In one embodiment, the present invention is a method for simulating fracture propagation in a brittle material of a structure, said method comprises: defining fracture energy release rate and critical separation of each of the three fracture modes; defining a normalized traction-separation curve; calculating the maximum traction for said each of the three fracture modes; selecting a fracture damage failure criterion; and calculating actual traction of said each of the three fracture modes based upon the criterion selected at each solution cycle.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention discloses a system, method and software product for simulating fracture propagation in brittle material. According to one aspect, fracture energy release rate and critical separation of each of the three fundamental modes is obtained and defined. A curve representing normalized traction-separation law (TSL) is also defined. The area under the normalized TSL curve is calculated. The maximum traction for each of the modes is calculated using the fracture energy release rate and the critical separation of each individual mode. Depending upon which traction-separation law is selected, a set of constitutive equations for traction and separation is derived accordingly for loading and unloading conditions. The actual traction can be calculated at each solution cycle. The loading and unloading condition are monitored and tracked in order to simulate the fracture propagation accurately.

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.

FEA stands for Finite Element Analysis.

Implicit FEA refers to Ku=F, where K is the effective stiffness matrix, u is the unknown displacement array and F is the effective loads array. F is a right hand side loads array while K is a left hand side stiffness matrix. The solution is performed at the global level with a factorization of the effective stiffness matrix, which is function of the stiffness, mass and damping. One exemplary solution method is the Newmark integration scheme.

Explicit FEA refers to Ma=F, where M is the diagonal mass array, ‘a’ is the unknown nodal acceleration array and F is the effective loads array. The solution can be carried out at element level without factorization of a matrix. One exemplary solution method is called the central difference method.

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

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.

Solution cycle and cycle are used interchangeably to refer to the numbered solution states starting with cycle0at time0.

The time step refers to an interval of time between consecutive solution cycles.

Total delamination, total debonding, or total decohesive refer to total damage or complete failure and are interchangeably used in this document.

Damage onset refers to the initiation of damage.

Separation and jump, used interchangeably, refer to the relative displacement between two cohesive points. At the critical separation, the material is in a state of total damage or complete failure, and fracture has occurred.

Referring to the drawings, in which like numerals refer to like parts throughout several views. The present invention may be implemented using hardware, software or a combination thereof and may be implemented in a computer system or other processing system. In fact, in one embodiment, the invention is directed towards one or more computer systems capable of carrying out the functionality described herein. An example of a computer system200is shown inFIG. 2. The computer system200includes one or more processors, such as processor204. The processor204is connected to a computer system internal communication bus202. Various software embodiments are described in terms of this exemplary computer system. After reading this description, it will become apparent to a person skilled in the relevant art(s) how to implement the invention using other computer systems and/or computer architectures.

Computer system200also includes a main memory208, preferably random access memory (RAM), and may also include a secondary memory210. The secondary memory210may include, for example, one or more hard disk drives212and/or one or more removable storage drives214, representing a floppy disk drive, a magnetic tape drive, an optical disk drive, etc. The removable storage drive214reads from and/or writes to a removable storage unit218in a well-known manner. Removable storage unit218, represents a floppy disk, magnetic tape, optical disk, etc. which is read by and written to by removable storage drive214. As will be appreciated, the removable storage unit218includes a computer usable storage medium having stored therein computer software and/or data.

In alternative embodiments, secondary memory210may include other similar means for allowing computer programs or other instructions to be loaded into computer system200. Such means may include, for example, a removable storage unit222and an interface220. Examples of such may include a program cartridge and cartridge interface (such as that found in video game devices), a removable memory chip (such as an EPROM, or PROM) and associated socket, and other removable storage units222and interfaces220which allow software and data to be transferred from the removable storage unit222to computer system200. In general, Computer system200is 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®.

There may also be a communications interface224connecting to the bus206. Communications interface224allows software and data to be transferred between computer system200and external devices. Examples of communications interface224may include a modem, a network interface (such as an Ethernet card), a communications port, a PCMCIA slot and card, etc. Software and data transferred via communications interface224are in the form of signals228which may be electronic, electromagnetic, optical, or other signals capable of being received by communications interface224. These signals228are provided to communications interface224via a communications path (i.e., channel)226. This channel226carries signals (or data flows)228and may be implemented using wire or cable, fiber optics, a phone line, a cellular phone link, an RF link and other communications channels.

The channel226facilitates a data flow228between a data network and the computer200and typically executes a special set of rules (i.e., a protocol) to send data back and forth. One of the common protocols is TCP/IP (Transmission Control Protocol/Internet Protocol) commonly used in the Internet. In general, the communication interface224manages the assembling of a data file into smaller packets that are transmitted over the data network or reassembles received packets into the original data file. In addition, the communication interface224handles the address part of each packet so that it gets to the right destination or intercepts packets destined for the computer200.

In this document, the terms “computer program medium” and “computer usable medium” are used to generally refer to media such as removable storage drive214, a hard disk installed in hard disk drive212, and signals228. These computer program products are means for providing software to computer system200. The invention is directed to such computer program products.

The computer system200may also include an I/O interface230, which provides the computer system200to access monitor, keyboard, mouse, printer, scanner, plotter, and a like.

Computer programs (also called computer control logic) are stored as applications modules206in main memory208and/or secondary memory210. Computer programs may also be received via communications interface224. Such computer programs, when executed, enable the computer system200to perform the features of the present invention as discussed herein. In particular, the computer programs, when executed, enable the processor204to perform the features of the present invention. Accordingly, such computer programs represent controllers of the computer system200.

In an embodiment where the invention is implemented using software, the software may be stored in a computer program product and loaded into computer system200using removable storage drive214, hard drive212, or communications interface224. The application module206, when executed by the processor204, causes the processor204to perform the functions of the invention as described herein.

The main memory208may be loaded with one or more application modules206that can be executed by one or more processors204with or without a user input through the I/O interface230to achieve desired tasks. In operation, when at least one processor204executes one of the application modules206, the results are computed and stored in the secondary memory210(i.e., the hard disk drive212). In one embodiment, the results and/or pictures of fracture propagation simulation based upon the general cohesive material laws are reported to the user via the I/O interface230either as a list or a graph.

In one embodiment, an application module is configured to facilitate defining a normalized traction-separation curve and the fracture energy release rate and the critical separation for each of the three modes. Depending upon what traction-separation law is selected, the actual traction of three modes is calculated at each time cycle. Thereby the fracture propagation is simulated in a time domain. In another embodiment, an application module is configured to facilitate to track the separation parameter, either dimensional or dimensionless, to determine whether the fracture is under loading or unloading. If the former, the actual traction is calculated using the formula for loading condition. Otherwise, the formula for unloading condition is used. One exemplary implementation of this technique is included in a well-known engineering computer program product, LS-DYNA®, offered by Livermore Software Technology Corporation.

FIGS. 3A-3Dshow a plurality of exemplary traction-separation curves300A-D in accordance with one embodiment of the present invention. In this embodiment, the traction σJ302A plotted in the vertical axis represents the traction of stress or force in one of the fracture modes100A-C. And the separation u,304A-D plotted in the horizontal axis represents the separation or jump or relative displacement of that mode. All of the plurality of the curves show the characteristics of an elastic (i.e., reversible) ascending portion from an undeformed state of zero traction and separation (i.e., original state) to the peak traction σJmax312A-D and an inelastic irreversible descending portion from the peak traction to zero traction at the critical separation δJ314A-D. The point at the peak traction σJmax312A-D indicates the damage onset of a particular mode, and the point at the maximum separation indicates the total failure criterion for fracture propagation.FIG. 3Ashows a bi-linear example300A with an intersection at the peak traction σJmax312A.FIG. 3Bshows a tri-linear example300B with a horizontal platform at the level of the maximum traction σJmax312B.FIG. 3Cshows a smooth curve example300C with peak traction occurring at σJmax312C.FIG. 3Dshows an example with multiple linear segments300D with peak traction occurring at σJmax312D. The subscript J represents one of the three modes I, II and III.

FIGS. 4A and 4Bshow a couple of exemplary normalized traction-separation curves400A-B in accordance with one embodiment of the present invention. The dimensionless or normalized separations are defined as:

λJ=uJδJJ⁢⁢ε⁢⁢I,II,III(1)
The horizontal axis represents normalized separation λJ404A-B with value from 0 to 1; the vertical axis represents normalized traction f(λJ)402A-B with maxλf(λJ)=1.

λ=λI2+λII2+λIII2=(uIδI)2+(uIIδII)2+(uIIIδIII)2(2)
where uI, uII, uIIIand δI, δII, δIIIare the separations at each solution cycle and critical separations of modes I100A, II100B and III100C, respectively. ESP λ=0 represents materials in the original undeformed state and the complete failure or total damage of material occurs at ESP λ=1. ESP λ404A-B is computed and its historical maximum λmax424A-B is recorded at each solution cycle to determine whether the material is under loading, unloading or reloading. Before ESP λ404A-B reaches λ0422A-B, the material behavior is elastic and follows the elastic portion of the user defined curve for all conditions. Once ESP λ404A-B reaches λ0422A-B, the damage onset starts. If ESP λ=λmaxand dλ/dt≧0, the material is under loading condition, the material continues to follow the user defined curve. Otherwise, the material is under unloading condition. A cleavage unloading law is assumed for brittle materials, i.e., the unloading follows a secant slope from the state at λmax424A-B to the original state. This implies an irreversible condition where the traction vanishes completely at the origin. At this stage, the material behavior follows the secant slope for both unloading and reloading. The fracture energy release rates or toughnesses, GI, GII, and GIII, the critical separations δI, δII, and δIII(e.g., δJ314A-D), and the normalized traction-separation curve are obtained from experimental data and defined in the input deck by the user. The area Γ under the normalized traction-separation curve is calculated as follows:

Γ=∫01⁢f⁡(λ)⁢⁢ⅆλ(3)
The maximum traction for each of the three modes is calculated as follows:

σI,max=GIΓδIσII,max=GIIΓδIIσIII,max=GIIIΓδIII(4)
The ratios between fracture energy release rates are defined as follows:

αII⁢GIIGI=σIImax⁢δIIσImax⁢δIαIII⁢GIIIGI=σIIImax⁢δIIIσImax⁢δI(5)
The actual tractions of the three modes for loading case are calculated as follows:

σI={f⁡(λ)λ⁢σImax⁢uIδIfor⁢⁢uI≥0(1+p)⁢k⁢⁢σImax⁢uIδIfor⁢⁢uI<0⁢⁢σII=f⁡(λ)λ⁢σIImax⁢uIIδII⁢⁢σIII=f⁡(λ)λ⁢σIIImax⁢uIIIδIII(6)
where p is the penalty stiffness multiplier. The loading condition is determined when ESP404A-B is increasing (i.e., dλ/dt≧0) and λ=λmax, where λmax424A-B is the maximum ESP during the time history up till now. The formulae for unloading case are as follows:

In the following embodiments, a plurality of exemplary functions f(λ) is listed:(1) Smith-Ferrante's universal binding law: f(λ)=nλe1-nλ, where n is chosen to satisfy ne1-n≦ε.(2) Cubic polynomial:

For the embodiments using dimensional effective separation.FIG. 5Ashows an exemplary traction-separation curve500used in the power law or the Benzeggagh-Kenane (B-K) law failure criterion in accordance with one embodiment of the present invention. In this embodiment, the vertical axis represents the effective traction σ502and the horizontal axis the effective separation δ504. The separation δ0512is the minimum separation that corresponds to the maximum effective traction σ0522and the total failure occurs at the critical separation δc524. Although a tri-linear curve is shown inFIG. 5A, other type of curves may also be used.

The effective opening separation is defined as the resultant separation δ550as shown inFIG. 5B. The actual separations uI, uII, uIII540A-C of the three fracture modes100A-C are the components of the resultant separation at each solution cycle. Their relationship in a three-dimensional space is depicted inFIG. 5Band can be expressed physically by the following formulae:
δ=√{square root over (uI2+uII2+uIII2)}uI=δcosθ |uII|=δsinθcosω |uIII|=δsinθsinω  (8)
where•is the McCauley bracket defined asx=max(0,x), θ552and ω554are angles defined inFIG. 5B, and δ550is total separation (i.e., resultant of the actual separations.
Before the onset of fracture damage, the tractions for the three fracture modes can be calculated as follows:

σI=σImax⁢f⁡(λI)σII=σIImax⁢f⁡(λII)σIII=σIIImax⁢f⁡(λIII)σI=σImax⁢f⁡(uIδI)σII=σIImax⁢f⁡(uIIδII)σIII=σIIImax⁢f⁡(uIIIδIII)(9)
The damage onset is predicted using a quadratic failure criterion as follows:

(σIσImax)2+(σIIσIImax)2+(σIIIσIIImax)2=1(10)
which can be rewritten as an equation in terms of effective separation:

f2⁡(δcosθδI)+f2⁡(δsinθcosωδII)+f2⁡(δsinθsinωδIII)=1(11)
The solution of the above equation δ0falls in the interval [0,δ0] and is unique, where:

Fracture propagation is predicted using either of the following failure criteria: a) the power law and b) the Benzeggagh-Kenane (B-K) law. In the former case, the criterion for predicting complete decohesion or total damage is expressed as follows:

F⁡(x)=∫0x⁢f⁡(λJ)⁢⁢ⅆλJ.
A unique solution δcfalls in the interval [0,δc], where:

σI=σImax⁢uIδIσII=σIImax⁢uIIδIIσIII=σIIImax⁢uIIIδIII(15)a={δcδ⁢f⁡(λoδo⁢δ)if⁢⁢0≤δmax<δoδcδmax⁢f⁡(δmax-δoδc-δo+λo⁢δc-δmaxδc-δo)if⁢⁢δo≤δmax<δc0if⁢⁢δmax≥δc
The above constitutive equation implies the following damage variable:

In the case of the B-K law, the constitutive equation is the same as that of the power law. The difference is that the value δc524is calculated from the following B-K failure criterion:

FIG. 6shows a flowchart or process600of simulating fracture propagation using general cohesive material laws in accordance with one embodiment of the present invention. The process600, which is preferably understood in conjunction with the previous figures, may be implemented in software, hardware, or a combination of both. In this embodiment, the process600starts by defining fracture energy release rates (i.e., toughnesses GI,GII,GIII) and the critical separation δI, δII, and δIIIof the three fracture modes100A-C and a normalized traction-separation law (TSL) (e.g.,FIGS. 4A-4D) at602. The process600calculates the area Γ under the normalized TSL curve at604using eq. (3). Then the maximum tractions σImax, σIImax, σIIImaxfor each of the three modes are calculated at606using equation (4). The process600moves to a test608, which determines whether the user has chosen to use the dimensionless traction-separation law. If the test608is true, the process600performs the following steps for calculating actual traction forces of each of the three modes. At610, the process600calculates the dimensionless effective separation parameter λ using eq. (2). It is noted that λ is calculated at each solution cycle, since the actual separations can be constantly changing as time progresses. Next at612, the process600keeps tracking the maximum value of λ in order to determine whether the structure is undergoing loading or unloading at614. When λ is the historical maximum and dλ/dt is greater than or equal to zero, then the test614is true. The process600follows the “yes” branch to616, in which the actual tractions for each of the three modes are calculated using the loading case eq. (6) before the process600ends for a particular solution cycle. Otherwise if the test614is false, the process600calculates the actual tractions using the unloading case eq. (7) at618before ending.

Referring back to test608, when the dimensional traction-separation law is chosen, the result of the test608is false. The process600follows the “no” branch to620, in which a dimensional effective separation δ is calculated as the resultant of the actual separations uI, uII, uIIIof each of the three fracture modes using eq. (8). As a result, the relationship between the dimensional effective separation δ and the actual separations uI, uII, uIIIis also expressed in eq. (8). At621, the process600calculates the effective separation at damage onset δ0by solving eq. (10) or (11). It is noted that δ0is solved every solution cycle to check the condition of damage onset. Then the process600moves to another test622, in which the failure criterion is determined. There are two options for a user to choose: a) the power law and b) the B-K law. When the power law is chosen, the test622is true; the process600follows the “yes” branch to624. The power law effective separation δc524at the complete failure (i.e., delamination, debonding, decohesive) is calculated in accordance with eq. (13). Otherwise, the process600goes to626, in which the B-K law effective separation δc524is calculated using eq. (17). Finally at628, the process600calculates the actual tractions σI, σII, σIIIfor each of the three modes using the values δ0522and δc524in accordance with eq. (15) before ending at each solution cycle.

Although an exemplary embodiment of the 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 the same function. The appended claims cover the present invention.