Impact fatigue finite element simulation method considering scale effect

Abstract

The invention relates to an impact fatigue finite element simulation method considering a scale effect, and belongs to the technical field of material damage failure simulation analysis. According to the invention, a constitutive model suitable for advanced materials and materials prepared by an advanced process is established, static characteristics, dynamic characteristics, scale effects and damage evolution processes of the materials are combined, and impact damage and impact fatigue processes of the advanced materials can be effectively simulated; a non-localization cross-scale constitutive method is adopted, the efficiency requirement of engineering finite element simulation can be met, meanwhile, localization non-uniform modeling is compatible so as to face more accurate scientific research calculation, and good balance is achieved between precision and efficiency; by means of the method, integrated calculation in commercial finite element software can be achieved, impact damage and fatigue simulation considering the material scale effect can be more efficiently conducted, and application of advanced materials in the aerospace field impact and reciprocating impact load environment is promoted.

Description

Technical Field

[0001] The present invention relates to an impact fatigue finite element simulation method considering scale effect, which is a simulation method for realizing cross-scale impact damage and fatigue simulation by using a finite element program and a plug-in, and belongs to the technical field of material damage failure simulation analysis. Background Art

[0002] Impact loads and impact fatigue loads are among the important working conditions faced by structures in the aerospace field. These structures are related to safety while also affecting their efficiency. Currently, many advanced manufacturing processes are used on these structures. Advanced materials prepared by manufacturing and processing processes such as additive manufacturing and surface strengthening have typical multi-scale microstructure characteristics, which in turn brings about significant material scale effects. Dislocation evolution often dominates the mechanical properties of metals in the plastic stage and the characteristics of their scale effects. Taylor's dislocation theory and the strain gradient plasticity theory based on it well describe the scale effects of metal materials and have been applied and verified in many scenarios and working conditions.

[0003] The unique microstructures in these advanced materials will produce higher strain gradients in the materials, which can often improve their strength properties such as static strength, plastic stiffness, and hardness. At the same time, the high stress level introduced by the microstructure and the increase in initial defects such as structural interfaces often reduce the toughness properties such as ductility and impact toughness and fatigue properties of the materials, and present a complex antagonistic relationship. Simulating the static, dynamic, damage and scale effects of materials at the same time is still a technical difficulty. The existing simulation methods for impact damage and fatigue have certain shortcomings when simulating advanced materials, and there are no reports on the cross-scale mechanical properties under impact fatigue loads. To this end, a finite element simulation method for impact fatigue considering scale effects is proposed to further study and utilize the impact fatigue resistance of advanced materials. Summary of the invention

[0004] The purpose of this invention is to propose a constitutive model based on cross-scale static and dynamic structure-activity relationship and damage evolution relationship, and to develop a finite element simulation method for impact fatigue considering scale effect, so as to realize a simulation calculation method suitable for the research and engineering application of cross-scale impact damage and fatigue of advanced materials, and to provide verification methods and data references for the application of advanced materials in fields such as aerospace.

[0005] The present invention achieves the above-mentioned purpose through the following technical scheme, which is a finite element simulation method of impact fatigue considering scale effect, and mainly includes the following steps:

[0006] S1: Establish and assemble the geometric model according to the geometric characteristics of the research object;

[0007] S2: Obtain material parameters in the constitutive relationship by consulting the material manual and conducting experiments, and write the constitutive relationship into the VUMAT program using FORTRAN language;

[0008] S3: Write the calculation procedures of strain gradient and element integration point shape function in VEXTERNALDB and VUSDFLD programs and merge them into VUMAT program;

[0009] S4: Boundary conditions and contact conditions are established according to the research questions of the simulation experiment, and the fatigue analysis defines the loading process of fatigue load through VUAMP program;

