CN112100820A - Excel-based anisotropic yield function parameter calibration method - Google Patents

Excel-based anisotropic yield function parameter calibration method Download PDF

Info

Publication number
CN112100820A
CN112100820A CN202010865152.3A CN202010865152A CN112100820A CN 112100820 A CN112100820 A CN 112100820A CN 202010865152 A CN202010865152 A CN 202010865152A CN 112100820 A CN112100820 A CN 112100820A
Authority
CN
China
Prior art keywords
calibration
yield function
value
yield
function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202010865152.3A
Other languages
Chinese (zh)
Inventor
丁伟洋
宋世光
张赛军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN202010865152.3A priority Critical patent/CN112100820A/en
Publication of CN112100820A publication Critical patent/CN112100820A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Control Of Metal Rolling (AREA)

Abstract

The invention belongs to the field of metal plastic forming, and provides a method for solving common yield function parameters of metal materials based on Excel, which comprises the following steps: for anisotropic materials, the yield function to be used is selected, being any of Barlat1989, Barlat 1991, Hill 1948, and BBC 2005; selecting a calibration method according to experimental data, wherein the calibration method is one of r value calibration, sigma value calibration and optimization calibration; obtaining each unknown parameter value in the selected yield function through the calibration method; finally, a complete expression of the selected yield function is obtained, and the finally obtained yield function is combined with a fracture criterion and a hardening model, so that a constitutive model of the material can be obtained through coupling to describe the mechanical property of the material. By the method, a proper yield function can be selected, and an unknown parameter value in the yield function is obtained by a calibration method, so that the complete expression of the yield function can be obtained, and the mechanical property can be conveniently and accurately described.

Description

