METHOD FOR ESTIMATING STRESS INTENSITY FACTORS AND METHOD FOR CALCULATING ASSOCIATED SERVICE LIFE

Abstract

The invention relates to a method for estimating (100) stress intensity factors (FIC) in a numerically modelled part, in the context of a fatigue crack propagation model, comprising the following steps: (E2): obtaining, from a numerical model of the part to be analysed (20), a plurality of values simulated at various points of the part to be analysed (20), (E3): determining a converted value (K) of the effective stress intensity amplitude corresponding to a straight front of a planar crack, as well as a converted length (a) of the crack corresponding to a planar crack with a straight front, said converted values being determined by equalisation of the energy dissipated in the three-dimensional crack of the numerical model and the energy dissipated in the crack of a standard model of a planar crack with a straight front, (E4): interpolating converted values of effective stress intensity factor amplitude between two successive converted lengths (a), (E5): storing the converted values of effective stress intensity factor amplitude interpolated in this way.

Claims

1. A method for estimating (100) the stress intensity factor in a numerically modeled part, within the context of a fatigue crack propagation model, implemented by a system (10) comprising a data processing unit (12), the method for estimating (100) comprising the following steps: (E2): obtaining from a numerical model of the part to be analyzed (20), a plurality of values simulated at different points of the part to be analyzed (20), said plurality of simulated values comprising, for different simulated propagation increments a three-dimensional crack which appears on the numerical model of the part, a set of simulated values of the stress intensity factor at associated points, the position of these points, and data relating to the cracked surface of the crack, (E3): determining, for each propagation increment and for the set of simulated values therefor, a converted amplitude value of the equivalent global effective stress intensity factor corresponding to that of a plane crack with a straight front and a converted length of the equivalent crack considered, said converted values being determined by equalizing the energy dissipated in the three-dimensional crack of the numerical model and the energy dissipated in the crack of a standard model of a plane crack with a straight front, the energies themselves being determined according to the stress intensity factors, (E4): Interpolating converted amplitude values of equivalent global effective stress intensity factor between two successive converted crack lengths, (E5): Storing the converted amplitude values of equivalent global effective stress intensity factor thus interpolated with the associated crack lengths.

2. The method for estimating as claimed in claim 1, wherein the interpolation step (E4) implements a piecewise linear interpolation.

3. The method for estimating as claimed in claim 1, wherein the interpolation step (E4) implements an interpolation by curvature energy minimization.

4. The method for estimating as claimed in claim 1, wherein the interpolation step (E4) comprises the following substeps: (E41) interpolation using a linear interpolation, (E42) interpolation using an interpolation by curvature energy minimization, the two interpolations being interchangeable, (E43) calculation of at least one magnitude representative of a difference in amplitude values of the equivalent global effective stress intensity factor between the two interpolations, said difference being calculated strictly between two values of lengths of the crack corresponding to two successive increments, (E44) comparison of said magnitude representative of the difference with a predetermined threshold, (E45) if the magnitude representative of the difference is above the threshold, generation of an instruction to calculate from the numerical simulation of the part to be analyzed (10) stress intensity factor values and values of the position of the crack, between the two successive increments.

5. The method for estimating as claimed in claim 1, comprising, prior to the step of obtaining (E2), a step of calculating (E1) implementing the finite element simulation at successive increments of the evolution of the crack in the part to be analyzed (10), said simulation being performed by the processing unit (12).

6. The method for estimating as claimed in claim 1, wherein the numerical data retrieved in the step of obtaining (E2) the data are obtained by a finite element or extended finite element simulation.

7. The method for estimating as claimed in claim 1, wherein, for the step of determining (E3), a first data subset makes it possible to find the effective amplitude of the energy release rate G via the following relationships and their equivalence in terms of amplitude: $G=\frac{1}{E^{*}}.Math.\left(K_I^2+K_{\mathrm{II}}^2\right)+\frac{1}{}.Math.K_{\mathrm{III}}^2$ $E^{*}=E.Math..Math.\mathrm{in}.Math..Math.\mathrm{plane}.Math..Math.\mathrm{stresses}$ $E^{*}=E\left(1-{}^2\right).Math..Math.\mathrm{in}.Math..Math.\mathrm{plane}.Math..Math.\mathrm{strains}$ ${\left(.Math..Math.K_{\mathrm{eff}}\left(s\right)\right)}^2=.Math..Math.G_{\mathrm{eff}}\left(s\right).E^{*}$ Where K.sub.I, K.sub.II and K.sub.III are the stress intensity factor SIF coefficients corresponding to the crack opening, plane shear and anti-plane shear modes, respectively, K.sub.eff(s) is the amplitude of the effective stress intensity factor, G.sub.eff is the amplitude of the effective energy release rate, s is the curvilinear abscissa along the crack front, E is Young's modulus, E* is the equivalent Young's modulus, is the shear modulus, v is Poisson's ratio, a second data subset makes it possible to find the cracked surface increment dSurf, these formulations make it possible to obtain the value of the dissipated energy:
Dissipated energy=.sub.crack frontG.sub.eff(S).Math.dSurf(s).Math.dl wherein this dissipated energy is equalized with the dissipated energy of the crack of a standard model of plane crack with a straight front, being expressed in the form:
Dissipated energy.sup.glob=F.sub.eff.sup.glob.Math.dSurf.sup.glob.Math.L.sub.front Where the notation glob relates to this standard model, and L.sub.front is the length of the crack front, wherein, thanks to modeling using a Paris law with the Elber crack closure parameters, the link is made between dSurf.sup.glob and G.sub.eff.sup.glob $\frac{{\mathrm{dSurf}}^{\mathrm{glob}}}{\mathrm{dN}}=C.{\left(\sqrt{.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}.E^{*}}\right)}^n$ where C and n are Paris coefficients, wherein, by equalizing the two dissipated energies, the following is obtained $.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}={\left(\frac{{}_{\mathrm{crack}.Math..Math.\mathrm{front}}.Math..Math..Math.G_{\mathrm{eff}}\left(s\right).\mathrm{dSurf}\left(s\right).\mathrm{dl}}{C.{\left(E^{*}\right)}^n.Math.L_{\mathrm{front}}}\right)}^{\frac{1}{n+1}}={\left(\frac{\mathrm{Dissipated}.Math..Math.\mathrm{energy}}{C.{\left(E^{*}\right)}^n.Math.L_{\mathrm{front}}}\right)}^{\frac{1}{n+1}}$ wherein, by means of the relationships between Geff and Keff, K.sub.eff.sup.glob is determined and wherein Surf.sup.glob is determined, i.e. the converted length of the equivalent crack, wherein the first data subset contains the stress intensity factor values at associated points on the crack front and their respective associated position, and the second data subset contains the data relating to the cracked surface.

