CN103955604B - A kind of containing crackle metal gradient material Prediction model for residual strength method - Google Patents

Abstract

The invention discloses a kind of containing crackle metal gradient material Prediction model for residual strength method.The method utilizes finite element model, Cracking behavior in the positional displacement interpolation function characterization model adopting expanding element to be interrupted, take crack length as variable, tracking is split sharp position and is dynamically judged the stress intensity factor of crack tip fracture toughness size corresponding to crack tip position, to determine whether crackle is expanded, realize the Prediction model for residual strength of functionally gradient material (FGM).The method without the need to cutting at geometry when setting up finite element model, is split sharp grid without the need to segmentation, is had higher precision.

Description

A kind of containing crackle metal gradient material Prediction model for residual strength method

Technical field

The invention belongs to the analog simulation of metal gradient material mechanical performance, is more particularly a kind of containing crackle metal gradient material Prediction model for residual strength method.

Background technology

Metal gradient material is with the modern development increasing material manufacture (3D printing) technique, function of reference functionally gradient material (FGM) changes material component to realize the design concept of specific function, the advanced composite material (ACM) system that the different metal component content that processing obtains changes along specific direction.This material not only remains traditional metal materials characteristic and stands good in structure-bearing position, decreases the weak link being mechanically connected between parts and causing simultaneously, and can need to realize specific function according to design as Functionally Graded Materials.

For metal material, because self carefully sees the effect of cyclic loading in defect and use procedure, often there is crackle in inside configuration, the quiet load-bearing capacity of Cracked structure is the residual intensity of this structure, for being against any misfortune property is destroyed, the prediction of residual intensity is very important to structural break safety evaluation.The domestic and international residual intensity for homogeneous metal material carries out extensive work for a long time, tries hard to obtain the Changing Pattern that increases with crackle of residual intensity, and just known possible crack size, predict whether meet residual intensity requirement.

Linear elastic fracture mechanics, through the research of long-term practice and theory, establishes the K>=K of Crack Extension _iCfailure criteria, wherein K is the stress intensity factor of crack tip under linear elasticity state, K _iCfor fracture toughness, be the performance parameters relevant with material.For conventional uniform metal material, K value has corresponding computing method, but for metal gradient material, existing forecasting techniques is no longer applicable, comprising:

1) the violent and r of crack tip change of stress field ^-1/2singularity, calculating the crack tip stress intensity factor by conventional finite elements method needs subdivided meshes.For metal gradient material, due to the change of material property, the stress distribution splitting point is more complicated, require higher to mesh quality, assess the cost excessive, there is larger cumulative errors, in the paper of " research of metal gradient material stress intensity factor " delivered as Shen Yiwen etc., then adopting the method for the subdivided meshes when building finite element model to calculate the crack tip stress intensity factor;

2) use expanding element can effectively reduce split sharp number of grid and reduce assess the cost.But for metal gradient material, existing expanding element can not realize the change of unit internal material performance, be the graded of simulation material performance, grid cell yardstick still has higher requirements;

3) due to material property graded, stress intensity factor calculates and no longer has nothing to do with path of integration, and Traditional calculating methods can cause maximum error, no longer applicable;

4) for functionally gradient material (FGM), fracture toughness K _iCbe no longer definite value, but the variable relevant to splitting sharp position place material component.

Summary of the invention

Instant invention overcomes in prior art and pass through finite element model subdivided meshes thus cause the defect that is excessive and the larger cumulative errors of existence that assesses the cost, Cracking behavior in the positional displacement interpolation function characterization model adopting expanding element to be interrupted, with crack length a for variable, tracking is split sharp position and is dynamically judged K>=K _iC, to determine whether crackle is expanded, realize the Prediction model for residual strength of functionally gradient material (FGM).

For solving the problems of the technologies described above, the present invention is open a kind of containing crackle metal gradient material Prediction model for residual strength method, specifically comprises the following steps:

Step 1, build the finite element model of this metal gradient material, using the node of crack tip place unit as splitting point strengthening node in this finite element model, node is strengthened using the node of crackle section place unit as crackle section, all the other nodes are regular node, wherein the node that crosses of crack tip place unit and crackle section place unit is classified as and splits point and strengthen node, and is numbered all nodes;

Step 2, include all point strengthening nodes that splits in crack tip strengthening set of node N _Λ, all crackle section strengthening node crackle section strengthening set of node N _Γ;