[0010] S5: Divide the finite element mesh according to the geometry and stress characteristics of the research object, and use the first-order fully integrated hexahedral unit or the second-order reduced integrated hexahedral unit to ensure that each unit contains eight integration points, i.e., material points, to realize the calculation of strain gradient within the unit; complete the impact fatigue finite element modeling considering the scale effect;

[0011] S6: Perform finite element simulation calculations and extract the required results from the output file.

[0012] Furthermore, in step S2, the material parameters in the constitutive relationship are obtained by querying the material manual and assisted by experiments, wherein the Young's modulus, yield strength, work hardening coefficient, yield limit, and rate hardening parameter of the material are obtained by querying the material manual, and the scale effect parameters and damage parameters of the material are obtained experimentally. The scale effect parameters and damage parameters can also be selected from literature reports, or all of the above parameters can be obtained through experiments.

[0013] Furthermore, in step S1, the geometric model is subjected to non-uniform modeling to simulate the specific geometric configuration of the microstructure, thereby reflecting the material microstructure effect.

[0014] Furthermore, the geometric model is homogenized and the microstructural effect of the material is reflected by applying different scale effect parameters.

[0015] Furthermore, in step S2, the VUMAT program can be converted from the corresponding UMAT program, or can be directly written according to the constitutive relationship.

[0016] Furthermore, the S5 divides the finite element mesh according to the geometry and stress characteristics of the research object, encrypts the mesh in the area with high stress level or geometric mutation or contact area, and determines the minimum mesh size according to the material wave velocity to ensure the model convergence speed.

[0017] Furthermore, in step S2, the specific implementation process of the VUMAT program is as follows:

[0018] S21: Call the VUMAT program to perform elastic calculations during initialization, and call the VUSDFLD program to obtain the number and coordinate x of each global unit integration point i ;

[0019] S22: Call VEXTERNALDB to calculate the shape function N in each global unit integration point i , Where ξ, ψ, η are the local coordinates of the integration point in the element, taking values ​​of ±1; and the plastic strain gradient η is calculated p , which is calculated as

[0020] where ε p is the plastic strain tensor; i, η, k are coordinate numbers.

[0021] S23: VUMAT obtains the above strain gradient η through global variables p , and update the equivalent stress σ according to the increment of displacement and time of the fatigue load input defined by the constitutive relationship and VUAMP eq , equivalent stress σ eq The update includes the equivalent strain rate Strain gradient η, plastic strain rate Damage degree D and its rate The common effects of In the above formulas, σ flow is the flow stress, σ 0 is the yield strength, E is Young's modulus, n is the power hardening coefficient, l is the scale effect parameter strain gradient characteristic length, μ is the shear modulus, b is the modulus of the Burgers vector, C is the rate hardening parameter, S and m are damage parameters, and α is the scale effect parameter;

[0023] S24: Perform Newton iteration until convergence, update stress, update damage value, traverse each integral point, and enter the next incremental step;

[0024] S25: Complete all design increment steps and analysis steps, and complete finite element calculations.

[0025] Furthermore, in step S23, for materials with complex microstructures and which are significantly capable of producing stronger scale effects in the material, the influence of back stress on flow stress is further considered, that is, The back stress σ b It is the long-range effect of the geometric dislocation GND, that is, the gradient of the GND density. The back stress σ is derived by the following method and formula b Calculation of

[0026] Near any material point, the GND density is expanded according to Taylor's formula and only the first-order linear term is retained, and the local coordinates are obtained. in is the GND density at the material point, is the Nye constant, b is the modulus of the Burgers vector; back stress σ b As a manifestation of the long-range effect of GND, each GND in the material will affect the flow stress at any material point. However, considering the attenuation of the force with distance and the convenience of calculation and solution, the influence of GND on the flow stress in a spherical domain with a radius of R is selected. Then the stress field component of a single z-direction edge dislocation at any material point is σ xz =0, σ yz =0, where v is Poisson's ratio; the integral of the stress field of a single edge dislocation in the z direction at any material point in the sphere with a radius of R is σ xz =0, σ yz =0, And for The item Then we have σ xx =σ yy =σ zz =0, that is, the amount There is only one non-zero term in There is only one non-zero term in the other two directions, which is The back stress can be related to the GND density gradient in the form of Mises equivalent stress: in The calculation is performed according to the method for calculating the strain gradient in step S22. After the back stress value is calculated, the flow stress can be continuously updated according to step S23.