8. (canceled)

9. A method for evaluating (200) the service life of a numerically modeled part (20), wherein a calculation is implemented of the service life of the part (20) involving the converted amplitude values of effective stress intensity factor interpolated according to the method of claim 1.

10. A system comprising a processing unit (12) comprising calculation means (14) and a memory (16), the unit being configured for implementing the method for estimating according to claim 1.

11. A computer program product configured for being implemented by a system, and comprising instructions for bringing about the implementation of the method for estimating as claimed in claim 1.

Description

DESCRIPTION OF THE FIGURES

[0051] Other features, objects and advantages of the invention will emerge from the following purely illustrative and non-restrictive description, which must be read with reference to the appended drawings in which:

[0052] FIG. 1 illustrates a crack front for different propagation increments,

[0053] FIG. 2 illustrates a system for implementing the invention,

[0054] FIG. 3 illustrates a part to be analyzed,

[0055] FIG. 4 illustrates two methods and associated embodiments according to the invention,

[0056] FIG. 5 illustrates points converted equivalently with a plane crack with a straight front, from a three-dimensional crack, according to an energy equivalence theory by thermodynamic analysis,

[0057] FIG. 6 illustrates these same points with two interpolation methods: piecewise linear and curvature energy minimization,

[0058] FIG. 7 illustrates the evolution of the polynomial interpolation with the degree of the polynomial,

[0059] FIG. 8 schematically illustrates a plane crack with a straight front,

[0060] FIGS. 9a to 9c illustrate 1D shape functions from a generalization of finite elements and used for calculating curvature energy,

[0061] FIGS. 10 to 10c illustrate a change in representation for smoothing the metrics.

DETAILED DESCRIPTION

[0062] Referring to FIG. 2, a system 10 for interpolating the value of stress intensity factors (SIF) K of a part to be analyzed 20, which is modeled numerically. As indicated in the introduction, the SIF K breaks down into three data KI, KII, KIII. The description will be given only for KI.

[0063] The part to be analyzed 20 is a part intended for aeronautics, for which it must be possible to estimate the service life. The part 10 is typically a compressor or turbine fan blade, an engine disk, a flange, an engine mounting, a housing as illustrated in FIG. 3. This list is illustrative, since for implementing the method, the type of part is not, however, involved.

[0064] The system 10 comprises a data processing unit 12, e.g. a computer or a server, having calculation means 14 configured for implementing a method that will be described in more detail below and, advantageously, having a memory 16. The calculation means 14 may, for example, be a processor, microprocessor, microcontroller, etc. type of computer. The memory 16 may be, for example, a hard disk, a flash memory or a delocalized, cloud type storage space.

[0065] The data processing unit 12 may also be suitable for implementing numerical simulations such as finite elements, of the part to be analyzed 20. Finite element simulation makes it possible to obtain, for each simulation increment, data relating to said part and to crack propagation.

[0066] In particular, it is possible to obtain for each propagation increment, the values of the SIF at associated points and the position of these points, generally along the crack front. The positions comprise, for example, the coordinates of the nodes of the mesh defining the crack front. It is possible to obtain data relating to the cracked surface: in addition to the node coordinates, the faces of the three-dimensional elements defining one among the faces of the crack. These data are stored in the form of a connectivity table, conventionally known as finite elements. Other data, such as the mesh of the cracked surface and of the front may also be obtained.

[0067] Referring to FIG. 4, a method 100 for interpolating the values of a SIF of a three-dimensional crack simulated by finite elements will be described.

[0068] This method 100 is also advantageously used for improving the interpolation of the SIF. This will be described later.

[0069] In a first step E1, a numerical simulation is performed on the part to be analyzed 10. This simulation makes it possible to obtain the aforementioned data.

[0070] In a preferred embodiment, this simulation is carried out by finite elements or extended finite elements.

[0071] This step E1 is performed either by the processing unit 12 of the system, or by another system.

[0072] In both cases, a step E2 of receiving said data (generated during step E1) by the processing unit 12 of the system is implemented. More precisely, for each simulated crack front a first data subset is defined containing the values of the SIF at the associated points on the crack front and their associated respective position, and a second data subset containing the data relating to the cracked surface (the node coordinates defining the crack front and the surface mesh of one of the faces of the crack, e.g. the connectivity table of the faces of the relevant three-dimensional elements) as they have been defined previously. Having a plurality of simulated fronts at different increments, a plurality of first and second data subsets is thereby obtained.

[0073] This is the first step of the interpolation method itself, insofar as the simulation step E1 is not necessarily implemented expressly for later interpolation.

[0074] In step E3, the problem of the geometry of the crack in three dimensions is reduced to a problem of a plane crack with a straight front. The plane crack with a straight front is a standard model known to the person skilled in the art. For this, a thermodynamic equivalence based on an equality of the dissipated energy between the two models is performed. The physical basis is explained below.

[0075] The first principle of thermodynamics, for a transformation of the infinitesimal crack front, expressed in energy density per unit length of crack front, gives:

[00004]$\underset{\begin{array}{c}\mathrm{internal}.Math..Math.\mathrm{energy}\\ \mathrm{density}\end{array}}{\underset{}{{\mathrm{dW}}_{\mathrm{el}}^i}}=\underset{\mathrm{heat}.Math..Math.\mathrm{density}}{\underset{}{.Math..Math.Q}}+\underset{\begin{array}{c}\mathrm{density}.Math..Math.\mathrm{of}.Math..Math.\mathrm{work}.Math..Math.\mathrm{of}\\ \mathrm{external}.Math..Math.\mathrm{forces}\end{array}}{\underset{}{.Math..Math.W_{\mathrm{ext}}}}$

[0076] The second principle of thermodynamics makes it possible to express the irreversibility of the transformation:

[00005]$T.\mathrm{dS}=.Math..Math.Q+\underset{\mathrm{dissipated}.Math..Math.\mathrm{energy}.Math..Math.\mathrm{density}}{\underset{}{G.\mathrm{dSurf}}}$

[0077] The dissipated energy density is always positive. The term dSurf represents the increment of cracked surface per unit length of crack front in the course of the infinitesimal transformation. This quantity is homogeneous with a length: it is a kind of increment of length of crack in the thermodynamic sense of the term. It thus appears a physical magnitude that may be considered for interpolating which has a physical meaning.

[0078] The preceding elements allow the initial three-dimensional cracking problem to be made equivalent to a problem of propagation of a plane crack with a straight front in plane elasticity. More generally, the case of adapted elasticity is considered, i.e. the material has already been able to undergo plastic deformation initially, but crack propagation takes place under cyclic elastic fatigue loading.

[0079] The temperature field may evolve over time, but the temperature field does not vary spatially. Consequently, a fatigue cycle is associated with a temperature.

[0080] One modeling hypothesis consists in considering that the law of crack propagation is a Paris law with Elber correction (references [5], [6]), given by the following equation:

[00006]$\frac{\mathrm{dSurf}}{\mathrm{dN}}=C\left(T\right).{\left(.Math..Math.K.U\left(R,T\right)\right)}^{n\left(T\right)}=C\left(T\right).{\left(.Math..Math.K_{\mathrm{eff}}\right)}^{n\left(T\right)}$$\mathrm{with}.Math..Math.U\left(R,T\right)=a\left(T\right).R+b\left(T\right)$

where C(T) and n(T) are the coefficients of the Paris law and a(T) and b(T) are the parameters of Elber's crack closure law, K is the amplitude of the SIF, which makes it possible to dispense with the load ratio (i.e. the ratio R=Kmin/Kmax) and K.sub.eff is the effective amplitude of the SIF, which takes into account the effect of the crack closure. The link between K and K.sub.eff has been established by Elber.

[0081] The cracked surface increment dSurf per cycle is now connected to the amplitude of the effective stress intensity factor K.sub.eff.

[0082] There is a relationship between K (or K.sub.eff, by equivalence) and the energy release rate G (or the amplitude of said effective rate G.sub.eff by equivalence), given by the following relationships:

[00007]$G=\frac{1}{E^{*}}.Math.\left(K_I^2+K_{\mathrm{II}}^2\right)+\frac{1}{}.Math.K_{\mathrm{III}}^2$$E^{*}=E.Math..Math.\mathrm{in}.Math..Math.\mathrm{plane}.Math..Math.\mathrm{stresses}$$E^{*}=E.Math.\text{/}.Math.\left(1-{}^2\right).Math..Math.\mathrm{in}.Math..Math.\mathrm{plane}.Math..Math.\mathrm{strains}$

where E is Young's modulus, E* is the equivalent Young's modulus, is the shear modulus, is Poisson's ratio.

[0083] The SIF is known at each point of the crack front and for each propagation increment using the numerical simulation of step E1 and retrieved in step E2.

[0084] By replacing the values with their amplitude in the preceding relationships, a relationship is obtained, by definition, between the energy release rate and the amplitude of the effective stress intensity factor (s being the curvilinear abscissa along the crack front):

(K.sub.eff(s)).sup.2=G.sup.eff(s).Math.E*

[0085] The energy dissipated in the course of the cracking process is obtained by the following relationship (by remaining in the infinitesimal transformation hypothesis, therefore considering only one loading cycle):

Dissipated energy=.sub.crack frontG.sub.eff(s).Math.dSurf(s).Math.dl

[0086] The dSurf is known using the data relating to the cracked surface calculated by the numerical simulation in step E1 and retrieved in step E2.

[0087] By postulating an equivalence of the 3D cracking problem with a problem of plane cracking with a straight front, the dissipated energies and the crack front lengths may be equalized. Consequently, it is legitimate to write:

Dissipated energy.sup.glob=G.sub.eff.sup.glob.Math.dSurf.sup.glob.Math.L.sub.front

[0088] The notation glob means that the datum is specific to the model in plane cracking with a straight front.

[0089] Plane cracking with a straight front also respects the Paris propagation law. Consequently, it is obtained that

[00008]$\frac{{\mathrm{dSurf}}^{\mathrm{glob}}}{\mathrm{dN}}=C.{\left(\sqrt{.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}.E^{*}}\right)}^n$

[0090] By equalizing the two dissipated energies and by rewriting the equation, the following relationship is obtained

[00009]$.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}={\left(\frac{\mathrm{Dissipated}.Math..Math.\mathrm{energy}}{C.{\left(E^{*}\right)}^n.Math.L_{\mathrm{front}}}\right)}^{\frac{1}{n+1}}={\left(\frac{{}_{\mathrm{crack}.Math..Math.\mathrm{front}}.Math..Math..Math.G_{\mathrm{eff}}\left(s\right).\mathrm{dSurf}\left(s\right).\mathrm{dl}}{C.{\left(E^{*}\right)}^n.Math.L_{\mathrm{front}}}\right)}^{\frac{1}{n+1}}$

[0091] However, dSurf.sup.glob/dN=dSurf.sup.glob (because dN=1 here, since only a single cycle is considered) is homogeneous with a crack length. Consequently, its integration in the course of the time between the first cycle and the last cycle makes it possible to obtain an equivalent crack length Surf.sup.glob, noted length a, in the thermodynamic sense. A length is actually obtained which has a physical meaning.

[0092] The previous relationship established between K.sub.eff.sup.glob and G.sub.eff.sup.glob makes it possible to inject the amplitude of the effective SIF into the preceding equation.

[0093] Thus a relationship is obtained connecting the amplitude of the equivalent global effective SIF in the modeling of the flat crack with a straight front as a function of the length of the crack Surf glob i.e. the length a. These data may then be interpolated.