Elastic modulus E, the Poisson ratio μ and fracture toughness K of all node respective material within the scope of step 3, acquisition gradient zones _iClinear changing relation's formula is in the y-direction as follows:

$\left\{\begin{array}{cc}E=E_1+(E_2-E_1)(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\\ \mathrm{\&mu;}={\mathrm{\&mu;}}_1+({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\\ K_{\mathrm{IC}}=K_{\mathrm{IC}1}+(K_{\mathrm{IC}2}-K_{\mathrm{IC}1})(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\end{array}\right.$

In formula, t is Thickness of Gradient Layer, E ₁, μ ₁, K _iC1be respectively the elastic modulus of the first component materials, Poisson ratio and fracture toughness, E ₂, μ ₂, K _iC2be respectively the elastic modulus of the second component materials, Poisson ratio and fracture toughness; With model geometric center for initial point, set up xy coordinate system, wherein y direction is parallel to material gradient direction;

Step 4, according to following formula Confirming model structural stiffness matrix K:

K=∑(∫ _vB(x,y) ^TCB(x,y)dv)

In formula, C is the constitutive matrix of material, and B (x, y) is cell geometry matrix, and subscript T is transpose of a matrix computing, and v represents that integral domain is individual unit scope;

Step 5, initialization external applied load F=0;

Step 6, solve following finite element discretization equation according to current external applied load F:

Ku _I=F

In formula, F is external applied load vector, u _ifor nodal displacement to be solved vector, its component comprises wherein u _ifor regular node displacement, a _ifor N _Γthe additional displacement that the crackle section that collection interior nodes causes due to crackle is discontinuous brought, for N _Λthe additional displacement that the crack tip that collection interior nodes causes due to crackle is discontinuous brought, wherein, α=1,2,3,4, subscript I is node serial number;

Step 7, based on nodal displacement required by step 6, strengthening displacement approximating function interpolation is utilized to obtain displacement field

N represents all node set, N _i(x, y), H _i(x, y) and be respectively standard element positional displacement interpolation function, crackle section interpolating function and crack tip interpolating function, expression is as follows:

$N_I(x,y)=\frac{1}{4}(1+{\mathrm{\&xi;}}_I\mathrm{\&xi;})(1+{\mathrm{\&eta;}}_I\mathrm{\&eta;})$

${\mathrm{\&Phi;}}_I^{\mathrm{\&alpha;}}(x,y)=\sqrt{r}\&CenterDot;\{\sin \frac{\mathrm{\&theta;}}{2},\cos \frac{\mathrm{\&theta;}}{2},\sin \mathrm{\&theta;}\sin \frac{\mathrm{\&theta;}}{2},\sin \mathrm{\&theta;}\cos \frac{\mathrm{\&theta;}}{2}\}$

(ξ, η) is Gauss coordinate, and (r, θ) is crack tip polar coordinates, (ξ _i, η _i) be node Gauss coordinate;

Step 8, the displacement field differentiate passed through step 7 is tried to achieve obtain strain field, by the constitutive matrix obtained in gained strain field and step 4, utilize physical equation solving model stress field;

Step 9, obtain the elasticity modulus of materials E of current crack tip position _tip, Poisson ratio μ _tipwith fracture toughness K _iCtip, and utilize the stress strength factor K of interactive integral and calculating crack tip _itip=(E _tip/ 2) I, wherein I is interactive integration;

The stress strength factor K of step 10, more current crack tip _itipfracture toughness K corresponding to crack tip position _iCtip; If K _itip<K _iCtip, then load F is increased Δ F, wherein Δ F is preset value, repeats step 6-10; Otherwise Crack Extension, export corresponding expansion load F, crack length is increased Δ a, wherein, Δ a gets the 1%-2% of Initial crack length, judges current crack length a and the maximum crack length a preset _fif: a<a _f, upgrade and split sharp current location, crack tip strengthening set of node N _Λwith crackle section strengthening set of node N _Γ, repeat step 5-10; If a>=a _f, calculate and terminate, obtain the variation tendency increasing required load F with crack length a, export maximum load F in expansion process _max, be the residual intensity of metal gradient material under Initial crack length.

Preferred version further, the present invention containing crackle metal gradient material Prediction model for residual strength method, in described step 3 constitutive matrix C to embody form as follows:

$C=\frac{E(x,y)}{1-\mathrm{\&mu;}{(x,y)}^2}\left[\begin{array}{ccc}1& \mathrm{\&mu;}(x,y)& 0\\ \mathrm{\&mu;}(x,y)& 1& 0\\ 0& 0& \frac{1-\mathrm{\&mu;}(x,y)}{2}\end{array}\right]$

In formula, E (x, y) is the elastic modulus set function of all node respective material, and μ (x, y) is the Poisson ratio set function of all node respective material.

Preferred version further, the present invention is containing crackle metal gradient material Prediction model for residual strength method, and in described step 10, load increment Δ F gets 0.5kN.

The present invention compared with prior art has following significant advantage: 1) use the interpolating function be interrupted to describe Cracking behavior, without the need to cutting at geometry when setting up finite element model, splitting sharp grid without the need to segmentation, having higher precision; 2), after Crack Extension, renewal is only needed to split sharp position, without the need to modeling or grid rezone again, can the propagation behavior of simulating crack easily, realize the dynamic similation of whole expansion process.

Accompanying drawing explanation

Fig. 1 is the process flow diagram that the present invention contains crackle metal gradient material Prediction model for residual strength method;

Fig. 2 is crack tip path of integration schematic diagram;

Fig. 3 is model structure schematic diagram in embodiment;

Fig. 4 is embodiment interior joint strengthening information schematic diagram;

Fig. 5 is prediction expansion load and the comparison of test results figure of embodiment.

Embodiment

As shown in Figure 1, the present invention's one, containing crackle metal gradient material Prediction model for residual strength method, specifically comprises the following steps:

Step 1, build the finite element model of this metal gradient material, using the node of crack tip place unit as splitting point strengthening node in this finite element model, node is strengthened using the node of crackle section place unit as crackle section, all the other nodes are regular node, wherein the node that crosses of crack tip place unit and crackle section place unit is classified as and splits point and strengthen node, and is numbered all nodes;

Step 2, include all point strengthening nodes that splits in crack tip strengthening set of node N _Λ, all crackle section strengthening nodes include crackle section strengthening set of node N in _Γ;

Elastic modulus E, the Poisson ratio μ and fracture toughness K of all node respective material within the scope of step 3, acquisition gradient zones _iClinear changing relation's formula is in the y-direction as follows:

$\left\{\begin{array}{cc}E=E_1+(E_2-E_1)(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\\ \mathrm{\&mu;}={\mathrm{\&mu;}}_1+({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\\ K_{\mathrm{IC}}=K_{\mathrm{IC}1}+(K_{\mathrm{IC}2}-K_{\mathrm{IC}1})(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\end{array}\right.$

In formula, t is Thickness of Gradient Layer, E ₁, μ ₁, K _iC1be respectively the elastic modulus of the first component materials, Poisson ratio and fracture toughness, E ₂, μ ₂, K _iC2be respectively the elastic modulus of the second component materials, Poisson ratio and fracture toughness; With model geometric center for initial point, set up xy coordinate system, wherein y direction is parallel to material gradient direction;

Step 4, according to following formula Confirming model structural stiffness matrix K:

K=∑(∫ _vB(x,y) ^TCB(x,y)dv)

In formula, C is the constitutive matrix of material, and B (x, y) is cell geometry matrix, is known quantity for the model determined, subscript T is transpose of a matrix computing, and v represents that integral domain is individual unit scope; Wherein, constitutive matrix C to embody form as follows:

$C=\frac{E(x,y)}{1-\mathrm{\&mu;}{(x,y)}^2}\left[\begin{array}{ccc}1& \mathrm{\&mu;}(x,y)& 0\\ \mathrm{\&mu;}(x,y)& 1& 0\\ 0& 0& \frac{1-\mathrm{\&mu;}(x,y)}{2}\end{array}\right]$

In formula, E (x, y) is the elastic modulus set function of all node respective material, and μ (x, y) is the Poisson ratio set function of all node respective material;

Step 5, initialization external applied load F=0;

Step 6, solve following finite element discretization equation according to current external applied load F: Ku _i=F

In formula, F is external applied load vector, u _ifor nodal displacement to be solved vector, its component comprises wherein u _ifor regular node displacement, a _ifor N _Γthe additional displacement that the crackle section that collection interior nodes causes due to crackle is discontinuous brought, for N _Λthe additional displacement that the crack tip that collection interior nodes causes due to crackle is discontinuous brought, wherein, α=1,2,3,4, subscript I is node serial number;

Step 7, based on nodal displacement required by step 6, strengthening displacement approximating function interpolation is utilized to obtain displacement field

N represents all node set, N _i(x, y), H _i(x, y) and be respectively standard element positional displacement interpolation function, crackle section interpolating function and crack tip interpolating function, expression is as follows:

$N_I(x,y)=\frac{1}{4}(1+{\mathrm{\&xi;}}_I\mathrm{\&xi;})(1+{\mathrm{\&eta;}}_I\mathrm{\&eta;}),$

${\mathrm{\&Phi;}}_I^{\mathrm{\&alpha;}}(x,y)=\sqrt{r}\&CenterDot;\{\sin \frac{\mathrm{\&theta;}}{2},\cos \frac{\mathrm{\&theta;}}{2},\sin \mathrm{\&theta;}\sin \frac{\mathrm{\&theta;}}{2},\sin \mathrm{\&theta;}\cos \frac{\mathrm{\&theta;}}{2}\}$

Above formula implication is: ${\mathrm{\&Phi;}}_I^1(x,y)=\sqrt{r}\&CenterDot;\sin \frac{\mathrm{\&theta;}}{2},{\mathrm{\&Phi;}}_I^2(x,y)=\sqrt{r}\&CenterDot;\cos \frac{\mathrm{\&theta;}}{2},{\mathrm{\&Phi;}}_I^3(x,y)=\sqrt{r}\&CenterDot;\sin \mathrm{\&theta;}\sin \frac{\mathrm{\&theta;}}{2}$

${\mathrm{\&Phi;}}_I^4(x,y)=\sqrt{r}\&CenterDot;\sin \mathrm{\&theta;}\cos \frac{\mathrm{\&theta;}}{2}$

The exponent number (quantity) got of interpolating function is higher in theory, and checkout result is also accurate, but calculated amount also can improve, generalized case is got 4 rank and is just enough met precision, and (ξ, η) is Gauss coordinate, (r, θ) is crack tip polar coordinates, (ξ _i, η _i) be node Gauss coordinate;

Step 8, the displacement field differentiate passed through step 7 is tried to achieve obtain strain field, by the constitutive matrix obtained in gained strain field and step 4, utilize physical equation solving model stress field, are specially:

Strain field ε _ijby geometric equation obtain, stress field σ _ijby physical equation σ _ij=C _ijklε _klobtain, in formula, tensor C _ijklfor the component of constitutive matrix C, subscript ijkl is tensor expression-form, i=x or i=y; J=x or j=y, the comma in tensor following table represents differentiate, is defined as the parameter before comma and carries out differentiate to the direction after comma, be i.e. u _i,jrepresent that the displacement in i direction is to the differentiate of j direction; Calculating about ess-strain belongs to general knowledge known in this field, and design parameter is shown in: Guo Xiu, Elasticity and tensor analysis [M], Higher Education Publishing House, 2003;

Step 9, according to the linear changing relation in step 3, obtain the elasticity modulus of materials E of current crack tip position _tip, Poisson ratio μ _tipwith fracture toughness K _iCtip, and utilize the stress strength factor K of interactive integral and calculating crack tip _itip=(E _tip/ 2) I, wherein I is interactive integration:

$I={\&Integral;}_{\mathrm{\&Omega;}}({\mathrm{\&sigma;}}_{\mathrm{ij}}{u}_{i,1}^{\mathrm{aux}}+{\mathrm{\&sigma;}}_{\mathrm{ij}}^{\mathrm{aux}}{u}_{i,1}-{\mathrm{\&sigma;}}_{\mathrm{ik}}{\mathrm{\&epsiv;}}_{\mathrm{ik}}^{\mathrm{aux}}{\mathrm{\&delta;}}_{1j})q_{,j}\mathrm{d\&Omega;}+{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&sigma;}}_{\mathrm{ij}}(u_{i,1j}^{\mathrm{aux}}-{\mathrm{\&epsiv;}}_{\mathrm{ij}.1}^{\mathrm{aux}})\mathrm{qd\&Omega;}-{\&Integral;}_{\mathrm{\&Omega;}}{C}_{\mathrm{ijkl},1}{\mathrm{\&epsiv;}}_{\mathrm{ij}}{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{aux}}\mathrm{qd\&Omega;}$

In formula, as shown in Figure 2, Ω splits the annular limit of integration that point is the center of circle, the desirable arbitrary value of internal-and external diameter of annular region, and q is the weight function relevant with path of integration, the peripheral Γ of limit of integration ₀on _q=0, Γ inside limit of integration _supper q=1, crackle section path of integration Γ ₊and Γ _-upper q is the continuously linear change from 0 to 1;

Be assist field containing target on aux, for auxiliary displacement field, for auxiliary stress field, for auxiliary strain field, its expression formula is as follows:

$u_i^{\mathrm{aux}}=\frac{1+{\mathrm{\&mu;}}_{\mathrm{tip}}}{E_{\mathrm{tip}}}\sqrt{\frac{r}{2\mathrm{\&pi;}}}{g}_i\left(\mathrm{\&theta;}\right)$

${\mathrm{\&sigma;}}_{\mathrm{ij}}^{\mathrm{aux}}=\frac{1}{\sqrt{2\mathrm{\&pi;r}}}{f}_{\mathrm{ij}}\left(\mathrm{\&theta;}\right)$

${\mathrm{\&epsiv;}}_{\mathrm{ij}}^{\mathrm{aux}}=S_{\mathrm{ijkl}}{\mathrm{\&sigma;}}_{\mathrm{kl}}^{\mathrm{aux}}$

Wherein, S _ijklfor Flexibility tensor, i.e. constitutive tensor C _ijklinvert, g _i(θ), f _ij(θ) be standard angle transforming function transformation function, represent respectively and split the Changing Pattern of sharp displacement and stress fields along angle, its expression formula is as follows:

$g_i\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}\cos \frac{\mathrm{\&theta;}}{2}(\frac{3-{\mathrm{\&mu;}}_{\mathrm{tip}}}{1+{\mathrm{\&mu;}}_{\mathrm{tip}}}-1+2\mathrm{si}{n}^2\frac{\mathrm{\&theta;}}{2})& i=x\\ \sin \frac{\mathrm{\&theta;}}{2}(\frac{3-{\mathrm{\&mu;}}_{\mathrm{tip}}}{1+{\mathrm{\&mu;}}_{\mathrm{tip}}}-1+\mathrm{co}{s}^2\frac{\mathrm{\&theta;}}{2})& i=y\end{array}\right.$

About utilizing interactive integration also to belong to general knowledge known in this field to calculate crack tip stress intensity factor, specifically see: " research of the metal gradient material crack tip stress intensity factor " paper that Shen Yiwen delivers for 2012, repeats no more herein! Preset about the tensor in computation process can see Guo Xiu simultaneously, Elasticity and tensor analysis [M], Higher Education Publishing House, 2003;

The stress strength factor K of step 10, more current crack tip _itipfracture toughness K corresponding to crack tip position _iCtip; If K _itip<K _iCtip, then load F is increased Δ F, wherein increment Delta F can be preset value, also can pass through fundamental strength theoretical calculation model breakdown strength magnitude, get the 1%-2% of this magnitude, repeats step 6-10; Otherwise Crack Extension, export corresponding expansion load F, crack length is increased Δ a, wherein, Δ a gets the 1%-2% of Initial crack length, judges the maximum crack length a of current crack length a and permission _fif: a<a _f, upgrade and split sharp current location, crack tip strengthening set of node N _Λwith crackle section strengthening set of node N _Γ, repeat step 5-10; If a>=a _f, calculate and terminate, the variation tendency increasing required load F with crack length a, exports maximum load F in expansion process _max, under being Initial crack length, the residual intensity of metal gradient material.

Embodiment

Step 1, set up finite element model as Fig. 3, model height H=24mm, material gradient region is positioned at high degree of symmetry place its thickness t=1mm, span S=96mm between fulcrum.The present embodiment process edge crack, therefore only have 1 crack tip, mark x place is respectively splits point and crackle starting point, and crackle starting point is positioned at border mid point on the downside of model, and splitting point with starting point spacing is Initial crack length a ₀=8.6mm, definition the present embodiment fracture-arrest length a _f=15mm.

Step 2, use rectangular element carry out discrete to model, and the grid length of side is 1mm, amounts to 2640 unit; Processing node strengthening information, as Fig. 4, is included the cell node relevant to splitting point in crack tip and is strengthened set of node N _Λ, include the cell node relevant to crackle section in crackle section and strengthen set of node N _Γ.

Step 3, model have 2 kinds of component materials, and its performance is as follows

By the model parameter of establishment of coordinate system shown in Fig. 3 field, material property linear change in gradient layer, because bi-material Poisson ratio is identical, in whole model, Poisson ratio is definite value:

$\left\{\begin{array}{cc}E=E_1& -12\&le;y<-0.5\\ E=E_1+(E_2-E_1)(y+0.5)& -0.5\&le;y<0.5\\ E=E_2& 0.5\&le;y\&le;12\end{array}\right.$

$\left\{\begin{array}{cc}{K}_{\mathrm{IC}}=K_{\mathrm{IC}1}& -12\&le;y<-0.5\\ K_{\mathrm{IC}}=K_{\mathrm{IC}1}+(K_{\mathrm{IC}2}-K_{\mathrm{IC}1})(y+0.5)& -0.5\&le;y<0.5\\ K_{\mathrm{IC}}=K_{\mathrm{IC}2}& 0.5\&le;y\&le;12\end{array}\right.$

Step 4, by $C=\frac{E(x,y)}{1-\mathrm{\&mu;}{(x,y)}^2}\left[\begin{array}{ccc}1& \mathrm{\&mu;}(x,y)& 0\\ \mathrm{\&mu;}(x,y)& 1& 0\\ 0& 0& \frac{1-\mathrm{\&mu;}(x,y)}{2}\end{array}\right]$ Calculate constitutive matrix, by K=∑ (∫ _vb (x, y) ^tcB (x, y) dv) computation model stiffness matrix, B is geometric matrix, and the present embodiment uses 4 degree Rectangular Elements, then B=[B ₁b ₂b ₃b ₄], for different node, submatrix expression formula is as follows:

Regular node: $B_i=\left[\begin{array}{cc}{N}_{i,x}& 0\\ 0& N_{i,y}\\ N_{i,y}& N_{i,x}\end{array}\right]$

Crackle section strengthening node: $B_i=\left[\begin{array}{cc}{\left(N_i{H}_i\right)}_{,x}& 0\\ 0& {\left(N_i{H}_i\right)}_{,y}\\ {\left(N_i{H}_i\right)}_{,y}& {\left(N_i{H}_i\right)}_{,x}\end{array}\right]$

Crack tip strengthening node: $B_i=\left[\begin{array}{cccc}{B}_i^1& B_i^1& B_i^3& B_i^4\end{array}\right]$

$B_i^{\mathrm{\&alpha;}}=\left[\begin{array}{cc}{\left(N_i{\mathrm{\&Phi;}}_i^{\mathrm{\&alpha;}}\right)}_{,x}& 0\\ 0& {\left(N_i{\mathrm{\&Phi;}}_i^{\mathrm{\&alpha;}}\right)}_{,y}\\ {\left(N_i{\mathrm{\&Phi;}}_i^{\mathrm{\&alpha;}}\right)}_{,y}& {\left(N_i{\mathrm{\&Phi;}}_i^{\mathrm{\&alpha;}}\right)}_{,x}\end{array}\right]\mathrm{\&alpha;}=\mathrm{1,2,3,4}$

In formula, for known models, i is four nodes of known units, N _i, H _iwith for positional displacement interpolation function, its expression is shown in step 7;

Step 5, initialization load F=0;

Step 6, according to current external applied load F solving finite element discrete equation Ku _i=F obtains nodal displacement vector u _i; Step 7, based on nodal displacement u required by step 6 _i, utilize strengthening displacement approximating function interpolation to obtain displacement field

Each interpolating function expression formula is as follows:

$N_I(x,y)=\frac{1}{4}(1+{\mathrm{\&xi;}}_I\mathrm{\&xi;})(1+{\mathrm{\&eta;}}_I\mathrm{\&eta;}),$

${\mathrm{\&Phi;}}_I^1(x,y)=\sqrt{r}\&CenterDot;\sin \frac{\mathrm{\&theta;}}{2},{\mathrm{\&Phi;}}_I^2(x,y)=\sqrt{r}\&CenterDot;\cos \frac{\mathrm{\&theta;}}{2}$

${\mathrm{\&Phi;}}_I^{3}x,y=\sqrt{r}\&CenterDot;\sin \mathrm{\&theta;}\sin \frac{\mathrm{\&theta;}}{2},{\mathrm{\&Phi;}}_I^4(x,y)=\sqrt{r}\&CenterDot;\sin \mathrm{\&theta;}\cos \frac{\mathrm{\&theta;}}{2}$

Step 8, by geometric equation obtain and solve strain field ε _ij, by physical equation σ _ij=C _ijklε _klsolve stress field σ _ij;

Step 9, obtain the elasticity modulus of materials E of current crack tip position _tip, Poisson ratio μ _tipwith fracture toughness K _iCtip, as when crack length is 12mm, splitting sharp coordinate is (0,0), now E _tip=117.5GPa, μ _tip=0.33, $K_{\mathrm{ICtip}}=77.5\mathrm{MPa}\sqrt{m};$

By formula K _itip=(E _tip/ 2) I solves the crack tip stress intensity factor, and wherein I is interactive integration:

$I={\&Integral;}_{\mathrm{\&Omega;}}({\mathrm{\&sigma;}}_{\mathrm{ij}}{u}_{i,1}^{\mathrm{aux}}+{\mathrm{\&sigma;}}_{\mathrm{ij}}^{\mathrm{aux}}{u}_{i,1}-{\mathrm{\&sigma;}}_{\mathrm{ik}}{\mathrm{\&epsiv;}}_{\mathrm{ik}}^{\mathrm{aux}}{\mathrm{\&delta;}}_{1j})q_{,j}\mathrm{d\&Omega;}+{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&sigma;}}_{\mathrm{ij}}(u_{i,1j}^{\mathrm{aux}}-{\mathrm{\&epsiv;}}_{\mathrm{ij}.1}^{\mathrm{aux}})\mathrm{qd\&Omega;}-{\&Integral;}_{\mathrm{\&Omega;}}{C}_{\mathrm{ijkl},1}{\mathrm{\&epsiv;}}_{\mathrm{ij}}{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{aux}}\mathrm{qd\&Omega;}$

For the present embodiment, limit of integration internal diameter gets 1/3 of crack length, and limit of integration external diameter gets 2/3 of crack length;

The stress strength factor K of step 10, more current crack tip _itipfracture toughness K corresponding to crack tip position _iCtip; If K _itip<K _iCtip, increase load, the present embodiment load increment gets 0.5kN, repeats step 6-10; For each LOAD FOR result of crack length 12mm as table 1, when load is less than 14kN, crackle is not expanded: table 1:

If K _itip>=K _iCtipcrack Extension, exports corresponding expansion load F, for upper table result of calculation, meets expansion criterion when load is 14kN, therefore when crack length is 12mm, its expansion load is 14kN.

Increase crack length, the present embodiment crack length increment gets 0.2mm, judges current crack length a and the maximum crack length a preset _fif: a<a _f, upgrade and split sharp current location, crack tip strengthening set of node N _Λwith crackle section strengthening set of node N _Γ, repeat step 5-10; If a>=a _f, calculate and terminate.

Calculate expansion load under each crack length test findings as 3 testpieces in returning in Fig. 5, figure, contrast finds that analog result and test findings are coincide better.

The maximum load of getting in crack propagation process is model residual intensity, and result of calculation and test findings are as table 2, and Relative Error, within 10%, meets the requirement to prediction of strength in engineering.

Table 2

Claims (3)

1., containing a crackle metal gradient material Prediction model for residual strength method, the method obtains current containing crackle metal gradient material residual intensity by carrying out analysis to the expansion process of I mode-Ⅲ crack, it is characterized in that, specifically comprises the following steps:

Step 1, build the finite element model of this metal gradient material, using the node of crack tip place unit as splitting point strengthening node in this finite element model, node is strengthened using the node of crackle section place unit as crackle section, all the other nodes are regular node, wherein the node that crosses of crack tip place unit and crackle section place unit is classified as and splits point and strengthen node, and is numbered all nodes;

Step 2, include all point strengthening nodes that splits in crack tip strengthening set of node N _Λ, all crackle section strengthening nodes include crackle section strengthening set of node N in _Γ;

Elastic modulus E, the Poisson ratio μ and fracture toughness K of all node respective material within the scope of step 3, acquisition gradient zones _iClinear changing relation's formula is in the y-direction as follows:

$\left\{\begin{array}{cc}E=E_1+(E_2-E_1)(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\\ \mathrm{\&mu;}={\mathrm{\&mu;}}_1+({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\\ K_{IC}=K_{IC1}+(K_{IC2}-K_{IC1})(y+\frac{t}{2})& -\frac{t}{2}\&le;y<\frac{t}{2}\end{array}\right.$

In formula, t is Thickness of Gradient Layer, E ₁, μ ₁, K _iC1be respectively the elastic modulus of the first component materials, Poisson ratio and fracture toughness, E ₂, μ ₂, K _iC2be respectively the elastic modulus of the second component materials, Poisson ratio and fracture toughness; With model geometric center for initial point, set up xy coordinate system, wherein y direction is parallel to material gradient direction;

Step 4, according to following formula Confirming model structural stiffness matrix K:

K＝∑(∫ _vB(x,y) ^TCB(x,y)dv)

In formula, C is the constitutive matrix of material, and B (x, y) is cell geometry matrix, and subscript T is transpose of a matrix computing, and v represents that integral domain is individual unit scope;

Step 5, initialization external applied load F=0;

Step 6, solve following finite element discretization equation according to current external applied load F:

Ku _I＝F

In formula, F is external applied load vector, u _ifor nodal displacement to be solved vector, its component comprises wherein u _ifor regular node displacement, a _ifor N _Γthe additional displacement that the crackle section that collection interior nodes causes due to crackle is discontinuous brought, for N _Λthe additional displacement that the crack tip that collection interior nodes causes due to crackle is discontinuous brought, wherein, α=1,2,3,4, subscript I is node serial number;

Step 7, based on nodal displacement required by step 6, strengthening displacement approximating function interpolation is utilized to obtain displacement field

N represents all node set, N _i(x, y), H _i(x, y) and be respectively standard element positional displacement interpolation function, crackle section interpolating function and crack tip interpolating function, expression is as follows:

$N_I(x,y)=\frac{1}{4}(1+{\mathrm{\&xi;}}_I\mathrm{\&xi;})(1+{\mathrm{\&eta;}}_I\mathrm{\&eta;})$

${\mathrm{\&Phi;}}_I^{\mathrm{\&alpha;}}(x,y)=\sqrt{r}\&CenterDot;\{sin\frac{\mathrm{\&theta;}}{2},cos\frac{\mathrm{\&theta;}}{2},sin\mathrm{\&theta;}sin\frac{\mathrm{\&theta;}}{2},sin\mathrm{\&theta;}cos\frac{\mathrm{\&theta;}}{2}\}$

(ξ, η) is Gauss coordinate, and (r, θ) is crack tip polar coordinates, (ξ _i, η _i) be node Gauss coordinate;

Step 8, the displacement field differentiate passed through step 7 is tried to achieve obtain strain field, then by the constitutive matrix obtained in gained strain field and step 4, utilize physical equation solving model stress field;

Step 9, obtain the elasticity modulus of materials E of current crack tip position _tip, Poisson ratio μ _tipwith fracture toughness K _iCtip, and utilize the stress strength factor K of interactive integral and calculating crack tip _itip=(E _tip/ 2) I, wherein I is interactive integration;

The stress strength factor K of step 10, more current crack tip _itipfracture toughness K corresponding to crack tip position _iCtipif: K _itip<K _iCtip, then load F is increased Δ F, wherein Δ F is preset value, repeats step 6-10; Otherwise Crack Extension, export corresponding expansion load F, crack length is increased Δ a, wherein, Δ a gets the 1%-2% of Initial crack length, judges current crack length a and the maximum crack length a preset _fif: a<a _f, upgrade and split sharp current location, crack tip strengthening set of node N _Λwith crackle section strengthening set of node N _Γ, repeat step 5-10; If a>=a _f, calculate and terminate, obtain the variation tendency increasing required load F with crack length a, maximum load F in the expansion process of output _max, be the residual intensity of metal gradient material under Initial crack length.

2., according to claim 1 containing crackle metal gradient material Prediction model for residual strength method, it is characterized in that, in described step 4 constitutive matrix C to embody form as follows:

$C=\frac{E(x,y)}{1-\mathrm{\&mu;}{(x,y)}^2}\left[\begin{array}{ccc}1& \mathrm{\&mu;}(x,y)& 0\\ \mathrm{\&mu;}(x,y)& 1& 0\\ 0& 0& \frac{1-\mathrm{\&mu;}(x,y)}{2}\end{array}\right]$

In formula, E (x, y) is the elastic modulus set function of all node respective material, and μ (x, y) is the Poisson ratio set function of all node respective material.

3., according to claim 1 containing crackle metal gradient material Prediction model for residual strength method, it is characterized in that, described Δ F=0.5kN.