Excel-based anisotropic yield function parameter calibration method
Technical Field
The invention belongs to the field of metal plastic forming, and particularly relates to a common anisotropic yield function parameter calibration method based on Excel.
Background
Stamping is an important process in modern manufacturing, and most of the processes are rolling metal plates, which have significant anisotropy. In order to be able to describe the anisotropy of the slab, the scholars have proposed various anisotropic yield functions. Under certain deformation conditions (deformation temperature, deformation speed, etc.), the mass point starts to enter the plastic state only when certain relations are met among the stress components, which is the mechanical condition that must be observed to describe that the stressed object enters the plastic state and the plastic deformation continues, and the mechanical condition can be generally expressed as a function of time, and the function is called as a yield function. The yield function is the basis of finite element simulation of metal plastic forming, and only if the mechanical property of the material in the forming process is accurately described, the qualified product can be obtained by improving the forming process in actual production. However, the accuracy of parameter calibration in the anisotropic yield function has a great influence on the prediction capability, and different parameter calibration methods have different relationships of the obtained anisotropic yield function, so that the yield function predicts the anisotropy of the material and obviously differs, and the finally described mechanical properties of the material are different.
In 1989 Barlat indicated that since Hosford's yield criterion does not contain a shear stress component, it is impossible to cope with the case where the principal axis of anisotropy does not coincide with the principal axis of stress, and proposed a yield criterion that takes into account in-plane anisotropy under planar stress conditions, in particular in the form of
Figure BDA0002649494780000011
In the formula:
Figure BDA0002649494780000012
Figure BDA0002649494780000013
wherein m is a non-quadratic yield function index; x, y and z are directions parallel to the rolling direction, perpendicular to the rolling direction and perpendicular to the plane of the plate, respectively; a. h and p are material parameters for characterizing anisotropy. Barlat1989 yield function, a and c are not independent, and m is related to the crystal structure. Therefore, four parameters of a, c, h and p need to be calibrated.
There are two calculation methods for solving the parameters a, c, h, p. One is obtained by a stress calculation method, and the assumption is that the sigma is90、τs1And τs2Yield stress when single-drawn in a direction 90 DEG to the rolling direction, shear stress is σ22＝-σ11＝τs2When, σ 220; when sigma is11＝σ22When equal to 0, σ12＝τs1At this time, it can be
Figure BDA0002649494780000021
Figure BDA0002649494780000022
Figure BDA0002649494780000023
Another calculation method is based on the thickness anisotropy index r0、r45、r90Is calculated to obtain
Figure BDA0002649494780000024
Figure BDA0002649494780000025
The p value cannot be resolved, but when a, h are known, forUniaxial stretching, rφIs univalent to p, and can be obtained by iterative method
Figure BDA0002649494780000026
In the formula, σ45Is the yield strength at 45 ° single draw from the rolling direction; for body centered cubic materials, m is 6, and for face centered cubic materials, m is 8. When m is 2, it is the yield criterion by Hill in 1948.
It can be seen from the above solving process that when solving the Barlart 1998 yield function calibration parameters, the manual calculation and drawing workload is large, a large amount of repetitive work needs to be done, the precision is low, and when the experimental data is incomplete, the parameter calibration cannot be performed. Some scholars calibrate the yield function by using MATLAB, but software programming is not easy and the operation is difficult; in addition, the software running configuration condition is higher, and the software copyright problem is solved, so that the parameter calibration is difficult and the efficiency is not high.
Disclosure of Invention
The invention aims to provide a method for solving common anisotropic yield function parameters based on Excel, which utilizes the function calculation function of Excel to automatically solve each unknown parameter in the yield function and provides three calibration methods: the yield function capable of well describing the mechanical property of the material can be obtained by selecting a proper yield function and a proper calibration method, the mechanical property of the material can be accurately described by combining a fracture criterion and a hardening model, and the forming process can be conveniently improved in production so as to obtain a qualified product.
In order to achieve the above object, the present invention provides a method for solving a common anisotropic yield function parameter based on Excel, comprising the following steps:
for anisotropic materials, selecting the yield function to be used;
selecting a calibration method according to experimental data, wherein the calibration method is one of r value calibration, sigma value calibration and optimization calibration;
obtaining each unknown parameter value in the selected yield function through the calibration method;
and finally, obtaining a complete expression of the selected yield function, and coupling the finally obtained yield function with a fracture criterion and a hardening model to obtain a constitutive model of the material so as to describe the mechanical property of the material.
Further, the yield function is any one of Barlat1989, Barlat 1991, Hill 1948, and BBC 2005.
Further, when the yield function selected is Barlat1989, the calibration method selected is any one of r-value calibration, σ -value calibration, and optimization calibration; when the selected yield function is Barlat 1991, the selected calibration method is any one of sigma value calibration and optimization calibration; when the selected yield function is Hill 1948, the selected calibration method is any one of r value calibration, sigma value calibration and optimization calibration; when the yield function is chosen to be BBC 2005, the calibration method chosen is optimized calibration.
Further, the selecting a calibration method according to the experimental data specifically includes: if r is obtained from the obtained experimental data0、r45、 r90And a fixed value m, then selecting the value r for calibration; if the obtained experimental data obtain sigma0、σ45、σ90、σbAnd a fixed value m, then selecting a value sigma for calibration; and if the obtained experimental data can not carry out r value calibration or sigma value calibration, selecting optimized calibration.
Further, when the yield function is Barlat1989, the unknown parameter values to be obtained are a, c, h and p; when the yield function is Barlat 1991, the unknown parameter values to be obtained are a, b, c, and h; when the yield function is Hill 1948, the unknown parameter values to be obtained are F, G, H, N, L and M, and when the yield function is BBC 2005, the unknown parameter values to be obtained are a, b, L, M, N, P, Q, and R.
Further, when the calibration method is r value calibration, the unknown parameter value is calculated as follows:
when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
Figure BDA0002649494780000031
Figure BDA0002649494780000032
a+c＝2；
Figure BDA0002649494780000041
wherein
Figure BDA0002649494780000042
C＝ma(A+B)|A+B|m-2+ma(A-B)|A-B|m-2， D＝ma(A+B)|A+B|m-2-ma(A-B)|A-B|m-2+4mcB|2B|m-2；
When the yield function is Hill 1948, the calculation formula for F, G, H, N, L and M is as follows:
Figure BDA0002649494780000043
Figure BDA0002649494780000044
Figure BDA0002649494780000045
Figure BDA0002649494780000046
the parameter L ═ M ═ 3, which indicates in-plane isotropy of the material;
wherein m is a constant value and is related to the crystal structure, r0、r45、r90Is along the rolling direction of the plate and in a direction of 45 degrees with the rolling directionThe anisotropy coefficient of the metal material perpendicular to the rolling direction can be obtained by experiment.
Further, when the calibration method is sigma value calibration, the unknown parameter value is calculated as follows:
when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
Figure BDA0002649494780000047
Figure BDA0002649494780000048
Figure BDA0002649494780000049
Figure BDA0002649494780000051
wherein sigmabIs equal biaxial tensile stress, σeIs a selected reference stress, typically taken as σ0Or σb。
When the yield function is Barlat 1991, the calculation formulas of a, b, c, and h are as follows:
Figure BDA0002649494780000052
Figure BDA0002649494780000053
Figure BDA0002649494780000054
Figure BDA0002649494780000055
when the yield function is Hill 1948, the unknown parameters F, G, H, N, L and M are calculated as follows:
Figure BDA0002649494780000056
Figure BDA0002649494780000057
Figure BDA0002649494780000058
Figure BDA0002649494780000059
the parameter L ═ M ═ 3, which indicates in-plane isotropy of the material;
in the formula, σ0、σ45、σ90Is the stress value, sigma, of the metal material along the rolling direction of the sheet, at 45 degrees to the rolling direction and perpendicular to the rolling directionbThe stress values measured for equal biaxial tension can be obtained experimentally.
Further, when the calibration method is an optimized calibration, the unknown parameter value is calculated as follows:
when the yield function is Barlat1989, the a, c, h and p values and the unknown experimental data are calculated as follows:
Figure BDA0002649494780000061
Figure BDA0002649494780000062
Figure BDA0002649494780000063
using the above 8 equations as constraints f (x) 0, a, c, h, p and unknown experimental data as unknowns, the objective function is MIN f1 2+f2 2+f3 2+f4 2+f5 2+f6 2+f7 2+f8 2(ii) a Solving a, c, h and p values and unknown experimental data by planning and solving under the condition of the minimum value of the objective function;
when the yield function is Barlat 1991, the calculation formulas of a, b, c and h are as follows:
the theoretical prediction formula of the thickness anisotropy coefficients of different oriented uniaxial tension is as follows
Figure BDA0002649494780000064
Wherein
Figure BDA0002649494780000065
Figure BDA0002649494780000066
E＝(2B|2B|m-2+(A-B)|A-B|m-2+2(A+B)|A+B|m-2) F＝(-2B|2B|m-2+2(A-B)|A-B|m-2+(A+B)|A+B|m-2) (ii) a When theta is 0 degrees, 45 degrees and 90 degrees, three values respectively corresponding to r can be obtained0、r45、r90The equation of (c); in addition, the
Figure BDA0002649494780000067
Wherein
Figure BDA0002649494780000068
Figure BDA0002649494780000069
To obtain rbCorresponding constraint, four formula structures calibrated with Barlat 1991 sigma valueAnd (4) forming 8 constraints, and planning and solving to obtain a, c, h and p values and unknown experimental data.
When the yield function is Hill 1948, the unknown parameters F, G, H, N, L and M are calculated as follows:
Figure BDA00026494947800000610
Figure BDA0002649494780000071
Figure BDA0002649494780000072
using the above 8 equations as constraints f (x) 0, F, G, H, N and unknown experimental data as unknowns, the objective function is MIN f1 2+f2 2+f3 2+f4 2+f5 2+f6 2+f7 2+f8 2(ii) a Solving F, G, H, N values and unknown experimental data by planning under the condition of the minimum value of the objective function;
when the yield function is BBC 2005, the unknown parameters a, b, L, M, N, P, Q, and R are calculated as follows:
for 0 degree direction uniaxial tension σ0The corresponding equation is
Figure BDA0002649494780000073
Wherein0＝L，Λ0＝|N|，Ψ0＝|Q|；
For uniaxial tension a in the 45 degree direction45The corresponding equation is
Figure BDA0002649494780000074
Wherein
Figure BDA0002649494780000075
For theUnidirectional stretching sigma in 90 degree direction90The corresponding equation is
Figure BDA0002649494780000076
Wherein90＝M，Λ90＝|P|，Ψ90＝|R|；
Equal biaxial tension sigmabThe corresponding equation is
Figure BDA0002649494780000077
Whereinb＝L+M，Λb＝|N-P|，Ψb＝|Q-R|；
For 0 degree direction uniaxial tension r0The corresponding equation is
Figure BDA0002649494780000078
Wherein C is10＝a(Λ0+0)2k+a(Λ0-0)2k+b(Λ0+Ψ0)2k+b(Λ0-Ψ0)2k,C20＝a(Λ0+0)m-1-a(Λ0-0)m-1, C30＝b(Λ0+Ψ0)m-1-b(Λ0-Ψ0)m-1，C40＝a(Λ0+0)m-1+a(Λ0-0)m-1+b(Λ0+Ψ0)m-1+b(Λ0-Ψ0)m-1，
Figure BDA0002649494780000079
For the unidirectional stretching r in the 45-degree direction45The corresponding equation is
Figure BDA00026494947800000710
Wherein C is145＝a(Λ45+45)2k+a(Λ45-45)2k+b(Λ45+Ψ45)2k+b(Λ45-Ψ45)2k,C245＝a(Λ45+45)m-1-a(Λ45-45)m-1 C345＝b(Λ45+Ψ45)m-1-b(Λ45-Ψ45)m-1，C445＝a(Λ45+45)m-1+a(Λ45-45)m-1+b(Λ45+Ψ45)m-1+b(Λ45-Ψ45)m-1
Figure BDA0002649494780000081
For the 90 degree direction uniaxial tension r90The corresponding equation is
Figure BDA0002649494780000082
Wherein C is190＝a(Λ90+90)2k+a(Λ90-90)2k+b(Λ90+Ψ90)2k+b(Λ90-Ψ90)2k,C290＝a(Λ90+90)m-1-a(Λ90-90)m-1 C390＝b(Λ90+Ψ90)m-1-b(Λ90-Ψ90)m-1，C490＝a(Λ90+90)m-1+a(Λ90-90)m-1+b(Λ90+Ψ90)m-1+b(Λ90-Ψ90)m-1
Figure BDA0002649494780000083
Equal biaxial tension rbThe corresponding equation is
Figure BDA0002649494780000084
Wherein C is2b＝a(Λb+b)m-1-a(Λb-b)m-1 C3b＝b(Λb+Ψb)m-1-b(Λb-Ψb)m-1，C4b＝a(Λb+b)m-1+a(Λb-b)m-1+b(Λb+Ψb)m-1+b(Λb-Ψb)m-1
Figure BDA0002649494780000085
The above 8 formulas are arranged into a form of f (x) equal to 0 as a constraint, a, b, L, M, N, P, Q and R are unknowns, and the objective function is
MIN＝f1 2+f2 2+f3 2+f4 2+f5 2+f6 2+f7 2+f8 2
And (4) solving the values of a, b, L, M, N, P, Q and R parameters by planning and solving under the condition of the minimum value of the objective function.
Furthermore, the rupture criterion, the hardening model and the yield function are coupled into a constitutive model of the material by a program language by using a Fortran language.
The yield function is combined with the fracture criterion and the hardening model, a constitutive model of the material can be obtained through coupling, and the mechanical property of the material can be predicted through the constitutive model. And the accurate expression of the yield function has great influence on the accuracy of the mechanical property of the material predicted by the constitutive model. The invention provides a plurality of calibration methods, and provides a method for solving unknown parameters in different yield functions under different calibration methods so as to obtain complete expression of the corresponding yield functions. During prediction, a yield function with the best prediction effect can be selected from the yield function, a proper calibration method is selected according to data obtained through experiments, so that an expression of the yield function with the best prediction effect is obtained, a constitutive model with the best prediction capability can be obtained by combining a fracture criterion and a hardening model, and the mechanical property of the material can be predicted more accurately, so that the forming process can be improved in actual production, and qualified products can be obtained.
Drawings
Fig. 1 is a schematic flow chart of the present embodiment.
FIG. 2 is a graph of r θ of the method for calibrating the r value of the yield function of Barlat1989 in accordance with an embodiment of the present invention.
FIG. 3 is a σ θ diagram of the method for calibrating the r value of the Barlat1989 yield function according to an embodiment of the invention.
FIG. 4 is a yield surface diagram of the calibration method for the r value of the Barlat1989 yield function according to the embodiment of the present invention.
FIG. 5 is a graph of r θ of the yield function σ value calibration method of Barlat1989 in accordance with an embodiment of the present invention.
FIG. 6 is a σ θ diagram of the method for calibrating the σ value of the yield function in Barlat1989 in accordance with an embodiment of the present invention.
FIG. 7 is a yield surface diagram of the method for calibrating the sigma value of the yield function in Barlat1989 according to an embodiment of the present invention.
FIG. 8 is a graph of r θ of the calibration method of the yield function OPT solution of Barlat1989 in accordance with an embodiment of the present invention.
FIG. 9 is a σ θ diagram of the calibration method of the yield function OPT solution of Barlat1989 in accordance with an embodiment of the present invention.
FIG. 10 is a yield surface diagram of the calibration method of the Barlat1989 yield function OPT solution according to the embodiment of the present invention.
FIG. 11 is a predicted load displacement curve for a 0 orientation of the yield function of Hill 1948, according to an embodiment of the invention.
FIG. 12 is a plot of predicted load displacement for a 45 orientation of the Hill 1948 yield function for an embodiment of the invention.
FIG. 13 is a predicted 90 orientation load displacement curve for the yield function of Hill 1948, according to an embodiment of the present invention.
FIG. 14 is a diagram of 5 fracture specimens in example of the present invention.
FIG. 15 is a graphical comparison of the results of the prediction and experiment performed on various samples using Hill 1948 in this example.
Detailed Description
The following detailed description of embodiments of the invention refers to the accompanying drawings. It should be understood that the detailed description and specific examples, while indicating embodiments of the invention, are given by way of illustration and explanation only, not limitation.
The method for calibrating the parameters of the common anisotropic yield function based on Excel provided by the embodiment comprises the following steps:
step 1: for anisotropic materials, the yield function to be used is selected, which is any of Barlat1989, Barlat 1991, Hill 1948, and BBC 2005. The yield function of the embodiment is suitable for anisotropic metal materials, when a specific yield function is selected for a specific material, if parameter calibration is performed on an existing material (a material which has been tested by others), a skilled person in the art can select a yield function capable of better and more accurately describing the mechanical property of the material according to experience, and for a new anisotropic material, various yield functions are generally adopted to perform parameter calibration respectively, and finally a yield function capable of more accurately describing the mechanical property of the material is selected according to a parameter calibration result.
(1) The expression for the yield function Barlat1989 is as follows:
Figure BDA0002649494780000091
Figure BDA0002649494780000092
Figure BDA0002649494780000093
a＝2-C
where m is a constant value, related to the crystal structure, obtained experimentally, σ11、σ12、σ22Is an unknown term in the yield function expression, and a, c, h and p are constant terms in the function expression, and the specific values of a, c, h and p can be obtained by the calibration method provided by the embodiment11And σ2The yield function Barlat1989 expression. SigmaeAll represent reference stress, and can be selected from0Or σb。
(2) The expression of the yield function Barlat 1991 is as follows:
Figure BDA0002649494780000101
Figure BDA0002649494780000102
wherein S isiIs the principal stress value of the bias stress tensor S, i is 1, 2, 3;
under in-plane stress conditions, the bias stress tensor S is:
Figure BDA0002649494780000103
wherein a, b, c and h are material parameters and unknown constant terms in the yield function Barlat 1991, and are obtained by the calibration method provided by the embodiment, and after the specific values of a, b, c and h are obtained by the calibration method, the specific values about sigma can be obtained11、σ22The yield function Barlat 1991 expression.
(3) The expression of the yield function Hill 1948 is as follows:
Figure BDA0002649494780000104
wherein σeFor reference stress, σ can be chosen0Or σbF, G, H, N, L and M are unknown constant terms in the functional expression, which is obtained by the calibration method provided by the embodiment, σ11、σ12、σ22、σ22、σ31And σ33Is an unknown term in the functional expression, and after obtaining F, G, H, N, L and M specific values by calibration method, the values for sigma can be obtained11、σ22And σ33The yield function Barlat1989 expression.
(4) The expression of the yield function BBC 2005 is as follows:
Figure BDA0002649494780000111
＝Lσ11+Mσ22
Figure BDA0002649494780000112
Figure BDA0002649494780000113
wherein a, b, L, M, N, P, Q and R are unknown constant terms in the function expression, and are obtained by the calibration method provided by the embodiment, and σ is11、σ12And σ22Is an unknown term in the function expression, and after the specific values of a, b, L, M, N, P, Q and R are obtained by a calibration method, the unknown term sigma can be obtained11And σ22Yield function BBC 2005.
Step 2: and selecting a calibration method according to the experimental data, wherein the calibration method is one of r value calibration, sigma value calibration and optimization calibration (also called an OPT solution).
When the selected yield function is Barlat1989, the selected calibration method is any one of r value calibration, sigma value calibration and optimization calibration; when the selected yield function is Barlat 1991, the selected calibration method is any one of sigma value calibration and optimization calibration; when the selected yield function is Hill 1948, the selected calibration method is any one of r value calibration, sigma value calibration and optimization calibration; when the yield function is chosen to be BBC 2005, the calibration method chosen is optimized calibration.
When the experimental data obtained by the experiment has r0、r45、r90When m is fixed, r is selected for calibration; when the experimental data obtained by the experiment has σ0、σ45、σ90、σbWhen m is fixed, the sigma value is selected for calibration; and when the experimental data obtained through the experiment can not meet the r value calibration or sigma value calibration, adopting optimized calibration.
And step 3: and obtaining each unknown parameter value in the selected yield function through a calibration method.
(1) When the calibration method is used for calibrating the r value, the unknown parameter value of each yield function is calculated in the following mode:
(1.1) when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
Figure BDA0002649494780000114
Figure BDA0002649494780000115
a+c＝2；
Figure BDA0002649494780000116
wherein
Figure BDA0002649494780000121
C＝ma(A+B)|A+B|m-2+ma(A-B)|A-B|m-2， D＝ma(A+B)|A+B|m-2-ma(A-B)|A-B|m-2+4mcB|2B|m-2
(1.2) when the yield function is Hill 1948, the calculation formula for F, G, H, N, L and M is as follows:
Figure BDA0002649494780000122
Figure BDA0002649494780000123
Figure BDA0002649494780000124
Figure BDA0002649494780000125
the parameter L ═ M ═ 3, which indicates in-plane isotropy of the material;
wherein m is a constant value and is related to the crystal structure, r0、r45、r90The anisotropy coefficient of the metal material along the rolling direction of the plate, 45 degrees to the rolling direction and vertical to the rolling direction can be obtained through experiments.
(2) When the calibration method is used for calibrating the sigma value, the unknown parameter value of each yield function is calculated in the following mode:
(2.1) when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
Figure BDA0002649494780000126
Figure BDA0002649494780000127
Figure BDA0002649494780000128
Figure BDA0002649494780000129
(2.2) when the yield function is Barlat 1991, the calculation formulas of a, b, c, and h are as follows:
Figure BDA00026494947800001210
Figure BDA0002649494780000131
Figure BDA0002649494780000132
Figure BDA0002649494780000133
(2.3) when the yield function is Hill 1948, the unknown parameters F, G, H, N, L and M are calculated as follows:
Figure BDA0002649494780000134
Figure BDA0002649494780000135
Figure BDA0002649494780000136
Figure BDA0002649494780000137
the parameter L ═ M ═ 3, which indicates in-plane isotropy of the material; in the formula, σ0、σ45、σ90The stress values are the stress values of the metal material along the rolling direction of the plate, in the direction of 45 degrees with the rolling direction and perpendicular to the rolling direction, and the sigma b is the stress value measured by equal biaxial tension and can be obtained through experiments.
(3) When the calibration method is optimized calibration, the unknown parameter values of each yield function are calculated in the following mode:
(3.1) when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
the 8 solving formulas calibrated by the Barlat1989 r value and the sigma value are arranged into the form of f (x) -0 as constraint, namely
Figure BDA0002649494780000138
Figure BDA0002649494780000139
Figure BDA00026494947800001310
Using the above 8 equations as constraints f (x) 0, a, c, h, p and unknown experimental data as unknowns, the objective function is MIN f1 2+f2 2+f3 2+f4 2+f5 2+f6 2+f7 2+f8 2(ii) a Solving a, c, h and p values and unknown experimental data by planning and solving under the condition of the minimum value of the objective function;
(3.2) when the yield function is Barlat 1991, the calculation formulas of a, c, h and p are as follows:
the theoretical prediction formula of the thickness anisotropy coefficients of different oriented uniaxial tension is as follows
Figure BDA0002649494780000141
Wherein
Figure BDA0002649494780000142
Figure BDA0002649494780000143
E＝(2B|2B|m-2+(A-B)|A-B|m-2+2(A+B)|A+B|m-2) F＝(-2B|2B|m-2+2(A-B)|A-B|m-2+(A+B)|A+B|m-2) (ii) a When theta is 0 degrees, 45 degrees and 90 degrees, three values respectively corresponding to r can be obtained0、r45、r90The equation of (c); in addition, the
Figure BDA0002649494780000144
Wherein
Figure BDA0002649494780000145
Figure BDA0002649494780000146
To obtain rbCorresponding constraints and four formulas calibrated with the Barlat 1991 sigma value form 8 constraints, and the values of a, c, h and p and unknown experimental data are obtained by planning and solving.
(3.3) when the yield function is Hill 1948, the unknown parameters F, G, H, N, L and M are calculated as follows:
the 8 solving formulas calibrated by Hill 1948r value and sigma value are arranged into a form of f (x) 0 as constraint, namely
Figure BDA0002649494780000147
Figure BDA0002649494780000148
Figure BDA0002649494780000149
Using the above 8 equations as constraints f (x) 0, F, G, H, N and unknown experimental data as unknowns, the objective function is MIN f1 2+f2 2+f3 2+f4 2+f5 2+f6 2+f7 2+f8 2(ii) a Solving F, G, H, N values and unknown experimental data by planning under the condition of the minimum value of the objective function;
(3.4) when the yield function is BBC 2005, the unknown parameters a, b, L, M, N, P, Q and R are calculated as follows: for 0 degree direction uniaxial tension σ0The corresponding equation is
Figure BDA0002649494780000151
Wherein0＝L，Λ0＝|N|，Ψ0＝|Q|；
For uniaxial tension a in the 45 degree direction45The corresponding equation is
Figure BDA0002649494780000152
Wherein
Figure BDA0002649494780000153
For a uniaxial extension of 90 degrees90The corresponding equation is
Figure BDA0002649494780000154
Wherein90＝M，Λ90＝|P|，Ψ90＝|R|；
Equal biaxial tension sigmabThe corresponding equation is
Figure BDA0002649494780000155
Whereinb＝L+M，Λb＝|N-P|，Ψb＝|Q-R|；
For 0 degree direction uniaxial tension r0The corresponding equation is
Figure BDA0002649494780000156
Wherein C is10＝a(Λ0+0)2k+a(Λ0-0)2k+b(Λ0+Ψ0)2k+b(Λ0-Ψ0)2k,C20＝a(Λ0+0)m-1-a(Λ0-0)m-1, C30＝b(Λ0+Ψ0)m-1-b(Λ0-Ψ0)m-1，C40＝a(Λ0+0)m-1+a(Λ0-0)m-1+b(Λ0+Ψ0)m-1+b(Λ0-Ψ0)m-1，
Figure BDA0002649494780000157
For the unidirectional stretching r in the 45-degree direction45The corresponding equation is
Figure BDA0002649494780000158
Wherein C is145＝a(Λ45+45)2k+a(Λ45-45)2k+b(Λ45+Ψ45)2k+b(Λ45-Ψ45)2k,C245＝a(Λ45+45)m-1-a(Λ45-45)m-1 C345＝b(Λ45+Ψ45)m-1-b(Λ45-Ψ45)m-1，C445＝a(Λ45+45)m-1+a(Λ45-45)m-1+b(Λ45+Ψ45)m-1+b(Λ45-Ψ45)m-1
Figure BDA0002649494780000159
For the 90 degree direction uniaxial tension r90The corresponding equation is
Figure BDA00026494947800001510
Wherein C is190＝a(Λ90+90)2k+a(Λ90-90)2k+b(Λ90+Ψ90)2k+b(Λ90-Ψ90)2k,C290＝a(Λ90+90)m-1-a(Λ90-90)m-1 C390＝b(Λ90+Ψ90)m-1-b(Λ90-Ψ90)m-1，C490＝a(Λ90+90)m-1+a(Λ90-90)m-1+b(Λ90+Ψ90)m-1+b(Λ90-Ψ90)m-1
Figure BDA0002649494780000161
Equal biaxial tension rbThe corresponding equation is
Figure BDA0002649494780000162
Wherein C is2b＝a(Λb+b)m-1-a(Λb-b)m-1 C3b＝b(Λb+Ψb)m-1-b(Λb-Ψb)m-1，C4b＝a(Λb+b)m-1+a(Λb-b)m-1+b(Λb+Ψb)m-1+b(Λb-Ψb)m-1
Figure BDA0002649494780000163
The above 8 formulas are arranged into a form of f (x) equal to 0 as a constraint, a, b, L, M, N, P, Q and R are unknowns, and the objective function is
MIN＝f1 2+f2 2+f3 2+f4 2+f5 2+f6 2+f7 2+f8 2
And (4) solving the values of a, b, L, M, N, P, Q and R parameters by planning and solving under the condition of the minimum value of the objective function.
And 4, step 4: finally, a complete expression of the selected yield function is obtained, and the finally obtained yield function is combined with a fracture criterion and a hardening model, so that a constitutive model of the material can be obtained through coupling to describe the mechanical property of the material.
In this embodiment, the yield function, the fracture criterion and the hardening model are coupled to obtain the constitutive model of the material through the VUAMT subprogram of abaqus, that is, the fracture criterion, the hardening model and the yield function are written together in a program language by using a Fortran language and coupled to form the constitutive model of the material. Breaking criteria and hardening models are well known in the art and will not be described in detail herein. The constitutive model of the material obtained by the Fortran language coupling is a common means in the prior art and will not be described in detail here.
The principle of the method for solving the parameters of the common anisotropic yield function by using Excel is described below.
Establishing an Excel file named as a common calibration method of anisotropic yield function parameters based on Excel, and respectively editing formulas of relevant calibration methods in yield functions of Barlat1989, Barlat 1991, Hill 1948 and BBC 2005 in 4 workbooks of Excel.
Experimental data for materials of examples, material name: A6181-T4; experimental data of r value and sigma value of material: r0 ═ 1.557, r45 ═ 1.342, r90 ═ 1.964, rb ═ 1.3, σ 0 ═ 160, σ 45 ═ 171, σ 90 ═ 167, σ b ═ 184; the relevant index of the Barlat series yield function is usually m-8, the relevant index of the BBC series yield function is usually k-4, and the value of the material for the Barlat series m is 8. In this example, the parameter is solved by using Excel, taking Barlat1989 yield function as an example. The EXCEL used in the present embodiment may be the version 2010 or above.
The method comprises the steps of firstly, opening an Excel file, selecting a 'Barlat 1989' workbook based on a common anisotropic yield function parameter calibration method of Excel, editing three calibration methods in the workbook, obviously editing a sigma value calibration formula and an optimized calibration formula in the 'Barlat 1991' workbook, editing the three calibration method formulas in the 'Hill 1948' workbook, and editing the optimized calibration formula in the 'BBC 2005' workbook. Of course, in other embodiments, a calibration method may be established by a workbook, according to the user's usage habits or preferences.
Secondly, as shown in fig. 2, selecting an r-value calibration method, inputting an m value of 8 into a corresponding cell C3, inputting r values of 1.557, 1.342, 1.964 into cells C5, D5, E5, and simultaneously specifying a reference stress value, and specifying a reference stress value σ f ═ σ 0 ═ 160 or σ b value into a cell G5, wherein the reference stress is specified and plays a role in normalization; at the moment, Excel automatically solves the values of the parameters a, c and h;
TABLE 1 values of a, c, h and p
Figure BDA0002649494780000171
There are two ways to solve the p value, one is: inputting a calculation formula of a parameter p in a target cell, finding a 'simulation analysis' in a pull-down menu of a 'data' column of Excel software, namely 'single-variable solving', and generating a 'single-variable solving' dialog box after clicking; setting target cell, inputting 0 in the 'target value', inputting variable cell, and finally solving p value according to determination. Of course, for the convenience of calculation, the calculation formula of the parameter P may be linked to the macro button "parameter solution 1" by using the VBA program, and in the actual calculation, the macro button is directly clicked to obtain the value P, that is, the parameter P is obtained through the following third step.
And thirdly, clicking a ' parameter solving 1 ' button, solving the implicit function of the single variable by using a ' simulation analysis ' -single variable solving ' function in Excel, and simultaneously realizing one-key solving of the parameter p by using a VBA program, as shown in the table.
Step four, after parameter solving is completed in the step two and the step three, Excel automatically calculates the r theta and sigma theta values corresponding to all angles; clicking a macro button to solve the yield surface 1, and solving the horizontal and vertical coordinates of the yield surface in batch by using a VBA program;
wherein, the calculation formula of r theta and sigma theta is as follows:
Figure BDA0002649494780000172
Figure BDA0002649494780000173
wherein
Figure BDA0002649494780000174
Figure BDA0002649494780000175
C＝ma(A+B)|A+B|m-2+ma(A-B)|A-B|m-2 D＝ma(A+B)|A+B|m-2-ma(A-B)|A-B|m-2+4mcB|2B|m-2
Wherein, solving the horizontal and vertical coordinates of the yielding surface, the function of which is an implicit function, specifically: a (ABS (sigma 1)) ^ m + a (ABS (h sigma 2)) ^ m + c (ABS (sigma 1-h sigma 2)) ^ m-2 sigma f Lambda m is 0, sigma 1 is regarded as a fixed value (since the yield function is an implicit function, direct drawing is not easy here, therefore, the initial value of sigma 1 is designated as-272, the interval is divided by 3.2, the interval can be smaller), sigma 2 is regarded as a variable, and the sigma 2 value corresponding to the sigma 1 value, namely the horizontal and vertical coordinates, is solved by using 'single variable solution'; and the fast solving is convenient, and the batch solving is carried out by using the for cycle in the VBA program. As shown in table 2, a series of values of r θ and σ θ are obtained by solving, the abscissa (σ 1) and the ordinate (σ 2) of the yield surface are finally obtained by solving, then θ (°) is used as the abscissa, and r θ is used as the ordinate, a r θ curve is drawn, as shown in fig. 2, θ (°) is used as the abscissa, and σ θ is used as the ordinate, a σ θ curve is drawn, as shown in fig. 2, the value of r is a thickness anisotropy coefficient, r θ represents the thickness anisotropy coefficient of the material in different rolling directions, and reflects the difference of strain capacity in the plate plane direction and the thickness direction, σ θ represents stress values in different directions, and actual values (square points in the graph) obtained by experiments are given on each graph, as can be seen from the graph, the yield function Barlat1989 and r value calibration are adopted to perform parameter calibration on the material a 1-T4, and the theoretical values and the actual values are basically approximate in the finally obtained r θ graph, the theoretical value and the actual value in the σ θ diagram are greatly different. A yield surface graph is drawn by taking σ 1 as an abscissa and σ 2 as an ordinate, as shown in fig. 4, the yield surface is a basis for judging that the material is plastically deformed, if a stress state of a point in the main stress space is located on the yield surface, the point is in a plastic state (is plastically deformed), and if the stress state is located on the yield surface, the point is in an elastic state.
TABLE 2 values of r θ, σ 1, and σ 2
θ(°) rθ σθ σ1 σ2 Equation of surface of yield
0 1.557 160 -272 -5.38804E+13 1.01E+110
1 1.556667 160.0096254 -268.8 -3.1895E+11 1.528E+92
2 1.555669 160.0384715 -265.6 999949.8082 1.428E+48
3 1.554013 160.0864479 -262.4 -142.1874674 1.571E+19
4 1.551706 160.1534043 -259.2 -709.1273554 3.453E+22
5 1.548763 160.2391301 -256 -54753550832 1.152E+86
6 1.545199 160.3433547 -252.8 -140.3574602 1.144E+19
7 1.541035 160.4657478 -249.6 -159.6246602 1.036E+19
。。。 。。。 。。。 。。。 。。。 。。。
For convenience of drawing, the abscissa and the ordinate can be respectively connected to the VBA program and respectively connected to the macro buttons of 'drawing r theta diagram', 'drawing sigma theta diagram' and 'drawing yield surface', and when drawing is needed, the corresponding macro button is directly clicked, so that the corresponding diagram can be drawn by the VBA program. And the actual data points are plotted on the respective graphs accordingly.
Fifthly, selecting a sigma value calibration method, inputting an m value of 8 in a corresponding cell, and respectively inputting sigma 0, sigma 45, sigma 90 and sigma b values in the corresponding cells: 160. 171, 167, 184, and at the same time, it is necessary to specify the value of the reference stress, where the reference stress value σ f ═ σ 0 ═ 160 or σ b may be specified in the cell U5, and the normalization function is performed by specifying the reference stress value; at this time, Excel automatically finds out that the parameter h is 0.958083832.
And sixthly, the functions f1, f2 and f3 are equal to 0, and the values of the parameters a, c and p can be obtained through solving.
Of course, for convenience of calculation, in the embodiment, the function formulas f1, f2 and f3 are linked to the macro button "plan solution" through the VBA program, and when calculation is needed, the "plan solution" is found in the "data" column of the Excel software, or the "plan solution" is loaded by the "loading macro"; clicking 'planning solution' to generate a 'planning solution parameter' dialog box; setting the target cell as one of three functions, and selecting the target value as 0 in the 'to' options; inputting variable cells, inputting constraint condition formulas with f1, f2 and f3 equal to zero one by one after 'obeying constraint' single click 'addition', wherein a 'selection solution method' is nonlinear GRG and is finally determined; after clicking the 'solution', a 'planning solution result' dialog box appears, the solution result is prompted, clicking 'determination' is carried out, the information closest to the planning solution result is filled into the target cell and the variable cell in the worksheet, and the storage of the planning solution result is completed; the a, c, p values are output in the cells.
And seventhly, obtaining a series of values of r theta and sigma theta in the same step as the fourth step. A series of values of sigma 1 and sigma 2 are obtained, as shown in the following table, an r theta diagram, a sigma theta diagram and a yield surface diagram are correspondingly drawn, as shown in FIGS. 5-6, when a yield function Barlat1989 and sigma value calibration are adopted to carry out parameter calibration on a material A6181-T4, in the finally obtained r theta diagram, the difference between a theoretical value and an actual value is large, and the theoretical value and the actual value in the sigma theta diagram are basically coincident, namely, the sigma theta diagram can better embody the mechanical property of the material in the method.
TABLE 3 values of r θ, σ 1, and σ 2 in the σ value calibration
θ(°) rθ σθ σ1 σ2 Equation of surface of yield
0 3.44288541 160 -272 -1.2E+15 5.043E+120
1 3.44242201 160.0102 -268.8 -1.7E+10 1.14662E+82
2 3.44103475 160.0407 -265.6 -2.1E+12 6.2172E+98
3 3.43873249 160.0914 -262.4 -7.9E+14 2.143E+119
4 3.43552995 160.1622 -259.2 -6.4E+13 3.7766E+110
5 3.43144768 160.2529 -256 -171.22 6.40294E+18
6 3.42651201 160.3632 -252.8 -187.058 5.92724E+18
7 3.42075492 160.4926 -249.6 -1.6E+15 4.8883E+121
8 3.41421402 160.641 -246.4 -4.3E+14 1.7115E+117
9 3.40693238 160.8076 -243.2 -141.521 3.89298E+18
10 3.39895841 160.9921 -240 -144.661 3.41991E+18
11 3.39034571 161.1939 -236.8 -146.772 2.99153E+18
。。。 。。。 。。。 。。。 。。。 。。。
Eighthly, selecting an OPT calibration method, inputting a value m of 8 into a corresponding cell, inputting data r0, r45, r90, 1.964, rb, 1.3, sigma 0, sigma 45, sigma 171, sigma 90, sigma 167 and sigma b 184 which are obtained in an experiment into the cell, and if some of the data are not obtained, taking the data which are not obtained and parameters a, c, h and p as unknown quantities; in the embodiment, the data is complete, the data is filled into corresponding cells, and reference stress is specified;
and ninthly, enabling all constraint functions f 1-f 8 to be equal to 0, and solving under the condition that the objective function MIN reaches the minimum value to obtain specific values of a, c, h and p.
Certainly, for convenience of calculation, in the embodiment, the function formulas f1 to f8 are linked to the macro button "plan solution 1" through the VBA program, and when calculation is needed, the "plan solution" is found in a pull-down menu of an Excel software "tool" column, or the "plan solution" is loaded by the "loading macro". Clicking the 'planning solution' and generating a 'planning solution parameter' dialog box. Setting a target cell: selecting a minimum value from a 'to' option, inputting four variables of a, c, h and p in a 'variable cell', inputting constraint condition formulas when each f 1-f 8 is equal to 0 after 'single click' addition 'in' constraint ', checking' enabling an unconstrained variable to be a non-negative number (K) ',' selecting a solving method 'to be' nonlinear GRG ', and finally' solving 'according to the formula'; after clicking 'solving', a 'planning solving result' dialog box appears, the solving result is prompted, clicking 'determining' is carried out, the objective function solving result closest to the planning can be filled into the objective cells in the worksheet, and the storage of the planning solving result is completed.
The parameters a, c, h and p are solved by the yield function Barlat1989 and by three different calibration methods.
And step ten, like the step four, obtaining a series of values of r theta and sigma theta. A series of values of sigma 1 and sigma 2 are obtained, as shown in the following table, an r theta diagram, a sigma theta diagram and a yield surface diagram are correspondingly drawn, as shown in FIGS. 8-9, a yield function Barlat1989 and an optimized calibration are adopted to carry out parameter calibration on the material A6181-T4, in the finally obtained r theta diagram, a theoretical value is basically close to an actual value, the difference between the theoretical value and the actual value in the sigma theta diagram is larger, and the mechanical property of the material can be better reflected by the r theta diagram in the method.
Similarly, by adopting the similar steps, other yield functions, namely Barlat 1991 or Hill 1948 or BBC 2005, are selected, the calibration method is correspondingly selected, the specific value of the unknown parameter of the corresponding yield function can be calculated and obtained through the calculation formula of the selected calibration method, and the r theta diagram, the sigma theta diagram and the yield surface diagram in the corresponding calibration method can be drawn to represent the mechanical characteristics of the material.
Taking four samples of samples with arc notch radiuses of 5mm and 20mm, samples with center holes and shear samples as examples, a load displacement curve of each sample is predicted by using Hill 1948, Hill 1948_ sigma (a yield function is Hill 1948, and a calibration method is sigma value calibration) is applied to a correlation flow criterion, and yield stress and an r value of each orientation are determined by Hill48_ sigma; hill48_ r was also applied to the correlation flow criteria, the yield stress and r-value for each orientation were determined by Hill 1948_ r (yield function is Hill 1948, calibration method is r-value calibration), and the actual load-displacement curve was experimentally determined, see FIGS. 11-12. By comparing the predicted load-displacement curve with the experimentally measured load-displacement curve, it is possible to obtain:
(1) on the prediction of the load-displacement curves of the notched specimens with 0 °, 45 ° and 90 ° orientations and the specimens with central holes, the yield criterion of Hill48 — σ applied to the correlation flow criterion was accurate and quite close to the experimental value.
(2) The predicted result of the yield stress of Hill 48-sigma is lower than the experimental value on the load displacement curve prediction of shear samples with 0 degrees, 45 degrees and 90 degrees, because the yield stress and the r value of Hill48 are determined by Hill 48-sigma when the criterion is applied to the correlation flow criterion; hill48 σ is only calibrated by yield stress, which does not accurately predict the direction of flow of the material, whereas shear specimens are more sensitive to the direction of flow of the material due to their particular stress state.
Therefore, the method is applied to the Hill48 yield function under the relevant flow criterion, and can accurately predict the mechanical property of the material in the yield stage, thereby laying a foundation for finite element simulation of material forming.
The description of the yield function, in combination with the constitutive model obtained by coupling the fracture criterion and the hardening model, on the mechanical properties of the material is verified through experimental results.
Taking Hill 1948 as an example, the influence of the shape of the yield surface on the fracture prediction accuracy was verified.
Taking materials A6181-T4 as an example, the following data are obtained by experiments:
TABLE 4 Experimental values of aluminum alloy materials
Figure BDA0002649494780000211
The following parameter values were obtained by r-value calibration (hereinafter referred to as Hill48_ r) and σ -value calibration (hereinafter referred to as Hill48_ s):
values of the respective parameters obtained in Table 4
Figure BDA0002649494780000212
And (3) coupling the DF2012 toughness fracture criterion, the Swift-Voce hardening model and the Hill 1948 yield criterion which is respectively calibrated by the r value and the sigma value into different material constitutive models through a user subprogram (VUMAT) so as to research the influence of different yield surface shapes on the fracture prediction accuracy. The fracture prediction results for the two yield surface shapes for the 5 fracture samples are shown in fig. 14. The 5 fracture samples are samples with arc notch radii of 5mm, 10mm and 20mm, a center hole sample and a shear sample, respectively, as shown in fig. 14. As can be seen from fig. 14, for the circular arc notch samples R5 and R10 and the center hole sample, the predicted fracture position of Hill48_ s is closer to the experimental result, and the fracture time is earlier than the experiment; for the circular arc notch sample R20, the predicted fracture position of Hill48_ R is closer to the experimental result, the fracture time is earlier than that of the experiment, and the predicted fracture time of Hill48_ s is later than that of the experiment; for the shear samples, the predicted results of the two yield surface shapes are greatly different from the experimental values, the predicted maximum load is obviously higher than the experimental values, and as can be seen from fig. 12, the predicted shear stress of Hill _ s in the shear region is higher than that of Hill _ r, so that the predicted load of Hill _ s is also higher. Therefore, unknown parameter values in the relevant yield functions can be obtained through different calibration methods, so that expressions of the yield functions under different calibration methods can be obtained, the obtained yield functions are coupled with the fracture criterion and the hardening model, so that constitutive models of materials under different calibration methods can be obtained, the mechanical properties of the materials under different constitutive models can be respectively researched, and the yield function with the most accurate description can be selected from the constitutive models. Therefore, in the subsequent actual production, the calibration method with the most accurate description can be adopted to obtain the corresponding yield function, and the most accurate constitutive model is obtained through coupling, so that the mechanical property of the material in the forming process can be described most accurately, and the qualified product can be obtained through improving the forming process.
The r value calibration and the sigma value calibration can utilize a simulation analysis function-single variable solving function in Excel to solve calibration parameters according to a formula of the selected yield function, and an optimized calibration method is adopted when the obtained experimental data is not enough to carry out r value calibration or sigma value calibration. And the Excel VBA can be used for programming during actual calculation, so that the operation is performed on a visual interface, and the method is simple to operate, accurate in calculation result, high in practicability, automatic in drawing and attractive in graph. The method can be stored as an application program and can be used repeatedly, the solved parameter value can be obtained as long as the numerical value obtained by inputting the test data into the corresponding cell, the complex calculation and drawing processing is avoided, a large amount of time is saved, the working efficiency is improved, the practicability is high, and the method can be widely applied to test teaching and scientific research.
It should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Other variations and modifications will be apparent to persons skilled in the art in light of the above description. And are neither required nor exhaustive of all embodiments. 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 claims of the present invention.