[0027] Beneficial effects of the present invention

[0028] 1) A set of constitutive models suitable for advanced materials and materials made by advanced processes has been established. The model is based on a variety of mature structure-activity relationships with real physical meanings and processes, and combines the static properties, dynamic properties, scale effects, and damage evolution processes of materials. It can effectively simulate the impact damage and impact fatigue processes of advanced materials;

[0029] 2) The non-localized cross-scale constitutive method is adopted, which can meet the efficiency requirements of engineering finite element simulation, and is also compatible with localized non-uniform modeling for more accurate scientific research calculations, achieving a good balance between accuracy and efficiency;

[0030] 3) A plug-in program suitable for mature commercial finite element software was written, making full use of traditional finite element software to simulate related problems of advanced materials.

[0031] In summary, the present invention establishes an impact damage constitutive model that takes scale effects into consideration based on mature theories, and forms a FORTRAN program that utilizes the constitutive model and fatigue analysis process. This can realize integrated calculations in commercial finite element software, and more efficiently simulate impact damage and fatigue that takes material scale effects into consideration, thereby promoting the application of advanced materials in impact and reciprocating impact load environments in the aerospace field. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is a schematic diagram of the process of the finite element simulation applicable to cross-scale impact fatigue of the present invention;

[0033] Figure 2 It is a flowchart of the program formed by the constitutive relationship of the present invention;

[0034] Figure 3 is a model diagram in an impact fatigue embodiment of the present invention;

[0035] Figure 4 is a stress cloud diagram in an impact fatigue embodiment of the present invention;

[0036] Figure 5 is a damage cloud diagram in an impact fatigue embodiment of the present invention;

[0037] Figure 6 It is a damage evolution relationship diagram in the impact fatigue embodiment of the present invention. DETAILED DESCRIPTION

[0038] The present invention is described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0039] like Figure 1-2 As shown, a finite element simulation method for impact fatigue considering scale effect mainly includes the following steps:

[0040] S1: Establish and assemble a geometric model based on the geometric features of the research object; the geometric model is subjected to non-homogeneous modeling to simulate the specific geometric configuration of the microstructure, thereby reflecting the microstructural effect of the material. The geometric model is subjected to homogenous modeling to reflect the microstructural effect of the material by applying different scale effect parameters.

[0041] S2: query the material manual and obtain the material parameters in the constitutive relationship through experiments, and write the constitutive relationship into the VUMAT program using FORTRAN language; the query material manual and obtain the material parameters in the constitutive relationship through experiments, wherein the Young's modulus, yield strength, work hardening coefficient, yield limit, rate hardening parameter of the material are obtained through querying the material manual, and the scale effect parameters and damage parameters of the material obtained through experiments can also be selected from literature reports, or all the above parameters can be obtained through experiments. The VUMAT program can be converted from the corresponding UMAT program, or it can be directly written according to the constitutive relationship.

[0042] The specific implementation process of the VUMAT program involved is as follows:

[0043] S21: Call the VUMAT program to perform elastic calculations during initialization, and call the VUSDFLD program to obtain the number and coordinate x of each global unit integration point i ;

[0044] S22: Call VEXTERNALDB to calculate the shape function n in each element integration point globally i , Where ξ, ψ, η are the local coordinates of the integration point in the element, taking values ​​of ±1; and the plastic strain gradient η is calculated p , which is calculated as

[0045] where ε p is the plastic strain tensor; i, j, k are the coordinate numbers.

[0046] S23: VUMAT obtains the above strain gradient η through global variables p , and update the equivalent stress σ according to the increment of displacement and time of the fatigue load input defined by the constitutive relationship and VUAMP eq , equivalent stress σ eq The update includes the equivalent strain rate Strain gradient η, plastic strain rate Damage degree D and its rate The common effects of In the above formulas, σ flow is the flow stress, σ 0is the yield strength, E is Young's modulus, n is the power hardening coefficient, l is the scale effect parameter strain gradient characteristic length, μ is the shear modulus, b is the modulus of the Burgers vector, C is the rate hardening parameter, S and m are damage parameters, and α is the scale effect parameter. Most of the above physical quantities have specific physical meanings and can be obtained by referring to manuals and literature. The physical meanings of the last four parameters are not strictly clear but have been recognized by the academic community and can be obtained by referring to literature or through experiments.

[0047] For materials with complex microstructures and significant ability to produce stronger scale effects in the material, the effect of back stress on flow stress is further considered, that is, The back stress σ b It is the long-range effect of the geometric dislocation GND, that is, the gradient of the GND density. The back stress σ is derived by the following method and formula b Calculation of

[0048] Near any material point, the GND density is expanded according to Taylor's formula and only the first-order linear term is retained, and the local coordinates are obtained. in is the GND density at the material point, is the Nye constant, b is the modulus of the Burgers vector; back stress σ b As a manifestation of the long-range effect of GND, each GND in the material will affect the flow stress at any material point. However, considering the attenuation of the force with distance and the convenience of calculation and solution, the influence of GND on the flow stress in a spherical domain with a radius of R is selected. Then the stress field component of a single z-direction edge dislocation at any material point is σ xz =0, σ yz =0, Where v is Poisson's ratio. Since the appearance of GND is to coordinate the uneven deformation in the material, it is inferred that its optimal arrangement is that the direction in which it appears is consistent with the uneven deformation that needs to be coordinated. The overall distribution of GND roughly follows this law, so that deformation coordination can be achieved at a smaller GND density.

[0049] Based on this, the integral of the stress field of a single edge dislocation in the z direction at any material point in the sphere with a radius of R is σ xz =0, σ yz =0, And for The item Then we have σ xx =σ yy =σ zz=0, that is, the amount There is only one non-zero term in There is only one non-zero term in the other two directions, which is The back stress can be related to the GND density gradient in the form of Mises equivalent stress: in The calculation is performed according to the method for calculating the strain gradient in step S22. After the back stress value is calculated, the flow stress can be continuously updated according to step S23.

[0050] S24: Perform Newton iteration until convergence, update stress, update damage value, traverse each integral point, and enter the next incremental step;

[0051] S25: Complete all design increment steps and analysis steps, and complete finite element calculations.

[0052] S3: Write the calculation procedures of strain gradient and element integration point shape function in VEXTERNALDB and VUSDFLD programs and merge them into VUMAT program;

[0053] S4: Boundary conditions and contact conditions are established according to the research questions of the simulation experiment, and the fatigue analysis defines the loading process of fatigue load through VUAMP program;

[0054] S5: Divide the finite element mesh according to the geometry and force characteristics of the research object, and use first-order fully integrated hexahedral units or second-order reduced integrated hexahedral units to ensure that each unit contains eight integration points, i.e., material points, to realize the calculation of strain gradients within the unit; complete the impact fatigue finite element modeling considering the scale effect; the S5 divides the finite element mesh according to the geometry and force characteristics of the research object, encrypts the mesh in areas with high stress levels or geometric mutations or contact areas, and determines the minimum mesh size according to the material wave velocity to ensure the model convergence speed.

[0055] S6: Perform finite element simulation calculations and extract the required results from the output file.

[0056] like Figure 3 As shown, the present invention is directed to an impact fatigue finite element simulation taking into account scale effect, and the implementation method is as follows: a specimen simulation model is designed according to the national Charpy impact standard (GB / T229-2007; GB / T12778-2008), the specimen is 55 mm long, 10 mm wide, and 10 mm high, is a 45° V-notch specimen, has a notch depth of 2 mm, and a notch chamfer radius of 0.25 mm, and the model adopts a symmetrical model in the width direction.

[0057] The material used in the embodiment is high-strength steel, with material parameters: Young's modulus 210GPa, Poisson's ratio 0.3, yield strength 1500MPa, power hardening coefficient 0.1, viscoplastic coefficient 20, strain gradient characteristic length 0.015mm, rate hardening parameter 0.017, Lemaitre damage model parameters 12 and 10.101, and these parameters are passed into programs such as VUMAT: Call VEXTERNALDB to calculate the shape function D in each unit integration point globally i , And calculate the strain gradient η, which is calculated as follows:

[0058] where ε p is the plastic strain tensor. VUMAT obtains the above strain gradients through global variables, and updates the stress σ according to the deformation and time increments passed in by the constitutive relationship and the load defined by VUAMP. eq , the stress update will include strain ε, strain gradient η, strain rate The common effects of damage degree D are as follows: For materials with complex microstructures, further consideration is given In the above formulas, σ flow is the flow stress, σ 0 is the yield strength, E is Young's modulus, n is the power hardening coefficient, l is the scale effect parameter strain gradient characteristic length, C is the rate hardening parameter, S and m are damage parameters, which have been given in this embodiment according to high strength steel.

[0059] The embodiment is a Charpy impact condition with state variable controlled damage deletion, so global contact needs to be established. Other boundary conditions include: the bottom support is rigidly fixed, all degrees of freedom outside the impact direction of the impact hammer are constrained and set to rigid, and symmetry constraints are imposed on the sample on the symmetry plane.

[0060] The embodiment is an impact fatigue example. A cyclic load is introduced through the VUAMP program. The load peak is 15000N. Under symmetric boundary conditions, a cyclic load with a peak value of 7500N is applied.

[0061] The finite element mesh divided in the embodiment has the highest mesh density in the notch area, the minimum mesh side length in the area is 0.05 mm, the material is high-strength steel, the wave speed is 5000 m / s, so the convergence step of the model is 10 -8 s, the element type is a first-order fully integrated hexahedral element with 8 integration points.

[0062] After completing the modeling of the embodiment, the relevant results can be calculated and obtained. The stress cloud diagram is as follows: Figure 4 As shown, the damage cloud diagram is as follows Figure 5 The damage evolution process is shown in Figure 6 shown.

[0063] It should be noted that, according to the needs of implementation, the steps described in this embodiment can be split and combined into new steps to achieve the purpose of the present invention. It is easy for those skilled in the art to understand that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modification, equivalent replacement and improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A finite element simulation method for impact fatigue considering scale effect, characterized in that: The following steps are involved: S1: Establish and assemble the geometric model according to the geometric characteristics of the research object; S2: Obtain material parameters in the constitutive relationship by consulting the material manual and conducting experiments, and write the constitutive relationship into the VUMAT program using FORTRAN language; S3: Write the calculation procedures of strain gradient and element integration point shape function in VEXTERNALDB and VUSDFLD programs and merge them into VUMAT program; S4: Boundary conditions and contact conditions are established according to the research questions of the simulation experiment, and the fatigue analysis defines the loading process of fatigue load through VUAMP program; S5: Divide the finite element mesh according to the geometry and stress characteristics of the research object, and use the first-order fully integrated hexahedral unit or the second-order reduced integrated hexahedral unit to ensure that each unit contains eight integration points, i.e., material points, to realize the calculation of strain gradient within the unit; complete the impact fatigue finite element modeling considering the scale effect; S6: Perform finite element simulation calculations and extract the required results from the output file.

2. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: In the step S2, the material parameters in the constitutive relationship are obtained by querying the material manual and assisted by experiments, wherein the material manual is queried to obtain the Young's modulus, yield strength, work hardening coefficient, yield limit, and rate hardening parameter of the material, and the scale effect parameter and damage parameter of the material are obtained experimentally. The scale effect parameter and damage parameter may be selected from literature reports, or all of the above parameters may be obtained through experiments.

3. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: In the step S1, the geometric model is subjected to non-uniform modeling to simulate the specific geometric configuration of the microstructure, thereby reflecting the material microstructure effect.

4. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: In the step S1, the geometric model is homogenized and modeled, and the microstructural effect of the material is reflected by applying different scale effect parameters.

5. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: In step S2, the VUMAT program can be converted from the corresponding UMAT program, or can be directly written according to the constitutive relationship.

6. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: The S5 divides the finite element grid according to the geometry and force characteristics of the research object, encrypts the grid in the area with high stress level or geometric mutation or contact area, and determines the minimum grid size according to the material wave velocity to ensure the model convergence speed.

7. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: In step S2, the specific process of implementing the VUMAT program involved is as follows: S21: Call the VUMAT program to perform elastic calculations during initialization, and call the VUSDFLD program to obtain the number and coordinate x of each global unit integration point i ; S22: Call VEXTERNALDB to calculate the shape function N in each global unit integration point i , i = 1-8, where ξ, ψ, η are the local coordinates of the integration point in the unit, taking values ​​of ±1; and calculate the plastic strain gradient η p , which is calculated as where ε p is the plastic strain tensor, i,η,k are the coordinate numbers; S23: VUMAT program obtains the strain gradient η through global variables p , and updates the equivalent stress σ according to the increment of displacement and time of the fatigue load defined by the constitutive relationship and VUAMP program eq , equivalent stress σ eq The update includes the equivalent strain rate Strain gradient η, plastic strain rate Damage degree D and its rate The common effects of In the above formulas, σ flow is the flow stress, σ0 is the yield strength, E is the Young's modulus, n is the power hardening coefficient, l is the scale effect parameter strain gradient characteristic length, μ is the shear modulus, b is the modulus of the Burgers vector, C is the rate hardening parameter, S and m are damage parameters, and α is the scale effect parameter; S24: Perform Newton iteration until convergence, update stress, update damage value, traverse each integral point, and enter the next incremental step; S25: Complete all design increment steps and analysis steps, and complete finite element calculations.

8. The impact fatigue finite element simulation method considering scale effect according to claim 1 is characterized in that: In step S23, for materials with complex microstructures and significantly stronger scale effects in the material, the influence of back stress on flow stress is further considered, that is, The back stress σ b It is the long-range effect of the geometric dislocation GND, that is, the gradient of the GND density. The back stress σ is derived by the following method and formula b Calculation of Near any material point, the GND density is expanded according to Taylor's formula and only the first-order linear term is retained, and the local coordinates are obtained. in is the GND density at the material point, is the Nye constant, b is the modulus of the Burgers vector; back stress σ b As a manifestation of the long-range effect of GND, each GND in the material will affect the flow stress at any material point. However, considering the attenuation of the force with distance and the convenience of calculation and solution, the influence of GND on the flow stress in a spherical domain with a radius of R is selected. Then the stress field component of a single z-direction edge dislocation at any material point is σ xz =0, σ yz =0, where ν is Poisson's ratio; the integral of the stress field of a single edge dislocation in the z direction at any material point in the sphere with a radius of R is σ xz =0, σ yz =0, And for The item Then we have σ xx =σ yy =σ zz =0, that is, the amount There is only one non-zero term in There is only one non-zero term in the other two directions, which is The back stress can be related to the GND density gradient in the form of Mises equivalent stress: in The calculation is performed according to the method for calculating the strain gradient in step S22; after the back stress value is calculated, the flow stress can be continuously updated according to step S23.

Images