[0094] This physical construction is valid in all generality. However, it requires knowing all the information of the problem for each loading cycle which is generally too costly in terms of calculation time and memory space for saving all the results (at least a thousand cycles should be simulated). The finite element calculations are often made with no fixed maximum propagation increment along the crack front, or pre-integrate over a given number of cycles (10, 100 or 1000, for example). This amounts to de-refining the discretization in propagation increment in order to reduce the calculation time.

[0095] Consequently, the propagation law is replaced by:

[00010]$\frac{{\mathrm{dSurf}}_{\mathrm{eq}}^{\mathrm{glob}}}{\mathrm{dN}}=C.{\mathrm{coeff}\left(\sqrt{.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}.E^{*}}\right)}^n=C_{\mathrm{eq}}.{\left(\sqrt{.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}.E^{*}}\right)}^{n_{\mathrm{eq}}}$

where coeff is the coefficient which allows the pre-integration of multiple cycles.

[0096] This amounts to modifying the coefficients of the Paris law C and n by the coefficients C.sub.eq=C.coeff and n.sub.eq=n. The link between dSurf.sub.eq.sup.glob and dSurf.sup.glob is therefore the following:

dSurf.sub.eq.sup.glob=coeff.Math.dSurf.sup.glob

[0097] In addition, the term dSurf(s) is replaced by the incremented surface between two simulated propagation increments dSurf.sub.eq(s) (therefore not infinitesimal).

[0098] In conclusion, the above reasoning applies approximately by performing the preceding substitutions, and the equivalent crack length is still deduced in the same way.

[0099] Thus, in step E3, from the plurality of first and second data subsets retrieved in step E2, i.e. the SIF data along the different fronts of simulated cracks at different increments and the data relating to the cracked surfaces, the means of which have been previously explained, the global effective amplitude of the SIF K.sub.eff.sup.glob and the associated length a of the plane crack with a straight front are determined. This determination is implemented by the processing unit 12.

[0100] For a given crack front, the simulated data of the SIF make it possible to find K (first data subset) and using the preceding equations make it possible to calculate the K.sub.eff.sup.glob and the cracked surface data dSurf.sub.eq (second data subset) make it possible to calculate the length a. Due to the plurality of data corresponding to each front, a plurality of pairs (a, K.sub.eff.sup.glob) is obtained.

[0101] The equivalence obtained may also be applied to a plane crack with a straight front. The following relationship shows that the system remains invariant:

[00011]$\begin{array}{c}.Math..Math.G_{\mathrm{eff}}^{\mathrm{glob}}=.Math.{\left(\frac{{}_{\mathrm{crack}.Math..Math.\mathrm{front}}.Math..Math..Math.G_{\mathrm{eff}}\left(s\right).{\mathrm{dSurf}}_{\mathrm{eq}}\left(s\right).\mathrm{dl}}{C_{\mathrm{eq}}.{\left(E^{*}\right)}^{n_{\mathrm{eq}}}.L_{\mathrm{front}}}\right)}^{\frac{1}{n_{\mathrm{eq}}+1}}\\ =.Math.{\left(\frac{.Math..Math.G_{\mathrm{eff}}.C_{\mathrm{eq}}.{\left(.Math..Math.G_{\mathrm{eff}}.E^{*}\right)}^{n_{\mathrm{eq}}}}{C_{\mathrm{eq}}.{\left(E^{*}\right)}^{n_{\mathrm{eq}}}}\right)}^{\frac{1}{n_{\mathrm{eq}}+1}}\\ =.Math..Math..Math.G_{\mathrm{eff}}\end{array}$

[0102] From this it follows that the cracked surfaces are the same, so that the final system is identical to the initial problem. This proves the consistency of the method based on a thermodynamic analysis of the problem.

[0103] FIG. 5 illustrates the representation of the pairs (a, K.sub.eff.sup.glob) in the plane. The references F1, F2, F3, F4 are given with reference to FIG. 1, for indicating to which front of the original crack the points of the graph correspond. This is a writing misuse.

[0104] It is recalled just now that the goal is to be able to determine SIFs for crack states that have not been numerically simulated, i.e. for states that occur strictly between two successive increments.

[0105] To this end, in a step E4, an interpolation of the K.sub.eff.sup.glob as a function of the length a is performed. The interpolation is performed by the calculation means 14 of the processing unit 12.

[0106] The interpolation may be performed according to several methods. Interpolations are preferred which do not create information for the problem or whereby a minimum of information is added. For this, a piecewise linear interpolation IL or interpolation by curvature energy minimization IMC would be suitable. Both interpolations may also be performed, insofar as they may have different applications.

[0107] FIG. 6 illustrates these two interpolations.

[0108] Interpolation by global curvature energy minimization may be seen as a method based on physical considerations specific to the cracking problem. Indeed, the scientific literature on the subject shows that if a crack is propagated without branching, then the evolution of the stress intensity factor with crack extension may be represented by an infinitely differentiable function. This means that calculating the curvature of such a function is meaningful.

[0109] In addition, minimizing a curvature energy takes place in some physical situations such as, for example, in the physics of soap bubbles or more generally in the physics of membranes.

[0110] FIG. 7 illustrates a linear interpolation between three points. This example shows us that it is possible to improve on approximating an interpolation by polynomials. The higher the degree of the polynomial, the better the approximation will be. It is also possible to demonstrate that the maximum difference between linear interpolation and polynomial approximation tends toward zero when the degree of the polynomial tends toward infinity. Thus, linear interpolation may be seen as the passage to the limit at infinity that would be applied to an infinitely differentiable function.

[0111] This approach confers a physical character on both piecewise linear interpolation and interpolation based on curvature energy minimization.

[0112] Performing an interpolation consists in completing the missing information using a principle which should at least respect the physics of the problem considered, i.e. not introducing any elements that would violate certain fundamental equations that the problem obeys.

[0113] For this reason, interpolation by curvature energy minimization is a physical method because it does not contradict the theory of fracture mechanics. The method consisting of piecewise linear interpolations is not initially physical but may be seen as the passage to the limit of a method of global polynomial interpolation, which would be physical and by extension would make the piecewise linear interpolation method almost physical too.