Claims (9)

1. A common anisotropic yield function parameter calibration method based on Excel is characterized by comprising the following steps:
for anisotropic materials, a common yield function is selected;
selecting a calibration method according to experimental data, wherein the calibration method is one of r value calibration, sigma value calibration and optimization calibration;
obtaining each unknown parameter value in the selected yield function through the selected calibration method;
and finally, obtaining a complete expression of the selected yield function, and coupling the finally obtained yield function with a fracture criterion and a hardening model to obtain a constitutive model of the material so as to describe the mechanical property of the material.
2. The Excel-based calibration method for common anisotropic yield function parameters according to claim 1, characterized in that: the yield function is any of Barlat1989, Barlat 1991, Hill 1948, and BBC 2005.
3. The Excel-based calibration method for common anisotropic yield function parameters according to claim 2, characterized in that: when the selected yield function is Barlat1989, the selected calibration method is any one of r value calibration, sigma value calibration and optimization calibration; when the selected yield function is Barlat 1991, the selected calibration method is any one of sigma value calibration and optimization calibration; when the selected yield function is Hill 1948, the selected calibration method is any one of r value calibration, sigma value calibration and optimization calibration; when the yield function is chosen to be BBC 2005, the calibration method chosen is optimized calibration.
4. The Excel-based calibration method for common anisotropic yield function parameters according to claim 1, characterized in that the calibration method is selected according to experimental data, specifically: if r is obtained from the obtained experimental data0、r45、r90And a fixed value m, then selecting the value r for calibration; if the obtained experimental data obtain sigma0、σ45、σ90、σbAnd a fixed value m, then selecting a value sigma for calibration; and if the obtained experimental data are not enough to carry out r value calibration or sigma value calibration, selecting the optimized calibration.
5. The Excel-based calibration method for common anisotropic yield function parameters according to any one of claims 1 to 4, characterized in that: when the yield function is Barlat1989, the unknown parameter values to be obtained are a, c, h and p; when the yield function is Barlat 1991, the unknown parameter values to be obtained are a, b, c, and h; when the yield function is Hill 1948, the unknown parameter values to be obtained are F, G, H, N, L and M, and when the yield function is BBC 2005, the unknown parameter values to be obtained are a, b, L, M, N, P, Q, and R.
6. The Excel-based calibration method for common anisotropic yield function parameters according to claim 5, wherein when the calibration method is r value calibration, the unknown parameter values are calculated as follows:
when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
Figure FDA0002649494770000011
Figure FDA0002649494770000012
a+c＝2；
Figure FDA0002649494770000021
wherein the content of the first and second substances,
Figure FDA0002649494770000022
C＝ma(A+B)|A+B|m-2+ma(A-B)|A-B|m-2，D＝ma(A+B)|A+B|m-2-ma(A-B)|A-B|m-2+4mcB|2B|m-2；
when the yield function is Hill 1948, the calculation formula for F, G, H, N, L and M is as follows:
Figure FDA0002649494770000023
Figure FDA0002649494770000024
Figure FDA0002649494770000025
Figure FDA0002649494770000026
the parameter L ═ M ═ 3, which indicates in-plane isotropy of the material;
where m is constant and is related to the crystal structure, r0、r45、r90Is along the rolling direction of the plate and at an angle of 45 degrees with the rolling directionAnd the anisotropy coefficient of the metal material perpendicular to the rolling direction can be obtained by experiment.
7. The Excel-based calibration method for common anisotropic yield function parameters according to claim 5, wherein when the calibration method is calibration of σ value, the unknown parameter value is calculated as follows:
when the yield function is Barlat1989, the calculation formulas for a, c, h, and p are as follows:
Figure FDA0002649494770000027
Figure FDA0002649494770000028
Figure FDA0002649494770000029
Figure FDA0002649494770000031
wherein σbIs equal biaxial tensile stress, σeIs a selected reference stress, typically taken as σ0Or σb；
When the yield function is Barlat 1991, the values of a, b, c, and h can be solved in combination with the following four equations:
Figure FDA0002649494770000032
Figure FDA0002649494770000033
Figure FDA0002649494770000034
Figure FDA0002649494770000035
when the yield function is Hill 1948, the unknown parameters F, G, H, N, L and M are calculated as follows:
Figure FDA0002649494770000036
Figure FDA0002649494770000037
Figure FDA0002649494770000038
Figure FDA0002649494770000039
the parameter L ═ M ═ 3, which indicates in-plane isotropy of the material;
in the formula, σ0、σ45、σ90The stress values of the metal material along the plate rolling direction, the direction at 45 degrees to the rolling direction and the direction perpendicular to the rolling direction are shown, and the sigma b is the stress value measured by equal biaxial tension and can be obtained through experiments.
8. The Excel-based calibration method for common anisotropic yield function parameters according to claim 5, wherein when the calibration method is optimized calibration, the unknown parameter values are calculated as follows:
when the yield function is Barlat1989, the a, c, h and p values and the unknown experimental data are calculated as follows:
Figure FDA0002649494770000041
Figure FDA0002649494770000042
f4＝a+c-2；
Figure FDA0002649494770000043
Figure FDA0002649494770000044
Figure FDA0002649494770000045
Figure FDA0002649494770000046
using the above 8 formulas as constraints f (x) 0, a, c, h, p and unknown experimental data as unknowns, the objective function is
Figure FDA0002649494770000047
Solving a, c, h and p values and unknown experimental data by planning and solving under the condition of the minimum value of the objective function;
when the yield function is Barlat 1991, the calculation methods of a, b, c, h and unknown experimental data are as follows:
the theoretical prediction formula of the thickness anisotropy coefficient of the unidirectional stretching with different orientations is
Figure FDA0002649494770000048
Wherein
Figure FDA0002649494770000049
Figure FDA00026494947700000410
E＝(2B|2B|m-2+(A-B)|A-B|m-2+2(A+B)|A+B|m-2)，F＝(-2B|2B|m-2+2(A-B)|A-B|m-2+(A+B)|A+B|m-2)；
When theta is 0 degrees, 45 degrees and 90 degrees, three values respectively corresponding to r can be obtained0、r45、r90The equation of (c);
in addition, the
Figure FDA00026494947700000411
Wherein
Figure FDA00026494947700000412
Figure FDA00026494947700000413
Figure FDA00026494947700000414
Thereby r can be obtainedbCorresponding constraints and 8 constraints formed by four formulas calibrated by the sigma value of Barlat 1991 are planned and solved to obtain a, b, c and h values and unknown experimental data;
when the yield function is Hill 1948, the unknown parameters F, G, H, N, L, M and the unknown experimental data are calculated as follows:
Figure FDA0002649494770000051
Figure FDA0002649494770000052
Figure FDA0002649494770000053
using the above 8 equations as constraints f (x) 0, F, G, H, N and unknown experimental data as unknowns, the objective function is
Figure FDA0002649494770000054
Solving F, G, H, N values and unknown experimental data by planning under the condition of the minimum value of the objective function;
when the yield function is BBC 2005, the unknown parameters a, b, L, M, N, P, Q, and R are calculated as follows:
for 0 degree direction uniaxial tension σ0The corresponding equation is
Figure FDA0002649494770000055
Wherein0＝L，Λ0＝|N|，Ψ0＝|Q|；
For uniaxial tension a in the 45 degree direction45The corresponding equation is
Figure FDA0002649494770000056
Wherein the content of the first and second substances,
Figure FDA0002649494770000057
for a uniaxial extension of 90 degrees90The corresponding equation is
Figure FDA0002649494770000058
Wherein the content of the first and second substances,90＝M，Λ90＝|P|，Ψ90＝|R|；
equal biaxial tension sigmabThe corresponding equation is
Figure FDA0002649494770000061
Wherein the content of the first and second substances,b＝L+M，Λb＝|N-P|，Ψb＝|Q-R|；
for 0 degree direction uniaxial tension r0The corresponding equation is
Figure FDA0002649494770000062
Wherein C is10＝a(Λ0+0)2k+a(Λ0-0)2k+b(Λ0+Ψ0)2k+b(Λ0-Ψ0)2k，C20＝a(Λ0+0)m-1-a(Λ0-0)m-1,C30＝b(Λ0+Ψ0)m-1-b(Λ0-Ψ0)m-1，C40＝a(Λ0+0)m-1+a(Λ0-0)m-1+b(Λ0+Ψ0)m-1+b(Λ0-Ψ0)m-1，
Figure FDA0002649494770000063
For the unidirectional stretching r in the 45-degree direction45The corresponding equation is
Figure FDA0002649494770000064
Wherein C is145＝a(Λ45+45)2k+a(Λ45-45)2k+b(Λ45+Ψ45)2k+b(Λ45-Ψ45)2k，C245＝a(Λ45+45)m-1-a(Λ45-45)m-1，C345＝b(Λ45+Ψ45)m-1-b(Λ45-Ψ45)m-1，C445＝a(Λ45+45)m-1+a(Λ45-45)m-1+b(Λ45+Ψ45)m-1+b(Λ45-Ψ45)m-1，
Figure FDA0002649494770000065
For the 90 degree direction uniaxial tension r90The corresponding equation is
Figure FDA0002649494770000066
Wherein C is190＝a(Λ90+90)2k+a(Λ90-90)2k+b(Λ90+Ψ90)2k+b(Λ90-Ψ90)2k，C290＝a(Λ90+90)m-1-a(Λ90-90)m-1，C390＝b(Λ90+Ψ90)m-1-b(Λ90-Ψ90)m-1，C490＝a(Λ90+90)m-1+a(Λ90-90)m-1+b(Λ90+Ψ90)m-1+b(Λ90-Ψ90)m-1，
Figure FDA0002649494770000067
Equal biaxial tension rbThe corresponding equation is
Figure FDA0002649494770000071
Wherein
C2b＝a(Λb+b)m-1-a(Λb-b)m-1，C3b＝b(Λb+Ψb)m-1-b(Λb-Ψb)m-1，C4b＝a(Λb+b)m-1+a(Λb-b)m-1+b(Λb+Ψb)m-1+b(Λb-Ψb)m-1，
Figure FDA0002649494770000072
The above 8 formulas are arranged into a form of f (x) equal to 0 as a constraint, a, b, L, M, N, P, Q and R are unknowns, and the objective function is
Figure FDA0002649494770000073
And (4) solving the values of a, b, L, M, N, P, Q and R parameters by planning and solving under the condition of the minimum value of the objective function.
9. The Excel-based calibration method for common anisotropic yield function parameters according to claim 1, characterized in that a Fortran language is used to couple a fracture criterion, a hardening model and a yield function into a material constitutive model.
CN202010865152.3A 2020-08-25 2020-08-25 Excel-based anisotropic yield function parameter calibration method Pending CN112100820A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010865152.3A CN112100820A (en) 2020-08-25 2020-08-25 Excel-based anisotropic yield function parameter calibration method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010865152.3A CN112100820A (en) 2020-08-25 2020-08-25 Excel-based anisotropic yield function parameter calibration method

