CN112861305A - Crack propagation direction prediction method and device and storage medium - Google Patents
Crack propagation direction prediction method and device and storage medium Download PDFInfo
- Publication number
- CN112861305A CN112861305A CN202011523316.0A CN202011523316A CN112861305A CN 112861305 A CN112861305 A CN 112861305A CN 202011523316 A CN202011523316 A CN 202011523316A CN 112861305 A CN112861305 A CN 112861305A
- Authority
- CN
- China
- Prior art keywords
- time
- displacement
- strain
- stress
- initial
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 93
- 238000006073 displacement reaction Methods 0.000 claims abstract description 92
- 238000004364 calculation method Methods 0.000 claims abstract description 31
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 31
- 239000000463 material Substances 0.000 claims abstract description 31
- 239000011159 matrix material Substances 0.000 claims description 82
- 230000001133 acceleration Effects 0.000 claims description 37
- 230000006870 function Effects 0.000 claims description 26
- 230000010354 integration Effects 0.000 claims description 14
- 238000004590 computer program Methods 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 10
- 230000002441 reversible effect Effects 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000000670 limiting effect Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
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
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
The application relates to a method, a device and a storage medium for predicting crack propagation direction, wherein the prediction method comprises the following steps: calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method; calculating a stress-strain field according to the calculation result of the displacement field; and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion. On the basis of a gridless method, the method provides the strain-rotation and decomposition theorem, forms a new crack propagation direction prediction method, can improve the accuracy of numerical analysis of the nonlinear deformation fracture problem, and provides new ideas and means for the numerical analysis of the fracture problem.
Description
Technical Field
The application relates to the technical field of fracture problem numerical analysis, in particular to a crack propagation direction prediction method and device and a storage medium.
Background
The problem of fracturing has long been one of the most important and complex problems in various engineering fields. Fracture analysis becomes particularly important, especially with the rapid development of advanced materials and structures. For the fracture problem, the analytical solution can be found only in a very simple case, and the approximate calculation of the fracture problem can be usually realized only by a numerical method.
Finite element is one of the numerical calculation methods commonly used in engineering calculation, and the finite element based expansion finite element method (XFEM) is a very common method for dealing with the crack problem at present. However, the finite element method needs to be processed by means of grid cells, and has the problems of cell distortion, complex preprocessing process, time consumption, large calculation amount and the like, so that the problem of complex crack propagation caused by the computational efficiency and accuracy of the discontinuous problem (such as the crack propagation problem) is greatly reduced; in addition, because the expansion path and the original unit interface can not be matched, the finite element is difficult to realize accurate calculation, and if the large deformation of materials and structures is accompanied before and after the fracture, the distortion of the unit can seriously influence the solving precision of the finite element.
The emerging meshless-Free Methods in the last 30 years can partially or completely eliminate meshes and directly utilize discrete nodes to analyze problems, and the unique characteristic can overcome partial problems existing in finite elements and ensure the precision of results, so the meshless method has great advantages for processing crack propagation problems. However, all the fracture problem numerical methods based on the meshless method are based on the small deformation theory basis at present, and although the linear elastic fracture mechanics is concise and perfect, the linear elastic fracture mechanics is strictly not suitable for fracture cusp field analysis due to the limitation of the linear elastic theory in the aspect of describing larger deformation, so that the accuracy of the numerical analysis of the nonlinear deformation fracture problem is limited from the theory itself.
Disclosure of Invention
In view of the above technical problems, the present application provides a method, an apparatus, and a storage medium for predicting a crack propagation direction, so as to improve accuracy of numerical analysis of a fracture problem due to nonlinear deformation, and provide a new idea and means for numerical analysis of the fracture problem.
In order to solve the above technical problem, the present application provides a method for predicting a crack propagation direction, including: calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method; calculating a stress-strain field according to the calculation result of the displacement field; and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion.
In one embodiment, the strain-rotation and resolution theorem includes:
F=S+R
wherein F is a reversible linear differential transformation function; s, S,Is the strain tensor; r, R,Is the tensor of rotation; u. ofi|jIs a displacement component uiTo coordinate xjA derivative of (a);is ui|jThe transposed matrix of (2);is a kronecker function;is a unit vector of the rotating shaft; theta is the average angle of full rotation.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
obtaining an incremental variational equation according to the strain-rotation and decomposition theorem:
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,is a stress tensor,Is an elastic matrix; Δ t is the time increment;is the strain rate;is the rate of rotation;t+ΔtWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
according to the gridless method, a moving least square interpolation function is adopted to obtain the displacement increment and the speed at the time t:
wherein,tψk(x) A moving least squares function for node k at time t;tΔukandtVkrespectively the displacement increment and the speed of the node k at the time t;tΔ u is the displacement increment at the time t,tAnd V is the speed at the time t.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
substituting the displacement increment of the time t and the speed of the time t into the increment variational equation to obtain a system discrete equation:
t+ΔtKe tΔu=t+ΔtFe
wherein,t+ΔtKe、t+ΔtFerespectively an equivalent stiffness matrix and an equivalent force vector at the moment of t + delta t;tKL、tKN、tM、tf is a linear stiffness matrix, a nonlinear stiffness matrix, a mass matrix and a stress integral vector at the moment t respectively; kαFor the penalty function stiffness matrix, FαIs a penalty function force vector;t+Δtr is an external force vector at the moment of t + delta t;ta is the acceleration at the time t; the parameter β is determined according to the accuracy and stability requirements of the integration.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
according to the Newmark method, the relationship among the displacement increment, the speed and the acceleration at different moments is obtained:
wherein,t+Δta is the acceleration at time t + Deltat,tA is the acceleration at time t,t+ΔtV is the speed at time t + Deltat,tV is the speed at time t,tAnd delta u is the displacement increment at the moment t, delta t is the time increment, and the parameters beta and gamma are determined according to the requirements of the integral precision and stability.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
obtaining an initial equivalent stiffness matrix, an initial mass matrix, an initial equivalent force vector and an initial displacement, and calculating an initial acceleration according to the following formula:
A=M-1(Fe-Keu0)
wherein A is the initial acceleration, KeIs the initial equivalent stiffness matrix, M is the initial mass matrix, FeIs the initial equivalent force vector, u0Is the initial displacement.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
obtaining a correction increment variational equation according to the strain-rotation and decomposition theorem:
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,is a stress tensor,Is an elastic matrix; Δ t is the time increment;is the strain rate;a rate of rotation;t+Δ tWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;is a penalty factor; s is a boundary; vi、VjIs the speed;is the velocity determined by the boundary conditions.
In one embodiment, the step of calculating a stress-strain field from the calculation of the displacement field includes:
and acquiring a stress tensor according to the strain tensor by adopting the following formula:
In one embodiment, the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and the crack propagation criterion includes:
and obtaining the main stress and the main direction thereof by adopting the following formula according to the weighted stress of the crack tip:
Wherein σtipA weighted stress of the crack tip; sigma1、σ2Is the principal stress, n1,n2Is the primary direction.
In one embodiment, the step of calculating the crack propagation direction according to the calculation result of the stress field and the crack propagation criterion includes:
the weighted stress of the crack tip is calculated according to the following formula:
wherein σ1Selecting the stress of a node near the tip of the crack; n isqThe number of integration points near the crack tip; w is aIIs the weight of the integral point; dmaxIs the crack tip andmaximum value of distance between selected nodes in the vicinity thereof, dIFor the corresponding integration point xiA distance from the tip of the crack, and dIAnd dmaxSatisfy the relationship 0 < dI≤2dmax。
In one embodiment, before the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and the crack propagation criterion, the method comprises: judging whether the crack length reaches the maximum value; judging whether the position of the crack tip reaches the boundary of the problem domain of the material; stopping the numerical calculation process if the crack length reaches the maximum value and/or the crack tip position reaches the problem domain boundary; and if the crack length does not reach the maximum value and the crack tip position does not reach the problem domain boundary, executing the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining the crack propagation criterion.
In one embodiment, the prediction method includes the following steps: acquiring an initial characteristic matrix, initial displacement and initial speed, and obtaining initial acceleration according to the initial characteristic matrix and the initial displacement; selecting a time step length and calculating an integral constant; acquiring a stiffness characteristic matrix at the time t, and acquiring an equivalent stiffness matrix at the time t + delta t according to the stiffness characteristic matrix at the time t and the integral constant; acquiring an equivalent force characteristic matrix at the time t, the speed at the time t and the acceleration at the time t, and acquiring an equivalent force vector at the time t + delta t according to the equivalent force characteristic matrix at the time t, the speed at the time t, the acceleration at the time t and the integral constant; obtaining a displacement increment at the t moment according to the equivalent stiffness matrix at the t + delta t moment and the equivalent force vector at the t + delta t moment; acquiring the displacement at the time t, and acquiring the displacement at the time t + delta t according to the displacement at the time t and the displacement increment at the time t; obtaining the speed at the t + delta t moment and the acceleration at the t + delta t moment according to the displacement increment at the t moment, the speed at the t moment and the acceleration at the t moment; and (5) circulating the time step, wherein t is t + delta t, and returning to execute the step three.
The application also provides a device for predicting the crack propagation direction, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the steps of the prediction method when executing the computer program.
The present application also provides a computer-readable storage medium, in which a computer program is stored, which computer program, when being executed by a processor, realizes the steps of the above-mentioned prediction method.
The method provides a strain-rotation and decomposition theorem on the basis of a gridless method, forms a new crack propagation direction prediction method, can effectively improve the accuracy of numerical analysis of the nonlinear deformation fracture problem, and provides new ideas and means for the numerical analysis of the fracture problem.
Drawings
FIG. 1 is a schematic flow chart of a prediction method according to an embodiment of the present disclosure;
FIG. 2 is a schematic diagram of a tow coordinate system provided in an embodiment of the present application;
FIG. 3 is a schematic view of a crack tip and its corresponding node according to an embodiment of the present disclosure;
fig. 4 is a schematic flowchart of a prediction method according to a second embodiment of the present application;
fig. 5 is a schematic flowchart of a prediction method provided in the third embodiment of the present application;
fig. 6 is a schematic structural diagram of a prediction apparatus according to a fourth embodiment of the present application.
Detailed Description
The following description of the embodiments of the present application is provided for illustrative purposes, and other advantages and capabilities of the present application will become apparent to those skilled in the art from the present disclosure.
In the following description, reference is made to the accompanying drawings that describe several embodiments of the application. It is to be understood that other embodiments may be utilized and that mechanical, structural, electrical, and operational changes may be made without departing from the spirit and scope of the present application. The following detailed description is not to be taken in a limiting sense, and the scope of embodiments of the present application is defined only by the claims of the issued patent. The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the application.
Although the terms first, second, etc. may be used herein to describe various elements in some instances, these elements should not be limited by these terms. These terms are only used to distinguish one element from another.
Also, as used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context indicates otherwise. It will be further understood that the terms "comprises," "comprising," "includes" and/or "including," when used in this specification, specify the presence of stated features, steps, operations, elements, components, items, species, and/or groups, but do not preclude the presence, or addition of one or more other features, steps, operations, elements, components, species, and/or groups thereof. The terms "or" and/or "as used herein are to be construed as inclusive or meaning any one or any combination. Thus, "A, B or C" or "A, B and/or C" means "any of the following: a; b; c; a and B; a and C; b and C; A. b and C ". An exception to this definition will occur only when a combination of elements, functions, steps or operations are inherently mutually exclusive in some way.
Fig. 1 is a schematic flowchart of a prediction method according to an embodiment of the present application. As shown in fig. 1, the prediction method provided by the present application includes the following steps:
step S101: calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method;
in one embodiment, the strain-rotation and resolution theorem includes:
F=S+R (1)
wherein F is a reversible linear differential transformation function; s, S,Is the strain tensor; r, R,Is the tensor of rotation; u. ofi|jIs a displacement component uiTo coordinate xjA derivative of (a);is ui|jThe transposed matrix of (2);is a kronecker function; theta is the average angle of full rotation.
Taking a two-dimensional problem as an example, in order to describe the motion of any deformable body in space, two reference frames can be selected, as shown in fig. 2: (1) global fixed reference system { Xi(i ═ 1, 2); (2) natural towing coordinate system { x embedded in deformable bodyiAnd (i is 1,2), the coordinate line of which deforms along with the deformation of the object. The motion deformation process of the object is completely carried out by an initial towing coordinate systemAnd instantaneous tow coordinate system (g)i) Are defined and derived. The strain-rotation and decomposition theorem states that a unique direct and decomposition exists for any one reversible linear differential transformation function F. Wherein r, tr、t+ΔtrIs relative to the O pointA position vector of (a); u is t0The displacement to the time t, Δ u, is the displacement increment from the time t to the time t + Δ t.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
1) obtaining an incremental variational equation according to the strain-rotation and decomposition theorem:
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,is a stress tensor,Is an elastic matrix; Δ t is the time increment;is the strain rate;is the rate of rotation;t+ΔtWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
2) according to the gridless method, a moving least square interpolation shape function is adopted for approximation, and displacement increment and speed at the time t are obtained:
wherein,tψk(x) A moving least squares function for node k at time t;tΔukandtVkrespectively the displacement increment and the speed of the node k at the time t;tΔ u is the displacement increment at time t,tV is the velocity at time t.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
substituting equations (7) and (8) into equation (6) to obtain a system discrete equation:
t+ΔtKe tΔu=t+ΔtFe (9)
wherein,t+ΔtKe、t+ΔtFerespectively an equivalent stiffness matrix and an equivalent force vector at the moment of t + delta t;tKL、tKN、tM、tf is a linear stiffness matrix, a nonlinear stiffness matrix, a mass matrix and a stress integral vector at the moment t respectively; kαFor the penalty function stiffness matrix, FαIs a penalty function force vector;t+Δtr is an external force vector at the moment of t + delta t;ta is the acceleration at the time t; the parameter β is determined according to the accuracy and stability requirements of the integration.
It is worth mentioning that for equation (9), if the stiffness matrix term is not considered, the solution of the statics problem can be achieved.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
according to the Newmark method, the relationship among the displacement increment, the speed and the acceleration at different moments is obtained:
wherein,t+Δta is the acceleration at time t + Deltat,tA is the acceleration at time t,t+ΔtV is the speed at time t + Deltat,tV is the speed at time t,tAnd delta u is the displacement increment at the moment t, delta t is the time increment, and the parameters beta and gamma are determined according to the requirements of the integral precision and stability.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method comprises:
obtaining an initial equivalent stiffness matrix, an initial mass matrix, an initial equivalent force vector and an initial displacement, and calculating an initial acceleration according to the following formula:
A=M-1(Fe-Keu0) (13)
wherein A is initial acceleration, KeIs an initial equivalent stiffness matrix, M is an initial mass matrix, FeIs an initial equivalent force vector, u0Is the initial displacement.
Since the least squares approximation does not have the property of the kronecker function, the application of the intrinsic boundary condition can adopt a penalty function method or a lagrange multiplier method, and if the penalty function method is adopted, the formula (6) needs to be modified appropriately.
In one embodiment, the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and the gridless method further comprises:
obtaining a correction increment variational equation according to the strain-rotation and decomposition theorem:
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,is a stress tensor,Is an elastic matrix; Δ t is the time increment;is the strain rate;a rate of rotation;t+ΔtWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;is a penalty factor; s is a boundary; vi、VjIs the speed;is the velocity determined by the boundary conditions.
Step S102: calculating a stress-strain field according to the calculation result of the displacement field;
in one embodiment, the step of calculating a stress-strain field from the calculation of the displacement field includes:
and acquiring a stress tensor according to the strain tensor by adopting the following formula:
Step S103: and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion.
The prediction method provided by the application can adapt to different crack propagation criteria, and the maximum shear stress criteria are taken as an example for illustration. To obtain a smoother crack propagation path, the crack propagation direction is calculated using the weighted stresses of all integration points near the crack tip. As shown in fig. 3, the integration point corresponding to the black solid node is the integration point considered in the crack propagation direction calculation.
In one embodiment, the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and the crack propagation criterion includes:
and obtaining the main stress and the main direction thereof by adopting the following formula according to the weighted stress of the crack tip:
Wherein σtipWeighted stress for crack tip; sigma1、σ2Is principal stress, n1,n2Is the main direction.
In one embodiment, the weighted stress of the crack tip is calculated according to the following formula:
wherein σISelecting the stress of a node near the tip of the crack; n isqThe number of integration points near the crack tip; w is aIIs the weight of the integral point; dmaxMaximum value of the distance between the crack tip and the selected node in the vicinity thereof, dIFor the corresponding integration point xiA distance from the tip of the crack, and dIAnd dmaxSatisfy the relationship 0 < dI≤2dmax。
Fig. 4 is a schematic flowchart of a prediction method according to a second embodiment of the present application. As shown in fig. 4, the prediction method provided by the present application includes the following steps:
step S201: performing node dispersion on the problem domain of the material;
step S202: creating a background integral grid and integral points according to the node discrete domain;
step S203: performing visual processing and display on the crack model;
step S204: calculating a displacement field;
step S205: calculating a stress-strain field by using the calculation result of the displacement field;
step S206: judging whether the crack length reaches the maximum value;
if the crack length reaches the maximum value, the process proceeds to step S207: stopping the numerical calculation process; if the crack length does not reach the maximum value, the step S208 is executed; optionally, whether the crack length reaches the maximum value is judged by using a crack propagation criterion, and if the crack length reaches the maximum value, the crack is judgedtip≤σ0The crack length reaches a maximum value, if σtip>σ0If so, the crack length does not reach the maximum value; wherein σtipThe weighted stress of the crack tip is calculated according to the formula (17); sigma0The fracture critical shear stress of the material is measured by experiments;
step S208: judging whether the position of the crack tip reaches the boundary of the problem domain;
if the crack tip position reaches the problem domain boundary, the process proceeds to step S207: if the crack tip position does not reach the problem domain boundary, the step S209 is carried out; optionally, whether the crack tip position reaches the problem domain boundary is judged by using a coordinate system, and if the coordinate of the crack tip position coincides with the coordinate of any point on the problem domain boundary, the crack tip position reaches the problem domain boundary.
Step S209: calculating the crack propagation direction by using a crack propagation criterion based on the calculation result of the stress-strain field;
step S210: the position of the new crack tip is determined and the process returns to step S202.
With the loop from step S202 to step S210, the displacement field in step S204 is also calculated with loop iteration, which is specifically described with reference to the third embodiment; step S205 and step S209 can refer to the first embodiment, and are not described herein again.
Fig. 5 is a schematic specific flowchart of a prediction method according to a third embodiment of the present application. As shown in fig. 5, the prediction method provided by the present application includes the following steps:
step S301: acquiring an initial characteristic matrix, initial displacement and initial speed, and obtaining initial acceleration according to the initial characteristic matrix and the initial displacement;
wherein the initial characteristic matrix comprises an initial equivalent stiffness matrix KeInitial quality matrix M, initial equivalent force vector Fe(ii) a In particular, according to an initial equivalent stiffness matrix KeInitial quality matrix M, initial equivalent force vector FeInitial displacement u0Calculating the initial acceleration A by adopting a formula (13) in the first embodiment; the initial velocity belongs to the initial condition, and normally the structure is stationary at the initial moment, so the initial velocity V is typically zero.
Step S302: selecting a time step length and calculating an integral constant;
specifically, a time step delta t is selected, a parameter beta is combined, and an integral constant A is calculated1、A2、A3;
the parameter β is determined according to the accuracy and stability requirements of the integration.
Step S303: acquiring a stiffness characteristic matrix at the time t, and acquiring an equivalent stiffness matrix at the time t + delta t according to the stiffness characteristic matrix and an integral constant at the time t; acquiring an equivalent force characteristic matrix at the time t, the speed at the time t and the acceleration at the time t, and acquiring an equivalent force vector at the time t + delta t according to the equivalent force characteristic matrix at the time t, the speed at the time t, the acceleration at the time t and an integral constant;
specifically, according to the stiffness characteristic matrix and the integral constant at the time t, the formula (10) in the first embodiment is adopted to obtain an equivalent stiffness matrix at the time t + Δ t; and obtaining an equivalent force vector at the t + delta t moment by adopting the formula (11) in the first embodiment according to the equivalent force characteristic matrix at the t moment, the speed at the t moment, the acceleration at the t moment and the integral constant. Wherein the stiffness characteristic matrix at the time t comprises a linear stiffness matrix at the time ttKLTime t non-linear stiffness matrixtKNPenalty function stiffness matrix KαAnd the quality matrix at time ttM; the isodynamic characteristic matrix at the time t comprises an external force vector at the time t + delta tt+ΔtR, t stress integration vectortF. Penalty function force vector FαAnd the quality matrix at time ttM。
Step S304: obtaining a displacement increment at the t moment according to the equivalent stiffness matrix at the t + delta t moment and the equivalent force vector at the t + delta t moment; acquiring the displacement at the time t, and acquiring the displacement at the time t + delta t according to the displacement at the time t and the displacement increment at the time t; obtaining the speed at the t + delta t moment and the acceleration at the t + delta t moment according to the displacement increment at the t moment, the speed at the t moment and the acceleration at the t moment;
specifically, according to the equivalent stiffness matrix and the equivalent force vector at the time t + Δ t, the displacement increment at the time t is obtained by adopting the formula (9) in the first embodimenttΔ u; displacement and displacement increment according to t timetΔ u, resulting in time t + Δ tWherein the displacement at time t + Δ tt+ΔtDisplacement at time u-ttDisplacement increment at time u + ttΔ u; increment of displacement according to t timetΔ u, speedtV and accelerationtA, adopting the formula (12) in the first embodiment, the acceleration at the time of t + delta t is obtainedt+ΔtA. Speed of rotationt+ΔtV。
Step S305: and (4) a time step loop, wherein t is t + delta t, and the step S303 is executed in a returning mode.
Specifically, the displacement at time t + Δ t is updatedt+Δtu, accelerationt+ΔtA. Speed of rotationt+ΔtThe stiffness characteristic matrix and the equivalent effect characteristic matrix at the time V, t + delta t are obtained by executing the step S303 and the step S304 to obtain the node displacement at the time t +2 delta tt+2Δtu, accelerationt+2ΔtA. Speed of rotationt+2ΔtAnd V, repeating the steps to realize the cyclic iterative computation of the displacement field.
The prediction method provided by the application is based on a gridless method, combines the strain-rotation and decomposition theorem to form a new prediction method of the crack propagation direction, effectively improves the accuracy of numerical analysis of the nonlinear deformation fracture problem, and provides a new idea and means for the numerical analysis of the fracture problem.
Fig. 6 is a schematic structural diagram of a prediction apparatus according to a fourth embodiment of the present application. As shown in fig. 6, the prediction apparatus of this embodiment includes: a processor 110, a memory 111 and a computer program 112 stored in said memory 111 and executable on said processor 110. The processor 110, when executing the computer program 112, implements the steps in the various prediction method embodiments described above, such as the steps S101 to S103 shown in fig. 1.
The prediction device may include, but is not limited to, a processor 110, a memory 111. Those skilled in the art will appreciate that fig. 6 is only an example of a terminal and is not intended to be limiting and may include more or fewer components than those shown, or some components may be combined, or different components, e.g., the terminal may also include input-output devices, network access devices, buses, etc.
The Processor 110 may be a Central Processing Unit (CPU), other general purpose Processor, a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), an off-the-shelf Programmable Gate Array (FPGA) or other Programmable logic device, discrete Gate or transistor logic device, discrete hardware component, etc. A general purpose processor may be a microprocessor or the processor may be any conventional processor or the like.
The storage 111 may be an internal storage unit of the terminal, such as a hard disk or a memory of the terminal. The memory 111 may also be an external storage device of the terminal, such as a plug-in hard disk, a Smart Media Card (SMC), a Secure Digital (SD) Card, a Flash memory Card (Flash Card), and the like provided on the terminal. Further, the memory 111 may also include both an internal storage unit and an external storage device of the terminal. The memory 111 is used for storing the computer program and other programs and data required by the terminal. The memory 111 may also be used to temporarily store data that has been output or is to be output.
The present application also provides a computer-readable storage medium having stored thereon a computer program which, when executed by a processor, implements the steps of the prediction method as described above.
It should be understood that, the sequence numbers of the steps in the foregoing embodiments do not imply an execution sequence, and the execution sequence of each process should be determined by its function and inherent logic, and should not constitute any limitation to the implementation process of the embodiments of the present application.
The above embodiments are merely illustrative of the principles and utilities of the present application and are not intended to limit the application. Any person skilled in the art can modify or change the above-described embodiments without departing from the spirit and scope of the present application. Accordingly, it is intended that all equivalent modifications or changes which can be made by those skilled in the art without departing from the spirit and technical concepts disclosed in the present application shall be covered by the claims of the present application.
Claims (14)
1. A method for predicting a crack propagation direction, comprising:
calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method;
calculating a stress-strain field according to the calculation result of the displacement field;
and calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining a crack propagation criterion.
2. The prediction method of claim 1, wherein the strain-rotation and decomposition theorem comprises:
F=S+R
wherein F is a reversible linear differential transformation function; s, S,Is the strain tensor; r, R,Is the tensor of rotation; u. ofi|jIs a displacement component uiTo coordinate xjA derivative of (a);is ui|jThe transposed matrix of (2);is a kronecker function;is a unit vector of the rotating shaft; theta is the average angle of full rotation.
3. The prediction method of claim 2, wherein the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method comprises:
obtaining an incremental variational equation according to the strain-rotation and decomposition theorem:
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,is a stress tensor,Is an elastic matrix; Δ t is the time increment;is the strain rate;is the rate of rotation;t+Δ tWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;
according to the gridless method, a moving least square interpolation function is adopted to obtain the displacement increment and the speed at the time t:
wherein,tψk(x) A moving least squares function for node k at time t;tΔukandtVkrespectively the displacement increment and the speed of the node k at the time t;tΔ u is the displacement increment at the time t,tAnd V is the speed at the time t.
4. The prediction method according to claim 3, wherein the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method comprises:
substituting the displacement increment of the time t and the speed of the time t into the increment variational equation to obtain a system discrete equation:
t+ΔtKe tΔu=t+ΔtFe
wherein,t+ΔtKe、t+ΔtFeat times t + Δ t, respectivelyAn effective stiffness matrix and an equivalent force vector;tKL、tKN、tM、tf is a linear stiffness matrix, a nonlinear stiffness matrix, a mass matrix and a stress integral vector at the moment t respectively; kαFor the penalty function stiffness matrix, FαIs a penalty function force vector;t+Δtr is an external force vector at the moment of t + delta t;ta is the acceleration at the time t; the parameter β is determined according to the accuracy and stability requirements of the integration.
5. The prediction method of claim 1, wherein the step of calculating the displacement field of the material according to strain-rotation and decomposition theorem and a gridless method further comprises:
according to the Newmark method, the relationship among the displacement increment, the speed and the acceleration at different moments is obtained:
wherein,t+Δta is the acceleration at time t + Deltat,tA is the acceleration at time t,t+ΔtV is the speed at time t + Deltat,tV is the speed at time t,tAnd delta u is the displacement increment at the moment t, delta t is the time increment, and the parameters beta and gamma are determined according to the requirements of the integral precision and stability.
6. The prediction method according to any one of claims 4 to 5, wherein the step of calculating the displacement field of the material according to the strain-rotation and decomposition theorem and a gridless method comprises:
obtaining an initial equivalent stiffness matrix, an initial mass matrix, an initial equivalent force vector and an initial displacement, and calculating an initial acceleration according to the following formula:
A=M-1(Fe-Keu0)
wherein A is the initial acceleration, KeIs the initial equivalent stiffness matrix, M is the initial mass matrix, FeIs the initial equivalent force vector, u0Is the initial displacement.
7. The prediction method of claim 1, wherein the step of calculating the displacement field of the material according to strain-rotation and decomposition theorem and a gridless method further comprises:
obtaining a correction increment variational equation according to the strain-rotation and decomposition theorem:
wherein the prime symbol "-" represents the physical quantity under the initial tow at time t in the tow coordinate system,is a stress tensor,Is an elastic matrix; Δ t is the time increment;is the strain rate;a rate of rotation;t+Δ tWineandt+ΔtWextrespectively is inertia force power and external force virtual power; delta is a variation sign; Ω is the area of the material;is a penalty factor; s is a boundary; vi、VjIs the speed;is the velocity determined by the boundary conditions.
8. The prediction method according to claim 2, wherein the step of calculating a stress-strain field from the calculation result of the displacement field comprises:
and acquiring a stress tensor according to the strain tensor by adopting the following formula:
9. The prediction method according to claim 1, wherein the step of calculating the crack propagation direction based on the calculation result of the stress-strain field in combination with the crack propagation criterion comprises:
and obtaining the main stress and the main direction thereof by adopting the following formula according to the weighted stress of the crack tip:
Wherein σtipA weighted stress of the crack tip; sigma1、σ2Is the principal stress, n1,n2Is the primary direction.
10. The prediction method according to claim 9, wherein the step of calculating the crack propagation direction based on the calculation result of the stress-strain field in combination with the crack propagation criterion comprises:
the weighted stress of the crack tip is calculated according to the following formula:
wherein σ1Selecting the stress of a node near the tip of the crack; n isqThe number of integration points near the crack tip; w is aIIs the weight of the integral point; dmaxMaximum value of the distance between the crack tip and the selected node in the vicinity thereof, dIFor the corresponding integration point xiA distance from the tip of the crack, and dIAnd dmaxSatisfy the relationship 0 < dI≤2dmax。
11. The prediction method according to claim 1, wherein before the step of calculating a crack propagation direction in combination with a crack propagation criterion based on the calculation of the stress-strain field, comprises:
judging whether the crack length reaches the maximum value;
judging whether the position of the crack tip reaches the boundary of the problem domain of the material;
stopping the numerical calculation process if the crack length reaches the maximum value and/or the crack tip position reaches the problem domain boundary;
and if the crack length does not reach the maximum value and the crack tip position does not reach the problem domain boundary, executing the step of calculating the crack propagation direction according to the calculation result of the stress-strain field and by combining the crack propagation criterion.
12. The prediction method according to claim 1, characterized in that the prediction method comprises the steps of:
acquiring an initial characteristic matrix, initial displacement and initial speed, and obtaining initial acceleration according to the initial characteristic matrix and the initial displacement;
selecting a time step length and calculating an integral constant;
acquiring a stiffness characteristic matrix at the time t, and acquiring an equivalent stiffness matrix at the time t + delta t according to the stiffness characteristic matrix at the time t and the integral constant; acquiring an equivalent force characteristic matrix at the time t, the speed at the time t and the acceleration at the time t, and acquiring an equivalent force vector at the time t + delta t according to the equivalent force characteristic matrix at the time t, the speed at the time t, the acceleration at the time t and the integral constant;
obtaining a displacement increment at the t moment according to the equivalent stiffness matrix at the t + delta t moment and the equivalent force vector at the t + delta t moment; acquiring the displacement at the time t, and acquiring the displacement at the time t + delta t according to the displacement at the time t and the displacement increment at the time t; obtaining the speed at the t + delta t moment and the acceleration at the t + delta t moment according to the displacement increment at the t moment, the speed at the t moment and the acceleration at the t moment;
and (5) circulating the time step, wherein t is t + delta t, and returning to execute the step three.
13. A prediction apparatus of crack propagation direction comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that the processor implements the steps of the prediction method according to any of claims 1 to 12 when executing the computer program.
14. A computer-readable storage medium, in which a computer program is stored which, when being executed by a processor, carries out the steps of the prediction method according to any one of claims 1 to 12.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011523316.0A CN112861305A (en) | 2020-12-21 | 2020-12-21 | Crack propagation direction prediction method and device and storage medium |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011523316.0A CN112861305A (en) | 2020-12-21 | 2020-12-21 | Crack propagation direction prediction method and device and storage medium |
Publications (1)
Publication Number | Publication Date |
---|---|
CN112861305A true CN112861305A (en) | 2021-05-28 |
Family
ID=75997827
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011523316.0A Pending CN112861305A (en) | 2020-12-21 | 2020-12-21 | Crack propagation direction prediction method and device and storage medium |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112861305A (en) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113984515A (en) * | 2021-10-18 | 2022-01-28 | 同济大学 | Method for observing initiation and expansion of hydraulic fracture by using confocal microscope |
CN115906583A (en) * | 2022-12-16 | 2023-04-04 | 中国人民解放军陆军工程大学 | Method and system for simulating and analyzing structural integrity of explosive column based on virtual unit method |
-
2020
- 2020-12-21 CN CN202011523316.0A patent/CN112861305A/en active Pending
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113984515A (en) * | 2021-10-18 | 2022-01-28 | 同济大学 | Method for observing initiation and expansion of hydraulic fracture by using confocal microscope |
CN115906583A (en) * | 2022-12-16 | 2023-04-04 | 中国人民解放军陆军工程大学 | Method and system for simulating and analyzing structural integrity of explosive column based on virtual unit method |
CN115906583B (en) * | 2022-12-16 | 2023-08-01 | 中国人民解放军陆军工程大学 | Grain structure integrity simulation analysis method and system based on virtual unit method |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Farhat et al. | Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy‐based mesh sampling and weighting for computational efficiency | |
CN111859766B (en) | Lagrange integral point finite element numerical simulation system and method of variable calculation domain | |
Witteveen | Explicit and robust inverse distance weighting mesh deformation for CFD | |
Biancolini et al. | Static aeroelastic analysis of an aircraft wind-tunnel model by means of modal RBF mesh updating | |
US20230102815A1 (en) | Turbulence field update method and apparatus, and related device thereof | |
CN112861305A (en) | Crack propagation direction prediction method and device and storage medium | |
Liu et al. | Impact force identification via sparse regularization with generalized minimax-concave penalty | |
Noventa et al. | A high-order discontinuous Galerkin solver for unsteady incompressible turbulent flows | |
Moxey et al. | Optimising the performance of the spectral/hp element method with collective linear algebra operations | |
CN116861822B (en) | Cartesian grid-based object plane boundary processing method and device | |
CN116702571B (en) | Numerical simulation method and device based on multiple smoothness measurement factors | |
Tamamidis | A new upwind scheme on triangular meshes using the finite volume method | |
Bauer et al. | Simulation of tropical-cyclone-like vortices in shallow-water ICON-hex using goal-oriented r-adaptivity | |
CN112733415B (en) | Non-grid processing method and device for thin-wall elastomer boundary, terminal equipment and computing medium | |
CN114692529B (en) | CFD high-dimensional response uncertainty quantification method and device, and computer equipment | |
CN113901552B (en) | Method and device for establishing weir body model | |
Sonar | Classical finite volume methods | |
JP3587827B2 (en) | Airfoil performance estimation method and airfoil performance estimation program | |
Gan et al. | Super-convergent second derivative recovery for lower-order strain gradient plasticity | |
CN115995278B (en) | Method, device, equipment and readable storage medium for evaluating thermodynamic characteristics of material | |
Hashemabadi et al. | Implicit second-order CUSP gridless method for unsteady moving boundary simulations | |
CN115995279B (en) | Material mechanical property evaluation method, device, equipment and readable storage medium | |
CN116611369B (en) | Interpolation method and device based on smoothness magnitude and candidate template point number | |
Van Der Helm et al. | Comparison of artificial dissipation and limited flux schemes in Arbitrary Lagrangian–Eulerian finite element formulations | |
Leifsson et al. | Fast multi-objective aerodynamic optimization using space-mapping-corrected multi-fidelity models and kriging interpolation |
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 |