[0114] Another point of view regarding the method of linear interpolation would be to see it as a method of minimizing the distance between points. This way of seeing things then places it in the category of optimal methods. Linear interpolation minimizes the distance while the other minimizes the curvature energy under the position constraints of the points to be connected. Minimizing distances may be seen as a means of connecting points without introducing additional information to the initial problem, and this results in a continuous field. In the case of the second method, at least one additional piece of information is introduced, which is that the field is regular. In this case, curvature energy minimization is an optimal method in the sense that the physical information to be introduced is minimized whereas the linear interpolation method is optimal in the sense that the purely mathematical information is minimized.

[0115] Finally, in a step E5, after the interpolation step, the interpolated data are stored, in order that other applications may have access thereto. Storage typically takes place in the memory 16. It is generally referred to as a form. The form is a table that gives the stress intensity factors as a function of the length of the crack.

[0116] Consequently, the interpolation step E4 generates a file, in the form of a text or a table containing the form, i.e. combining the converted data and the interpolated converted data. Step E5 consists in storing this file.

[0117] Indeed, as indicated in the introduction, obtaining the form is generally used for calculating the service life of a part. Consequently, with reference to FIG. 4, a method 200 is described for determining the service life of the part to be analyzed 20.

[0118] In a step E1, the processing unit 12 receives the form calculated in step E4 and/or stored in step E5 of the preceding method. This retrieval step may simply consist in having access to the memory 16.

[0119] In a step E2, a method is implemented for calculating the fatigue crack propagation service life of the part. Document [7] describes such a method.

[0120] It is sufficient to use an existing code and replace the SIF form with that provided in this invention for determining the crack propagation service life.

[0121] This method 200 may be implemented by a system other than that described previously.

Error Refinement

[0122] When all of the information is available concerning the problem to be resolved, the result of interpolation will always be the same regardless of the interpolation method used and which meets all physical constraints of the problem. Typically, over the intervals I1 in FIG. 6, the interpolation intuitively appears to be of good quality; over the I2 intervals, it intuitively appears to be of lower quality.

[0123] It is thus possible to see the differences between two physical and or mathematical optimal interpolations as a means of identifying a lack of information that would be useful for improving the interpolation or the degree of confidence in the interpolation.

[0124] This consideration helps set up a step of checking relevance, illustrated in FIG. 6.

[0125] In two interchangeable substeps E41 and E42 of interpolation, a linear piecewise interpolation and curvature energy interpolation are performed.

[0126] It is noted that strictly between two length values corresponding to two simulated successive increments, i.e. strictly in an interpolated area, there are differences in value K.sub.eff.sup.glob between the two interpolations.

[0127] In a substep E43, at least one of these differences is calculated, which is compared, in a substep E44, to a predetermined threshold value VSP.

[0128] The predetermined threshold VSP is selected according to the desired quality of the interpolation.

[0129] It is also possible to compare a value of these differences for each area between two successive increments, or to compare an average or a maximum value, etc. Thus a magnitude r representative of the difference will be referred to more generally. This magnitude, r indicates that there is a difference that could be quantified by the functions mentioned previously.

[0130] Finally, a substep of comparison E45 between the representative magnitude r and the threshold VSP is performed.

[0131] If the representative magnitude r is below the threshold VSP, it may be considered that the two interpolations are of good quality and that the data, if they had been simulated, would be close to the two interpolation values.

[0132] If the representative magnitude r is above the threshold VSP, then it may be considered that there is too much uncertainty regarding the interpolation, and obtaining new simulated values is useful, or even necessary. To this end, a calculation instruction is generated from the numerical simulation, typically by finite elements (or extended finite elements) of the part to be analyzed 20 of at least a value of the stress intensity factor SIF and the position of the crack front between the two successive increments. In other words, the processing unit 12 sends an instruction to perform a part of step E1, with a reduced increment for at least one simulation.

[0133] This is referred to as refinement of the discretization in crack length.

[0134] By way of example, a predetermined threshold value VSP of 2% may be suitable. The criterion depends on the specifications.

[0135] It may be decided, for example, that r=5 and that as soon as r/max

(K.sub.eff.sup.glob among linear and curvature interpolation)>2% for a given value of length a, an instruction for recalculating the simulation is generated.

[0136] Finally, the step of storing E5 may comprise the storing of both the interpolations and the representative magnitudes of the differences r and/or the differences .

Appendix 1: Example of Interpolation by Curvature Energy Minimization

[0137] Interpolation by curvature energy minimization only gives consistent physical predictions within the context of a plane crack with a straight front in plane elasticity, as illustrated in FIG. 8, where represents the stress, and a the length of the crack. However, the method described previously makes it possible to thermodynamically consider any three-dimensional crack as a crack satisfying these properties.

[0138] To perform this interpolation, the curvature energy has to be determined.

[0139] The function f is defined in the following way (x is the virtual crack length, previously referred to as a):

[00012]$\frac{{\mathrm{dK}}_I}{\mathrm{dx}}=f\left(x\right)$

[0140] The need to introduce virtual crack lengths is explained in APPENDIX 2.

[0141] The function f is written as a linear combination of the piecewise reference polynomial functions (of degree 5 for ensuring a regularity C2 over the whole evolution range of the virtual crack length, which involves a generalization of the conventional finite elements) given below:

[00013]$\{\begin{array}{c}{C}_{1,1}\left(z\right)=0,5-0.9375.z+0,625..Math.z^3-0,1875..Math.z^5\\ C_{1,2}\left(z\right)=2..Math.\left(\begin{array}{c}0,3125-0.4375.z-0.375.z^2+\\ 0,625..Math.z^3+0,0625..Math.z^4-0,1875..Math.z^5\end{array}\right)\\ C_{1,3}\left(z\right)=4.Math.\left(\begin{array}{c}0,0625-0.0625.z-0.125.z^2+\\ 0,125.Math.z^3+0,0625..Math.z^4-0,0625..Math.z^5\end{array}\right)\\ C_{2,1}\left(z\right)=0,5+0.9375.z-0,625..Math.z^3+0,1875..Math.z^5\\ C_{2,2}\left(z\right)=2..Math.\left(\begin{array}{c}-0,3125-0.4375.z+0.375.z^2+\\ 0,625..Math.z^3-0,025..Math.z^4=0,1875..Math.z^5\end{array}\right)\\ C_{2,3}\left(z\right)=4..Math.\left(\begin{array}{c}0,0625+0.0625.z-0.125.z^2-\\ 0,125..Math.z^3+0,0625..Math.z^4+0,0625..Math.z^5\end{array}\right)\\ z=\frac{x-x^{j+1}}{x^{j+1}-x^j}+\frac{x-x^j}{x^{j+1}-x^j}\end{array}$

[0142] The first index indicates the number of the node on which the nodal value is taken.

[0143] The second index indicates the nature of the nodal value:

1: Nodal value of the stress intensity factor (SIF)
2: Nodal value of the derivative with respect to the variable z of the SIF
3: Nodal value of the second derivative with respect to the variable z of the SIF

[0144] The one-dimensional element associated with these shape functions is defined so that the first node is located at z=1 and the second at z=1.

[0145] The curvilinear abscissa associated with the curve giving the evolution of the stress intensity factor as functions of the virtual crack length is given by:

[00014]$\begin{array}{c}s\left(x\right)=.Math.{}_{x_1}^x.Math.\sqrt{\left(1+{\left(f\left(u\right)\right)}^2\right)}.Math.\mathrm{du}\\ =.Math.\underset{i=1}{\overset{j-1}{.Math.}}.Math.{}_{x_1}^{x_{i+1}}.Math.\sqrt{\left(1+{\left(f\left(u\right)\right)}^2\right)}.Math.\mathrm{du}+\\ .Math.{}_{x_j}^x.Math.\sqrt{\left(1+{\left(f\left(u\right)\right)}^2\right)}.Math.\mathrm{du}\end{array}$$\mathrm{It}.Math..Math.\mathrm{will}.Math..Math.\mathrm{be}.Math..Math.\mathrm{noted}.Math.\text{:}$$\{\begin{array}{c}z=\frac{x-x^{j+1}}{x^{j+1}-x^j}+\frac{x-x^j}{x^{j+1}-x^j}\\ f\left(x\right)=f^j\left(z\right)=f_{11}^j.Math.C_{1,1}\left(z\right)+f_{12}^j.Math.C_{1,2}\left(z\right)+\\ f_{13}^j.Math.C_{1,3}\left(z\right)+f_{11}^{j+1}.Math.C_{2,1}\left(z\right)+f_{12}^{j+1}+C_{2,2}\left(z\right)+\\ f_{13}^{j+1}.Math.C_{2,3}\left(z\right)\end{array}$

[0146] It will be noted that the degrees of freedom f.sub.11.sup.j and f.sub.11.sup.j+1=f.sub.21.sup.j are known. Optimization will focus on the other degrees of freedom.

[0147] Whence the following is obtained (where L.sub.x is the length of the virtual element):

[00015]$s_x\left(x\right)=\left(\underset{i=1}{\overset{j-1}{.Math.}}.Math.{}_{-1}^1.Math.\sqrt{\left(1+{\left(f^i\left(z\right)\right)}^2\right)}.Math.\mathrm{dz}+{}_{-1}^z.Math.\sqrt{\left(1+{\left(f^j\left(z\right)\right)}^2\right)}.Math.\mathrm{du}\right).Math.\frac{L_x}{2}$

[0148] The unit tangent vector is obtained as follows:

[00016]$t\left(s_x\right)=\frac{{\mathrm{dK}}_I}{{\mathrm{ds}}_x}.Math.\left(s_x\right)$

[0149] FIGS. 9a to 9c illustrate the values of the shape functions C11, C12, C13, C21, C22 and C23 and their derivative and second derivative. It is noted that these curves satisfy their specific constraints (0 or 1) at the nodes.

[0150] The unit tangent vector to the curve K.sub.I may now be calculated. The formulae given previously lead to:

[00017]$\begin{array}{c}{t}_x=.Math.\frac{d.Math..Math.{}_x}{{\mathrm{ds}}_x}=\frac{d.Math..Math.{}_x}{\mathrm{dx}}.Math.\frac{\mathrm{dx}}{{\mathrm{ds}}_x}\\ =.Math.\left(\begin{array}{c}1\\ f\left(x\right)\end{array}\right).Math.\frac{1}{\frac{{\mathrm{ds}}_x}{\mathrm{dz}}.Math.\frac{\mathrm{dz}}{\mathrm{dx}}}\\ =.Math.\left(\begin{array}{c}1\\ f\left(x\right)\end{array}\right).Math.{\left[1+{\left(f^j\left(\frac{x-x^{j+1}}{x^{j+1}-x^j}+\frac{x-x^j}{x^{j+1}-x^j}\right)\right)}^2\right]}^{-\frac{1}{2}}\end{array}$

[0151] All that remains is the calculation of the curvature. For this, the derivative of t.sub.x is calculated with respect to the curvilinear abscissa:

[00018]$\begin{array}{c}\frac{{\mathrm{dt}}_x}{{\mathrm{ds}}_x}=.Math.\frac{d}{{\mathrm{ds}}_x}.Math.\left(\frac{d.Math..Math.{}_x}{\mathrm{dx}}.Math.\frac{\mathrm{dx}}{{\mathrm{ds}}_x}\right)\\ =.Math.\frac{{d.Math.}^2.Math.{}_x}{{\mathrm{dx}}^2}.Math.{\left(\frac{\mathrm{dx}}{{\mathrm{ds}}_x}\right)}^2+\frac{d.Math..Math.{}_x}{\mathrm{dx}}.Math.\frac{d^2.Math.x}{{d\left(s_x\right)}^2}\\ =.Math.\left(\begin{array}{c}0\\ f^{}\left(x\right)\end{array}\right).Math.{\left(\frac{\mathrm{dx}}{{\mathrm{ds}}_x}\right)}^2+\left(\begin{array}{c}1\\ f\left(x\right)\end{array}\right).Math.\frac{d^2.Math.x}{{d\left(s_x\right)}^2}\end{array}$

[0152] For determining the derivative terms with respect to the curvilinear abscissa, the formulae for differential passage from the curvilinear coordinate to the virtual coordinate are expressed:

[00019]$\left(\begin{array}{c}\frac{d}{{\mathrm{ds}}_x}\\ \frac{d^2}{{d\left(s_x\right)}^2}\end{array}\right)=\underset{{\mathrm{Diff}}_{s_x}}{\underset{}{\left[\begin{array}{cc}\frac{\mathrm{dx}}{{\mathrm{ds}}_x}& 0\\ \frac{d^2.Math.x}{{d\left(s_x\right)}^2}& {\left(\frac{\mathrm{dx}}{{\mathrm{ds}}_x}\right)}^2\end{array}\right]}}.Math.\left(\begin{array}{c}\frac{d}{\mathrm{dx}}\\ \frac{d^2}{{\mathrm{dx}}^2}\end{array}\right)$

[0153] This matrix equation is invertible, which makes it possible to obtain:

[00020]$\left(\begin{array}{c}\frac{d}{\mathrm{dx}}\\ \frac{d^2}{{\mathrm{dx}}^2}\end{array}\right)=\underset{{\mathrm{Diff}}_{s_x}}{\underset{}{\left[\begin{array}{cc}\frac{{\mathrm{ds}}_x}{\mathrm{dx}}& 0\\ \frac{d^2.Math.s_x}{{\mathrm{dx}}^2}& {\left(\frac{{\mathrm{ds}}_x}{\mathrm{dx}}\right)}^2\end{array}\right]}}.Math.\left(\begin{array}{c}\frac{d}{{\mathrm{ds}}_x}\\ \frac{d^2}{{d\left(s_x\right)}^2}\end{array}\right)$

[0154] It is noted that Diff.sub.x.sup.1=Diff.sub.s.sub.x, which makes it possible to easily determine

[00021]${\left(\frac{{\mathrm{ds}}_x}{\mathrm{dx}}\right)}^2.Math..Math.\mathrm{and}.Math..Math.\frac{{\mathrm{ds}}_x}{\mathrm{dx}}.$

Finally, it is found that:

[00022]$\frac{{\mathrm{dt}}_x}{{\mathrm{ds}}_x}=\left(\begin{array}{c}0\\ f^{}\left(x\right)\end{array}\right).Math.{\mathrm{Diff}}_x^{-1}\left(2,2\right)+\left(\begin{array}{c}1\\ f\left(x\right)\end{array}\right).Math.{\mathrm{Diff}}_x^{-1}\left(2,1\right)$$\mathrm{With}.Math.\text{:}$${\mathrm{Diff}}_x=\left[\begin{array}{cc}\frac{{\mathrm{ds}}_x}{\mathrm{dx}}& 0\\ \frac{d^2.Math.s_x}{{\mathrm{dx}}^2}& {\left(\frac{{\mathrm{ds}}_x}{\mathrm{dx}}\right)}^2\end{array}\right]$

[0155] The curvature is automatically deduced therefrom:

[00023]$\left(x\right)=\sqrt{{\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}^2+{\left(f^{}\left(x\right).{\mathrm{Diff}}_x^{-1}\left(2,2\right)+f\left(x\right).{\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}^2}$

[0156] The matrix Diff.sub.x.sup.1 must now be explained. For this, the following conventional 22 matrix inversion formula is used:

[00024]$A=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right).Math.A^{-1}=\frac{1}{\det \left(A\right)}.Math.\left(\begin{array}{cc}d& -b\\ -c& a\end{array}\right).Math..Math.\mathrm{with}$$\det \left(A\right)=a.d-b.c$$\mathrm{Hence}.Math.\text{:}$$\{\begin{array}{c}z=\frac{x-x^{j+1}}{x^{j+1}-x^j}+\frac{x-x^j}{x^{j+1}-x^j}\\ \det \left({\mathrm{Diff}}_x\right)={\left(\frac{{\mathrm{ds}}_x}{\mathrm{dx}}\right)}^3={\left[1+{\left(f^j\left(z\right)\right)}^2\right]}^3\\ {\mathrm{Diff}}_x^{-1}\left(x\right)=\frac{1}{\det \left({\mathrm{Diff}}_x\right)}.Math.\left[\begin{array}{cc}{\left[1+{\left(f^j\left(z\right)\right)}^2\right]}^2& 0\\ -\frac{4}{L_x}.Math.f^j\left(z\right).Math.\frac{{\mathrm{df}}^j}{\mathrm{dz}}.Math.\left(z\right)& 1+{\left(f^j\left(z\right)\right)}^2\end{array}\right]\end{array}$

[0157] The functional that it is sought to minimize is the following (if the mesh consists of n1 element):

W.sub.x=.sub.x.sub.1.sup.x.sup.n{square root over ((1+(f(u)).sup.2))}.Math.((u)).sup.2.Math.du

[0158] Preferably a Newton algorithm is then used for determining the parameters f.sub.ij.sup.k that minimize this functional. It is thus necessary to know the gradient vectors and Hessian matrices (with respect to the unknown variables) of the magnitudes involved in the optimization problem.

[00025]$\{\begin{array}{c}{\mathrm{Grad}\left({\left(u\right)}^2\right)}_i=\left[\begin{array}{c}2..Math.\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_i+\\ 2..Math.\left(f^{}\left(x\right).{\mathrm{Diff}}_x^{-1}\left(2,2\right)+f\left(x\right).{\mathrm{Diff}}_x^{-1}\left(2,1\right)\right).Math.\\ \left(\begin{array}{c}{\mathrm{Grad}\left(f^{}\left(x\right)\right)}_i.{\mathrm{Diff}}_x^{-1}\left(2,2\right)+f^{}\left(x\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_i+\\ {\mathrm{Grad}\left(f\left(x\right)\right)}_i.{\mathrm{Diff}}_x^{-1}\left(2,1\right)+f\left(x\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_i\end{array}\right)\end{array}\right]\\ {H\left({\left(u\right)}^2\right)}_i=\left[\begin{array}{c}2..Math.\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right).Math.{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_i+\\ 2..Math.\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right).{H\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_{\mathrm{ij}}+\\ 2..Math.\left(\begin{array}{c}{\mathrm{Grad}\left(f^{}\left(x\right)\right)}_j.{\mathrm{Diff}}_x^{-1}\left(2,2\right)+f^{}\left(x\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_j\\ {\mathrm{Grad}\left(f\left(x\right)\right)}_j.{\mathrm{Diff}}_x^{-1}\left(2,1\right)+f\left(x\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_j\end{array}\right).Math.\\ \left(\begin{array}{c}{\mathrm{Grad}\left(f^{}\left(x\right)\right)}_i.{\mathrm{Diff}}_x^{-1}\left(2,2\right)+f^{}\left(x\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_i+\\ {\mathrm{Grad}\left(f\left(x\right)\right)}_i.{\mathrm{Diff}}_x^{-1}\left(2,1\right)+f\left(x\right).{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_i\end{array}\right)+\\ 2..Math.\left(f^{}\left(x\right).{\mathrm{Diff}}_x^{-1}\left(2,2\right)+f\left(x\right).{\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)\\ \left(\begin{array}{c}{H\left(f^{}\left(x\right)\right)}_{\mathrm{ij}}.Math.{\mathrm{Diff}}_x^{-1}\left(2,2\right)+{\mathrm{Grad}\left(f^{}\left(x\right)\right)}_i.{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_j+\\ {\mathrm{Grad}\left(f^{}\left(x\right)\right)}_j.{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_i+f^{}\left(x\right).{H\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_{\mathrm{ij}}\\ {H\left(f\left(x\right)\right)}_{\mathrm{ij}}.Math.{\mathrm{Diff}}_x^{-1}\left(2,1\right)+{\mathrm{Grad}\left(f\left(x\right)\right)}_i.{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_j+\\ {\mathrm{Grad}\left(f\left(x\right)\right)}_j.{\mathrm{Grad}\left({\mathrm{Diff}}_x^{-1}\left(2,1\right)\right)}_i+f\left(x\right).{H\left({\mathrm{Diff}}_x^{-1}\left(2,2\right)\right)}_{\mathrm{ij}}\end{array}\right)\end{array}\right]\end{array}.Math..Math.\{\begin{array}{c}{\mathrm{Grad}\left(W_x\right)}_i=\left(\begin{array}{c}{}_{x_1}^{x_n}.Math.{\left(1+{\left(f\left(u\right)\right)}^2\right)}^{\frac{1}{2}}.{\mathrm{Grad}\left(f\left(u\right)\right)}_i.{\left(\left(u\right)\right)}^2.\mathrm{du}+\\ 2..Math.{}_{x_1}^{x_n}.Math.\sqrt{\left(1+{\left(f\left(u\right)\right)}^2\right)}.{\mathrm{Grad}\left({\left(u\right)}^2\right)}_i.Math.\mathrm{du}\end{array}\right)\\ {H\left(W_x\right)}_{\mathrm{ij}}=\left(\begin{array}{c}-{}_{x_1}^{x_n}.Math.{\left(1+{\left(f\left(u\right)\right)}^2\right)}^{-\frac{3}{2}}.{\mathrm{Grad}\left(f\left(u\right)\right)}_j.{\mathrm{Grad}\left(f\left(u\right)\right)}_i.{\left(\left(u\right)\right)}^2.\mathrm{du}+\\ {}_{x_1}^{x_n}.Math.{\left(1+{\left(f\left(u\right)\right)}^2\right)}^{-\frac{1}{2}}.{H\left(f\left(u\right)\right)}_{\mathrm{ij}}.{\left(\left(u\right)\right)}^2.\mathrm{du}+\\ {}_{x_1}^{x_n}.Math.\sqrt{\left(1+{\left(f\left(u\right)\right)}^2\right)}.{H\left({\left(u\right)}^2\right)}_{\mathrm{ij}}.Math.\mathrm{du}\end{array}\right)\end{array}$

[0159] This algorithm may be implemented by the processing unit 12.

Appendix 2: Variable Metric Method

[0160] FIG. 10a illustrates the points of the front with different propagation increments. It is noted that the distance between the same two points of the front is not constant between two increments, which means that the length scale associated with each element is different. This is then referred to as different metrics.

[0161] The fact of having elements with different metrics (therefore discontinuous at the point of intersection) while imposing a regularity connection C2 between the elements leads to numerical instabilities.

[0162] Another point of view consists in saying that the elements do not have the same contributions to the functional (as the result of their different dimensions) while they all convey the same amount of information. Whatever the point of view, this may induce a numerical bias has to be suitably addressed.

[0163] To resolve this problem, it is sufficient to utilize a virtual representation space in which the interpolation elements all have the same dimensions. For example, if the discretized crack lengths are x{0,1; 0,2; 0,4; 0,6; 0,8}, then x.sub.virtual{1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11} will be taken.

[0164] For passing from virtual representation to real representation, a bijection is used between the two reference frames (see FIG. 10b).

[0165] The same interpolation functions are used as before. The correspondence between the two representation spaces is implemented by imposing the degrees of freedom f.sub.11.sup.j and f.sub.11.sup.j+1=f.sub.21.sup.j on each node. Then, in order to have a metric that evolves in a soft way between each element, the method described previously is implemented. This allows a smooth correspondence to be obtained between the two representation spaces as represented in FIG. 10c.

REFERENCES

[0166] [1] Bueckner H G Z (1970), A novel principle for the computation of stress intensity factors, Angew, Math. Mech., Vol. 50, p. 529-546. [0167] [2] Rice J. R. (1972), Some remarks on elastic crack-tip stress fields, Int. J. Solids and Structures, 8, 751-758. [0168] [3] Rice J. R. (1989), Weight function theory for three-dimensional elastic crack analysis, Fracture Mechanics: Perspectives and Directions (Twentieth Symposium), ASTM STP 1020, R. P. Wei and R. P. Gangloff, Eds., American Society for Testing and Materials, Philadelphia, pp. 29-57. [0169] [4] V. Lazarus (1997), Some three-dimensional problems of the mechanics of brittle fracture, Thesis, University of Paris 6. [0170] [5] W. Elber. The significance of fatigue crack closure. ASTM STP, 486: 230-242, 1971. [0171] [6] PC. Paris, F. Erdogan, 1963, A critical analysis of crack propagation laws. J Basic Eng 85, pp 528-534. [0172] [7] NASGROFracture Mechanics and Fatigue Crack Growth Analysis SoftwareReference Manual, Version 6.2, 2011.