Publications (1)

Publication Number Publication Date
CN112100820A true CN112100820A (en) 2020-12-18

Family

ID=73753338

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010865152.3A Pending CN112100820A (en) 2020-08-25 2020-08-25 Excel-based anisotropic yield function parameter calibration method

Country Status (1)

Country Link
CN (1) CN112100820A (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010061249A (en) * 2008-09-02 2010-03-18 Jfe Steel Corp Method for calculating material characteristic parameter
JP2014079767A (en) * 2012-10-15 2014-05-08 Jfe Steel Corp Press forming analysis method, press forming analysis program, and storage medium
CN108256211A (en) * 2018-01-08 2018-07-06 大连理工大学 Timber constitutive relation method for numerical simulation based on ABAQUS
CN110763566A (en) * 2019-11-28 2020-02-07 大连理工大学 Method for determining circumferential thickness anisotropy coefficient of anisotropic pipe

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010061249A (en) * 2008-09-02 2010-03-18 Jfe Steel Corp Method for calculating material characteristic parameter
JP2014079767A (en) * 2012-10-15 2014-05-08 Jfe Steel Corp Press forming analysis method, press forming analysis program, and storage medium
CN108256211A (en) * 2018-01-08 2018-07-06 大连理工大学 Timber constitutive relation method for numerical simulation based on ABAQUS
CN110763566A (en) * 2019-11-28 2020-02-07 大连理工大学 Method for determining circumferential thickness anisotropy coefficient of anisotropic pipe

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
BINGTAO TANG ET AL: "Effect of Anisotropic Yield Functions on the Accuracy of Material Flow and its Experimental Verification", ACTA MECHANICA SOLIDA SINICA, vol. 32, 9 August 2018 (2018-08-09), pages 50 - 68 *
DOREL BANABIC ET AL: "INFLUENCE OF CONSTITUTIVE EQUATIONS ON THE ACCURACY OF PREDICTION IN SHEET METAL FORMING SIMULATION", NUMISHEET, 5 September 2008 (2008-09-05), pages 1 - 6 *
KAUSHIK BANDYOPADHYAY ET AL: "obust multi objective optimization of anisotropic yield function coefficients", MATERIALS AND DESIGN, vol. 156, 19 June 2018 (2018-06-19), pages 184 - 197, XP085440277, DOI: 10.1016/j.matdes.2018.06.033 *
倪向贵, 吴恒安, 王宇, 王秀喜: "各向异性本构关系在板料成形数值模拟中的应用", 计算力学学报, vol. 20, no. 02, 30 April 2003 (2003-04-30), pages 231 - 235 *
展学鹏: "基于球形压痕的材料塑性各向异性性能识别研究", 中国优秀硕士学位论文全文数据库工程科技Ⅰ辑, 15 June 2020 (2020-06-15), pages 022 - 7 *
李健强;张赛军;龚小龙;黄珍媛;袁宁;: "基于优化方法的复杂各向异性屈服函数参数标定", 塑性工程学报, vol. 24, no. 01, 28 February 2007 (2007-02-28), pages 160 - 167 *

Similar Documents

Publication Publication Date Title
Wen et al. New model for ductile fracture of metal alloys. I: Monotonic loading
Shlyannikov et al. Characterization of crack tip stress fields in test specimens using mode mixity parameters
Dæhli et al. A Lode-dependent Gurson model motivated by unit cell analyses
CN103886125B (en) A kind of titanium alloy hot combined shaping method for numerical simulation
CN103745114B (en) Method for computing stress relaxation numerical values and resilience of titanium alloy
Nesládek et al. An Abaqus plugin for fatigue predictions
CN114462124A (en) Method for establishing and numerically simulating concrete three-dimensional multiphase mesoscopic model
Gusel et al. Application of genetic programming for modelling of material characteristics
WO2019062727A1 (en) Method and device for calculating equivalent mechanical parameters of film layer etching region
Voyiadjis et al. Effects of stress invariants and reverse loading on ductile fracture initiation
Mu et al. Anisotropic hardening and evolution of r-values for sheet metal based on evolving non-associated Hill48 model
Hamzehkolaei et al. Reliability-based design optimization of rotating FGM cylindrical shells with temperature-dependent probabilistic frequency constraints
Wang et al. An investigation on the anisotropic plastic behavior and forming limits of an Al-Mg-Li alloy sheet
Smith Photoelasticity in fracture mechanics: Use of photoelasticity to obtain stress-intensity factors in three-dimensional cracked-body problems
CN112100820A (en) Excel-based anisotropic yield function parameter calibration method
Safikhani et al. The strain gradient approach for determination of forming limit stress and strain diagrams
Zhang et al. Prediction of ductile fracture for DP590 high-strength steel with a new semi-coupled ductile fracture criterion
Xie et al. Asymmetric yield effect evolving with internal variables during continuous large deformations and its FEM validation
CN106682328B (en) Vibration deformation measurement calculation method for vertical high-rise structure vibration isolation system
Chatziioannou et al. An implicit numerical scheme for cyclic elastoplasticity and ratcheting under plane stress conditions
CN114974478A (en) Crystal metal material right-angle micro-cutting modeling method and system considering strain rate
Zhang et al. Constitutive modeling based on non-associated flow rule for anisotropic sheet metals forming
Mu et al. Applicability of Hill48 Yield Model and Effect of Anisotropic Parameter Determination Methods on Anisotropic Prediction
JP2002530197A (en) A Modeling Method for Anisotropic Metal Sheet Forming
CN116050225B (en) Method and device for determining nonuniform pressure in plate forming process

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination